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