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