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