0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[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 #include <Geom2d_Curve.hxx>
81 #endif
82 //static Standard_Integer compteur = 0;
83 #endif
84
85 struct aFuncStruct
86 {
87   aFuncStruct() // Empty constructor.
88   : mySqProjOrtTol(0.0),
89     myTolU(0.0),
90     myTolV(0.0)
91   {
92     memset(myPeriod, 0, sizeof (myPeriod));
93   }
94
95   Handle(Adaptor3d_HSurface) mySurf; // Surface where to project.
96   Handle(Adaptor3d_HCurve)   myCurve; // Curve to project.
97   Handle(Adaptor2d_HCurve2d) myInitCurve2d; // Initial 2dcurve projection.
98   Standard_Real mySqProjOrtTol; // Used to filter non-orthogonal projected point.
99   Standard_Real myTolU;
100   Standard_Real myTolV;
101   Standard_Real myPeriod[2]; // U and V period correspondingly.
102 };
103
104 //=======================================================================
105 //function : computePeriodicity
106 //purpose  : Compute period information on adaptor.
107 //=======================================================================
108 static void computePeriodicity(const Handle(Adaptor3d_HSurface)& theSurf,
109                                Standard_Real &theUPeriod,
110                                Standard_Real &theVPeriod)
111 {
112   theUPeriod = 0.0;
113   theVPeriod = 0.0;
114
115   // Compute once information about periodicity.
116   // Param space may be reduced in case of rectangular trimmed surface,
117   // in this case really trimmed bounds should be set as unperiodic.
118   Standard_Real aTrimF, aTrimL, aBaseF, aBaseL, aDummyF, aDummyL;
119   Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theSurf->Surface(), Standard_False); // Not trim.
120   // U param space.
121   if (theSurf->IsUPeriodic())
122   {
123     theUPeriod = theSurf->UPeriod();
124   }
125   else if(theSurf->IsUClosed())
126   {
127     theUPeriod = theSurf->LastUParameter() - theSurf->FirstUParameter();
128   }
129   if (theUPeriod != 0.0)
130   {
131     aTrimF = theSurf->FirstUParameter(); // Trimmed first
132     aTrimL = theSurf->LastUParameter(); // Trimmed last
133     aS->Bounds(aBaseF, aBaseL, aDummyF, aDummyL); // Non-trimmed values.
134     if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
135     {
136       // Param space reduced.
137       theUPeriod = 0.0;
138     }
139   }
140
141   // V param space.
142   if (theSurf->IsVPeriodic())
143   {
144     theVPeriod = theSurf->VPeriod();
145   }
146   else if(theSurf->IsVClosed())
147   {
148     theVPeriod = theSurf->LastVParameter() - theSurf->FirstVParameter();
149   }
150   if (theVPeriod != 0.0)
151   {
152     aTrimF = theSurf->FirstVParameter(); // Trimmed first
153     aTrimL = theSurf->LastVParameter(); // Trimmed last
154     aS->Bounds(aDummyF, aDummyL, aBaseF, aBaseL); // Non-trimmed values.
155     if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
156     {
157       // Param space reduced.
158       theVPeriod = 0.0;
159     }
160   }
161 }
162
163 //=======================================================================
164 //function : aFuncValue
165 //purpose  : compute functional value in (theU,theV) point
166 //=======================================================================
167 static Standard_Real anOrthogSqValue(const gp_Pnt& aBasePnt,
168                                      const Handle(Adaptor3d_HSurface)& Surf,
169                                      const Standard_Real theU,
170                                      const Standard_Real theV)
171 {
172   // Since find projection, formula is:
173   // F1 = Dot(S_U, Vec(aBasePnt, aProjPnt))
174   // F2 = Dot(S_V, Vec(aBasePnt, aProjPnt))
175
176   gp_Pnt aProjPnt;
177   gp_Vec aSu, aSv;
178
179   Surf->D1(theU, theV, aProjPnt, aSu, aSv);
180   gp_Vec aBaseVec(aBasePnt, aProjPnt);
181
182   if (aSu.SquareMagnitude() > Precision::SquareConfusion())
183     aSu.Normalize();
184
185   if (aSv.SquareMagnitude() > Precision::SquareConfusion())
186     aSv.Normalize();
187
188   Standard_Real aFirstPart = aSu.Dot(aBaseVec);
189   Standard_Real aSecondPart = aSv.Dot(aBaseVec);
190   return (aFirstPart * aFirstPart + aSecondPart * aSecondPart);
191 }
192
193 //=======================================================================
194 //function : Value
195 //purpose  : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
196 //            <InitCurve2d> use for calculate start 2D point.
197 //=======================================================================
198 static gp_Pnt2d Function_Value(const Standard_Real theU,
199                                const aFuncStruct& theData)
200 {
201   gp_Pnt2d p2d = theData.myInitCurve2d->Value(theU) ;
202   gp_Pnt p = theData.myCurve->Value(theU);
203   gp_Pnt aSurfPnt = theData.mySurf->Value(p2d.X(), p2d.Y());
204   Standard_Real aSurfPntDist = aSurfPnt.SquareDistance(p);
205
206   Standard_Real Uinf, Usup, Vinf, Vsup;
207   Uinf = theData.mySurf->Surface().FirstUParameter();
208   Usup = theData.mySurf->Surface().LastUParameter();
209   Vinf = theData.mySurf->Surface().FirstVParameter();
210   Vsup = theData.mySurf->Surface().LastVParameter();
211
212   // Check case when curve is close to co-parametrized isoline on surf.
213   if (Abs (p2d.X() - Uinf) < Precision::PConfusion() ||
214       Abs (p2d.X() - Usup) < Precision::PConfusion() )
215   {
216     // V isoline.
217     gp_Pnt aPnt;
218     theData.mySurf->D0(p2d.X(), theU, aPnt);
219     if (aPnt.SquareDistance(p) < aSurfPntDist)
220       p2d.SetY(theU);
221   }
222
223   if (Abs (p2d.Y() - Vinf) < Precision::PConfusion() ||
224       Abs (p2d.Y() - Vsup) < Precision::PConfusion() )
225   {
226     // U isoline.
227     gp_Pnt aPnt;
228     theData.mySurf->D0(theU, p2d.Y(), aPnt);
229     if (aPnt.SquareDistance(p) < aSurfPntDist)
230       p2d.SetX(theU);
231   }
232
233   Standard_Integer decalU = 0, decalV = 0;
234   Standard_Real U0 = p2d.X(), V0 = p2d.Y();
235
236   GeomAbs_SurfaceType Type = theData.mySurf->GetType();
237   if((Type != GeomAbs_BSplineSurface) && 
238      (Type != GeomAbs_BezierSurface)  &&
239      (Type != GeomAbs_OffsetSurface)    )
240   {
241     // Analytical cases.
242     Standard_Real S = 0., T = 0.;
243     switch (Type)
244     {
245     case GeomAbs_Cylinder:
246       {
247         gp_Cylinder Cylinder = theData.mySurf->Cylinder();
248         ElSLib::Parameters( Cylinder, p, S, T);
249         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
250         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
251         S += decalU*2*M_PI;
252         break;
253       }
254     case GeomAbs_Cone:
255       {
256         gp_Cone Cone = theData.mySurf->Cone();
257         ElSLib::Parameters( Cone, p, S, T);
258         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
259         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
260         S += decalU*2*M_PI;
261         break;
262       }
263     case GeomAbs_Sphere:
264       {
265         gp_Sphere Sphere = theData.mySurf->Sphere();
266         ElSLib::Parameters( Sphere, p, S, T);
267         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
268         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
269         S += decalU*2*M_PI;
270         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
271         if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
272         T += decalV*2*M_PI;
273         if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI)
274         {
275           T = M_PI - T;
276           if(U0 < S)
277             S -= M_PI;
278           else
279             S += M_PI;
280         }
281         break;
282       }
283     case GeomAbs_Torus:
284       {
285         gp_Torus Torus = theData.mySurf->Torus();
286         ElSLib::Parameters( Torus, p, S, T);
287         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
288         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
289         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
290         if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
291         S += decalU*2*M_PI; T += decalV*2*M_PI;
292         break;
293       }
294     default:
295       throw Standard_NoSuchObject("ProjLib_ComputeApproxOnPolarSurface::Value");
296     }
297     return gp_Pnt2d(S, T);
298   }
299
300   // Non-analytical case.
301   Standard_Real Dist2Min = RealLast();
302   Standard_Real uperiod = theData.myPeriod[0],
303                 vperiod = theData.myPeriod[1],
304                 u, v;
305
306   // U0 and V0 are the points within the initialized period.
307   if(U0 < Uinf)
308   {
309     if(!uperiod)
310       U0 = Uinf;
311     else
312     {
313       decalU = int((Uinf - U0)/uperiod)+1;
314       U0 += decalU*uperiod;
315     }
316   }
317   if(U0 > Usup)
318   {
319     if(!uperiod)
320       U0 = Usup;
321     else
322     {
323       decalU = -(int((U0 - Usup)/uperiod)+1);
324       U0 += decalU*uperiod;
325     }
326   }
327   if(V0 < Vinf)
328   {
329     if(!vperiod)
330       V0 = Vinf;
331     else
332     {
333       decalV = int((Vinf - V0)/vperiod)+1;
334       V0 += decalV*vperiod;
335     }
336   }
337   if(V0 > Vsup)
338   {
339     if(!vperiod)
340       V0 = Vsup;
341     else
342     {
343       decalV = -int((V0 - Vsup)/vperiod)-1;
344       V0 += decalV*vperiod;
345     }
346   }
347
348   // The surface around (U0,V0) is reduced.
349   Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
350   Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
351   if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
352   if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
353   if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
354   if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
355
356   GeomAdaptor_Surface SurfLittle;
357   if (Type == GeomAbs_BSplineSurface)
358   {
359     Handle(Geom_Surface) GBSS(theData.mySurf->Surface().BSpline());
360     SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
361   }
362   else if (Type == GeomAbs_BezierSurface)
363   {
364     Handle(Geom_Surface) GS(theData.mySurf->Surface().Bezier());
365     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
366   }
367   else if (Type == GeomAbs_OffsetSurface)
368   {
369     Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(theData.mySurf->Surface());
370     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
371   }
372   else
373   {
374     throw Standard_NoSuchObject ("ProjLib_ComputeApproxOnPolarSurface::ProjectUsingInitialCurve2d() - unknown surface type");
375   }
376
377   // Try to run simple search with initial point (U0, V0).
378   Extrema_GenLocateExtPS  locext(SurfLittle, theData.myTolU, theData.myTolV);
379   locext.Perform(p, U0, V0);
380   if (locext.IsDone()) 
381   {
382     locext.Point().Parameter(u, v);
383     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
384     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
385         locext.SquareDistance() < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
386     {
387       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
388       return pnt;
389     }
390   }
391
392   // Perform whole param space search.
393   Extrema_ExtPS  ext(p, SurfLittle, theData.myTolU, theData.myTolV);
394   if (ext.IsDone() && ext.NbExt() >= 1)
395   {
396     Dist2Min = ext.SquareDistance(1);
397     Standard_Integer GoodValue = 1;
398     for (Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
399     {
400       if( Dist2Min > ext.SquareDistance(i))
401       {
402         Dist2Min = ext.SquareDistance(i);
403         GoodValue = i;
404       }
405     }
406     ext.Point(GoodValue).Parameter(u, v);
407     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
408     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
409         ext.SquareDistance(GoodValue) < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
410     {
411       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
412       return pnt;
413     }
414   }
415
416   // Both searches return bad values, use point from initial 2dcurve.
417   return p2d;
418 }
419
420
421 //=======================================================================
422 //function : ProjLib_PolarFunction
423 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
424 //=======================================================================
425
426 class ProjLib_PolarFunction : public AppCont_Function
427 {
428   aFuncStruct myStruct;
429
430   public :
431
432   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C, 
433                         const Handle(Adaptor3d_HSurface)& Surf,
434                         const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
435                         const Standard_Real Tol3d)
436   {
437     myNbPnt = 0;
438     myNbPnt2d = 1;
439
440     computePeriodicity(Surf, myStruct.myPeriod[0], myStruct.myPeriod[1]);
441
442     myStruct.myCurve = C;
443     myStruct.myInitCurve2d = InitialCurve2d;
444     myStruct.mySurf = Surf;
445     myStruct.mySqProjOrtTol = 10000.0 * Tol3d * Tol3d;
446     myStruct.myTolU = Surf->UResolution(Tol3d);
447     myStruct.myTolV = Surf->VResolution(Tol3d);
448   }
449
450   ~ProjLib_PolarFunction() {}
451
452   Standard_Real FirstParameter() const
453   {
454     return myStruct.myCurve->FirstParameter();
455   }
456
457   Standard_Real LastParameter() const
458   {
459     return myStruct.myCurve->LastParameter();
460   }
461
462   gp_Pnt2d Value(const Standard_Real t) const
463   {
464     return Function_Value(t, myStruct);
465   }
466
467   Standard_Boolean Value(const Standard_Real   theT,
468                          NCollection_Array1<gp_Pnt2d>& thePnt2d,
469                          NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
470   {
471     thePnt2d(1) = Function_Value(theT, myStruct);
472     return Standard_True;
473   }
474
475   Standard_Boolean D1(const Standard_Real   /*theT*/,
476                       NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
477                       NCollection_Array1<gp_Vec>&   /*theVec*/) const
478     {return Standard_False;}
479 };
480
481 //=======================================================================
482 //function : ProjLib_ComputeApproxOnPolarSurface
483 //purpose  : 
484 //=======================================================================
485
486 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
487 : myProjIsDone(Standard_False),
488   myTolerance(Precision::Approximation()),
489   myTolReached(-1.0),
490   myDegMin(-1), myDegMax(-1),
491   myMaxSegments(-1),
492   myMaxDist(-1.),
493   myBndPnt(AppParCurves_TangencyPoint),
494   myDist(0.)
495 {
496 }
497
498 //=======================================================================
499 //function : ProjLib_ComputeApproxOnPolarSurface
500 //purpose  : 
501 //=======================================================================
502
503 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
504                     (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
505                      const Handle(Adaptor3d_HCurve)&   theCurve,
506                      const Handle(Adaptor3d_HSurface)& theSurface,
507                      const Standard_Real               theTolerance3D)
508 : myProjIsDone(Standard_False),
509   myTolerance(theTolerance3D),
510   myTolReached(-1.0),
511   myDegMin(-1), myDegMax(-1),
512   myMaxSegments(-1),
513   myMaxDist(-1.),
514   myBndPnt(AppParCurves_TangencyPoint),
515   myDist(0.)
516 {
517   myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
518 }
519
520 //=======================================================================
521 //function : ProjLib_ComputeApproxOnPolarSurface
522 //purpose  : case without curve of initialization
523 //=======================================================================
524
525 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
526                       (const Handle(Adaptor3d_HCurve)&   theCurve,
527                        const Handle(Adaptor3d_HSurface)& theSurface,
528                        const Standard_Real               theTolerance3D)
529 : myProjIsDone(Standard_False),
530   myTolerance(theTolerance3D),
531   myTolReached(-1.0),
532   myDegMin(-1), myDegMax(-1),
533   myMaxSegments(-1),
534   myMaxDist(-1.),
535   myBndPnt(AppParCurves_TangencyPoint),
536   myDist(0.)
537 {
538   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
539   myBSpline = Perform(anInitCurve2d, theCurve, theSurface);  
540
541
542 //=======================================================================
543 //function : ProjLib_ComputeApproxOnPolarSurface
544 //purpose  : Process the case of sewing
545 //=======================================================================
546
547 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
548                 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
549                  const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
550                  const Handle(Adaptor3d_HCurve)&   theCurve,
551                  const Handle(Adaptor3d_HSurface)& theSurface,
552                  const Standard_Real               theTolerance3D)
553 : myProjIsDone(Standard_False),
554   myTolerance(theTolerance3D),
555   myTolReached(-1.0),
556   myDegMin(-1), myDegMax(-1),
557   myMaxSegments(-1),
558   myMaxDist(-1.),
559   myBndPnt(AppParCurves_TangencyPoint),
560   myDist(0.)
561 {
562   // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing 
563   Handle(Geom2d_BSplineCurve) bsc =
564     Perform(theInitialCurve2d, theCurve, theSurface);
565
566   if(myProjIsDone) {
567     gp_Pnt2d P2dproj, P2d, P2dBis;
568     P2dproj = bsc->StartPoint();
569     P2d    = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
570     P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
571
572     Standard_Real Dist, DistBis;
573     Dist    = P2dproj.Distance(P2d);
574     DistBis = P2dproj.Distance(P2dBis);
575     if( Dist < DistBis)  {
576       // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
577       myBSpline = bsc;
578       Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
579       my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
580     }
581     else {
582       my2ndCurve = bsc;
583       Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
584       myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
585     }
586   }
587 }
588
589
590 //=======================================================================
591 //function : Concat
592 //purpose  : 
593 //=======================================================================
594
595 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
596                                           Handle(Geom2d_BSplineCurve) C2,
597                                           Standard_Real theUJump,
598                                           Standard_Real theVJump)
599 {
600   Standard_Integer deg, deg1, deg2;
601   deg1 = C1->Degree();
602   deg2 = C2->Degree();
603   
604   if ( deg1 < deg2) {
605     C1->IncreaseDegree(deg2);
606     deg = deg2;
607   }
608   else if ( deg2 < deg1) {
609     C2->IncreaseDegree(deg1);
610     deg = deg1;
611   }
612   else deg = deg1;
613
614   Standard_Integer np1,np2,nk1,nk2,np,nk;
615   np1 = C1->NbPoles();
616   nk1 = C1->NbKnots();
617   np2 = C2->NbPoles();
618   nk2 = C2->NbKnots();
619   nk = nk1 + nk2 -1;
620   np = np1 + np2 -1;
621
622   TColStd_Array1OfReal    K1(1,nk1); C1->Knots(K1);
623   TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
624   TColgp_Array1OfPnt2d    P1(1,np1); C1->Poles(P1);
625   TColStd_Array1OfReal    K2(1,nk2); C2->Knots(K2);
626   TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
627   TColgp_Array1OfPnt2d    P2(1,np2); C2->Poles(P2);
628
629   // Compute the new BSplineCurve
630   TColStd_Array1OfReal    K(1,nk);
631   TColStd_Array1OfInteger M(1,nk);
632   TColgp_Array1OfPnt2d    P(1,np);
633
634   Standard_Integer i, count = 0;
635   // Set Knots and Mults
636   for ( i = 1; i <= nk1; i++) {
637     count++;
638     K(count) = K1(i);
639     M(count) = M1(i);
640   }
641   M(count) = deg;
642   for ( i = 2; i <= nk2; i++) {
643     count++;
644     K(count) = K2(i);
645     M(count) = M2(i);
646   }
647   // Set the Poles
648   count = 0;
649   for (i = 1; i <= np1; i++) {
650     count++;
651     P(count) = P1(i);
652   }
653   for (i = 2; i <= np2; i++) {
654     count++;
655     P(count).SetX(P2(i).X() + theUJump);
656     P(count).SetY(P2(i).Y() + theVJump);
657   }
658
659   Handle(Geom2d_BSplineCurve) BS = 
660     new Geom2d_BSplineCurve(P,K,M,deg);
661   return BS;
662 }
663
664 //=======================================================================
665 //function : Perform
666 //purpose  : 
667 //=======================================================================
668
669 void ProjLib_ComputeApproxOnPolarSurface::Perform
670 (const Handle(Adaptor3d_HCurve)& Curve, const Handle(Adaptor3d_HSurface)& S)
671 {
672   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
673   myBSpline = Perform(anInitCurve2d, Curve, S);  
674 }
675
676 //=======================================================================
677 //function : Perform
678 //purpose  : 
679 //=======================================================================
680
681 Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
682 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
683  const Handle(Adaptor3d_HCurve)& Curve,
684  const Handle(Adaptor3d_HSurface)& S)
685 {
686   //OCC217
687   Standard_Real Tol3d = myTolerance; 
688   Standard_Real ParamTol = Precision::PApproximation();
689
690   Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
691   Handle(Adaptor3d_HCurve) AHC = Curve;
692   
693 // if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
694 // -> bug cts18237
695   GeomAbs_CurveType typeCurve = Curve->GetType();
696   if(typeCurve == GeomAbs_BSplineCurve) {
697     TColStd_ListOfTransient LOfBSpline2d;
698     Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
699     Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
700     if(nbInter > 1) {
701       Standard_Integer i, j;
702       Handle(Geom_TrimmedCurve) GTC;
703       Handle(Geom2d_TrimmedCurve) G2dTC;
704       TColStd_Array1OfReal Inter(1,nbInter+1);
705       Curve->Intervals(Inter,GeomAbs_C1);
706       Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
707       // initialization 3d
708       GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
709       AHC = new GeomAdaptor_HCurve(GTC);
710       
711       // if there is an initialization curve: 
712       // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
713       // and the sections C1 are taken
714       // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
715       
716       // initialization 2d
717       Standard_Integer nbInter2d;
718       Standard_Boolean C2dIsToCompute;
719       C2dIsToCompute = InitialCurve2d.IsNull();
720       Handle(Geom2d_BSplineCurve) BSC2d;
721       Handle(Geom2d_Curve) G2dC;
722       
723       if(!C2dIsToCompute) {
724         nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
725         TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
726         InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
727         j = 1;
728         for(i = 1,j = 1;i <= nbInter;i++)
729           if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
730           //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
731             if (j > nbInter2d) break;
732             j++;
733           }
734         if(j != (nbInter2d+1)) {
735           C2dIsToCompute = Standard_True;
736         }
737       }
738       
739       if(C2dIsToCompute) {
740         AHC2d = BuildInitialCurve2d(AHC, S);
741       }
742       else {
743         typeCurve = InitialCurve2d->GetType();
744         switch (typeCurve) {
745         case GeomAbs_Line: {
746           G2dC = new Geom2d_Line(InitialCurve2d->Line());
747           break;
748         }
749         case GeomAbs_Circle: {
750           G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
751           break;
752         }
753         case GeomAbs_Ellipse: {
754           G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
755           break;
756         }
757         case GeomAbs_Hyperbola: {
758           G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
759           break;
760         }
761         case GeomAbs_Parabola: {
762           G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
763           break;
764         }
765         case GeomAbs_BezierCurve: {
766           G2dC = InitialCurve2d->Bezier();
767           break;
768         }
769         case GeomAbs_BSplineCurve: {
770           G2dC = InitialCurve2d->BSpline();
771           break;
772         }
773         case GeomAbs_OtherCurve:
774         default:
775           break;
776         }
777         gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
778         gp_Pnt fps, lps, fpc, lpc;
779         S->D0(fp2d.X(), fp2d.Y(), fps);
780         S->D0(lp2d.X(), lp2d.Y(), lps);
781         Curve->D0(firstinter, fpc);
782         Curve->D0(secondinter, lpc);
783         //OCC217
784         if((fps.IsEqual(fpc, Tol3d)) &&
785            (lps.IsEqual(lpc, Tol3d))) {
786         //if((fps.IsEqual(fpc, myTolerance)) &&
787         //   (lps.IsEqual(lpc, myTolerance))) {
788           G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
789           Geom2dAdaptor_Curve G2dAC(G2dTC);
790           AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
791           myProjIsDone = Standard_True;
792         }
793         else {
794           AHC2d = BuildInitialCurve2d(AHC, S);
795           C2dIsToCompute = Standard_True;
796         }
797       }
798         
799       if(myProjIsDone) {
800         BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
801         if(BSC2d.IsNull()) 
802         {
803             return Handle(Geom2d_BSplineCurve)();
804         }
805         LOfBSpline2d.Append(BSC2d);
806       }
807       else {
808         return Handle(Geom2d_BSplineCurve)();
809       }
810       
811
812
813       Standard_Real iinter, ip1inter;
814       Standard_Integer nbK2d, deg;
815       nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
816
817       for(i = 2;i <= nbInter;i++) {
818         iinter = Inter.Value(i);
819         ip1inter = Inter.Value(i+1);
820         // general case 3d
821         GTC->SetTrim(iinter, ip1inter);
822         AHC = new GeomAdaptor_HCurve(GTC);
823         
824         // general case 2d
825         if(C2dIsToCompute) {
826           AHC2d = BuildInitialCurve2d(AHC, S);
827         }
828         else {
829           gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
830           gp_Pnt fps, lps, fpc, lpc;
831           S->D0(fp2d.X(), fp2d.Y(), fps);
832           S->D0(lp2d.X(), lp2d.Y(), lps);
833           Curve->D0(iinter, fpc);
834           Curve->D0(ip1inter, lpc);
835           //OCC217
836           if((fps.IsEqual(fpc, Tol3d)) &&
837              (lps.IsEqual(lpc, Tol3d))) {
838           //if((fps.IsEqual(fpc, myTolerance)) &&
839           //   (lps.IsEqual(lpc, myTolerance))) {
840             G2dTC->SetTrim(iinter, ip1inter);
841             Geom2dAdaptor_Curve G2dAC(G2dTC);
842             AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
843             myProjIsDone = Standard_True;
844           }
845           else {
846             AHC2d = BuildInitialCurve2d(AHC, S);
847           }
848         }
849         if(myProjIsDone) {
850           BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
851           if(BSC2d.IsNull()) {
852             return Handle(Geom2d_BSplineCurve)();
853           }
854           LOfBSpline2d.Append(BSC2d);
855           nbK2d += BSC2d->NbKnots() - 1;
856           deg = Max(deg, BSC2d->Degree());
857         }
858         else {
859           return Handle(Geom2d_BSplineCurve)();
860         }
861       }
862
863       Standard_Real anUPeriod, anVPeriod;
864       computePeriodicity(S, anUPeriod, anVPeriod);
865       Standard_Integer NbC = LOfBSpline2d.Extent();
866       Handle(Geom2d_BSplineCurve) CurBS;
867       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
868       LOfBSpline2d.RemoveFirst();
869       for (Standard_Integer ii = 2; ii <= NbC; ii++)
870       {
871         Handle(Geom2d_BSplineCurve) BS = 
872           Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
873
874         //Check for period jump in point of contact.
875         gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
876         gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
877         Standard_Real anUJump = 0.0, anVJump = 0.0;
878
879         if (anUPeriod > 0.0 &&
880             Abs (aC1End.X() - aC2Beg.X()) > (anUPeriod ) / 2.01)
881         {
882           Standard_Real aMultCoeff =  aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
883           anUJump = (anUPeriod) * aMultCoeff;
884         }
885
886         if (anVPeriod &&
887             Abs (aC1End.Y() - aC2Beg.Y()) > (anVPeriod) / 2.01)
888         {
889           Standard_Real aMultCoeff =  aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
890           anVJump = (anVPeriod) * aMultCoeff;
891         }
892
893         CurBS = Concat(CurBS,BS, anUJump, anVJump);
894         LOfBSpline2d.RemoveFirst();
895       }
896       return CurBS;
897     }
898   }
899   
900   if(InitialCurve2d.IsNull()) {
901     AHC2d = BuildInitialCurve2d(Curve, S);
902     if(!myProjIsDone) 
903       return Handle(Geom2d_BSplineCurve)(); 
904   }
905   return ProjectUsingInitialCurve2d(AHC, S, AHC2d);     
906 }
907
908
909 //=======================================================================
910 //function : ProjLib_BuildInitialCurve2d
911 //purpose  : 
912 //=======================================================================
913
914 Handle(Adaptor2d_HCurve2d) 
915      ProjLib_ComputeApproxOnPolarSurface::
916      BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)&   Curve,
917                          const Handle(Adaptor3d_HSurface)& Surf)
918 {
919   //  discretize the Curve with quasiuniform deflection
920   //  density at least NbOfPnts points
921   myProjIsDone = Standard_False;
922   
923   //OCC217
924   Standard_Real Tol3d = myTolerance; 
925   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
926   Standard_Real DistTol3d = 100.0*Tol3d;
927   if(myMaxDist > 0.)
928   {
929     DistTol3d = myMaxDist;
930   }
931   Standard_Real DistTol3d2 = DistTol3d * DistTol3d;
932   Standard_Real uperiod = 0.0, vperiod = 0.0;
933   computePeriodicity(Surf, uperiod, vperiod);
934
935   // NO myTol is Tol2d !!!!
936   //Standard_Real TolU = myTolerance, TolV = myTolerance;
937   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
938
939   Standard_Integer NbOfPnts = 61; 
940   GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
941   NbOfPnts = QUA.NbPoints();
942   TColgp_Array1OfPnt Pts(1,NbOfPnts);
943   TColStd_Array1OfReal Param(1,NbOfPnts);
944   Standard_Integer i, j;
945   for( i = 1; i <= NbOfPnts ; i++ ) { 
946     Param(i) = QUA.Parameter(i);
947     Pts(i) = Curve->Value(Param(i));
948   }
949   
950   TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
951   TColStd_Array1OfInteger Mult(1,NbOfPnts);
952   Mult.Init(1);
953   Mult(1) = Mult(NbOfPnts) = 2;
954   
955   Standard_Real Uinf, Usup, Vinf, Vsup;
956   Uinf = Surf->Surface().FirstUParameter();
957   Usup = Surf->Surface().LastUParameter();
958   Vinf = Surf->Surface().FirstVParameter();
959   Vsup = Surf->Surface().LastVParameter();
960   GeomAbs_SurfaceType Type = Surf->GetType();
961   if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
962      (Type != GeomAbs_OffsetSurface)) {
963     Standard_Real S, T;
964 //    Standard_Integer usens = 0, vsens = 0; 
965     // to know the position relatively to the period
966     switch (Type) {
967 //    case GeomAbs_Plane:
968 //      {
969 //      gp_Pln Plane = Surf->Plane();
970 //      for ( i = 1 ; i <= NbOfPnts ; i++) { 
971 //        ElSLib::Parameters( Plane, Pts(i), S, T);
972 //        Pts2d(i).SetCoord(S,T);
973 //      }
974 //      myProjIsDone = Standard_True;
975 //      break;
976 //      }
977     case GeomAbs_Cylinder:
978       {
979 //      Standard_Real Sloc, Tloc;
980         Standard_Real Sloc;
981         Standard_Integer usens = 0;
982         gp_Cylinder Cylinder = Surf->Cylinder();
983         ElSLib::Parameters( Cylinder, Pts(1), S, T);
984         Pts2d(1).SetCoord(S,T);
985         for ( i = 2 ; i <= NbOfPnts ; i++) { 
986           Sloc = S;
987           ElSLib::Parameters( Cylinder, Pts(i), S, T);
988           if(Abs(Sloc - S) > M_PI) {
989             if(Sloc > S)
990               usens++;
991             else
992               usens--;
993           }
994           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
995         }
996         myProjIsDone = Standard_True;
997         break;
998       }
999     case GeomAbs_Cone:
1000       {
1001 //      Standard_Real Sloc, Tloc;
1002         Standard_Real Sloc;
1003         Standard_Integer usens = 0;
1004         gp_Cone Cone = Surf->Cone();
1005         ElSLib::Parameters( Cone, Pts(1), S, T);
1006         Pts2d(1).SetCoord(S,T);
1007         for ( i = 2 ; i <= NbOfPnts ; i++) { 
1008           Sloc = S;
1009           ElSLib::Parameters( Cone, Pts(i), S, T);
1010           if(Abs(Sloc - S) > M_PI) {
1011             if(Sloc > S)
1012               usens++;
1013             else
1014               usens--;
1015           }
1016           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
1017         }
1018         myProjIsDone = Standard_True;
1019         break;
1020       }
1021     case GeomAbs_Sphere:
1022       {
1023         Standard_Real Sloc, Tloc;
1024         Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
1025         Standard_Boolean vparit = Standard_False;
1026         gp_Sphere Sphere = Surf->Sphere();
1027         ElSLib::Parameters( Sphere, Pts(1), S, T);
1028         Pts2d(1).SetCoord(S,T);
1029         for ( i = 2 ; i <= NbOfPnts ; i++) { 
1030           Sloc = S;Tloc = T;
1031           ElSLib::Parameters( Sphere, Pts(i), S, T);
1032           if(1.6*M_PI < Abs(Sloc - S)) {
1033             if(Sloc > S)
1034               usens += 2;
1035             else
1036               usens -= 2;
1037     }
1038           if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
1039             vparit = !vparit;
1040             if(Sloc > S)
1041               usens++;
1042             else
1043               usens--;
1044             if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
1045               vsens++;
1046             else
1047               vsens--;
1048           }
1049           if(vparit) {
1050             Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
1051           }       
1052           else {
1053             Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
1054             
1055           }
1056         }
1057         myProjIsDone = Standard_True;
1058         break;
1059       }
1060     case GeomAbs_Torus:
1061       {
1062         Standard_Real Sloc, Tloc;
1063         Standard_Integer usens = 0, vsens = 0;
1064         gp_Torus Torus = Surf->Torus();
1065         ElSLib::Parameters( Torus, Pts(1), S, T);
1066         Pts2d(1).SetCoord(S,T);
1067         for ( i = 2 ; i <= NbOfPnts ; i++) { 
1068           Sloc = S; Tloc = T;
1069           ElSLib::Parameters( Torus, Pts(i), S, T);
1070           if(Abs(Sloc - S) > M_PI) {
1071             if(Sloc > S)
1072               usens++;
1073             else
1074               usens--;
1075     }
1076           if(Abs(Tloc - T) > M_PI) {
1077             if(Tloc > T)
1078               vsens++;
1079             else
1080               vsens--;
1081     }
1082           Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
1083         }
1084         myProjIsDone = Standard_True;
1085         break;
1086       }
1087     default:
1088       throw Standard_NoSuchObject("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
1089     }
1090   }
1091   else {
1092     myProjIsDone = Standard_False;
1093     Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
1094     myDist = 0.;
1095     gp_Pnt pntproj;
1096
1097     TColgp_SequenceOfPnt2d Sols;
1098     Standard_Boolean areManyZeros = Standard_False;
1099     
1100     Curve->D0(Param.Value(1), pntproj) ;
1101     Extrema_ExtPS  aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
1102     Standard_Real aMinSqDist = RealLast();
1103     if (aExtPS.IsDone())
1104     {
1105       for (i = 1; i <= aExtPS.NbExt(); i++)
1106       {
1107         Standard_Real aSqDist = aExtPS.SquareDistance(i);
1108         if (aSqDist < aMinSqDist)
1109           aMinSqDist = aSqDist;
1110       }
1111     }
1112     if (aMinSqDist > DistTol3d2) //try to project with less tolerance
1113     {
1114       TolU = Min(TolU, Precision::PConfusion());
1115       TolV = Min(TolV, Precision::PConfusion());
1116       aExtPS.Initialize(Surf->Surface(),
1117                         Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(), 
1118                         Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
1119                         TolU, TolV);
1120       aExtPS.Perform(pntproj);
1121     }
1122
1123     if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
1124
1125       Standard_Integer GoodValue = 1;
1126
1127       for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
1128         if( aExtPS.SquareDistance(i) < DistTol3d2 ) {
1129           if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
1130             aExtPS.Point(i).Parameter(u,v);
1131             gp_Pnt2d p2d(u,v);
1132             Standard_Boolean isSame = Standard_False;
1133             for( j = 1; j <= Sols.Length(); j++ ) {
1134               if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
1135                 isSame = Standard_True;
1136                 break;
1137               }
1138             }
1139             if( !isSame ) Sols.Append( p2d );
1140           }
1141           if( Dist2Min > aExtPS.SquareDistance(i) ) {
1142             Dist2Min = aExtPS.SquareDistance(i);
1143             GoodValue = i;
1144           }
1145         }
1146       }
1147
1148       if( Sols.Length() > 1 ) areManyZeros = Standard_True;
1149
1150       if( Dist2Min <= DistTol3d2) {
1151         if( !areManyZeros ) {
1152           aExtPS.Point(GoodValue).Parameter(u,v);
1153           Pts2d(1).SetCoord(u,v);
1154           myProjIsDone = Standard_True;
1155         }
1156         else {
1157           Standard_Integer nbSols = Sols.Length();
1158           Standard_Real Dist2Max = -1.e+200;
1159           for( i = 1; i <= nbSols; i++ ) {
1160             const gp_Pnt2d& aP1 = Sols.Value(i);
1161             for( j = i+1; j <= nbSols; j++ ) {
1162               const gp_Pnt2d& aP2 = Sols.Value(j);
1163               Standard_Real aDist2 = aP1.SquareDistance(aP2);
1164               if( aDist2 > Dist2Max ) Dist2Max = aDist2;
1165             }
1166           }
1167       Standard_Real aMaxT2 = Max(TolU,TolV);
1168       aMaxT2 *= aMaxT2;
1169           if( Dist2Max > aMaxT2 ) {
1170             Standard_Integer tPp = 0;
1171             for( i = 1; i <= 5; i++ ) {
1172               Standard_Integer nbExtOk = 0;
1173               Standard_Integer indExt = 0;
1174               Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
1175               Curve->D0( Param.Value(iT), pntproj );
1176               Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1177               Dist2Min = 1.e+200;
1178               if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1179                 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
1180                   if( aTPS.SquareDistance(j) < DistTol3d2 ) {
1181                     nbExtOk++;
1182                     if( aTPS.SquareDistance(j) < Dist2Min ) {
1183                       Dist2Min = aTPS.SquareDistance(j);
1184                       indExt = j;
1185                     }
1186                   }
1187                 }
1188               }
1189               if( nbExtOk == 1 ) {
1190                 tPp = iT;
1191                 aTPS.Point(indExt).Parameter(u,v);
1192                 break;
1193               }
1194             }
1195
1196             if (tPp != 0 && tPp != NbOfPnts) {
1197               gp_Pnt2d aPp = gp_Pnt2d(u,v);
1198               gp_Pnt2d aPn;
1199               Standard_Boolean isFound = Standard_False;
1200               for (j = tPp + 1; j <= NbOfPnts; ++j)
1201               {
1202                 Curve->D0( Param.Value(j), pntproj );
1203                 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1204                 Dist2Min = RealLast();
1205                 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1206                   Standard_Integer indExt = 0;
1207                   for (i = 1; i <= aTPS.NbExt(); i++) {
1208                     if (aTPS.SquareDistance(i) < DistTol3d2 && aTPS.SquareDistance(i) < Dist2Min) {
1209                       Dist2Min = aTPS.SquareDistance(i);
1210                       indExt = i;
1211                     }
1212                   }
1213                   if (indExt > 0) {
1214                     aTPS.Point(indExt).Parameter(u,v);
1215                     aPn = gp_Pnt2d(u,v);
1216                     isFound = Standard_True;
1217                     break;
1218                   }     
1219                 }
1220               }
1221
1222               if( isFound ) {
1223                 gp_Vec2d atV(aPp,aPn);
1224                 Standard_Boolean isChosen = Standard_False;
1225                 for( i = 1; i <= nbSols; i++ ) {
1226                   const gp_Pnt2d& aP1 = Sols.Value(i);
1227                   gp_Vec2d asV(aP1,aPp);
1228                   if( asV.Dot(atV) > 0. ) {
1229                     isChosen = Standard_True;
1230                     Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1231                     myProjIsDone = Standard_True;
1232                     break;
1233                   }
1234                 }
1235                 if( !isChosen ) {
1236                   aExtPS.Point(GoodValue).Parameter(u,v);
1237                   Pts2d(1).SetCoord(u,v);
1238                   myProjIsDone = Standard_True;
1239                 }
1240               }
1241               else {
1242                 aExtPS.Point(GoodValue).Parameter(u,v);
1243                 Pts2d(1).SetCoord(u,v);
1244                 myProjIsDone = Standard_True;
1245               }
1246             }
1247             else {
1248               aExtPS.Point(GoodValue).Parameter(u,v);
1249               Pts2d(1).SetCoord(u,v);
1250               myProjIsDone = Standard_True;
1251             }
1252           }
1253           else {
1254             aExtPS.Point(GoodValue).Parameter(u,v);
1255             Pts2d(1).SetCoord(u,v);
1256             myProjIsDone = Standard_True;
1257           }
1258         }
1259       }
1260       
1261       //  calculate the following points with GenLocate_ExtPS
1262       // (and store the result and each parameter in a sequence)
1263       Standard_Integer usens = 0, vsens = 0; 
1264       // to know the position relatively to the period
1265       Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1266       // U0 and V0 are the points in the initialized period 
1267       // (period with u and v),
1268       // U1 and V1 are the points for construction of poles
1269       myDist = Dist2Min;
1270       for ( i = 2 ; i <= NbOfPnts ; i++) 
1271         if(myProjIsDone) {
1272           myProjIsDone = Standard_False;
1273           Dist2Min = RealLast();
1274           Curve->D0(Param.Value(i), pntproj);
1275           Extrema_GenLocateExtPS  aLocateExtPS(Surf->Surface(), TolU, TolV);
1276           aLocateExtPS.Perform(pntproj, U0, V0);
1277
1278           if (aLocateExtPS.IsDone())
1279           {
1280             if (aLocateExtPS.SquareDistance() < DistTol3d2)
1281             {  //OCC217
1282               //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1283               if (aLocateExtPS.SquareDistance() > myDist)
1284               {
1285                 myDist = aLocateExtPS.SquareDistance();
1286               }
1287               (aLocateExtPS.Point()).Parameter(U0,V0);
1288               U1 = U0 + usens*uperiod;
1289               V1 = V0 + vsens*vperiod;
1290               Pts2d(i).SetCoord(U1,V1);
1291               myProjIsDone = Standard_True;
1292             }
1293             else
1294             {
1295               Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1296               if (aGlobalExtr.IsDone())
1297               {
1298                 Standard_Real LocalMinSqDist = RealLast();
1299                 Standard_Integer imin = 0;
1300                 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1301                 {
1302                   Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1303                   if (aSqDist < LocalMinSqDist)
1304                   {
1305                     LocalMinSqDist = aSqDist;
1306                     imin = isol;
1307                   }
1308                 }
1309                 if (LocalMinSqDist < DistTol3d2)
1310                 {
1311                   if (LocalMinSqDist > myDist)
1312                   {
1313                     myDist = LocalMinSqDist;
1314                   }
1315                   Standard_Real LocalU, LocalV;
1316                   aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1317                   if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1318                   {
1319                     if (LocalU > U0)
1320                       usens = -1;
1321                     else
1322                       usens = 1;
1323                   }
1324                   if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1325                   {
1326                     if (LocalV > V0)
1327                       vsens = -1;
1328                     else
1329                       vsens = 1;
1330                   }
1331                   U0 = LocalU; V0 = LocalV;
1332                   U1 = U0 + usens*uperiod;
1333                   V1 = V0 + vsens*vperiod;
1334                   Pts2d(i).SetCoord(U1,V1);
1335                   myProjIsDone = Standard_True;
1336
1337                   if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
1338                   {//Make 1st point more precise for periodic surfaces
1339                     const Standard_Integer aSize = 3;
1340                     const gp_Pnt2d aP(Pts2d(2)); 
1341                     Standard_Real aUpar[aSize], aVpar[aSize];
1342                     Pts2d(1).Coord(aUpar[1], aVpar[1]);
1343                     aUpar[0] = aUpar[1] - uperiod;
1344                     aUpar[2] = aUpar[1] + uperiod;
1345                     aVpar[0] = aVpar[1] - vperiod;
1346                     aVpar[2] = aVpar[1] + vperiod;
1347
1348                     Standard_Real aSQdistMin = RealLast();
1349                     Standard_Integer aBestUInd = 1, aBestVInd = 1;
1350                     const Standard_Integer  aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
1351                                             aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
1352                     for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
1353                     {
1354                       for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
1355                       {
1356                         Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
1357                         if(aSQdist < aSQdistMin)
1358                         {
1359                           aSQdistMin = aSQdist;
1360                           aBestUInd = uInd;
1361                           aBestVInd = vInd;
1362                         }
1363                       }
1364                     }
1365
1366                     Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
1367                   }//if(i == 2) condition
1368                 }
1369               }
1370             }
1371           }
1372           if(!myProjIsDone && uperiod) {
1373             Standard_Real aUinf, aUsup, Uaux;
1374             aUinf = Surf->Surface().FirstUParameter();
1375             aUsup = Surf->Surface().LastUParameter();
1376             if((aUsup - U0) > (U0 - aUinf)) 
1377               Uaux = 2*aUinf - U0 + uperiod;
1378             else 
1379               Uaux = 2*aUsup - U0 - uperiod;
1380
1381             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1382             locext.Perform(pntproj, Uaux, V0);
1383
1384             if (locext.IsDone())
1385               if (locext.SquareDistance() < DistTol3d2) {  //OCC217
1386               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1387                 if (locext.SquareDistance() > myDist)
1388                 {
1389                   myDist = locext.SquareDistance();
1390                 }
1391                 (locext.Point()).Parameter(u,v);
1392                 if((aUsup - U0) > (U0 - aUinf)) 
1393                   usens--;
1394                 else 
1395                   usens++;
1396                 U0 = u; V0 = v;
1397                 U1 = U0 + usens*uperiod;
1398                 V1 = V0 + vsens*vperiod;
1399                 Pts2d(i).SetCoord(U1,V1);
1400                 myProjIsDone = Standard_True;
1401               }
1402           }
1403           if(!myProjIsDone && vperiod) {
1404             Standard_Real aVinf, aVsup, Vaux;
1405             aVinf = Surf->Surface().FirstVParameter();
1406             aVsup = Surf->Surface().LastVParameter();
1407             if((aVsup - V0) > (V0 - aVinf)) 
1408               Vaux = 2*aVinf - V0 + vperiod;
1409             else 
1410               Vaux = 2*aVsup - V0 - vperiod;
1411
1412             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1413             locext.Perform(pntproj, U0, Vaux);
1414
1415             if (locext.IsDone())
1416               if (locext.SquareDistance() < DistTol3d2) {  //OCC217
1417               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1418                 if (locext.SquareDistance() > myDist)
1419                 {
1420                   myDist = locext.SquareDistance();
1421                 }
1422                 (locext.Point()).Parameter(u, v);
1423                 if((aVsup - V0) > (V0 - aVinf)) 
1424                   vsens--;
1425                 else 
1426                   vsens++;
1427                 U0 = u; V0 = v;
1428                 U1 = U0 + usens*uperiod;
1429                 V1 = V0 + vsens*vperiod;
1430                 Pts2d(i).SetCoord(U1,V1);
1431                 myProjIsDone = Standard_True;
1432               }
1433           }     
1434           if(!myProjIsDone && uperiod && vperiod) {
1435             Standard_Real Uaux, Vaux;
1436             if((Usup - U0) > (U0 - Uinf)) 
1437               Uaux = 2*Uinf - U0 + uperiod;
1438             else 
1439               Uaux = 2*Usup - U0 - uperiod;
1440             if((Vsup - V0) > (V0 - Vinf)) 
1441               Vaux = 2*Vinf - V0 + vperiod;
1442             else 
1443               Vaux = 2*Vsup - V0 - vperiod;
1444
1445             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1446             locext.Perform(pntproj, Uaux, Vaux);
1447
1448             if (locext.IsDone())
1449               if (locext.SquareDistance() < DistTol3d2) {
1450               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1451                 if (locext.SquareDistance() > myDist)
1452                 {
1453                   myDist = locext.SquareDistance();
1454                 }
1455                 (locext.Point()).Parameter(u, v);
1456                 if((Usup - U0) > (U0 - Uinf)) 
1457                   usens--;
1458                 else 
1459                   usens++;
1460                 if((Vsup - V0) > (V0 - Vinf)) 
1461                   vsens--;
1462                 else 
1463                   vsens++;
1464                 U0 = u; V0 = v;
1465                 U1 = U0 + usens*uperiod;
1466                 V1 = V0 + vsens*vperiod;
1467                 Pts2d(i).SetCoord(U1,V1);
1468                 myProjIsDone = Standard_True;
1469               }
1470           }
1471           if(!myProjIsDone) {
1472             Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1473             if (ext.IsDone()) {
1474               Dist2Min = ext.SquareDistance(1);
1475               Standard_Integer aGoodValue = 1;
1476               for ( j = 2 ; j <= ext.NbExt() ; j++ )
1477                 if( Dist2Min > ext.SquareDistance(j)) {
1478                   Dist2Min = ext.SquareDistance(j);
1479                   aGoodValue = j;
1480                 }
1481               if (Dist2Min < DistTol3d2) {
1482               //if (Dist2Min < Tol3d * Tol3d) {
1483                 if (Dist2Min > myDist)
1484                 {
1485                   myDist = Dist2Min;
1486                 }
1487                 (ext.Point(aGoodValue)).Parameter(u, v);
1488                 if(uperiod) {
1489                   if((U0 - u) > (2*uperiod/3)) {
1490                     usens++;
1491                   }
1492                   else
1493                     if((u - U0) > (2*uperiod/3)) {
1494                       usens--;
1495                     }
1496     }
1497                 if(vperiod) {
1498                   if((V0 - v) > (vperiod/2)) {
1499                     vsens++;
1500                   }
1501                   else
1502                     if((v - V0) > (vperiod/2)) {
1503                       vsens--;
1504                     }
1505       }
1506                 U0 = u; V0 = v;
1507                 U1 = U0 + usens*uperiod;
1508                 V1 = V0 + vsens*vperiod;
1509                 Pts2d(i).SetCoord(U1,V1);
1510                 myProjIsDone = Standard_True;
1511               }
1512             }
1513           }
1514         }
1515         else break;
1516     }    
1517   }
1518   // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1519   if(myProjIsDone) {
1520     myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1521     //jgv: put the curve into parametric range
1522     gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1523     Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1524     Standard_Real sense = 0.;
1525     if (uperiod)
1526     {
1527       if (TestU < Uinf - TolU)
1528         sense = 1.;
1529       else if (TestU > Usup + TolU)
1530         sense = -1;
1531       while (TestU < Uinf - TolU || TestU > Usup + TolU)
1532         TestU += sense * uperiod;
1533     }
1534     if (vperiod)
1535     {
1536       sense = 0.;
1537       if (TestV < Vinf - TolV)
1538         sense = 1.;
1539       else if (TestV > Vsup + TolV)
1540         sense = -1.;
1541       while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1542         TestV += sense * vperiod;
1543     }
1544     gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1545     if (Abs(Offset.X()) > gp::Resolution() ||
1546         Abs(Offset.Y()) > gp::Resolution())
1547       myBSpline->Translate(Offset);
1548     //////////////////////////////////////////
1549     Geom2dAdaptor_Curve GAC(myBSpline);
1550     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1551 #ifdef OCCT_DEBUG
1552 //    char name [100];
1553 //    sprintf(name,"%s_%d","build",compteur++);
1554 //    DrawTrSurf::Set(name,myBSpline);
1555 #endif
1556     return IC2d;
1557   }
1558   else {
1559 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1560 //     Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1561 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1562     return Handle(Adaptor2d_HCurve2d)();
1563   }
1564 //  myProjIsDone = Standard_False;
1565 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1566 //   Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1567 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1568 }
1569
1570 //=======================================================================
1571 //function : ProjLib_ProjectUsingInitialCurve2d
1572 //purpose  : 
1573 //=======================================================================
1574 Handle(Geom2d_BSplineCurve) 
1575      ProjLib_ComputeApproxOnPolarSurface::
1576      ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1577                                 const Handle(Adaptor3d_HSurface)& Surf,
1578                                 const Handle(Adaptor2d_HCurve2d)& InitCurve2d)
1579 {  
1580   //OCC217
1581   Standard_Real Tol3d = myTolerance;
1582   Standard_Real DistTol3d = 100.0*Tol3d;
1583   if(myMaxDist > 0.)
1584   {
1585     DistTol3d = myMaxDist;
1586   }
1587   Standard_Real DistTol3d2 = DistTol3d * DistTol3d;
1588   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1589   Standard_Real Tol2d = Max(Sqrt(TolU*TolU + TolV*TolV), Precision::PConfusion());
1590
1591   Standard_Integer i;
1592   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1593   GeomAbs_CurveType TheTypeC = Curve->GetType();
1594   if(TheTypeS == GeomAbs_Plane) {
1595     Standard_Real S, T;
1596     gp_Pln Plane = Surf->Plane();
1597     if(TheTypeC == GeomAbs_BSplineCurve) {
1598       myTolReached = Precision::Confusion();
1599       Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1600       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1601       for(i = 1;i <= Curve->NbPoles();i++) {
1602         ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1603         Poles2d(i).SetCoord(S,T);
1604       }
1605       TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1606       BSC->Knots(Knots);
1607       TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1608       BSC->Multiplicities(Mults);
1609       if(BSC->IsRational()) {
1610         TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1611         BSC->Weights(Weights); 
1612         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1613                                        BSC->Degree(), BSC->IsPeriodic()) ;
1614       }
1615       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1616                                      BSC->Degree(), BSC->IsPeriodic()) ;
1617       
1618     }
1619     if(TheTypeC == GeomAbs_BezierCurve) {
1620       myTolReached = Precision::Confusion();
1621       Handle(Geom_BezierCurve) BC = Curve->Bezier();
1622       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1623       for(i = 1;i <= Curve->NbPoles();i++) {
1624         ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1625         Poles2d(i).SetCoord(S,T);
1626       }
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   if(TheTypeS == GeomAbs_BSplineSurface) {
1643     Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1644     if((BSS->MaxDegree() == 1) &&
1645        (BSS->NbUPoles() == 2) &&
1646        (BSS->NbVPoles() == 2)) {
1647       gp_Pnt p11 = BSS->Pole(1,1);
1648       gp_Pnt p12 = BSS->Pole(1,2);
1649       gp_Pnt p21 = BSS->Pole(2,1);
1650       gp_Pnt p22 = BSS->Pole(2,2);
1651       gp_Vec V1(p11,p12);
1652       gp_Vec V2(p21,p22);
1653       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){  
1654         Standard_Integer Dist2Min = IntegerLast();
1655         Standard_Real u,v;
1656         if(TheTypeC == GeomAbs_BSplineCurve) {
1657           myTolReached = Tol3d;
1658           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1659           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1660           for(i = 1;i <= Curve->NbPoles();i++) {
1661             myProjIsDone = Standard_False;
1662             Dist2Min = IntegerLast();
1663
1664             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1665             extrloc.Perform(BSC->Pole(i), (p11.X()+p22.X())/2, (p11.Y()+p22.Y())/2);
1666
1667             if (extrloc.IsDone()) {
1668               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1669               if (Dist2Min < DistTol3d2) {  
1670                 (extrloc.Point()).Parameter(u,v);
1671                 Poles2d(i).SetCoord(u,v);
1672                 myProjIsDone = Standard_True;
1673               }
1674               else break;
1675             }
1676             else break;
1677             if(!myProjIsDone) 
1678               break;
1679           }
1680           if(myProjIsDone) {
1681             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1682             BSC->Knots(Knots);
1683             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1684             BSC->Multiplicities(Mults);
1685             if(BSC->IsRational()) {
1686               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1687               BSC->Weights(Weights); 
1688               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1689                                              BSC->Degree(), BSC->IsPeriodic()) ;
1690             }
1691             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1692                                            BSC->Degree(), BSC->IsPeriodic()) ;
1693             
1694             
1695           }
1696         } 
1697         if(TheTypeC == GeomAbs_BezierCurve) {
1698           myTolReached = Tol3d;
1699           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1700           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1701           for(i = 1;i <= Curve->NbPoles();i++) {
1702             Dist2Min = IntegerLast();
1703
1704             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1705             extrloc.Perform(BC->Pole(i), 0.5, 0.5);
1706
1707             if (extrloc.IsDone()) {
1708               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1709               if (Dist2Min < DistTol3d2) {  
1710                 (extrloc.Point()).Parameter(u,v);
1711                 Poles2d(i).SetCoord(u,v);
1712                 myProjIsDone = Standard_True;
1713               }
1714               else break;
1715             }
1716             else break;
1717             if(myProjIsDone) 
1718               myProjIsDone = Standard_False;
1719             else break;
1720           }
1721           if(myProjIsDone) {
1722             TColStd_Array1OfReal Knots(1, 2);
1723             Knots.SetValue(1,0.0);
1724             Knots.SetValue(2,1.0);
1725             TColStd_Array1OfInteger Mults(1, 2);
1726             Mults.Init(BC->NbPoles());
1727             if(BC->IsRational()) {
1728               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1729               BC->Weights(Weights); 
1730               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1731                                                     BC->Degree(), BC->IsPeriodic()) ;
1732             }
1733             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1734                                                   BC->Degree(), BC->IsPeriodic()) ;
1735           }
1736         } 
1737       }
1738     }
1739   }
1740   else if(TheTypeS == GeomAbs_BezierSurface) {
1741     Handle(Geom_BezierSurface) BS = Surf->Bezier();
1742     if((BS->MaxDegree() == 1) &&
1743        (BS->NbUPoles() == 2) &&
1744        (BS->NbVPoles() == 2)) {
1745       gp_Pnt p11 = BS->Pole(1,1);
1746       gp_Pnt p12 = BS->Pole(1,2);
1747       gp_Pnt p21 = BS->Pole(2,1);
1748       gp_Pnt p22 = BS->Pole(2,2);
1749       gp_Vec V1(p11,p12);
1750       gp_Vec V2(p21,p22);
1751       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ 
1752         Standard_Integer Dist2Min = IntegerLast();
1753         Standard_Real u,v;
1754  
1755 //      gp_Pnt pntproj;
1756         if(TheTypeC == GeomAbs_BSplineCurve) {
1757           myTolReached = Tol3d;
1758           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1759           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1760           for(i = 1;i <= Curve->NbPoles();i++) {
1761             myProjIsDone = Standard_False;
1762             Dist2Min = IntegerLast();
1763
1764             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1765             extrloc.Perform(BSC->Pole(i), (p11.X()+p22.X())/2, (p11.Y()+p22.Y())/2);
1766
1767             if (extrloc.IsDone()) {
1768               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1769               if (Dist2Min < DistTol3d2) {  
1770                 (extrloc.Point()).Parameter(u,v);
1771                 Poles2d(i).SetCoord(u,v);
1772                 myProjIsDone = Standard_True;
1773               }
1774               else break;
1775             }
1776             else break;
1777             if(!myProjIsDone) 
1778               break;
1779           }
1780           if(myProjIsDone) {
1781             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1782             BSC->Knots(Knots);
1783             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1784             BSC->Multiplicities(Mults);
1785             if(BSC->IsRational()) {
1786               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1787               BSC->Weights(Weights); 
1788               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1789                                                     BSC->Degree(), BSC->IsPeriodic()) ;
1790             }
1791             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1792                                                   BSC->Degree(), BSC->IsPeriodic()) ;
1793             
1794             
1795           }
1796         } 
1797         if(TheTypeC == GeomAbs_BezierCurve) {
1798           myTolReached = Tol3d;
1799           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1800           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1801           for(i = 1;i <= Curve->NbPoles();i++) {
1802             Dist2Min = IntegerLast();
1803
1804             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1805             extrloc.Perform(BC->Pole(i), 0.5, 0.5);
1806
1807             if (extrloc.IsDone()) {
1808               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1809               if (Dist2Min < DistTol3d2) {  
1810                 (extrloc.Point()).Parameter(u,v);
1811                 Poles2d(i).SetCoord(u,v);
1812                 myProjIsDone = Standard_True;
1813               }
1814               else break;
1815             }
1816             else break;
1817             if(myProjIsDone) 
1818               myProjIsDone = Standard_False;
1819             else break;
1820           }
1821           if(myProjIsDone) {
1822             TColStd_Array1OfReal Knots(1, 2);
1823             Knots.SetValue(1,0.0);
1824             Knots.SetValue(2,1.0);
1825             TColStd_Array1OfInteger Mults(1, 2);
1826             Mults.Init(BC->NbPoles());
1827             if(BC->IsRational()) {
1828               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1829               BC->Weights(Weights); 
1830               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1831                                                     BC->Degree(), BC->IsPeriodic()) ;
1832             }
1833             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1834                                                   BC->Degree(), BC->IsPeriodic()) ;
1835           }
1836         } 
1837       }
1838     }
1839   }
1840
1841   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;  
1842
1843 #ifdef OCCT_DEBUG
1844   Standard_Integer Nb = 50;
1845   
1846   Standard_Real U, U1, U2;
1847   U1 = F.FirstParameter();
1848   U2 = F.LastParameter();
1849   
1850   TColgp_Array1OfPnt2d    DummyPoles(1,Nb+1);
1851   TColStd_Array1OfReal    DummyKnots(1,Nb+1);
1852   TColStd_Array1OfInteger DummyMults(1,Nb+1);
1853   DummyMults.Init(1);
1854   DummyMults(1) = 2;
1855   DummyMults(Nb+1) = 2;
1856   for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1857     U = (Nb-ij)*U1 + ij*U2;
1858     U /= Nb;
1859     DummyPoles(ij+1) = F.Value(U);
1860     DummyKnots(ij+1) = ij;
1861   }
1862   Handle(Geom2d_BSplineCurve) DummyC2d =
1863     new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1864 #ifdef DRAW
1865   Standard_CString Temp = "bs2d";
1866   DrawTrSurf::Set(Temp,DummyC2d);
1867 #endif
1868 //  DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1869   Handle(Geom2dAdaptor_HCurve) DDD = 
1870     Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1871   
1872 #ifdef DRAW
1873   Temp = "initc2d";
1874   DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1875 #endif
1876 //  DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1877 #endif
1878
1879   Standard_Integer Deg1,Deg2;
1880   Deg1 = 2; 
1881   Deg2 = 8;  
1882   if(myDegMin > 0)
1883   {
1884     Deg1 = myDegMin;
1885   }
1886   if(myDegMax > 0)
1887   {
1888     Deg2 = myDegMax;
1889   }
1890   Standard_Integer aMaxSegments = 1000;
1891   if(myMaxSegments > 0)
1892   {
1893     aMaxSegments = myMaxSegments;
1894   }
1895   AppParCurves_Constraint aFistC = AppParCurves_TangencyPoint, aLastC = AppParCurves_TangencyPoint;
1896   if(myBndPnt != AppParCurves_TangencyPoint)
1897   {
1898     aFistC = myBndPnt; 
1899     aLastC = myBndPnt;
1900   }
1901
1902   if (myDist > 10.*Tol3d)
1903   {
1904     aFistC = AppParCurves_PassPoint; 
1905     aLastC = AppParCurves_PassPoint;
1906   }
1907
1908   Approx_FitAndDivide2d Fit(Deg1, Deg2, Tol3d, Tol2d, Standard_True, aFistC, aLastC);
1909   Fit.SetMaxSegments(aMaxSegments);
1910   if (InitCurve2d->GetType() == GeomAbs_Line)
1911   {
1912     Fit.SetInvOrder(Standard_False);
1913   }
1914   Fit.Perform(F);
1915
1916   Standard_Real anOldTol2d = Tol2d;
1917   Standard_Real aNewTol2d = 0;
1918   if(Fit.IsAllApproximated()) {
1919     Standard_Integer j;
1920     Standard_Integer NbCurves = Fit.NbMultiCurves();
1921     Standard_Integer MaxDeg = 0;
1922     // To transform the MultiCurve into BSpline, it is required that all  
1923     // Bezier constituing it have the same degree -> Calculation of MaxDeg
1924     Standard_Integer NbPoles  = 1;
1925     for (j = 1; j <= NbCurves; j++) {
1926       Standard_Integer Deg = Fit.Value(j).Degree();
1927       MaxDeg = Max ( MaxDeg, Deg);
1928       Fit.Error(j,Tol3d, Tol2d);
1929       aNewTol2d = Max(aNewTol2d, Tol2d);
1930     }
1931     //
1932     myTolReached = Max(myTolReached, myTolerance * (aNewTol2d / anOldTol2d));
1933     //
1934     NbPoles = MaxDeg * NbCurves + 1;               //Tops on the BSpline
1935     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
1936       
1937     TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1938     
1939     TColStd_Array1OfReal Knots( 1, NbCurves + 1);  //Nodes of the BSpline
1940     
1941     Standard_Integer Compt = 1;
1942     for (i = 1; i <= NbCurves; i++) {
1943       Fit.Parameters(i, Knots(i), Knots(i+1)); 
1944       AppParCurves_MultiCurve MC = Fit.Value( i);       //Load the Ith Curve
1945       TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1946       MC.Curve(1, Poles2d);
1947       
1948       //Eventual augmentation of the degree
1949       Standard_Integer Inc = MaxDeg - MC.Degree();
1950       if ( Inc > 0) {
1951 //      BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(), 
1952         BSplCLib::IncreaseDegree( MaxDeg, Poles2d, BSplCLib::NoWeights(), 
1953                          TempPoles, BSplCLib::NoWeights());
1954         //update of tops of the PCurve
1955         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1956           Poles.SetValue( Compt, TempPoles( k));
1957           Compt++;
1958         }
1959       }
1960       else {
1961         //update of tops of the PCurve
1962         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1963           Poles.SetValue( Compt, Poles2d( k));
1964           Compt++;
1965         }
1966       } 
1967       
1968       Compt--;
1969     }
1970     
1971     //update of fields of  ProjLib_Approx
1972     Standard_Integer NbKnots = NbCurves + 1;
1973
1974     TColStd_Array1OfInteger   Mults( 1, NbKnots);
1975     Mults.Init(MaxDeg);
1976     Mults.SetValue( 1, MaxDeg + 1);
1977     Mults.SetValue(NbKnots, MaxDeg + 1);
1978     myProjIsDone = Standard_True;
1979     Handle(Geom2d_BSplineCurve) Dummy =
1980       new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1981     
1982     // try to smoother the Curve GeomAbs_C1.
1983
1984     Standard_Boolean OK = Standard_True;
1985     Standard_Real aSmoothTol = Max(Precision::Confusion(), aNewTol2d);
1986     if (myBndPnt == AppParCurves_PassPoint)
1987     {
1988       aSmoothTol *= 10.;
1989     }
1990     for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1991       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1, aSmoothTol);  
1992     }
1993 #ifdef OCCT_DEBUG
1994     if (!OK) {
1995       std::cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<std::endl;
1996     }
1997 #endif 
1998     return Dummy;
1999   }
2000   return Handle(Geom2d_BSplineCurve)();
2001 }
2002
2003 //=======================================================================
2004 //function : BSpline
2005 //purpose  : 
2006 //=======================================================================
2007
2008 Handle(Geom2d_BSplineCurve)  
2009      ProjLib_ComputeApproxOnPolarSurface::BSpline() const 
2010      
2011 {
2012   return myBSpline ;
2013 }
2014
2015 //=======================================================================
2016 //function : Curve2d
2017 //purpose  : 
2018 //=======================================================================
2019
2020 Handle(Geom2d_Curve)  
2021      ProjLib_ComputeApproxOnPolarSurface::Curve2d() const 
2022      
2023 {
2024   Standard_NoSuchObject_Raise_if
2025     (!myProjIsDone,
2026      "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
2027   return my2ndCurve ;
2028 }
2029
2030
2031 //=======================================================================
2032 //function : IsDone
2033 //purpose  : 
2034 //=======================================================================
2035
2036 Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const 
2037      
2038 {
2039   return myProjIsDone;
2040 }
2041 //=======================================================================
2042 //function : SetTolerance
2043 //purpose  : 
2044 //=======================================================================
2045
2046 void ProjLib_ComputeApproxOnPolarSurface::SetTolerance(const Standard_Real theTol) 
2047      
2048 {
2049   myTolerance = theTol;
2050 }
2051 //=======================================================================
2052 //function : SetDegree
2053 //purpose  : 
2054 //=======================================================================
2055 void ProjLib_ComputeApproxOnPolarSurface::SetDegree(
2056                                        const Standard_Integer theDegMin, 
2057                                        const Standard_Integer theDegMax)
2058 {
2059   myDegMin = theDegMin;
2060   myDegMax = theDegMax;
2061 }
2062 //=======================================================================
2063 //function : SetMaxSegments
2064 //purpose  : 
2065 //=======================================================================
2066 void ProjLib_ComputeApproxOnPolarSurface::SetMaxSegments(
2067                                    const Standard_Integer theMaxSegments)
2068 {
2069   myMaxSegments = theMaxSegments;
2070 }
2071
2072 //=======================================================================
2073 //function : SetBndPnt
2074 //purpose  : 
2075 //=======================================================================
2076 void ProjLib_ComputeApproxOnPolarSurface::SetBndPnt(
2077                                  const AppParCurves_Constraint theBndPnt)
2078 {
2079   myBndPnt = theBndPnt;
2080 }
2081
2082 //=======================================================================
2083 //function : SetMaxDist
2084 //purpose  : 
2085 //=======================================================================
2086 void ProjLib_ComputeApproxOnPolarSurface::SetMaxDist(
2087                                           const Standard_Real theMaxDist)
2088 {
2089   myMaxDist = theMaxDist;
2090 }
2091
2092 //=======================================================================
2093 //function : Tolerance
2094 //purpose  : 
2095 //=======================================================================
2096
2097 Standard_Real ProjLib_ComputeApproxOnPolarSurface::Tolerance() const 
2098      
2099 {
2100   return myTolReached;
2101 }
2102