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