0025621: CAST analysis - Avoid constructors not supplying an initial value for all...
[occt.git] / src / ProjLib / ProjLib_ComputeApproxOnPolarSurface.cxx
1 // Created by: Bruno DUMORTIER
2 // Copyright (c) 1995-1999 Matra Datavision
3 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <ProjLib_ComputeApproxOnPolarSurface.hxx>
17 #include <AppCont_Function.hxx>
18 #include <ElSLib.hxx>
19 #include <ElCLib.hxx>
20 #include <BSplCLib.hxx>
21 #include <PLib.hxx>
22 #include <Standard_NoSuchObject.hxx>
23 #include <Geom_UndefinedDerivative.hxx>
24 #include <gp_Trsf.hxx>
25 #include <Precision.hxx>
26 #include <Approx_FitAndDivide2d.hxx>
27 #include <math.hxx>
28 #include <AppParCurves_MultiCurve.hxx>
29 #include <Geom_Surface.hxx>
30 #include <Geom2d_BSplineCurve.hxx>
31 #include <Geom2d_BezierCurve.hxx>
32 #include <Geom2d_Line.hxx>
33 #include <Geom2d_Circle.hxx>
34 #include <Geom2d_Ellipse.hxx>
35 #include <Geom2d_Hyperbola.hxx>
36 #include <Geom2d_Parabola.hxx>
37 #include <Geom2d_TrimmedCurve.hxx>
38 #include <Geom_BSplineSurface.hxx>
39 #include <Geom_BezierSurface.hxx>
40 #include <Geom_BSplineCurve.hxx>
41 #include <Geom_BezierCurve.hxx>
42 #include <Geom_TrimmedCurve.hxx>
43
44 #include <TColgp_Array1OfPnt2d.hxx>
45 #include <TColgp_Array2OfPnt2d.hxx>
46 #include <TColgp_Array1OfPnt.hxx>
47 #include <TColgp_SequenceOfPnt2d.hxx>
48 #include <TColStd_Array1OfReal.hxx>
49 #include <TColStd_Array1OfInteger.hxx>
50 #include <TColStd_SequenceOfReal.hxx>
51 #include <TColStd_ListOfTransient.hxx>
52
53 #include <GeomAbs_SurfaceType.hxx>
54 #include <GeomAbs_CurveType.hxx>
55 #include <Adaptor3d_Surface.hxx>
56 #include <Adaptor3d_Curve.hxx>
57 #include <Adaptor3d_HSurface.hxx>
58 #include <Adaptor3d_HCurve.hxx>
59 #include <Adaptor2d_HCurve2d.hxx>
60 #include <Geom2dAdaptor_Curve.hxx>
61 #include <Geom2dAdaptor_HCurve.hxx>
62 #include <GeomAdaptor_HCurve.hxx>
63 #include <GeomAdaptor.hxx>
64 #include <GeomAdaptor_Surface.hxx>
65 #include <TColgp_SequenceOfPnt.hxx>
66
67 #include <gp_Pnt.hxx>
68 #include <gp_Pnt2d.hxx>
69 #include <gp_Vec2d.hxx>
70 #include <Extrema_GenLocateExtPS.hxx>
71 #include <Extrema_ExtPS.hxx>
72 #include <GCPnts_QuasiUniformAbscissa.hxx>
73 #include <Standard_DomainError.hxx>
74 //#include <GeomLib_IsIso.hxx>
75 //#include <GeomLib_CheckSameParameter.hxx>
76
77 #ifdef OCCT_DEBUG
78 #ifdef DRAW
79 #include <DrawTrSurf.hxx>
80 #endif
81 //static Standard_Integer compteur = 0;
82 #endif
83
84 //=======================================================================
85 //function : Value
86 //purpose  : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
87 //            <InitCurve2d> use for calculate start 2D point.
88 //=======================================================================
89
90 static gp_Pnt2d Function_Value(const Standard_Real U,
91                                const Handle(Adaptor3d_HSurface)& Surf,    
92                                const Handle(Adaptor3d_HCurve)& Curve,
93                                const Handle(Adaptor2d_HCurve2d)& InitCurve2d,
94                                //OCC217
95                                const Standard_Real DistTol3d, const Standard_Real tolU, const Standard_Real tolV)
96                                //const Standard_Real Tolerance)
97 {
98   //OCC217
99   //Standard_Real Tol3d = 100*Tolerance;
100
101   gp_Pnt2d p2d = InitCurve2d->Value(U) ;
102   gp_Pnt p = Curve->Value(U);
103 //  Curve->D0(U,p) ;
104   Standard_Real Uinf, Usup, Vinf, Vsup;
105   Uinf = Surf->Surface().FirstUParameter();
106   Usup = Surf->Surface().LastUParameter();
107   Vinf = Surf->Surface().FirstVParameter();
108   Vsup = Surf->Surface().LastVParameter();
109   Standard_Integer decalU = 0, decalV = 0;
110   Standard_Real U0 = p2d.X(), V0 = p2d.Y();
111   
112   GeomAbs_SurfaceType Type = Surf->GetType();
113   if((Type != GeomAbs_BSplineSurface) && 
114      (Type != GeomAbs_BezierSurface)  &&
115      (Type != GeomAbs_OffsetSurface)    ) {
116     Standard_Real S = 0., T = 0.;
117     switch (Type) {
118 //    case GeomAbs_Plane:
119 //      {
120 //      gp_Pln Plane = Surf->Plane();
121 //      ElSLib::Parameters( Plane, p, S, T);
122 //      break;
123 //      }
124     case GeomAbs_Cylinder:
125       {
126         gp_Cylinder Cylinder = Surf->Cylinder();
127         ElSLib::Parameters( Cylinder, p, S, T);
128         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
129         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
130         S += decalU*2*M_PI;
131         break;
132       }
133     case GeomAbs_Cone:
134       {
135         gp_Cone Cone = Surf->Cone();
136         ElSLib::Parameters( Cone, p, S, T);
137         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
138         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
139         S += decalU*2*M_PI;
140         break;
141       }
142     case GeomAbs_Sphere:
143       {
144         gp_Sphere Sphere = Surf->Sphere();
145         ElSLib::Parameters( Sphere, p, S, T);
146         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
147         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
148         S += decalU*2*M_PI;
149         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
150         if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
151         T += decalV*2*M_PI;
152         if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI) {
153           T = M_PI - T;
154           if(U0 < S)
155             S -= M_PI;
156           else
157             S += M_PI;
158         }
159         break;
160       }
161     case GeomAbs_Torus:
162       {
163         gp_Torus Torus = Surf->Torus();
164         ElSLib::Parameters( Torus, p, S, T);
165         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
166         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
167         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
168         if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
169         S += decalU*2*M_PI; T += decalV*2*M_PI;
170         break;
171       }
172     default:
173       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::Value");
174     }
175     return gp_Pnt2d(S, T);
176   }
177   
178   //////////////////
179   Standard_Real Dist2Min = RealLast();
180   //OCC217
181   //Standard_Real tolU,tolV ;
182   //tolU = Tolerance;
183   //tolV = Tolerance;
184   
185   Standard_Real uperiod =0, vperiod = 0, u, v;
186   // U0 and V0 are the points within the initialized period 
187   // (periode with u and v),
188   // U1 and V1 are the points for construction of tops
189   
190   if(Surf->IsUPeriodic() || Surf->IsUClosed()) {
191     uperiod =  Surf->LastUParameter() - Surf->FirstUParameter();
192   }
193   if(Surf->IsVPeriodic() || Surf->IsVClosed()) {
194     vperiod = Surf->LastVParameter() - Surf->FirstVParameter();
195   } 
196   if(U0 < Uinf) {
197     if(!uperiod)
198       U0 = Uinf;
199     else {
200       decalU = int((Uinf - U0)/uperiod)+1;
201       U0 += decalU*uperiod;
202     }
203   }
204   if(U0 > Usup) {
205     if(!uperiod)
206       U0 = Usup;
207     else {
208       decalU = -(int((U0 - Usup)/uperiod)+1);
209       U0 += decalU*uperiod;
210     }
211   }
212   if(V0 < Vinf) {
213     if(!vperiod)
214       V0 = Vinf;
215     else {
216       decalV = int((Vinf - V0)/vperiod)+1;
217       V0 += decalV*vperiod;
218     }
219   }
220   if(V0 > Vsup) {
221     if(!vperiod)
222       V0 = Vsup;
223     else {
224       decalV = -int((V0 - Vsup)/vperiod)-1;
225       V0 += decalV*vperiod;
226     }
227   }
228   
229   // The surface around U0 is reduced
230   Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
231   Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
232   if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
233   if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
234   if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
235   if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
236   
237   //  const Adaptor3d_Surface GAS = Surf->Surface();
238
239   GeomAdaptor_Surface SurfLittle;
240   if (Type == GeomAbs_BSplineSurface) {
241     Handle(Geom_Surface) GBSS(Surf->Surface().BSpline());
242     SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
243   }
244   else if (Type == GeomAbs_BezierSurface) {
245     Handle(Geom_Surface) GS(Surf->Surface().Bezier());
246     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
247   }
248   else if (Type == GeomAbs_OffsetSurface) {
249     Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(Surf->Surface());
250     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
251   }
252   else {
253     Standard_NoSuchObject::Raise("");
254   }
255
256   Extrema_GenLocateExtPS  locext(p, SurfLittle, U0, V0, tolU, tolV);
257   if (locext.IsDone()) {
258     Dist2Min = locext.SquareDistance();
259     if (Dist2Min < DistTol3d * DistTol3d) {
260       (locext.Point()).Parameter(u, v);
261       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
262       return pnt;
263     }
264   } 
265   
266   Extrema_ExtPS  ext(p, SurfLittle, tolU, tolV) ;
267   if (ext.IsDone() && ext.NbExt()>=1 ) {
268     Dist2Min = ext.SquareDistance(1);
269     Standard_Integer GoodValue = 1;
270     for ( Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ ) 
271       if( Dist2Min > ext.SquareDistance(i)) {
272         Dist2Min = ext.SquareDistance(i);
273         GoodValue = i;
274       }
275     if (Dist2Min < DistTol3d * DistTol3d) {
276       (ext.Point(GoodValue)).Parameter(u,v);
277       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
278       return pnt;
279     }
280   }
281   
282   return p2d;
283 }
284
285
286 //=======================================================================
287 //function : ProjLib_PolarFunction
288 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
289 //=======================================================================
290
291 class ProjLib_PolarFunction : public AppCont_Function
292 {
293   Handle(Adaptor3d_HCurve)           myCurve;
294   Handle(Adaptor2d_HCurve2d)         myInitialCurve2d;
295   Handle(Adaptor3d_HSurface)         mySurface;
296   Standard_Real                      myTolU, myTolV;
297   Standard_Real                      myDistTol3d;
298
299   public :
300
301   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C, 
302                         const Handle(Adaptor3d_HSurface)& Surf,
303                         const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
304                         const Standard_Real Tol3d)
305 : myCurve(C),
306   myInitialCurve2d(InitialCurve2d),
307   mySurface(Surf),
308   myTolU(Surf->UResolution(Tol3d)),
309   myTolV(Surf->VResolution(Tol3d)),
310   myDistTol3d(100.0*Tol3d)
311   {
312     myNbPnt = 0;
313     myNbPnt2d = 1;
314   }
315
316   ~ProjLib_PolarFunction() {}
317
318   Standard_Real FirstParameter() const
319   {
320     return myCurve->FirstParameter();
321   }
322
323   Standard_Real LastParameter() const
324   {
325     return myCurve->LastParameter();
326   }
327
328   gp_Pnt2d Value(const Standard_Real t) const 
329   {
330   return Function_Value
331     (t,mySurface,myCurve,myInitialCurve2d,myDistTol3d,myTolU,myTolV);
332   }
333
334   Standard_Boolean Value(const Standard_Real   theT,
335                          NCollection_Array1<gp_Pnt2d>& thePnt2d,
336                          NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
337   {
338     thePnt2d(1) = Function_Value(theT, mySurface, myCurve, myInitialCurve2d, myDistTol3d, myTolU, myTolV);
339     return Standard_True;
340   }
341
342   Standard_Boolean D1(const Standard_Real   /*theT*/,
343                       NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
344                       NCollection_Array1<gp_Vec>&   /*theVec*/) const
345     {return Standard_False;}
346 };
347
348 //=======================================================================
349 //function : ProjLib_ComputeApproxOnPolarSurface
350 //purpose  : 
351 //=======================================================================
352
353 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
354 : myProjIsDone(Standard_False),
355   myTolerance (-1.0)
356 {
357 }
358
359
360 //=======================================================================
361 //function : ProjLib_ComputeApproxOnPolarSurface
362 //purpose  : 
363 //=======================================================================
364
365 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
366                     (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
367                      const Handle(Adaptor3d_HCurve)&   theCurve,
368                      const Handle(Adaptor3d_HSurface)& theSurface,
369                      const Standard_Real               theTolerance3D)
370 : myProjIsDone(Standard_False),
371   myTolerance (theTolerance3D)
372 {
373   myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
374 }
375 //=======================================================================
376 //function : ProjLib_ComputeApproxOnPolarSurface
377 //purpose  : Process the case of sewing
378 //=======================================================================
379
380 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
381                 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
382                  const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
383                  const Handle(Adaptor3d_HCurve)&   theCurve,
384                  const Handle(Adaptor3d_HSurface)& theSurface,
385                  const Standard_Real               theTolerance3D)
386 : myProjIsDone(Standard_False),
387   myTolerance (theTolerance3D)
388 {
389   // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing 
390   Handle(Geom2d_BSplineCurve) bsc =
391     Perform(theInitialCurve2d, theCurve, theSurface);
392
393   if(myProjIsDone) {
394     gp_Pnt2d P2dproj, P2d, P2dBis;
395     P2dproj = bsc->StartPoint();
396     P2d    = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
397     P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
398
399     Standard_Real Dist, DistBis;
400     Dist    = P2dproj.Distance(P2d);
401     DistBis = P2dproj.Distance(P2dBis);
402     if( Dist < DistBis)  {
403       // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
404       myBSpline = bsc;
405       Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
406       my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
407     }
408     else {
409       my2ndCurve = bsc;
410       Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
411       myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
412     }
413   }
414 }
415
416 //=======================================================================
417 //function : ProjLib_ComputeApproxOnPolarSurface
418 //purpose  : case without curve of initialization
419 //=======================================================================
420
421 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
422                       (const Handle(Adaptor3d_HCurve)&   theCurve,
423                        const Handle(Adaptor3d_HSurface)& theSurface,
424                        const Standard_Real               theTolerance3D)
425 : myProjIsDone(Standard_False),
426   myTolerance (theTolerance3D)
427 {
428   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
429   myBSpline = Perform(anInitCurve2d, theCurve, theSurface);  
430
431
432 //=======================================================================
433 //function : Concat
434 //purpose  : 
435 //=======================================================================
436
437 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
438                                           Handle(Geom2d_BSplineCurve) C2,
439                                           Standard_Real theUJump,
440                                           Standard_Real theVJump)
441 {
442   Standard_Integer deg, deg1, deg2;
443   deg1 = C1->Degree();
444   deg2 = C2->Degree();
445   
446   if ( deg1 < deg2) {
447     C1->IncreaseDegree(deg2);
448     deg = deg2;
449   }
450   else if ( deg2 < deg1) {
451     C2->IncreaseDegree(deg1);
452     deg = deg1;
453   }
454   else deg = deg1;
455
456   Standard_Integer np1,np2,nk1,nk2,np,nk;
457   np1 = C1->NbPoles();
458   nk1 = C1->NbKnots();
459   np2 = C2->NbPoles();
460   nk2 = C2->NbKnots();
461   nk = nk1 + nk2 -1;
462   np = np1 + np2 -1;
463
464   TColStd_Array1OfReal    K1(1,nk1); C1->Knots(K1);
465   TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
466   TColgp_Array1OfPnt2d    P1(1,np1); C1->Poles(P1);
467   TColStd_Array1OfReal    K2(1,nk2); C2->Knots(K2);
468   TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
469   TColgp_Array1OfPnt2d    P2(1,np2); C2->Poles(P2);
470
471   // Compute the new BSplineCurve
472   TColStd_Array1OfReal    K(1,nk);
473   TColStd_Array1OfInteger M(1,nk);
474   TColgp_Array1OfPnt2d    P(1,np);
475
476   Standard_Integer i, count = 0;
477   // Set Knots and Mults
478   for ( i = 1; i <= nk1; i++) {
479     count++;
480     K(count) = K1(i);
481     M(count) = M1(i);
482   }
483   M(count) = deg;
484   for ( i = 2; i <= nk2; i++) {
485     count++;
486     K(count) = K2(i);
487     M(count) = M2(i);
488   }
489   // Set the Poles
490   count = 0;
491   for (i = 1; i <= np1; i++) {
492     count++;
493     P(count) = P1(i);
494   }
495   for (i = 2; i <= np2; i++) {
496     count++;
497     P(count).SetX(P2(i).X() + theUJump);
498     P(count).SetY(P2(i).Y() + theVJump);
499   }
500
501   Handle(Geom2d_BSplineCurve) BS = 
502     new Geom2d_BSplineCurve(P,K,M,deg);
503   return BS;
504 }
505
506
507 //=======================================================================
508 //function : Perform
509 //purpose  : 
510 //=======================================================================
511 Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
512 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
513  const Handle(Adaptor3d_HCurve)& Curve,
514  const Handle(Adaptor3d_HSurface)& S) 
515 {
516   //OCC217
517   Standard_Real Tol3d = myTolerance; 
518   Standard_Real ParamTol = Precision::PApproximation();
519
520   Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
521   Handle(Adaptor3d_HCurve) AHC = Curve;
522   
523 // if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
524 // -> bug cts18237
525   GeomAbs_CurveType typeCurve = Curve->GetType();
526   if(typeCurve == GeomAbs_BSplineCurve) {
527     TColStd_ListOfTransient LOfBSpline2d;
528     Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
529     Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
530     if(nbInter > 1) {
531       Standard_Integer i, j;
532       Handle(Geom_TrimmedCurve) GTC;
533       Handle(Geom2d_TrimmedCurve) G2dTC;
534       TColStd_Array1OfReal Inter(1,nbInter+1);
535       Curve->Intervals(Inter,GeomAbs_C1);
536       Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
537       // initialization 3d
538       GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
539       AHC = new GeomAdaptor_HCurve(GTC);
540       
541       // if there is an initialization curve: 
542       // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
543       // and the sections C1 are taken
544       // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
545       
546       // initialization 2d
547       Standard_Integer nbInter2d;
548       Standard_Boolean C2dIsToCompute;
549       C2dIsToCompute = InitialCurve2d.IsNull();
550       Handle(Geom2d_BSplineCurve) BSC2d;
551       Handle(Geom2d_Curve) G2dC;
552       
553       if(!C2dIsToCompute) {
554         nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
555         TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
556         InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
557         j = 1;
558         for(i = 1,j = 1;i <= nbInter;i++)
559           if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
560           //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
561             if (j > nbInter2d) break;
562             j++;
563           }
564         if(j != (nbInter2d+1)) {
565           C2dIsToCompute = Standard_True;
566         }
567       }
568       
569       if(C2dIsToCompute) {
570         AHC2d = BuildInitialCurve2d(AHC, S);
571       }
572       else {
573         typeCurve = InitialCurve2d->GetType();
574         switch (typeCurve) {
575         case GeomAbs_Line: {
576           G2dC = new Geom2d_Line(InitialCurve2d->Line());
577           break;
578         }
579         case GeomAbs_Circle: {
580           G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
581           break;
582         }
583         case GeomAbs_Ellipse: {
584           G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
585           break;
586         }
587         case GeomAbs_Hyperbola: {
588           G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
589           break;
590         }
591         case GeomAbs_Parabola: {
592           G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
593           break;
594         }
595         case GeomAbs_BezierCurve: {
596           G2dC = InitialCurve2d->Bezier();
597           break;
598         }
599         case GeomAbs_BSplineCurve: {
600           G2dC = InitialCurve2d->BSpline();
601           break;
602         }
603         case GeomAbs_OtherCurve:
604         default:
605           break;
606         }
607         gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
608         gp_Pnt fps, lps, fpc, lpc;
609         S->D0(fp2d.X(), fp2d.Y(), fps);
610         S->D0(lp2d.X(), lp2d.Y(), lps);
611         Curve->D0(firstinter, fpc);
612         Curve->D0(secondinter, lpc);
613         //OCC217
614         if((fps.IsEqual(fpc, Tol3d)) &&
615            (lps.IsEqual(lpc, Tol3d))) {
616         //if((fps.IsEqual(fpc, myTolerance)) &&
617         //   (lps.IsEqual(lpc, myTolerance))) {
618           G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
619           Geom2dAdaptor_Curve G2dAC(G2dTC);
620           AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
621           myProjIsDone = Standard_True;
622         }
623         else {
624           AHC2d = BuildInitialCurve2d(AHC, S);
625           C2dIsToCompute = Standard_True;
626         }
627       }
628         
629       if(myProjIsDone) {
630         BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);  
631         if(BSC2d.IsNull()) return Handle(Geom2d_BSplineCurve)(); //IFV
632         LOfBSpline2d.Append(BSC2d);
633       }
634       else {
635         return Handle(Geom2d_BSplineCurve)();
636       }
637       
638
639
640       Standard_Real iinter, ip1inter;
641       Standard_Integer nbK2d, deg;
642       nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
643
644       for(i = 2;i <= nbInter;i++) {
645         iinter = Inter.Value(i);
646         ip1inter = Inter.Value(i+1);
647         // general case 3d
648         GTC->SetTrim(iinter, ip1inter);
649         AHC = new GeomAdaptor_HCurve(GTC);
650         
651         // general case 2d
652         if(C2dIsToCompute) {
653           AHC2d = BuildInitialCurve2d(AHC, S);
654         }
655         else {
656           gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
657           gp_Pnt fps, lps, fpc, lpc;
658           S->D0(fp2d.X(), fp2d.Y(), fps);
659           S->D0(lp2d.X(), lp2d.Y(), lps);
660           Curve->D0(iinter, fpc);
661           Curve->D0(ip1inter, lpc);
662           //OCC217
663           if((fps.IsEqual(fpc, Tol3d)) &&
664              (lps.IsEqual(lpc, Tol3d))) {
665           //if((fps.IsEqual(fpc, myTolerance)) &&
666           //   (lps.IsEqual(lpc, myTolerance))) {
667             G2dTC->SetTrim(iinter, ip1inter);
668             Geom2dAdaptor_Curve G2dAC(G2dTC);
669             AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
670             myProjIsDone = Standard_True;
671           }
672           else {
673             AHC2d = BuildInitialCurve2d(AHC, S);
674           }
675         }
676         if(myProjIsDone) {
677           BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);  
678           if(BSC2d.IsNull()) {
679             return Handle(Geom2d_BSplineCurve)();
680           }
681           LOfBSpline2d.Append(BSC2d);
682           nbK2d += BSC2d->NbKnots() - 1;
683           deg = Max(deg, BSC2d->Degree());
684         }
685         else {
686           return Handle(Geom2d_BSplineCurve)();
687         }
688       }
689
690       Standard_Integer NbC = LOfBSpline2d.Extent();
691       Handle(Geom2d_BSplineCurve) CurBS;
692       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
693       LOfBSpline2d.RemoveFirst();
694       for (Standard_Integer ii = 2; ii <= NbC; ii++)
695       {
696         Handle(Geom2d_BSplineCurve) BS = 
697           Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
698
699         //Check for period jump in point of contact.
700         gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
701         gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
702         Standard_Real anUJump = 0.0, anVJump = 0.0;
703
704         if (S->IsUPeriodic() || S->IsUClosed())
705         {
706           if (Abs (aC1End.X() - aC2Beg.X()) > (S->LastUParameter() - S->FirstUParameter() ) / 2.01)
707           {
708             Standard_Real aMultCoeff =  aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
709             anUJump = (S->LastUParameter() - S->FirstUParameter() ) * aMultCoeff;
710           }
711         }
712
713         if (S->IsVPeriodic() || S->IsVClosed())
714         {
715           if (Abs (aC1End.Y() - aC2Beg.Y()) > (S->LastVParameter() - S->FirstVParameter() ) / 2.01)
716           {
717             Standard_Real aMultCoeff =  aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
718             anVJump = (S->LastVParameter() - S->FirstVParameter() ) * aMultCoeff;
719           }
720         }
721
722         CurBS = Concat(CurBS,BS, anUJump, anVJump);
723         LOfBSpline2d.RemoveFirst();
724       }
725       return CurBS;
726     }
727   }
728   
729   if(InitialCurve2d.IsNull()) {
730     AHC2d = BuildInitialCurve2d(Curve, S);
731     if(!myProjIsDone) 
732       return Handle(Geom2d_BSplineCurve)(); 
733   }
734   return ProjectUsingInitialCurve2d(AHC, S, AHC2d);    
735 }
736
737 //=======================================================================
738 //function : ProjLib_BuildInitialCurve2d
739 //purpose  : 
740 //=======================================================================
741
742 Handle(Adaptor2d_HCurve2d) 
743      ProjLib_ComputeApproxOnPolarSurface::
744      BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)&   Curve,
745                          const Handle(Adaptor3d_HSurface)& Surf)
746 {
747   //  discretize the Curve with quasiuniform deflection
748   //  density at least NbOfPnts points
749   myProjIsDone = Standard_False;
750   
751   //OCC217
752   Standard_Real Tol3d = myTolerance; 
753   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
754   Standard_Real DistTol3d = 100.0*Tol3d;
755
756   Standard_Real uperiod = 0., vperiod = 0.;
757   if(Surf->IsUPeriodic() || Surf->IsUClosed())
758     uperiod = Surf->LastUParameter() - Surf->FirstUParameter(); 
759   
760   if(Surf->IsVPeriodic() || Surf->IsVClosed())
761     vperiod = Surf->LastVParameter() - Surf->FirstVParameter(); 
762
763   
764   // NO myTol is Tol2d !!!!
765   //Standard_Real TolU = myTolerance, TolV = myTolerance;
766   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
767
768   Standard_Integer NbOfPnts = 61; 
769   GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
770   TColgp_Array1OfPnt Pts(1,NbOfPnts);
771   TColStd_Array1OfReal Param(1,NbOfPnts);
772   Standard_Integer i, j;
773   for( i = 1; i <= NbOfPnts ; i++ ) { 
774     Param(i) = QUA.Parameter(i);
775     Pts(i) = Curve->Value(Param(i));
776   }
777   
778   TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
779   TColStd_Array1OfInteger Mult(1,NbOfPnts);
780   Mult.Init(1);
781   Mult(1) = Mult(NbOfPnts) = 2;
782   
783   Standard_Real Uinf, Usup, Vinf, Vsup;
784   Uinf = Surf->Surface().FirstUParameter();
785   Usup = Surf->Surface().LastUParameter();
786   Vinf = Surf->Surface().FirstVParameter();
787   Vsup = Surf->Surface().LastVParameter();
788   GeomAbs_SurfaceType Type = Surf->GetType();
789   if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
790      (Type != GeomAbs_OffsetSurface)) {
791     Standard_Real S, T;
792 //    Standard_Integer usens = 0, vsens = 0; 
793     // to know the position relatively to the period
794     switch (Type) {
795 //    case GeomAbs_Plane:
796 //      {
797 //      gp_Pln Plane = Surf->Plane();
798 //      for ( i = 1 ; i <= NbOfPnts ; i++) { 
799 //        ElSLib::Parameters( Plane, Pts(i), S, T);
800 //        Pts2d(i).SetCoord(S,T);
801 //      }
802 //      myProjIsDone = Standard_True;
803 //      break;
804 //      }
805     case GeomAbs_Cylinder:
806       {
807 //      Standard_Real Sloc, Tloc;
808         Standard_Real Sloc;
809         Standard_Integer usens = 0;
810         gp_Cylinder Cylinder = Surf->Cylinder();
811         ElSLib::Parameters( Cylinder, Pts(1), S, T);
812         Pts2d(1).SetCoord(S,T);
813         for ( i = 2 ; i <= NbOfPnts ; i++) { 
814           Sloc = S;
815           ElSLib::Parameters( Cylinder, Pts(i), S, T);
816           if(Abs(Sloc - S) > M_PI) {
817             if(Sloc > S)
818               usens++;
819             else
820               usens--;
821           }
822           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
823         }
824         myProjIsDone = Standard_True;
825         break;
826       }
827     case GeomAbs_Cone:
828       {
829 //      Standard_Real Sloc, Tloc;
830         Standard_Real Sloc;
831         Standard_Integer usens = 0;
832         gp_Cone Cone = Surf->Cone();
833         ElSLib::Parameters( Cone, Pts(1), S, T);
834         Pts2d(1).SetCoord(S,T);
835         for ( i = 2 ; i <= NbOfPnts ; i++) { 
836           Sloc = S;
837           ElSLib::Parameters( Cone, Pts(i), S, T);
838           if(Abs(Sloc - S) > M_PI) {
839             if(Sloc > S)
840               usens++;
841             else
842               usens--;
843           }
844           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
845         }
846         myProjIsDone = Standard_True;
847         break;
848       }
849     case GeomAbs_Sphere:
850       {
851         Standard_Real Sloc, Tloc;
852         Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
853         Standard_Boolean vparit = Standard_False;
854         gp_Sphere Sphere = Surf->Sphere();
855         ElSLib::Parameters( Sphere, Pts(1), S, T);
856         Pts2d(1).SetCoord(S,T);
857         for ( i = 2 ; i <= NbOfPnts ; i++) { 
858           Sloc = S;Tloc = T;
859           ElSLib::Parameters( Sphere, Pts(i), S, T);
860           if(1.6*M_PI < Abs(Sloc - S)) {
861             if(Sloc > S)
862               usens += 2;
863             else
864               usens -= 2;
865     }
866           if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
867             vparit = !vparit;
868             if(Sloc > S)
869               usens++;
870             else
871               usens--;
872             if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
873               vsens++;
874             else
875               vsens--;
876           }
877           if(vparit) {
878             Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
879           }       
880           else {
881             Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
882             
883           }
884         }
885         myProjIsDone = Standard_True;
886         break;
887       }
888     case GeomAbs_Torus:
889       {
890         Standard_Real Sloc, Tloc;
891         Standard_Integer usens = 0, vsens = 0;
892         gp_Torus Torus = Surf->Torus();
893         ElSLib::Parameters( Torus, Pts(1), S, T);
894         Pts2d(1).SetCoord(S,T);
895         for ( i = 2 ; i <= NbOfPnts ; i++) { 
896           Sloc = S; Tloc = T;
897           ElSLib::Parameters( Torus, Pts(i), S, T);
898           if(Abs(Sloc - S) > M_PI) {
899             if(Sloc > S)
900               usens++;
901             else
902               usens--;
903     }
904           if(Abs(Tloc - T) > M_PI) {
905             if(Tloc > T)
906               vsens++;
907             else
908               vsens--;
909     }
910           Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
911         }
912         myProjIsDone = Standard_True;
913         break;
914       }
915     default:
916       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
917     }
918   }
919   else {
920     myProjIsDone = Standard_False;
921     Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
922     gp_Pnt pntproj;
923
924     TColgp_SequenceOfPnt2d Sols;
925     Standard_Boolean areManyZeros = Standard_False;
926     
927     Curve->D0(Param.Value(1), pntproj) ;
928     Extrema_ExtPS  aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
929     Standard_Real aMinSqDist = RealLast();
930     if (aExtPS.IsDone())
931     {
932       for (i = 1; i <= aExtPS.NbExt(); i++)
933       {
934         Standard_Real aSqDist = aExtPS.SquareDistance(i);
935         if (aSqDist < aMinSqDist)
936           aMinSqDist = aSqDist;
937       }
938     }
939     if (aMinSqDist > DistTol3d * DistTol3d) //try to project with less tolerance
940     {
941       TolU = Min(TolU, Precision::PConfusion());
942       TolV = Min(TolV, Precision::PConfusion());
943       aExtPS.Initialize(Surf->Surface(),
944                         Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(), 
945                         Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
946                         TolU, TolV);
947       aExtPS.Perform(pntproj);
948     }
949
950     if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
951
952       Standard_Integer GoodValue = 1;
953
954       for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
955         if( aExtPS.SquareDistance(i) < DistTol3d * DistTol3d ) {
956           if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
957             aExtPS.Point(i).Parameter(u,v);
958             gp_Pnt2d p2d(u,v);
959             Standard_Boolean isSame = Standard_False;
960             for( j = 1; j <= Sols.Length(); j++ ) {
961               if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
962                 isSame = Standard_True;
963                 break;
964               }
965             }
966             if( !isSame ) Sols.Append( p2d );
967           }
968           if( Dist2Min > aExtPS.SquareDistance(i) ) {
969             Dist2Min = aExtPS.SquareDistance(i);
970             GoodValue = i;
971           }
972         }
973       }
974
975       if( Sols.Length() > 1 ) areManyZeros = Standard_True;
976
977       if( Dist2Min <= DistTol3d * DistTol3d) {
978         if( !areManyZeros ) {
979           aExtPS.Point(GoodValue).Parameter(u,v);
980           Pts2d(1).SetCoord(u,v);
981           myProjIsDone = Standard_True;
982         }
983         else {
984           Standard_Integer nbSols = Sols.Length();
985           Standard_Real Dist2Max = -1.e+200;
986           for( i = 1; i <= nbSols; i++ ) {
987             const gp_Pnt2d& aP1 = Sols.Value(i);
988             for( j = i+1; j <= nbSols; j++ ) {
989               const gp_Pnt2d& aP2 = Sols.Value(j);
990               Standard_Real aDist2 = aP1.SquareDistance(aP2);
991               if( aDist2 > Dist2Max ) Dist2Max = aDist2;
992             }
993           }
994       Standard_Real aMaxT2 = Max(TolU,TolV);
995       aMaxT2 *= aMaxT2;
996           if( Dist2Max > aMaxT2 ) {
997             Standard_Integer tPp = 0;
998             for( i = 1; i <= 5; i++ ) {
999               Standard_Integer nbExtOk = 0;
1000               Standard_Integer indExt = 0;
1001               Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
1002               Curve->D0( Param.Value(iT), pntproj );
1003               Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1004               Dist2Min = 1.e+200;
1005               if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1006                 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
1007                   if( aTPS.SquareDistance(j) < DistTol3d * DistTol3d ) {
1008                     nbExtOk++;
1009                     if( aTPS.SquareDistance(j) < Dist2Min ) {
1010                       Dist2Min = aTPS.SquareDistance(j);
1011                       indExt = j;
1012                     }
1013                   }
1014                 }
1015               }
1016               if( nbExtOk == 1 ) {
1017                 tPp = iT;
1018                 aTPS.Point(indExt).Parameter(u,v);
1019                 break;
1020               }
1021             }
1022
1023             if( tPp != 0 ) {
1024               gp_Pnt2d aPp = gp_Pnt2d(u,v);
1025               gp_Pnt2d aPn;
1026               j = 1;
1027               Standard_Boolean isFound = Standard_False;
1028               while( !isFound ) { 
1029                 Curve->D0( Param.Value(tPp+j), pntproj );
1030                 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1031                 Dist2Min = 1.e+200;
1032                 Standard_Integer indExt = 0;
1033                 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1034                   for( i = 1 ; i <= aTPS.NbExt() ; i++ ) {
1035                     if( aTPS.SquareDistance(i) < DistTol3d * DistTol3d && aTPS.SquareDistance(i) < Dist2Min ) {
1036                       Dist2Min = aTPS.SquareDistance(i);
1037                       indExt = i;
1038                       isFound = Standard_True;
1039                     }
1040                   }
1041                 }
1042                 if( isFound ) {
1043                   aTPS.Point(indExt).Parameter(u,v);
1044                   aPn = gp_Pnt2d(u,v);
1045                   break;
1046                 }
1047                 j++;
1048                 if( (tPp+j) > NbOfPnts ) break;
1049               }
1050
1051               if( isFound ) {
1052                 gp_Vec2d atV(aPp,aPn);
1053                 Standard_Boolean isChosen = Standard_False;
1054                 for( i = 1; i <= nbSols; i++ ) {
1055                   const gp_Pnt2d& aP1 = Sols.Value(i);
1056                   gp_Vec2d asV(aP1,aPp);
1057                   if( asV.Dot(atV) > 0. ) {
1058                     isChosen = Standard_True;
1059                     Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1060                     myProjIsDone = Standard_True;
1061                     break;
1062                   }
1063                 }
1064                 if( !isChosen ) {
1065                   aExtPS.Point(GoodValue).Parameter(u,v);
1066                   Pts2d(1).SetCoord(u,v);
1067                   myProjIsDone = Standard_True;
1068                 }
1069               }
1070               else {
1071                 aExtPS.Point(GoodValue).Parameter(u,v);
1072                 Pts2d(1).SetCoord(u,v);
1073                 myProjIsDone = Standard_True;
1074               }
1075             }
1076             else {
1077               aExtPS.Point(GoodValue).Parameter(u,v);
1078               Pts2d(1).SetCoord(u,v);
1079               myProjIsDone = Standard_True;
1080             }
1081           }
1082           else {
1083             aExtPS.Point(GoodValue).Parameter(u,v);
1084             Pts2d(1).SetCoord(u,v);
1085             myProjIsDone = Standard_True;
1086           }
1087         }
1088       }
1089       
1090       //  calculate the following points with GenLocate_ExtPS
1091       // (and store the result and each parameter in a sequence)
1092       Standard_Integer usens = 0, vsens = 0; 
1093       // to know the position relatively to the period
1094       Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1095       // U0 and V0 are the points in the initialized period 
1096       // (period with u and v),
1097       // U1 and V1 are the points for construction of poles
1098
1099       for ( i = 2 ; i <= NbOfPnts ; i++) 
1100         if(myProjIsDone) {
1101           myProjIsDone = Standard_False;
1102           Dist2Min = RealLast();
1103           Curve->D0(Param.Value(i), pntproj);
1104           Extrema_GenLocateExtPS  aLocateExtPS
1105             (pntproj, Surf->Surface(), U0, V0, TolU, TolV) ;
1106           
1107           if (aLocateExtPS.IsDone())
1108           {
1109             if (aLocateExtPS.SquareDistance() < DistTol3d * DistTol3d)
1110             {  //OCC217
1111               //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1112               (aLocateExtPS.Point()).Parameter(U0,V0);
1113               U1 = U0 + usens*uperiod;
1114               V1 = V0 + vsens*vperiod;
1115               Pts2d(i).SetCoord(U1,V1);
1116               myProjIsDone = Standard_True;
1117             }
1118             else
1119             {
1120               Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1121               if (aGlobalExtr.IsDone())
1122               {
1123                 Standard_Real LocalMinSqDist = RealLast();
1124                 Standard_Integer imin = 0;
1125                 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1126                 {
1127                   Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1128                   if (aSqDist < LocalMinSqDist)
1129                   {
1130                     LocalMinSqDist = aSqDist;
1131                     imin = isol;
1132                   }
1133                 }
1134                 if (LocalMinSqDist < DistTol3d * DistTol3d)
1135                 {
1136                   Standard_Real LocalU, LocalV;
1137                   aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1138                   if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1139                   {
1140                     if (LocalU > U0)
1141                       usens = -1;
1142                     else
1143                       usens = 1;
1144                   }
1145                   if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1146                   {
1147                     if (LocalV > V0)
1148                       vsens = -1;
1149                     else
1150                       vsens = 1;
1151                   }
1152                   U0 = LocalU; V0 = LocalV;
1153                   U1 = U0 + usens*uperiod;
1154                   V1 = V0 + vsens*vperiod;
1155                   Pts2d(i).SetCoord(U1,V1);
1156                   myProjIsDone = Standard_True;
1157
1158                   if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
1159                   {//Make 1st point more precise for periodic surfaces
1160                     const Standard_Integer aSize = 3;
1161                     const gp_Pnt2d aP(Pts2d(2)); 
1162                     Standard_Real aUpar[aSize], aVpar[aSize];
1163                     Pts2d(1).Coord(aUpar[1], aVpar[1]);
1164                     aUpar[0] = aUpar[1] - uperiod;
1165                     aUpar[2] = aUpar[1] + uperiod;
1166                     aVpar[0] = aVpar[1] - vperiod;
1167                     aVpar[2] = aVpar[1] + vperiod;
1168
1169                     Standard_Real aSQdistMin = RealLast();
1170                     Standard_Integer aBestUInd = 1, aBestVInd = 1;
1171                     const Standard_Integer  aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
1172                                             aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
1173                     for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
1174                     {
1175                       for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
1176                       {
1177                         Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
1178                         if(aSQdist < aSQdistMin)
1179                         {
1180                           aSQdistMin = aSQdist;
1181                           aBestUInd = uInd;
1182                           aBestVInd = vInd;
1183                         }
1184                       }
1185                     }
1186
1187                     Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
1188                   }//if(i == 2) condition
1189                 }
1190               }
1191             }
1192           }
1193           if(!myProjIsDone && uperiod) {
1194             Standard_Real Uinf, Usup, Uaux;
1195             Uinf = Surf->Surface().FirstUParameter();
1196             Usup = Surf->Surface().LastUParameter();
1197             if((Usup - U0) > (U0 - Uinf)) 
1198               Uaux = 2*Uinf - U0 + uperiod;
1199             else 
1200               Uaux = 2*Usup - U0 - uperiod;
1201             Extrema_GenLocateExtPS  locext(pntproj, 
1202                                            Surf->Surface(), 
1203                                            Uaux, V0, TolU, TolV);
1204             if (locext.IsDone())
1205               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1206               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1207                 (locext.Point()).Parameter(u,v);
1208                 if((Usup - U0) > (U0 - Uinf)) 
1209                   usens--;
1210                 else 
1211                   usens++;
1212                 U0 = u; V0 = v;
1213                 U1 = U0 + usens*uperiod;
1214                 V1 = V0 + vsens*vperiod;
1215                 Pts2d(i).SetCoord(U1,V1);
1216                 myProjIsDone = Standard_True;
1217               }
1218           }
1219           if(!myProjIsDone && vperiod) {
1220             Standard_Real Vinf, Vsup, Vaux;
1221             Vinf = Surf->Surface().FirstVParameter();
1222             Vsup = Surf->Surface().LastVParameter();
1223             if((Vsup - V0) > (V0 - Vinf)) 
1224               Vaux = 2*Vinf - V0 + vperiod;
1225             else 
1226               Vaux = 2*Vsup - V0 - vperiod;
1227             Extrema_GenLocateExtPS  locext(pntproj, 
1228                                            Surf->Surface(), 
1229                                            U0, Vaux, TolU, TolV) ;
1230             if (locext.IsDone())
1231               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1232               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1233                 (locext.Point()).Parameter(u,v);
1234                 if((Vsup - V0) > (V0 - Vinf)) 
1235                   vsens--;
1236                 else 
1237                   vsens++;
1238                 U0 = u; V0 = v;
1239                 U1 = U0 + usens*uperiod;
1240                 V1 = V0 + vsens*vperiod;
1241                 Pts2d(i).SetCoord(U1,V1);
1242                 myProjIsDone = Standard_True;
1243               }
1244           }     
1245           if(!myProjIsDone && uperiod && vperiod) {
1246             Standard_Real Uaux, Vaux;
1247             if((Usup - U0) > (U0 - Uinf)) 
1248               Uaux = 2*Uinf - U0 + uperiod;
1249             else 
1250               Uaux = 2*Usup - U0 - uperiod;
1251             if((Vsup - V0) > (V0 - Vinf)) 
1252               Vaux = 2*Vinf - V0 + vperiod;
1253             else 
1254               Vaux = 2*Vsup - V0 - vperiod;
1255             Extrema_GenLocateExtPS  locext(pntproj, 
1256                                            Surf->Surface(), 
1257                                            Uaux, Vaux, TolU, TolV);
1258             if (locext.IsDone())
1259               if (locext.SquareDistance() < DistTol3d * DistTol3d) {
1260               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1261                 (locext.Point()).Parameter(u,v);
1262                 if((Usup - U0) > (U0 - Uinf)) 
1263                   usens--;
1264                 else 
1265                   usens++;
1266                 if((Vsup - V0) > (V0 - Vinf)) 
1267                   vsens--;
1268                 else 
1269                   vsens++;
1270                 U0 = u; V0 = v;
1271                 U1 = U0 + usens*uperiod;
1272                 V1 = V0 + vsens*vperiod;
1273                 Pts2d(i).SetCoord(U1,V1);
1274                 myProjIsDone = Standard_True;
1275               }
1276           }
1277           if(!myProjIsDone) {
1278             Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1279             if (ext.IsDone()) {
1280               Dist2Min = ext.SquareDistance(1);
1281               Standard_Integer GoodValue = 1;
1282               for ( j = 2 ; j <= ext.NbExt() ; j++ )
1283                 if( Dist2Min > ext.SquareDistance(j)) {
1284                   Dist2Min = ext.SquareDistance(j);
1285                   GoodValue = j;
1286                 }
1287               if (Dist2Min < DistTol3d * DistTol3d) {
1288               //if (Dist2Min < Tol3d * Tol3d) {
1289                 (ext.Point(GoodValue)).Parameter(u,v);
1290                 if(uperiod) {
1291                   if((U0 - u) > (2*uperiod/3)) {
1292                     usens++;
1293                   }
1294                   else
1295                     if((u - U0) > (2*uperiod/3)) {
1296                       usens--;
1297                     }
1298     }
1299                 if(vperiod) {
1300                   if((V0 - v) > (vperiod/2)) {
1301                     vsens++;
1302                   }
1303                   else
1304                     if((v - V0) > (vperiod/2)) {
1305                       vsens--;
1306                     }
1307       }
1308                 U0 = u; V0 = v;
1309                 U1 = U0 + usens*uperiod;
1310                 V1 = V0 + vsens*vperiod;
1311                 Pts2d(i).SetCoord(U1,V1);
1312                 myProjIsDone = Standard_True;
1313               }
1314             }
1315           }
1316         }
1317         else break;
1318     }    
1319   }
1320   // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1321   if(myProjIsDone) {
1322     myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1323     //jgv: put the curve into parametric range
1324     gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1325     Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1326     Standard_Real sense = 0.;
1327     if (uperiod)
1328     {
1329       if (TestU < Uinf - TolU)
1330         sense = 1.;
1331       else if (TestU > Usup + TolU)
1332         sense = -1;
1333       while (TestU < Uinf - TolU || TestU > Usup + TolU)
1334         TestU += sense * uperiod;
1335     }
1336     if (vperiod)
1337     {
1338       sense = 0.;
1339       if (TestV < Vinf - TolV)
1340         sense = 1.;
1341       else if (TestV > Vsup + TolV)
1342         sense = -1.;
1343       while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1344         TestV += sense * vperiod;
1345     }
1346     gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1347     if (Abs(Offset.X()) > gp::Resolution() ||
1348         Abs(Offset.Y()) > gp::Resolution())
1349       myBSpline->Translate(Offset);
1350     //////////////////////////////////////////
1351     Geom2dAdaptor_Curve GAC(myBSpline);
1352     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1353 #ifdef OCCT_DEBUG
1354 //    char name [100];
1355 //    sprintf(name,"%s_%d","build",compteur++);
1356 //    DrawTrSurf::Set(name,myBSpline);
1357 #endif
1358     return IC2d;
1359   }
1360   else {
1361 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1362 //     Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1363 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1364     return Handle(Adaptor2d_HCurve2d)();
1365   }
1366 //  myProjIsDone = Standard_False;
1367 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1368 //   Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1369 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1370 }
1371
1372
1373
1374
1375 //=======================================================================
1376 //function : ProjLib_ProjectUsingInitialCurve2d
1377 //purpose  : 
1378 //=======================================================================
1379 Handle(Geom2d_BSplineCurve) 
1380      ProjLib_ComputeApproxOnPolarSurface::
1381      ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1382                                 const Handle(Adaptor3d_HSurface)& Surf,
1383                                 const Handle(Adaptor2d_HCurve2d)& InitCurve2d) 
1384 {  
1385   //OCC217
1386   Standard_Real Tol3d = myTolerance;
1387   Standard_Real DistTol3d = 1.0*Tol3d;
1388   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1389   Standard_Real Tol2d = Sqrt(TolU*TolU + TolV*TolV);
1390
1391   Standard_Integer i;
1392   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1393   GeomAbs_CurveType TheTypeC = Curve->GetType();
1394 //  Handle(Standard_Type) TheTypeS = Surf->DynamicType();
1395 //  Handle(Standard_Type) TheTypeC = Curve->DynamicType();   // si on a :
1396 //  if(TheTypeS == STANDARD_TYPE(Geom_BSplineSurface)) {
1397   if(TheTypeS == GeomAbs_Plane) {
1398     Standard_Real S, T;
1399     gp_Pln Plane = Surf->Plane();
1400     if(TheTypeC == GeomAbs_BSplineCurve) {
1401       Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1402       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1403       for(i = 1;i <= Curve->NbPoles();i++) {
1404         ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1405         Poles2d(i).SetCoord(S,T);
1406       }
1407       TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1408       BSC->Knots(Knots);
1409       TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1410       BSC->Multiplicities(Mults);
1411       if(BSC->IsRational()) {
1412         TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1413         BSC->Weights(Weights); 
1414         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1415                                        BSC->Degree(), BSC->IsPeriodic()) ;
1416       }
1417       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1418                                      BSC->Degree(), BSC->IsPeriodic()) ;
1419       
1420     }
1421     if(TheTypeC == GeomAbs_BezierCurve) {
1422       Handle(Geom_BezierCurve) BC = Curve->Bezier();
1423       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1424       for(i = 1;i <= Curve->NbPoles();i++) {
1425         ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1426         Poles2d(i).SetCoord(S,T);
1427       }
1428       TColStd_Array1OfReal Knots(1, 2);
1429       Knots.SetValue(1,0.0);
1430       Knots.SetValue(2,1.0);
1431       TColStd_Array1OfInteger Mults(1, 2);
1432       Mults.Init(BC->NbPoles());
1433       if(BC->IsRational()) {
1434         TColStd_Array1OfReal Weights(1, BC->NbPoles());
1435         BC->Weights(Weights); 
1436         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1437                                        BC->Degree(), BC->IsPeriodic()) ;
1438       }
1439       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1440                                      BC->Degree(), BC->IsPeriodic()) ;
1441     }
1442   }
1443   if(TheTypeS == GeomAbs_BSplineSurface) {
1444     Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1445     if((BSS->MaxDegree() == 1) &&
1446        (BSS->NbUPoles() == 2) &&
1447        (BSS->NbVPoles() == 2)) {
1448       gp_Pnt p11 = BSS->Pole(1,1);
1449       gp_Pnt p12 = BSS->Pole(1,2);
1450       gp_Pnt p21 = BSS->Pole(2,1);
1451       gp_Pnt p22 = BSS->Pole(2,2);
1452       gp_Vec V1(p11,p12);
1453       gp_Vec V2(p21,p22);
1454       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){  //OCC217
1455       //if(V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1456         //  so the polar surface is plane
1457         //  and if it is enough to projet the  poles of Curve
1458         Standard_Integer Dist2Min = IntegerLast();
1459         Standard_Real u,v;
1460         //OCC217
1461         //Standard_Real TolU = Surf->UResolution(myTolerance)
1462         //  , TolV = Surf->VResolution(myTolerance);
1463 //      gp_Pnt pntproj;
1464         if(TheTypeC == GeomAbs_BSplineCurve) {
1465           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1466           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1467           for(i = 1;i <= Curve->NbPoles();i++) {
1468             myProjIsDone = Standard_False;
1469             Dist2Min = IntegerLast();
1470             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1471                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1472             if (extrloc.IsDone()) {
1473               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1474               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1475               //if (Dist2Min < myTolerance * myTolerance) {
1476                 (extrloc.Point()).Parameter(u,v);
1477                 Poles2d(i).SetCoord(u,v);
1478                 myProjIsDone = Standard_True;
1479               }
1480               else break;
1481             }
1482             else break;
1483             if(!myProjIsDone) 
1484               break;
1485           }
1486           if(myProjIsDone) {
1487             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1488             BSC->Knots(Knots);
1489             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1490             BSC->Multiplicities(Mults);
1491             if(BSC->IsRational()) {
1492               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1493               BSC->Weights(Weights); 
1494               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1495                                              BSC->Degree(), BSC->IsPeriodic()) ;
1496             }
1497             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1498                                            BSC->Degree(), BSC->IsPeriodic()) ;
1499             
1500             
1501           }
1502         } 
1503         if(TheTypeC == GeomAbs_BezierCurve) {
1504           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1505           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1506           for(i = 1;i <= Curve->NbPoles();i++) {
1507             Dist2Min = IntegerLast();
1508             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1509                                             0.5,TolU,TolV) ;
1510             if (extrloc.IsDone()) {
1511               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1512               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1513               //if (Dist2Min < myTolerance * myTolerance) {
1514                 (extrloc.Point()).Parameter(u,v);
1515                 Poles2d(i).SetCoord(u,v);
1516                 myProjIsDone = Standard_True;
1517               }
1518               else break;
1519             }
1520             else break;
1521             if(myProjIsDone) 
1522               myProjIsDone = Standard_False;
1523             else break;
1524           }
1525           if(myProjIsDone) {
1526             TColStd_Array1OfReal Knots(1, 2);
1527             Knots.SetValue(1,0.0);
1528             Knots.SetValue(2,1.0);
1529             TColStd_Array1OfInteger Mults(1, 2);
1530             Mults.Init(BC->NbPoles());
1531             if(BC->IsRational()) {
1532               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1533               BC->Weights(Weights); 
1534               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1535                                                     BC->Degree(), BC->IsPeriodic()) ;
1536             }
1537             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1538                                                   BC->Degree(), BC->IsPeriodic()) ;
1539           }
1540         } 
1541       }
1542     }
1543   }
1544   else if(TheTypeS == GeomAbs_BezierSurface) {
1545     Handle(Geom_BezierSurface) BS = Surf->Bezier();
1546     if((BS->MaxDegree() == 1) &&
1547        (BS->NbUPoles() == 2) &&
1548        (BS->NbVPoles() == 2)) {
1549       gp_Pnt p11 = BS->Pole(1,1);
1550       gp_Pnt p12 = BS->Pole(1,2);
1551       gp_Pnt p21 = BS->Pole(2,1);
1552       gp_Pnt p22 = BS->Pole(2,2);
1553       gp_Vec V1(p11,p12);
1554       gp_Vec V2(p21,p22);
1555       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ //OCC217
1556       //if (V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){ 
1557         //    and if it is enough to project the poles of Curve
1558         Standard_Integer Dist2Min = IntegerLast();
1559         Standard_Real u,v;
1560         //OCC217
1561         //Standard_Real TolU = Surf->UResolution(myTolerance)
1562         //  , TolV = Surf->VResolution(myTolerance);
1563
1564 //      gp_Pnt pntproj;
1565         if(TheTypeC == GeomAbs_BSplineCurve) {
1566           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1567           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1568           for(i = 1;i <= Curve->NbPoles();i++) {
1569             myProjIsDone = Standard_False;
1570             Dist2Min = IntegerLast();
1571             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1572                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1573             if (extrloc.IsDone()) {
1574               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1575               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1576               //if (Dist2Min < myTolerance * myTolerance) {
1577                 (extrloc.Point()).Parameter(u,v);
1578                 Poles2d(i).SetCoord(u,v);
1579                 myProjIsDone = Standard_True;
1580               }
1581               else break;
1582             }
1583             else break;
1584             if(!myProjIsDone) 
1585               break;
1586           }
1587           if(myProjIsDone) {
1588             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1589             BSC->Knots(Knots);
1590             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1591             BSC->Multiplicities(Mults);
1592             if(BSC->IsRational()) {
1593               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1594               BSC->Weights(Weights); 
1595               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1596                                                     BSC->Degree(), BSC->IsPeriodic()) ;
1597             }
1598             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1599                                                   BSC->Degree(), BSC->IsPeriodic()) ;
1600             
1601             
1602           }
1603         } 
1604         if(TheTypeC == GeomAbs_BezierCurve) {
1605           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1606           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1607           for(i = 1;i <= Curve->NbPoles();i++) {
1608             Dist2Min = IntegerLast();
1609             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1610                                             0.5,TolU,TolV) ;
1611             if (extrloc.IsDone()) {
1612               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1613               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1614               //if (Dist2Min < myTolerance * myTolerance) {
1615                 (extrloc.Point()).Parameter(u,v);
1616                 Poles2d(i).SetCoord(u,v);
1617                 myProjIsDone = Standard_True;
1618               }
1619               else break;
1620             }
1621             else break;
1622             if(myProjIsDone) 
1623               myProjIsDone = Standard_False;
1624             else break;
1625           }
1626           if(myProjIsDone) {
1627             TColStd_Array1OfReal Knots(1, 2);
1628             Knots.SetValue(1,0.0);
1629             Knots.SetValue(2,1.0);
1630             TColStd_Array1OfInteger Mults(1, 2);
1631             Mults.Init(BC->NbPoles());
1632             if(BC->IsRational()) {
1633               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1634               BC->Weights(Weights); 
1635               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1636                                                     BC->Degree(), BC->IsPeriodic()) ;
1637             }
1638             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1639                                                   BC->Degree(), BC->IsPeriodic()) ;
1640           }
1641         } 
1642       }
1643     }
1644   }
1645
1646   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;  //OCC217
1647   //ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, myTolerance) ;
1648
1649 #ifdef OCCT_DEBUG
1650   Standard_Integer Nb = 50;
1651   
1652   Standard_Real U, U1, U2;
1653   U1 = F.FirstParameter();
1654   U2 = F.LastParameter();
1655   
1656   TColgp_Array1OfPnt2d    DummyPoles(1,Nb+1);
1657   TColStd_Array1OfReal    DummyKnots(1,Nb+1);
1658   TColStd_Array1OfInteger DummyMults(1,Nb+1);
1659   DummyMults.Init(1);
1660   DummyMults(1) = 2;
1661   DummyMults(Nb+1) = 2;
1662   for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1663     U = (Nb-ij)*U1 + ij*U2;
1664     U /= Nb;
1665     DummyPoles(ij+1) = F.Value(U);
1666     DummyKnots(ij+1) = ij;
1667   }
1668   Handle(Geom2d_BSplineCurve) DummyC2d =
1669     new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1670 #ifdef DRAW
1671   Standard_CString Temp = "bs2d";
1672   DrawTrSurf::Set(Temp,DummyC2d);
1673 #endif
1674 //  DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1675   Handle(Geom2dAdaptor_HCurve) DDD = 
1676     Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1677   
1678 #ifdef DRAW
1679   Temp = "initc2d";
1680   DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1681 #endif
1682 //  DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1683 #endif
1684
1685   Standard_Integer Deg1,Deg2;
1686 //  Deg1 = 8;
1687 //  Deg2 = 8;
1688   Deg1 = 2; //IFV
1689   Deg2 = 8; //IFV 
1690
1691   Approx_FitAndDivide2d Fit(F,Deg1,Deg2,Tol3d,Tol2d, //OCC217
1692   //Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
1693                             Standard_True);
1694
1695   if(Fit.IsAllApproximated()) {
1696     Standard_Integer i;
1697     Standard_Integer NbCurves = Fit.NbMultiCurves();
1698     Standard_Integer MaxDeg = 0;
1699     // To transform the MultiCurve into BSpline, it is required that all  
1700     // Bezier constituing it have the same degree -> Calculation of MaxDeg
1701     Standard_Integer NbPoles  = 1;
1702     for (i = 1; i <= NbCurves; i++) {
1703       Standard_Integer Deg = Fit.Value(i).Degree();
1704       MaxDeg = Max ( MaxDeg, Deg);
1705     }
1706
1707     NbPoles = MaxDeg * NbCurves + 1;               //Tops on the BSpline
1708     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
1709       
1710     TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1711     
1712     TColStd_Array1OfReal Knots( 1, NbCurves + 1);  //Nodes of the BSpline
1713     
1714     Standard_Integer Compt = 1;
1715     for (i = 1; i <= NbCurves; i++) {
1716       Fit.Parameters(i, Knots(i), Knots(i+1)); 
1717       AppParCurves_MultiCurve MC = Fit.Value( i);       //Load the Ith Curve
1718       TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1719       MC.Curve(1, Poles2d);
1720       
1721       //Eventual augmentation of the degree
1722       Standard_Integer Inc = MaxDeg - MC.Degree();
1723       if ( Inc > 0) {
1724 //      BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(), 
1725         BSplCLib::IncreaseDegree( MaxDeg, Poles2d, PLib::NoWeights(), 
1726                          TempPoles, PLib::NoWeights());
1727         //update of tops of the PCurve
1728         for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
1729           Poles.SetValue( Compt, TempPoles( j));
1730           Compt++;
1731         }
1732       }
1733       else {
1734         //update of tops of the PCurve
1735         for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
1736           Poles.SetValue( Compt, Poles2d( j));
1737           Compt++;
1738         }
1739       } 
1740       
1741       Compt--;
1742     }
1743     
1744     //update of fields of  ProjLib_Approx
1745     Standard_Integer NbKnots = NbCurves + 1;
1746
1747     TColStd_Array1OfInteger   Mults( 1, NbKnots);
1748     Mults.Init(MaxDeg);
1749     Mults.SetValue( 1, MaxDeg + 1);
1750     Mults.SetValue(NbKnots, MaxDeg + 1);
1751     myProjIsDone = Standard_True;
1752     Handle(Geom2d_BSplineCurve) Dummy =
1753       new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1754     
1755     // try to smoother the Curve GeomAbs_C1.
1756
1757     Standard_Boolean OK = Standard_True;
1758
1759     for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1760       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,Tol3d);  //OCC217
1761       //OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,myTolerance);
1762     }
1763 #ifdef OCCT_DEBUG
1764     if (!OK) {
1765       cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<endl;
1766     }
1767 #endif 
1768     return Dummy;
1769   }
1770   return Handle(Geom2d_BSplineCurve)();
1771 }
1772
1773 //=======================================================================
1774 //function : BSpline
1775 //purpose  : 
1776 //=======================================================================
1777
1778 Handle(Geom2d_BSplineCurve)  
1779      ProjLib_ComputeApproxOnPolarSurface::BSpline() const 
1780      
1781 {
1782 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:46 2002 End
1783 //   Standard_NoSuchObject_Raise_if
1784 //     (!myProjIsDone,
1785 //      "ProjLib_ComputeApproxOnPolarSurface:BSpline");
1786 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:47 2002 End
1787   return myBSpline ;
1788 }
1789
1790 //=======================================================================
1791 //function : Curve2d
1792 //purpose  : 
1793 //=======================================================================
1794
1795 Handle(Geom2d_Curve)  
1796      ProjLib_ComputeApproxOnPolarSurface::Curve2d() const 
1797      
1798 {
1799   Standard_NoSuchObject_Raise_if
1800     (!myProjIsDone,
1801      "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
1802   return my2ndCurve ;
1803 }
1804
1805
1806 //=======================================================================
1807 //function : IsDone
1808 //purpose  : 
1809 //=======================================================================
1810
1811 Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const 
1812      
1813 {
1814   return myProjIsDone;
1815 }