0024510: Remove unused local variables
[occt.git] / src / ProjLib / ProjLib_ComputeApprox.cxx
1 // Created on: 1993-09-07
2 // Created by: Bruno DUMORTIER
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and / or modify it
9 // under the terms of the GNU Lesser General Public version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // modified by NIZHNY-OFV  Thu Jan 20 11:04:19 2005
18
19 #include <ProjLib_ComputeApprox.hxx>
20
21 #include <GeomAbs_SurfaceType.hxx>
22 #include <GeomAbs_CurveType.hxx>
23 #include <AppCont_Function2d.hxx>
24 #include <Convert_CompBezierCurves2dToBSplineCurve2d.hxx>
25 #include <ElSLib.hxx>
26 #include <ElCLib.hxx>
27 #include <BSplCLib.hxx>
28 #include <Standard_NoSuchObject.hxx>
29 #include <Geom_UndefinedDerivative.hxx>
30 #include <gp.hxx>
31 #include <gp_Trsf.hxx>
32 #include <Precision.hxx>
33 #include <Approx_FitAndDivide2d.hxx>
34 #include <AppParCurves_MultiCurve.hxx>
35 #include <Handle_Adaptor3d_HCurve.hxx>
36 #include <Adaptor3d_HCurve.hxx>
37 #include <Handle_Adaptor3d_HSurface.hxx>
38 #include <Adaptor3d_HSurface.hxx>
39 #include <TColgp_Array1OfPnt2d.hxx>
40 #include <TColgp_Array1OfPnt.hxx>
41 #include <TColStd_Array1OfReal.hxx>
42 #include <TColStd_Array1OfInteger.hxx>
43 #include <Geom_BSplineCurve.hxx>
44 #include <Geom_BezierCurve.hxx>
45 #include <Geom2d_BSplineCurve.hxx>
46 #include <Geom2d_BezierCurve.hxx>
47
48 #ifdef DRAW
49 #include <DrawTrSurf.hxx>
50 #endif
51 #ifdef DEB
52 static Standard_Boolean AffichValue = Standard_False;
53 #endif    
54
55 static 
56   void Parameters(const Handle(Adaptor3d_HCurve)&   myCurve,
57                   const Handle(Adaptor3d_HSurface)& mySurface,
58                   const gp_Pnt& aP1, 
59                   const Standard_Integer iFirst,
60                   const Standard_Real aTolU, 
61                   Standard_Real& aU, 
62                   Standard_Real& aV);
63
64 //=======================================================================
65 //function : IsEqual
66 //purpose  : 
67 //=======================================================================
68 // OFV:
69 static inline Standard_Boolean IsEqual(Standard_Real Check,Standard_Real With,Standard_Real Toler)
70 {
71   return ((Abs(Check - With) < Toler) ? Standard_True : Standard_False);
72 }
73
74
75 //=======================================================================
76 //function : Value
77 //purpose  : 
78 //=======================================================================
79
80 static gp_Pnt2d Function_Value(const Standard_Real U,
81                                const Handle(Adaptor3d_HCurve)&   myCurve,
82                                const Handle(Adaptor3d_HSurface)& mySurface,
83                                const Standard_Real U1,
84                                const Standard_Real U2, 
85                                const Standard_Real V1,
86                                const Standard_Real V2,
87                                const Standard_Boolean UCouture,
88                                const Standard_Boolean VCouture ) 
89 {
90   Standard_Real S = 0., T = 0.;
91
92   gp_Pnt P3d = myCurve->Value(U);
93   GeomAbs_SurfaceType SType = mySurface->GetType();
94
95   switch ( SType ) {
96     
97   case GeomAbs_Plane:
98     {
99       gp_Pln Plane = mySurface->Plane();
100       ElSLib::Parameters( Plane, P3d, S, T);
101       break;
102     }
103   case GeomAbs_Cylinder:
104     {
105       gp_Cylinder Cylinder = mySurface->Cylinder();
106       ElSLib::Parameters( Cylinder, P3d, S, T);
107       break;
108     }
109   case GeomAbs_Cone:
110     {
111       gp_Cone Cone = mySurface->Cone();
112       ElSLib::Parameters( Cone, P3d, S, T);
113       break;
114     }
115   case GeomAbs_Sphere:
116     {
117       gp_Sphere Sphere = mySurface->Sphere();
118       ElSLib::Parameters(Sphere, P3d, S, T);
119       break;
120     }
121   case GeomAbs_Torus:
122     {
123       gp_Torus Torus = mySurface->Torus();
124       ElSLib::Parameters( Torus, P3d, S, T);
125       break;
126     }
127   default:
128     Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
129   }
130
131   if ( UCouture) {
132     if(S < U1 || S > U2)
133       S = ElCLib::InPeriod(S, U1, U2);
134   }
135  
136   if ( VCouture) {
137     if(SType == GeomAbs_Sphere) {
138       if ( Abs( S - U1 ) > M_PI ) {
139               T = M_PI - T;
140               S = M_PI + S;
141       }
142       if(S > U1 || S < U2)
143         S = ElCLib::InPeriod(S, U1, U2);
144     }
145     if(T < V1 || T > V2)
146       T = ElCLib::InPeriod(T, V1, V2);
147   }
148   
149   return gp_Pnt2d(S, T);
150 }
151 //=======================================================================
152 //function : D1
153 //purpose  : 
154 //=======================================================================
155 static Standard_Boolean Function_D1( const Standard_Real U, 
156                                     gp_Pnt2d&            P,
157                                     gp_Vec2d&            D,
158                                     const Handle(Adaptor3d_HCurve)&   myCurve,
159                                     const Handle(Adaptor3d_HSurface)& mySurface,
160                                     const Standard_Real U1,
161                                     const Standard_Real U2, 
162                                     const Standard_Real V1,
163                                     const Standard_Real V2,
164                                     const Standard_Boolean UCouture,
165                                     const Standard_Boolean VCouture )
166 {
167   gp_Pnt P3d;
168   Standard_Real dU, dV;
169   
170   P = Function_Value(U,myCurve,mySurface,U1,U2,V1,V2,UCouture,VCouture);
171
172   GeomAbs_SurfaceType Type = mySurface->GetType();
173   
174   switch ( Type) {
175   case GeomAbs_Plane:
176   case GeomAbs_Cone:
177   case GeomAbs_Cylinder:
178   case GeomAbs_Sphere:
179   case GeomAbs_Torus: 
180     {
181       gp_Vec D1U, D1V;
182       gp_Vec T;
183       myCurve->D1(U,P3d,T);
184       mySurface->D1(P.X(),P.Y(),P3d,D1U,D1V);
185       
186       dU = T.Dot(D1U);
187       dV = T.Dot(D1V);
188       Standard_Real Nu = D1U.SquareMagnitude();
189       Standard_Real Nv = D1V.SquareMagnitude(); 
190       
191       if ( Nu < Epsilon(1.) || Nv < Epsilon(1.))
192               return Standard_False;
193       
194       dU /= Nu;
195       dV /= Nv;
196       D = gp_Vec2d( dU, dV);   
197     }
198     break;
199     
200   default:
201     return Standard_False;
202   }
203   
204   return Standard_True;
205 }
206
207 //=======================================================================
208 //function : Function_SetUVBounds
209 //purpose  : 
210 //=======================================================================
211 static void Function_SetUVBounds(Standard_Real& myU1, 
212                                  Standard_Real& myU2,
213                                  Standard_Real& myV1,
214                                  Standard_Real& myV2,
215                                  Standard_Boolean& UCouture,
216                                  Standard_Boolean& VCouture,
217                                  const Handle(Adaptor3d_HCurve)&   myCurve,
218                                  const Handle(Adaptor3d_HSurface)& mySurface) 
219 {
220   Standard_Real W1, W2, W;
221   gp_Pnt P1, P2, P;
222   //
223   W1 = myCurve->FirstParameter();
224   W2 = myCurve->LastParameter ();
225   W  = 0.5*(W1+W2);
226   // on ouvre l`intervalle
227   // W1 += 1.0e-9;
228   // W2 -= 1.0e-9;
229   P1 = myCurve->Value(W1);
230   P2 = myCurve->Value(W2);
231   P  = myCurve->Value(W);
232
233   switch ( mySurface->GetType()) {
234
235     case GeomAbs_Cone:    {
236       gp_Cone Cone = mySurface->Cone();
237       VCouture = Standard_False;
238
239       switch( myCurve->GetType() ){
240       case GeomAbs_Parabola:
241       case GeomAbs_Hyperbola:
242       case GeomAbs_Ellipse:{
243         Standard_Real U1, U2, V1, V2, U , V;
244         ElSLib::Parameters( Cone, P1, U1, V1);
245         ElSLib::Parameters( Cone, P2, U2, V2);
246         ElSLib::Parameters( Cone, P , U , V );
247         myU1 = Min(U1,U2);
248         myU2 = Max(U1,U2);
249         if (  ( U1 < U  &&  U < U2 ) && !myCurve->IsClosed() ) {
250           UCouture = Standard_False;
251         }      
252         else {
253           UCouture = Standard_True;
254           myU2 = myU1 + 2*M_PI;
255         }
256
257       }
258       break;
259       default:    { 
260         Standard_Real U1, V1, U , V, Delta = 0., d = 0., pmin = W1, pmax = W1, dmax = 0., Uf, Ul;
261         ElSLib::Parameters( Cone, P1, U1, V1);
262         ElSLib::Parameters( Cone, P2, Ul, V1);
263         myU1 = U1; myU2 = U1; Uf = U1; 
264         Standard_Real Step = .1;
265         Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
266         nbp = Max(nbp, 3);
267         Step = (W2 - W1) / (nbp - 1);
268         Standard_Boolean isclandper = (!(myCurve->IsClosed()) && !(myCurve->IsPeriodic()));
269         for(Standard_Real par = W1 + Step; par <= W2; par += Step) {
270           if(!isclandper) par += Step;
271           P = myCurve->Value(par);
272           ElSLib::Parameters( Cone, P, U, V);
273           U += Delta;
274           d = U - U1;
275           if(d > M_PI)  {
276             if( ( (IsEqual(U,(2*M_PI),1.e-10) && (U1 >= 0. && U1 <= M_PI)) && 
277               (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,0.,1.e-10)) ) && isclandper ) U = 0.;
278             else Delta -= 2*M_PI;
279             U += Delta;
280             d = U - U1;
281           }
282           else if(d < -M_PI)      {
283             if( ( (IsEqual(U,0.,1.e-10) && (U1 >= M_PI && U1 <= (2*M_PI))) &&
284               (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,(2*M_PI),1.e-10)) ) && isclandper ) U = 2*M_PI;
285             else Delta += 2*M_PI;
286             U += Delta;
287             d = U - U1;
288           }
289           dmax = Max(dmax, Abs(d));
290           if(U < myU1) {myU1 = U; pmin = par;}
291           if(U > myU2) {myU2 = U; pmax = par;}
292           U1 = U;
293         }
294
295         if(!(Abs(pmin - W1) <= Precision::PConfusion() || Abs(pmin - W2) <= Precision::PConfusion()) ) myU1 -= dmax*.5;
296         if(!(Abs(pmax - W1) <= Precision::PConfusion() || Abs(pmax - W2) <= Precision::PConfusion()) ) myU2 += dmax*.5;
297
298         if((myU1 >=0. && myU1 <= 2*M_PI) && (myU2 >=0. && myU2 <= 2*M_PI) ) UCouture = Standard_False;
299         else{
300           U = ( myU1 + myU2 ) /2.;
301           myU1 = U - M_PI;
302           myU2 = U + M_PI;
303           UCouture = Standard_True;
304         }
305       }
306       break;
307     }// switch curve type
308   }// case Cone
309   break;
310
311   case GeomAbs_Cylinder:    {
312     gp_Cylinder Cylinder = mySurface->Cylinder();
313     VCouture = Standard_False;
314
315     if (myCurve->GetType() == GeomAbs_Ellipse) {
316
317       Standard_Real U1, U2, V1, V2, U , V;
318       ElSLib::Parameters( Cylinder, P1, U1, V1);
319       ElSLib::Parameters( Cylinder, P2, U2, V2);
320       ElSLib::Parameters( Cylinder, P , U , V );
321       myU1 = Min(U1,U2);
322       myU2 = Max(U1,U2);
323
324       if ( !myCurve->IsClosed()) {
325         if ( myU1 < U && U < myU2) {
326           U = ( myU1 + myU2 ) /2.;
327           myU1 = U - M_PI;
328           myU2 = U + M_PI;
329         }
330         else {
331           U = ( myU1 + myU2 ) /2.;
332           if ( myU1 < U) {
333             myU1 = U - 2*M_PI;
334             myU2 = U;
335           }
336           else {
337             myU1 = U;
338             myU2 = U + 2*M_PI;
339           }
340         }
341         UCouture = Standard_True;
342       }
343       else {
344         gp_Vec D1U, D1V;
345         gp_Vec T;
346         gp_Pnt P3d;
347         myCurve->D1(W1,P3d,T);
348         mySurface->D1(U1,U2,P3d,D1U,D1V);
349         Standard_Real dU = T.Dot(D1U);
350
351         UCouture = Standard_True;
352         if ( dU > 0.) {
353           myU2 = myU1 + 2*M_PI;
354         }
355         else {
356           myU2 = myU1;
357           myU1 -= 2*M_PI;
358         }
359       }
360     }
361     else {
362       Standard_Real U1, V1, U , V;
363       ElSLib::Parameters( Cylinder, P1, U1, V1);
364       Standard_Real Step = .1, Delta = 0.;
365       Standard_Real eps = M_PI, dmax = 0., d = 0.;
366       Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
367       nbp = Max(nbp, 3);
368       Step = (W2 - W1) / (nbp - 1);
369       myU1 = U1; myU2 = U1;
370       Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
371       for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
372         P = myCurve->Value(par);
373         ElSLib::Parameters( Cylinder, P, U, V);
374         U += Delta;
375         d = U - U1;
376         if(d > eps) {
377           U -= Delta;
378           Delta -= 2*M_PI;
379           U += Delta;
380           d = U - U1;
381         }
382         else if(d < -eps) {
383           U -= Delta;
384           Delta += 2*M_PI;
385           U += Delta;
386           d = U - U1;
387         }
388         dmax = Max(dmax, Abs(d));
389         if(U < myU1) {myU1 = U; pmin = par;}
390         if(U > myU2) {myU2 = U; pmax = par;}
391         U1 = U;
392       }
393
394       if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
395         Abs(pmin - W2) <= Precision::PConfusion())  ) myU1 -= dmax*.5;
396       if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
397         Abs(pmax - W2) <= Precision::PConfusion())  ) myU2 += dmax*.5;
398
399       if((myU1 >=0. && myU1 <= 2*M_PI) &&
400         (myU2 >=0. && myU2 <= 2*M_PI)    ) {
401           UCouture = Standard_False;
402       }
403       else {
404         U = ( myU1 + myU2 ) /2.;
405         myU1 = U - M_PI;
406         myU2 = U + M_PI;
407         UCouture = Standard_True;
408       }
409     }
410   }
411   break;
412   //    
413   case GeomAbs_Sphere:{
414     VCouture = Standard_False;
415     gp_Sphere SP = mySurface->Sphere();
416     if ( myCurve->GetType() == GeomAbs_Circle) {
417       UCouture = Standard_True;
418
419       // on cherche a savoir le nombre de fois que la couture est
420       // traversee.
421       // si 0 ou 2 fois : la PCurve est fermee et dans l`intervalle 
422       //                  [Uc-PI, Uc+PI] (Uc: U du centre du cercle)
423       // si 1 fois      : la PCurve est ouverte et dans l`intervalle
424       //                  [U1, U1 +/- 2*PI]
425
426       // pour determiner le nombre de solution, on resoud le systeme
427       // x^2 + y^2 + z^2     = R^2  (1)
428       // A x + B y + C z + D = 0    (2)
429       // x > 0                      (3)
430       // y = 0                      (4)
431       // REM : (1) (2)     : equation du cercle
432       //       (1) (3) (4) : equation de la couture.
433       Standard_Integer NbSolutions = 0;
434       Standard_Real A, B, C, D, R, Tol = 1.e-10;
435       Standard_Real U1, U2, V1, V2;
436       gp_Trsf Trsf;
437       //
438       gp_Circ Circle = myCurve->Circle();
439       Trsf.SetTransformation(SP.Position());
440       Circle.Transform(Trsf);
441       //
442       R = SP.Radius();
443       gp_Pln Plane( gp_Ax3(Circle.Position()));
444       Plane.Coefficients(A,B,C,D);
445       //
446       if ( Abs(C) < Tol) {
447         if ( Abs(A) > Tol) {
448           if ( (D/A) < 0.) {
449             if      ( ( R - Abs(D/A))  > Tol)  NbSolutions = 2;
450             else if ( Abs(R - Abs(D/A))< Tol)  NbSolutions = 1;
451             else                               NbSolutions = 0;
452           }
453         }
454       }
455       else {
456         Standard_Real delta = R*R*(A*A+C*C) - D*D;
457         delta *= C*C;
458         if ( Abs(delta) < Tol*Tol) {
459           if ( A*D > 0.) NbSolutions = 1;
460         }
461         else if ( delta > 0) {
462           Standard_Real xx;
463           delta = Sqrt(delta);
464           xx = -A*D+delta;
465           //      
466           if ( xx > Tol) NbSolutions++;
467           xx = -A*D-delta;
468           //    
469           if ( xx > Tol) NbSolutions++;
470         }
471       }
472       //
473
474       // box+sphere >>
475       Standard_Real UU = 0.;
476       ElSLib::Parameters(SP, P1, U1, V1);
477       Standard_Real eps = 2.*Epsilon(1.);
478       Standard_Real dt = Max(Precision::PConfusion(), 0.01*(W2-W1)); 
479       if(Abs(U1) < eps)
480       {
481         //May be U1 must be equal 2*PI?
482         gp_Pnt Pd  = myCurve->Value(W1+dt);
483         Standard_Real ud, vd;
484         ElSLib::Parameters(SP, Pd, ud, vd);
485         if(Abs(U1 - ud) > M_PI)
486         {
487           U1 = 2.*M_PI;
488         }
489       }
490       else if(Abs(2.*M_PI - U1) < eps)
491       {
492         //maybe U1 = 0.?
493         gp_Pnt Pd  = myCurve->Value(W1+dt);
494         Standard_Real ud, vd;
495         ElSLib::Parameters(SP, Pd, ud, vd);
496         if(Abs(U1 - ud) > M_PI)
497         {
498           U1 = 0.;
499         }
500       }
501       //
502       ElSLib::Parameters(SP, P2, U2, V1);
503       if(Abs(U2) < eps)
504       {
505         //May be U2 must be equal 2*PI?
506         gp_Pnt Pd  = myCurve->Value(W2-dt);
507         Standard_Real ud, vd;
508         ElSLib::Parameters(SP, Pd, ud, vd);
509         if(Abs(U2 - ud) > M_PI)
510         {
511           U2 = 2.*M_PI;
512         }
513       }
514       else if(Abs(2.*M_PI - U2) < eps)
515       {
516         //maybe U2 = 0.?
517         gp_Pnt Pd  = myCurve->Value(W2-dt);
518         Standard_Real ud, vd;
519         ElSLib::Parameters(SP, Pd, ud, vd);
520         if(Abs(U2 - ud) > M_PI)
521         {
522           U2 = 0.;
523         }
524       }
525       //
526       ElSLib::Parameters(SP, P, UU, V1);
527       Standard_Real UUmi = Min(Min(U1,UU),Min(UU,U2));
528       Standard_Real UUma = Max(Max(U1,UU),Max(UU,U2));
529       Standard_Boolean reCalc = ((UUmi >= 0. && UUmi <= M_PI) && (UUma >= 0. && UUma <= M_PI));
530       // box+sphere <<
531       P2 = myCurve->Value(W1+M_PI/8);
532       ElSLib::Parameters(SP,P2,U2,V2);
533       //
534       if ( NbSolutions == 1) {
535         if ( Abs(U1-U2) > M_PI) { // on traverse la couture
536           if ( U1 > M_PI) {
537             myU1 = U1;
538             myU2 = U1+2*M_PI;
539           }
540           else {
541             myU2 = U1;
542             myU1 = U1-2*M_PI;
543           }
544         }
545         else { // on ne traverse pas la couture
546           if ( U1 > U2) {
547             myU2 = U1;
548             myU1 = U1-2*M_PI;
549           }
550           else {
551             myU1 = U1;
552             myU2 = U1+2*M_PI;
553           }
554         }
555       }
556       else { // 0 ou 2 solutions
557         gp_Pnt Center = Circle.Location();
558         Standard_Real U,V;
559         ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
560         myU1 = U-M_PI;
561         myU2 = U+M_PI;
562       }
563       //
564       // eval the VCouture.
565       if ( (C==0) || Abs(Abs(D/C)-R) > 1.e-10) {
566         VCouture = Standard_False;
567       }
568       else {
569         VCouture = Standard_True;
570         UCouture = Standard_True;
571
572         if ( D/C < 0.) {
573           myV1 =   - M_PI / 2.;
574           myV2 = 3 * M_PI / 2.;
575         }
576         else {
577           myV1 = -3 * M_PI / 2.;
578           myV2 =      M_PI / 2.;
579         }
580
581         // si P1.Z() vaut +/- R on est sur le sommet : pas significatif.
582         gp_Pnt pp = P1.Transformed(Trsf);
583
584         if ( Abs( Abs(pp.Z()) - R) < Tol) {
585           gp_Pnt Center = Circle.Location();
586           Standard_Real U,V;
587           ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
588           myU1 = U-M_PI;
589           myU2 = U+M_PI;
590           VCouture = Standard_False;
591         }
592       }
593
594       // box+sphere >>
595       myV1 = -1.e+100; myV2 = 1.e+100;
596       Standard_Real UU1 = myU1, UU2 = myU2;
597       if((Abs(UU1) <= (2.*M_PI) && Abs(UU2) <= (2.*M_PI)) && NbSolutions == 1 && reCalc) {
598         gp_Pnt Center = Circle.Location();
599         Standard_Real U,V;
600         ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
601         myU1 = U-M_PI;
602         myU1 = Min(UU1,myU1);
603         myU2 = myU1 + 2.*M_PI;
604       }
605       // box+sphere <<
606
607     }//if ( myCurve->GetType() == GeomAbs_Circle)
608
609     else {
610       Standard_Real U1, V1, U , V;
611       ElSLib::Parameters( SP, P1, U1, V1);
612       Standard_Real Step = .1, Delta = 0.;
613       Standard_Real eps = M_PI, dmax = 0., d = 0.;
614       Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
615       nbp = Max(nbp, 3);
616       Step = (W2 - W1) / (nbp - 1);
617       myU1 = U1; myU2 = U1;
618       Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
619       for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
620         P = myCurve->Value(par);
621         ElSLib::Parameters( SP, P, U, V);
622         U += Delta;
623         d = U - U1;
624         if(d > eps) {
625           U -= Delta;
626           Delta -= 2*M_PI;
627           U += Delta;
628           d = U - U1;
629         }
630         else if(d < -eps) {
631           U -= Delta;
632           Delta += 2*M_PI;
633           U += Delta;
634           d = U - U1;
635         }
636         dmax = Max(dmax, Abs(d));
637         if(U < myU1) {myU1 = U; pmin = par;}
638         if(U > myU2) {myU2 = U; pmax = par;}
639         U1 = U;
640       }
641
642       if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
643         Abs(pmin - W2) <= Precision::PConfusion())  ) myU1 -= dmax*.5;
644       if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
645         Abs(pmax - W2) <= Precision::PConfusion())  ) myU2 += dmax*.5;
646
647       if((myU1 >=0. && myU1 <= 2*M_PI) &&
648         (myU2 >=0. && myU2 <= 2*M_PI)    ) {
649           myU1 = 0.;
650           myU2 = 2.*M_PI;
651           UCouture = Standard_False;
652       }
653       else {
654         U = ( myU1 + myU2 ) /2.;
655         myU1 = U - M_PI;
656         myU2 = U + M_PI;
657         UCouture = Standard_True;
658       }
659
660       VCouture = Standard_False;
661     }
662   }
663   break;
664   //     
665   case GeomAbs_Torus:{
666     gp_Torus TR = mySurface->Torus();
667     Standard_Real U1, V1, U , V;
668     ElSLib::Parameters( TR, P1, U1, V1);
669     Standard_Real Step = .1, DeltaU = 0., DeltaV = 0.;
670     Standard_Real eps = M_PI, dmaxU = 0., dU = 0., dmaxV = 0., dV = 0.;
671     Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
672     nbp = Max(nbp, 3);
673     Step = (W2 - W1) / (nbp - 1);
674     myU1 = U1; myU2 = U1;
675     myV1 = V1; myV2 = V1;
676     Standard_Real pminU = W1, pmaxU = W1, pminV = W1, pmaxV = W1,
677       plim = W2+.1*Step;
678     for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
679       P = myCurve->Value(par);
680       ElSLib::Parameters( TR, P, U, V);
681       U += DeltaU;
682       V += DeltaV;
683       dU = U - U1;
684       dV = V - V1;
685       if(dU > eps) {
686         U -= DeltaU;
687         DeltaU -= 2*M_PI;
688         U += DeltaU;
689         dU = U - U1;
690       }
691       else if(dU < -eps) {
692         U -= DeltaU;
693         DeltaU += 2*M_PI;
694         U += DeltaU;
695         dU = U - U1;
696       }
697       if(dV > eps) {
698         V -= DeltaV;
699         DeltaV -= 2*M_PI;
700         V += DeltaV;
701         dV = V - V1;
702       }
703       else if(dV < -eps) {
704         V -= DeltaV;
705         DeltaV += 2*M_PI;
706         V += DeltaV;
707         dV = V - V1;
708       }
709       dmaxU = Max(dmaxU, Abs(dU));
710       dmaxV = Max(dmaxV, Abs(dV));
711       if(U < myU1) {myU1 = U; pminU = par;}
712       if(U > myU2) {myU2 = U; pmaxU = par;}
713       if(V < myV1) {myV1 = V; pminV = par;}
714       if(V > myV2) {myV2 = V; pmaxV = par;}
715       U1 = U;
716       V1 = V;
717     }
718
719     if(!(Abs(pminU - W1) <= Precision::PConfusion() ||
720       Abs(pminU - W2) <= Precision::PConfusion())  ) myU1 -= dmaxU*.5;
721     if(!(Abs(pmaxU - W1) <= Precision::PConfusion() ||
722       Abs(pmaxU - W2) <= Precision::PConfusion())  ) myU2 += dmaxU*.5;
723     if(!(Abs(pminV - W1) <= Precision::PConfusion() ||
724       Abs(pminV - W2) <= Precision::PConfusion())  ) myV1 -= dmaxV*.5;
725     if(!(Abs(pmaxV - W1) <= Precision::PConfusion() ||
726       Abs(pmaxV - W2) <= Precision::PConfusion())  ) myV2 += dmaxV*.5;
727
728     if((myU1 >=0. && myU1 <= 2*M_PI) &&
729       (myU2 >=0. && myU2 <= 2*M_PI)    ) {
730         myU1 = 0.;
731         myU2 = 2.*M_PI;
732         UCouture = Standard_False;
733     }
734     else {
735       U = ( myU1 + myU2 ) /2.;
736       myU1 = U - M_PI;
737       myU2 = U + M_PI;
738       UCouture = Standard_True;
739     }
740     if((myV1 >=0. && myV1 <= 2*M_PI) &&
741       (myV2 >=0. && myV2 <= 2*M_PI)    ) {
742         VCouture = Standard_False;
743     }
744     else {
745       V = ( myV1 + myV2 ) /2.;
746       myV1 = V - M_PI;
747       myV2 = V + M_PI;
748       VCouture = Standard_True;
749     }
750
751   }
752   break;
753
754   default:
755     {
756       UCouture = Standard_False;
757       VCouture = Standard_False;
758     }
759     break;
760   }
761 }
762 //
763 //
764 //=======================================================================
765 //classn : ProjLib_Function
766 //purpose  : 
767 //=======================================================================
768 class ProjLib_Function : public AppCont_Function2d
769 {
770   Handle(Adaptor3d_HCurve)   myCurve;
771   Handle(Adaptor3d_HSurface) mySurface;
772
773   public :
774
775   Standard_Real    myU1,myU2,myV1,myV2;
776   Standard_Boolean UCouture,VCouture;
777   
778   ProjLib_Function(const Handle(Adaptor3d_HCurve)&   C, 
779                    const Handle(Adaptor3d_HSurface)& S) :
780   myCurve(C), mySurface(S),
781   myU1(0.0),
782   myU2(0.0),
783   myV1(0.0),
784   myV2(0.0),
785   UCouture(Standard_False),
786   VCouture(Standard_False)
787     {Function_SetUVBounds(myU1,myU2,myV1,myV2,UCouture,VCouture,myCurve,mySurface);}
788   
789   Standard_Real FirstParameter() const
790     {return (myCurve->FirstParameter() + 1.e-9);}
791   
792   Standard_Real LastParameter() const
793     {return (myCurve->LastParameter() -1.e-9);}
794   
795   
796   gp_Pnt2d Value(const Standard_Real t) const
797     {return Function_Value(t,myCurve,mySurface,myU1,myU2,myV1,myV2,UCouture,VCouture);}
798   
799   Standard_Boolean D1(const Standard_Real t, gp_Pnt2d& P, gp_Vec2d& V) const
800     {return Function_D1(t,P,V,myCurve,mySurface,myU1,myU2,myV1,myV2,UCouture,VCouture);}
801 };
802
803 //=======================================================================
804 //function : ProjLib_ComputeApprox
805 //purpose  : 
806 //=======================================================================
807
808 ProjLib_ComputeApprox::ProjLib_ComputeApprox
809   (const Handle(Adaptor3d_HCurve)   & C,
810    const Handle(Adaptor3d_HSurface) & S,
811    const Standard_Real              Tol )
812 {
813   // if the surface is a plane and the curve a BSpline or a BezierCurve,
814   // don`t make an Approx but only the projection of the poles.
815
816   myTolerance = Max(Precision::PApproximation(),Tol);
817   Standard_Integer NbKnots, NbPoles ;
818   GeomAbs_CurveType   CType = C->GetType();
819   GeomAbs_SurfaceType SType = S->GetType();
820
821   Standard_Boolean SurfIsAnal = (SType != GeomAbs_BSplineSurface) &&
822                                 (SType != GeomAbs_BezierSurface)  &&
823                                 (SType != GeomAbs_OtherSurface)     ;
824
825   Standard_Boolean CurvIsAnal = (CType != GeomAbs_BSplineCurve) &&
826                                 (CType != GeomAbs_BezierCurve)  &&
827                                 (CType != GeomAbs_OtherCurve)     ;
828
829   Standard_Boolean simplecase = SurfIsAnal && CurvIsAnal;
830
831   if (CType == GeomAbs_BSplineCurve &&
832       SType == GeomAbs_Plane ) {
833     
834     // get the poles and eventually the weights
835     Handle(Geom_BSplineCurve) BS = C->BSpline();
836     NbPoles = BS->NbPoles();
837     TColgp_Array1OfPnt P3d( 1, NbPoles);
838     TColgp_Array1OfPnt2d Poles( 1, NbPoles);
839     TColStd_Array1OfReal Weights( 1, NbPoles);
840     if ( BS->IsRational()) BS->Weights(Weights);
841     BS->Poles( P3d);
842     gp_Pln Plane = S->Plane();
843     Standard_Real U,V;
844     for ( Standard_Integer i = 1; i <= NbPoles; i++) {
845       ElSLib::Parameters( Plane, P3d(i), U, V);
846       Poles.SetValue(i,gp_Pnt2d(U,V));
847     }
848     NbKnots = BS->NbKnots();
849     TColStd_Array1OfReal     Knots(1,NbKnots);
850     TColStd_Array1OfInteger  Mults(1,NbKnots);
851     BS->Knots(Knots) ;
852     BS->Multiplicities(Mults) ; 
853     // get the knots and mults if BSplineCurve
854     if ( BS->IsRational()) {
855       myBSpline = new Geom2d_BSplineCurve(Poles,
856                                           Weights,
857                                           Knots,
858                                           Mults,
859                                           BS->Degree(),
860                                           BS->IsPeriodic());
861     }
862     else {
863       myBSpline = new Geom2d_BSplineCurve(Poles,
864                                           Knots,
865                                           Mults,
866                                           BS->Degree(),
867                                           BS->IsPeriodic());
868     }
869   }
870   else if (CType == GeomAbs_BezierCurve &&
871            SType == GeomAbs_Plane ) {
872     
873     // get the poles and eventually the weights
874     Handle(Geom_BezierCurve) BezierCurvePtr = C->Bezier() ;
875     NbPoles = BezierCurvePtr->NbPoles();
876     TColgp_Array1OfPnt P3d( 1, NbPoles);
877     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
878     TColStd_Array1OfReal   Weights( 1, NbPoles);
879     if ( BezierCurvePtr->IsRational()) {
880       BezierCurvePtr->Weights(Weights);
881     }
882     BezierCurvePtr->Poles( P3d);  
883     
884     // project the 3D-Poles on the plane
885     
886     gp_Pln Plane = S->Plane();
887     Standard_Real U,V;
888     for ( Standard_Integer i = 1; i <= NbPoles; i++) {
889       ElSLib::Parameters( Plane, P3d(i), U, V);
890       Poles.SetValue(i,gp_Pnt2d(U,V));
891     }
892     if (  BezierCurvePtr->IsRational()) {
893       myBezier =  new Geom2d_BezierCurve(Poles, Weights);
894     }
895     else {
896       myBezier =  new Geom2d_BezierCurve(Poles);
897     }
898   }
899   else {
900     ProjLib_Function F( C, S);
901
902 #ifdef DEB
903     if ( AffichValue) {
904       Standard_Integer Nb = 20;
905       Standard_Real U1, U2, dU, U;
906       U1 = F.FirstParameter();
907       U2 = F.LastParameter();
908       dU = ( U2 - U1) / Nb;
909       TColStd_Array1OfInteger Mults(1,Nb+1);
910       TColStd_Array1OfReal    Knots(1,Nb+1);
911       TColgp_Array1OfPnt2d    Poles(1,Nb+1);
912       for ( Standard_Integer i = 1; i <= Nb+1; i++) {
913               U = U1 + (i-1)*dU;
914               Poles(i) = F.Value(U);
915               Knots(i) = i;
916               Mults(i) = 1;
917       }
918       Mults(1)    = 2;
919       Mults(Nb+1) = 2;
920 #ifdef DRAW
921 // POP pour NT
922       char* ResultName = "Result";
923       DrawTrSurf::Set(ResultName,new Geom2d_BSplineCurve(Poles,Knots,Mults,1));
924 //      DrawTrSurf::Set("Result",new Geom2d_BSplineCurve(Poles,Knots,Mults,1));
925 #endif
926     }
927 #endif    
928
929 //-----------
930     Standard_Integer Deg1, Deg2;
931     if(simplecase) {
932       Deg1 = 8; 
933       Deg2 = 10; 
934     }
935     else {
936       Deg1 = 8; 
937       Deg2 = 12; 
938     }
939 //-------------
940     Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
941                               Standard_True);
942     if(Fit.IsAllApproximated()) {
943       Standard_Integer i;
944       Standard_Integer NbCurves = Fit.NbMultiCurves();
945     
946     // on essaie de rendre la courbe au moins C1
947       Convert_CompBezierCurves2dToBSplineCurve2d Conv;
948
949       myTolerance = 0;
950       Standard_Real Tol3d,Tol2d;
951       for (i = 1; i <= NbCurves; i++) {
952               Fit.Error(i,Tol3d, Tol2d);
953               myTolerance = Max(myTolerance, Tol2d);
954               AppParCurves_MultiCurve MC = Fit.Value( i);       //Charge la Ieme Curve
955               TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Recupere les poles
956               MC.Curve(1, Poles2d);
957             
958               Conv.AddCurve(Poles2d);
959       }
960     
961     //mise a jour des fields de ProjLib_Approx
962       Conv.Perform();
963     
964       NbPoles    = Conv.NbPoles();
965       NbKnots    = Conv.NbKnots();
966
967       //7626
968       if(NbPoles <= 0 || NbPoles > 100000)
969               return;
970       if(NbKnots <= 0 || NbKnots > 100000)
971               return;
972
973       TColgp_Array1OfPnt2d    NewPoles(1,NbPoles);
974       TColStd_Array1OfReal    NewKnots(1,NbKnots);
975       TColStd_Array1OfInteger NewMults(1,NbKnots);
976     
977       Conv.KnotsAndMults(NewKnots,NewMults);
978       Conv.Poles(NewPoles);
979     
980       BSplCLib::Reparametrize(C->FirstParameter(),
981                               C->LastParameter(),
982                               NewKnots);
983     
984       // il faut recadrer les poles de debut et de fin:
985       // ( Car pour les problemes de couture, on a du ouvrir l`intervalle
986       // de definition de la courbe.)
987       // On choisit de calculer ces poles par prolongement de la courbe
988       // approximee.
989     
990       gp_Pnt2d P;
991       Standard_Real U;
992     
993       U = C->FirstParameter() - 1.e-9;
994       BSplCLib::D0(U, 
995                    0, 
996                    Conv.Degree(), 
997                    Standard_False,
998                    NewPoles, 
999                    BSplCLib::NoWeights(), 
1000                    NewKnots, 
1001                    NewMults,
1002                    P);
1003       NewPoles.SetValue(1,P);
1004       U = C->LastParameter() + 1.e-9;
1005       BSplCLib::D0(U, 
1006                    0, 
1007                    Conv.Degree(), 
1008                    Standard_False,
1009                    NewPoles, 
1010                    BSplCLib::NoWeights(), 
1011                    NewKnots, 
1012                    NewMults,
1013                    P);
1014       NewPoles.SetValue(NbPoles,P);
1015       myBSpline = new Geom2d_BSplineCurve (NewPoles,
1016                                            NewKnots,
1017                                            NewMults,
1018                                            Conv.Degree());
1019     }
1020     else {
1021       Standard_Integer NbCurves = Fit.NbMultiCurves();
1022       if(NbCurves != 0) {
1023               Standard_Real Tol3d,Tol2d;
1024               Fit.Error(NbCurves,Tol3d, Tol2d);
1025               myTolerance = Tol2d;
1026       }
1027     }
1028
1029     //Return curve home
1030     Standard_Real UFirst = F.FirstParameter();
1031     gp_Pnt P3d = C->Value( UFirst );
1032     Standard_Real u = 0., v = 0.;
1033     switch (SType)
1034     {
1035     case GeomAbs_Plane:
1036       {
1037         gp_Pln Plane = S->Plane();
1038         ElSLib::Parameters( Plane, P3d, u, v );
1039         break;
1040       }
1041     case GeomAbs_Cylinder:
1042       {
1043         gp_Cylinder Cylinder = S->Cylinder();
1044         ElSLib::Parameters( Cylinder, P3d, u, v );
1045         break;
1046       }
1047     case GeomAbs_Cone:
1048       {
1049         gp_Cone Cone = S->Cone();
1050         ElSLib::Parameters( Cone, P3d, u, v );
1051         break;
1052       }
1053     case GeomAbs_Sphere:
1054       {
1055         gp_Sphere Sphere = S->Sphere();
1056         ElSLib::Parameters( Sphere, P3d, u, v );
1057         break;
1058       }
1059     case GeomAbs_Torus:
1060       {
1061         gp_Torus Torus = S->Torus();
1062         ElSLib::Parameters( Torus, P3d, u, v );
1063         break;
1064       }
1065     default:
1066       Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
1067     }
1068     Standard_Boolean ToMirror = Standard_False;
1069     Standard_Real du = 0., dv = 0.;
1070     Standard_Integer number;
1071     if (F.VCouture)
1072     { 
1073       if (SType == GeomAbs_Sphere && Abs(u-F.myU1) > M_PI)
1074       {
1075         ToMirror = Standard_True;
1076         dv = -M_PI;
1077         v = M_PI - v;
1078       }
1079       Standard_Real newV = ElCLib::InPeriod( v, F.myV1, F.myV2 );
1080       number = (Standard_Integer) (Floor((newV-v)/(F.myV2-F.myV1)));
1081       dv -= number*(F.myV2-F.myV1);
1082     }
1083     if (F.UCouture || (F.VCouture && SType == GeomAbs_Sphere))
1084     {
1085       gp_Pnt2d P2d = F.Value( UFirst );
1086       number = (Standard_Integer) (Floor((P2d.X()-u)/M_PI + Epsilon(M_PI)));
1087       du = -number*M_PI;
1088     }
1089
1090     if (!myBSpline.IsNull())
1091     {
1092       if (du != 0. || dv != 0.)
1093         myBSpline->Translate( gp_Vec2d(du,dv) );
1094       if (ToMirror)
1095       {
1096         gp_Ax2d Axe( gp_Pnt2d(0.,0.), gp_Dir2d(1.,0.) );
1097         myBSpline->Mirror( Axe );
1098       }
1099     }
1100   }
1101 }
1102
1103 //=======================================================================
1104 //function : BSpline
1105 //purpose  : 
1106 //=======================================================================
1107
1108 Handle(Geom2d_BSplineCurve)  ProjLib_ComputeApprox::BSpline() const 
1109
1110 {
1111   return myBSpline ;
1112 }
1113
1114 //=======================================================================
1115 //function : Bezier
1116 //purpose  : 
1117 //=======================================================================
1118
1119 Handle(Geom2d_BezierCurve)  ProjLib_ComputeApprox::Bezier()  const 
1120
1121 {
1122   return myBezier ;
1123 }
1124
1125
1126 //=======================================================================
1127 //function : Tolerance
1128 //purpose  : 
1129 //=======================================================================
1130
1131 Standard_Real ProjLib_ComputeApprox::Tolerance() const 
1132 {
1133   return myTolerance;
1134 }
1135
1136