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