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