b4ac5189fa8f7c38ab2a7a3778f91ee4611a825e
[occt.git] / src / ProjLib / ProjLib_ProjectedCurve.cxx
1 // Created on: 1993-08-25
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
22
23 //  Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272
24
25 #include <GeomAbs_SurfaceType.hxx>
26 #include <Standard_NoSuchObject.hxx>
27 #include <Standard_NotImplemented.hxx>
28 #include <ProjLib_ProjectedCurve.hxx>
29 #include <ProjLib_CompProjectedCurve.hxx>
30 #include <ProjLib_HCompProjectedCurve.hxx>
31 #include <ProjLib_ComputeApproxOnPolarSurface.hxx>
32 #include <ProjLib_ComputeApprox.hxx>
33 #include <ProjLib_Projector.hxx>
34 #include <Handle_Adaptor3d_HCurve.hxx>
35 #include <Handle_Adaptor3d_HSurface.hxx>
36 #include <Adaptor3d_HCurve.hxx>
37 #include <Adaptor3d_HSurface.hxx>
38 #include <Approx_CurveOnSurface.hxx>
39 #include <ProjLib_Plane.hxx>
40 #include <ProjLib_Cylinder.hxx>
41 #include <ProjLib_Cone.hxx>
42 #include <ProjLib_Sphere.hxx>
43 #include <ProjLib_Torus.hxx>
44 #include <Precision.hxx>
45 #include <Handle_Geom_BSplineCurve.hxx>
46 #include <Geom2d_BSplineCurve.hxx>
47 #include <Handle_Geom2d_BSplineCurve.hxx>
48 #include <Geom2d_BezierCurve.hxx>
49 #include <Handle_Geom2d_BezierCurve.hxx>
50 #include <Handle_Adaptor2d_HCurve2d.hxx>
51 #include <gp_Vec2d.hxx>
52 #include <StdFail_NotDone.hxx>
53 #include <gp_XY.hxx>
54 #include <TColgp_HArray1OfPnt2d.hxx>
55 #include <TColStd_HArray1OfReal.hxx>
56 #include <Geom2dConvert_CompCurveToBSplineCurve.hxx>
57 #include <TColStd_Array1OfReal.hxx>
58 #include <TColStd_Array1OfInteger.hxx>
59 #include <TColgp_Array1OfPnt2d.hxx>
60 #include <TColgp_HArray1OfVec2d.hxx>
61 #include <TColStd_HArray1OfBoolean.hxx>
62 #include <BSplCLib.hxx>
63 #include <GeomAbs_IsoType.hxx>
64
65 //=======================================================================
66 //function : IsoIsDeg
67 //purpose  : 
68 //=======================================================================
69
70 static Standard_Boolean IsoIsDeg  (const Adaptor3d_Surface& S,
71                                    const Standard_Real      Param,
72                                    const GeomAbs_IsoType    IT,
73                                    const Standard_Real      TolMin,
74                                    const Standard_Real      TolMax) 
75 {
76     Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T;
77     Standard_Boolean Along = Standard_True;
78     U1 = S.FirstUParameter();
79     U2 = S.LastUParameter();
80     V1 = S.FirstVParameter();
81     V2 = S.LastVParameter();
82     gp_Vec D1U,D1V;
83     gp_Pnt P;
84     Standard_Real Step,D1NormMax;
85     if (IT == GeomAbs_IsoV) 
86     {
87       Step = (U2 - U1)/10;
88       D1NormMax=0.;
89       for (T=U1;T<=U2;T=T+Step) 
90       {
91         S.D1(T,Param,P,D1U,D1V);
92         D1NormMax=Max(D1NormMax,D1U.Magnitude());
93       }
94
95       if (D1NormMax >TolMax || D1NormMax < TolMin ) 
96            Along = Standard_False;
97     }
98     else 
99     {
100       Step = (V2 - V1)/10;
101       D1NormMax=0.;
102       for (T=V1;T<=V2;T=T+Step) 
103       {
104         S.D1(Param,T,P,D1U,D1V);
105         D1NormMax=Max(D1NormMax,D1V.Magnitude());
106       }
107
108       if (D1NormMax >TolMax || D1NormMax < TolMin ) 
109            Along = Standard_False;
110
111
112     }
113     return Along;
114 }
115
116 //=======================================================================
117 //function : Interpolate
118 //purpose  : 
119 //=======================================================================
120
121 static Handle(Geom2d_BSplineCurve) Interpolate(const Handle(TColgp_HArray1OfPnt2d)& myPoints,
122                                                const Handle(TColStd_HArray1OfReal)& myParameters,
123                                                const gp_Vec2d& InitialTangent, 
124                                                const gp_Vec2d& FinalTangent)
125 {
126   Handle(Geom2d_BSplineCurve) myCurve = NULL;
127
128 // This code is extraction from Geom2dAPI_Interpolate with small correction
129 // This is done to avoid of cyclic dependency if Geom2dAPI is used in ProjLib 
130
131   Standard_Integer degree,
132   ii,
133   jj,
134   index,
135   index1,
136   index2,
137   index3,
138   mult_index,
139   inversion_problem,
140   num_points,
141   num_distinct_knots,
142   num_poles  ;
143   
144   gp_Pnt2d a_point ;
145
146   Standard_Boolean myTangentRequest = Standard_True;
147
148   Handle(TColgp_HArray1OfVec2d) myTangents = 
149      new TColgp_HArray1OfVec2d(myPoints->Lower(),
150                               myPoints->Upper()) ;
151   Handle(TColStd_HArray1OfBoolean) myTangentFlags =
152       new TColStd_HArray1OfBoolean(myPoints->Lower(),
153                                    myPoints->Upper()) ;
154   myTangentFlags->Init(Standard_False);
155   
156   myTangentFlags->SetValue(1,Standard_True) ;
157   myTangentFlags->SetValue(myPoints->Length(),Standard_True) ;
158   myTangents->SetValue(1,InitialTangent) ;
159   myTangents->SetValue(myPoints->Length(),FinalTangent);
160
161   num_points =
162   num_distinct_knots =
163   num_poles = myPoints->Length() ;
164   if (num_poles == 2 &&   !myTangentRequest)  {
165     degree = 1 ;
166   } 
167   else if (num_poles == 3 && !myTangentRequest) {
168     degree = 2 ;
169     num_distinct_knots = 2 ;
170   }
171   else {
172     degree = 3 ;
173     num_poles += 2 ;
174     //if (myTangentRequest) 
175       //for (ii = myTangentFlags->Lower() + 1 ; 
176         //   ii < myTangentFlags->Upper() ; ii++) {
177         //if (myTangentFlags->Value(ii)) {
178           //num_poles += 1 ;
179         //}
180       //}
181     }
182   
183   
184   TColStd_Array1OfReal     parameters(1,num_poles) ;  
185   TColStd_Array1OfReal     flatknots(1,num_poles + degree + 1) ;
186   TColStd_Array1OfInteger  mults(1,num_distinct_knots) ;
187   TColStd_Array1OfReal     knots(1,num_distinct_knots) ;
188   TColStd_Array1OfInteger  contact_order_array(1, num_poles) ;
189   TColgp_Array1OfPnt2d       poles(1,num_poles) ;
190
191   for (ii = 1 ; ii <= degree + 1 ; ii++) {
192     flatknots.SetValue(ii,myParameters->Value(1)) ;
193     flatknots.SetValue(ii + num_poles,
194                        myParameters->Value(num_points)) ;
195   }
196   for (ii = 1 ; ii <= num_poles ; ii++) {
197     contact_order_array.SetValue(ii,0)  ;
198   }
199   for (ii = 2 ; ii < num_distinct_knots ; ii++) {
200     mults.SetValue(ii,1) ; 
201   }
202   mults.SetValue(1,degree + 1) ;
203   mults.SetValue(num_distinct_knots ,degree + 1) ;
204   
205   switch (degree) {
206   case 1:
207     for (ii = 1 ; ii <= num_poles ; ii++) {
208       poles.SetValue(ii ,myPoints->Value(ii)) ;
209     }
210     myCurve =
211       new Geom2d_BSplineCurve(poles,
212                             myParameters->Array1(),
213                             mults,
214                             degree) ;
215     //myIsDone = Standard_True ;
216     break ;
217   case 2:
218     knots.SetValue(1,myParameters->Value(1)) ;
219     knots.SetValue(2,myParameters->Value(3)) ;
220     for (ii = 1 ; ii <= num_poles ; ii++) {
221       poles.SetValue(ii,myPoints->Value(ii)) ;
222       
223     }
224     BSplCLib::Interpolate(degree,
225                           flatknots,
226                           myParameters->Array1(),
227                           contact_order_array,
228                           poles,
229                           inversion_problem) ;
230     if (!inversion_problem) {
231       myCurve =
232         new Geom2d_BSplineCurve(poles,
233                               knots,
234                               mults,
235                               degree) ;
236       //myIsDone = Standard_True ;
237     }
238     break ;
239   case 3:
240 //
241 // check if the boundary conditions are set
242 //
243     //if (num_points >= 3) {
244 //
245 // cannot build the tangents with degree 3 with only 2 points
246 // if those where not given in advance 
247 //
248       //BuildTangents(myPoints->Array1(),
249                     //myTangents->ChangeArray1(),
250                     //myTangentFlags->ChangeArray1(),
251                     //myParameters->Array1()) ;
252     //}
253     contact_order_array.SetValue(2,1)  ;
254     parameters.SetValue(1,myParameters->Value(1)) ;
255     parameters.SetValue(2,myParameters->Value(1)) ;
256     poles.SetValue(1,myPoints->Value(1)) ;
257     for (jj = 1 ; jj <= 2 ; jj++) {
258       a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
259
260     }
261     poles.SetValue(2,a_point) ;
262     mult_index = 2 ;
263     index = 3 ;
264     index1 = 2 ;
265     index2 = myPoints->Lower() + 1 ;
266     index3 = degree + 2 ;
267     if (myTangentRequest) {
268       for (ii = myParameters->Lower() + 1 ; 
269            ii < myParameters->Upper() ; ii++) {
270         parameters.SetValue(index,myParameters->Value(ii)) ;
271         poles.SetValue(index,myPoints->Value(index2)) ;
272         flatknots.SetValue(index3,myParameters->Value(ii)) ; 
273         index += 1 ;
274         index3 += 1 ;
275         if (myTangentFlags->Value(index1)) {
276 //
277 // set the multiplicities, the order of the contact, the 
278 // the flatknots,
279 //
280           mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
281           contact_order_array(index) = 1 ;
282           flatknots.SetValue(index3, myParameters->Value(ii)) ;
283           parameters.SetValue(index,
284                               myParameters->Value(ii)) ;
285           for (jj = 1 ; jj <= 2 ; jj++) {
286             a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
287           }
288           poles.SetValue(index,a_point) ;
289           index += 1  ;
290           index3 += 1 ;
291         }
292         mult_index += 1 ;
293         index1 += 1 ;
294         index2 += 1 ;
295
296       }
297     }
298     else {
299       index1 = 2 ;
300       for(ii = myParameters->Lower()  ; ii <= myParameters->Upper()  ; ii++) {
301         parameters.SetValue(index1, 
302                             myParameters->Value(ii)) ;
303         index1 += 1 ;
304       }
305       index = 3 ;
306       for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() - 1 ; ii++) {
307         poles.SetValue(index,
308                        myPoints->Value(ii)) ;
309         index += 1 ;
310       }
311    
312    
313       index = degree + 1 ;
314       for(ii = myParameters->Lower()  ; ii <= myParameters->Upper()  ; ii++) {
315         flatknots.SetValue(index,
316                            myParameters->Value(ii)) ;
317         index += 1 ;
318       }
319     }
320     for (jj = 1 ; jj <= 2 ; jj++) {
321       a_point.SetCoord(jj,
322                        myTangents->Value(num_points).Coord(jj)) ;
323     }
324     poles.SetValue(num_poles-1 ,a_point) ;
325
326     contact_order_array.SetValue(num_poles - 1,1) ;
327     parameters.SetValue(num_poles,
328                         myParameters->Value(myParameters->Upper())) ;
329     parameters.SetValue(num_poles -1,
330                         myParameters->Value(myParameters->Upper())) ;
331
332     poles.SetValue(num_poles,
333                    myPoints->Value(num_points)) ;
334
335     BSplCLib::Interpolate(degree,
336                           flatknots,
337                           parameters,
338                           contact_order_array,
339                           poles,
340                           inversion_problem) ;
341     if (!inversion_problem) {
342       myCurve =
343         new Geom2d_BSplineCurve(poles,
344                               myParameters->Array1(),
345                               mults,
346                               degree) ;
347       //myIsDone = Standard_True ;
348     }
349     break ;
350  
351   }
352
353
354   return myCurve;
355 }
356
357 //=======================================================================
358 //function : TrimC3d
359 //purpose  : 
360 //=======================================================================
361
362 static void TrimC3d(Handle(Adaptor3d_HCurve)& myCurve,
363                     Standard_Boolean* IsTrimmed,
364                     const Standard_Real dt,
365                     const gp_Pnt& Pole)
366 {
367   Standard_Real f = myCurve->FirstParameter();
368   Standard_Real l = myCurve->LastParameter();
369
370   gp_Pnt P = myCurve->Value(f);
371
372   if(P.Distance(Pole) < Precision::Confusion()) {
373     IsTrimmed[0] = Standard_True;
374     f = f+dt;
375     myCurve = myCurve->Trim(f, l, Precision::Confusion());
376   }
377   
378   P = myCurve->Value(l);
379   if(P.Distance(Pole) < Precision::Confusion()) {
380     IsTrimmed[1] = Standard_True;
381     l = l-dt;
382     myCurve = myCurve->Trim(f, l, Precision::Confusion());
383   }
384 }
385
386 //=======================================================================
387 //function : ExtendC2d
388 //purpose  : 
389 //=======================================================================
390
391 static void ExtendC2d(Handle(Geom2d_BSplineCurve)& aRes,
392                       const Standard_Real t,
393                       const Standard_Real dt,
394                       const Standard_Real u1,
395                       const Standard_Real u2,
396                       const Standard_Real v1,
397                       const Standard_Real v2)
398 {
399   gp_Pnt2d P0;
400   gp_Vec2d V01, V02;
401   aRes->D2(t, P0, V01, V02);
402   
403   gp_XY XYP1 = P0.XY() + V01.XY()*dt + .5*V02.XY()*dt*dt;
404   
405   gp_Vec2d V11 = V01 + V02*dt;
406   
407   if(XYP1.X() < u1) XYP1.SetX(u1);
408   if(XYP1.X() > u2) XYP1.SetX(u2);
409   if(XYP1.Y() < v1) XYP1.SetY(v1);
410   if(XYP1.Y() > v2) XYP1.SetY(v2);
411   
412   Handle(TColgp_HArray1OfPnt2d) aPnts = new TColgp_HArray1OfPnt2d(1, 2);
413   Handle(TColStd_HArray1OfReal) aPars = new TColStd_HArray1OfReal(1, 2);
414   
415   if(dt < 0.) {
416     aPnts->SetValue(1, gp_Pnt2d(XYP1));
417     aPnts->SetValue(2, P0);
418     aPars->SetValue(1, t + dt);
419     aPars->SetValue(2, t);
420   }
421   else {
422     aPnts->SetValue(2, gp_Pnt2d(XYP1));
423     aPnts->SetValue(1, P0);
424     aPars->SetValue(2, t + dt);
425     aPars->SetValue(1, t);
426   }
427   
428   Handle(Geom2d_BSplineCurve) aC;  
429   
430   if(dt < 0.) {
431     aC = Interpolate(aPnts, aPars, V11, V01); 
432   }
433   else {
434     aC = Interpolate(aPnts, aPars, V01, V11); 
435   }
436   
437   
438   Geom2dConvert_CompCurveToBSplineCurve aConcat(aRes);
439   aConcat.Add(aC, Precision::PConfusion());
440   
441   aRes = aConcat.BSplineCurve();
442   
443 }  
444
445 //=======================================================================
446 //function : Project
447 //purpose  : 
448 //=======================================================================
449
450 static void Project(ProjLib_Projector& P, Handle(Adaptor3d_HCurve)& C)
451 {
452   GeomAbs_CurveType CType = C->GetType();
453   switch (CType) {
454     case GeomAbs_Line:
455       P.Project(C->Line());
456       break;
457     case GeomAbs_Circle:
458       P.Project(C->Circle());
459       break;
460     case GeomAbs_Ellipse:
461       P.Project(C->Ellipse());
462       break;
463     case GeomAbs_Hyperbola:
464       P.Project(C->Hyperbola());
465       break;
466     case GeomAbs_Parabola:
467       P.Project(C->Parabola());
468       break;
469     case GeomAbs_BSplineCurve:
470     case GeomAbs_BezierCurve:
471     case GeomAbs_OtherCurve:    // try the approximation
472       break;
473     default:
474       Standard_NoSuchObject::Raise(" ");
475   }
476 }
477
478 //=======================================================================
479 //function : ProjLib_ProjectedCurve
480 //purpose  : 
481 //=======================================================================
482
483 ProjLib_ProjectedCurve::ProjLib_ProjectedCurve()
484
485 {
486   myTolerance = Precision::PApproximation();
487 }
488
489
490 //=======================================================================
491 //function : ProjLib_ProjectedCurve
492 //purpose  : 
493 //=======================================================================
494
495 ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
496 (const Handle(Adaptor3d_HSurface)& S)
497 {
498   myTolerance = Precision::PApproximation();
499   Load(S);
500 }
501
502
503 //=======================================================================
504 //function : ProjLib_ProjectedCurve
505 //purpose  : 
506 //=======================================================================
507
508 ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
509 (const Handle(Adaptor3d_HSurface)& S,
510  const Handle(Adaptor3d_HCurve)& C)
511 {
512   myTolerance = Precision::PApproximation();
513   Load(S);
514   Load(C);
515 }
516
517
518 //=======================================================================
519 //function : ProjLib_ProjectedCurve
520 //purpose  : 
521 //=======================================================================
522
523 ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
524 (const Handle(Adaptor3d_HSurface)& S,
525  const Handle(Adaptor3d_HCurve)&   C,
526  const Standard_Real             Tol)
527 {
528   myTolerance = Max(Tol, Precision::PApproximation());
529   Load(S);
530   Load(C);
531 }
532
533
534 //=======================================================================
535 //function : Load
536 //purpose  : 
537 //=======================================================================
538
539 void ProjLib_ProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
540 {
541   mySurface = S ;
542 }
543
544
545 //=======================================================================
546 //function : Load
547 //purpose  : 
548 //=======================================================================
549
550 void ProjLib_ProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
551 {
552   myTolerance = Max(myTolerance,Precision::PApproximation());
553   myCurve = C;
554   GeomAbs_SurfaceType SType = mySurface->GetType();    
555   GeomAbs_CurveType   CType = myCurve->GetType();
556
557   switch (SType) {
558
559     case GeomAbs_Plane:
560       {
561         ProjLib_Plane P(mySurface->Plane());
562         Project(P,myCurve);
563         myResult = P;
564       }
565       break;
566
567     case GeomAbs_Cylinder:
568       {
569         ProjLib_Cylinder P(mySurface->Cylinder());
570         Project(P,myCurve);
571         myResult = P;
572       }
573       break;
574
575     case GeomAbs_Cone:
576       {
577         ProjLib_Cone P(mySurface->Cone());
578         Project(P,myCurve);
579         myResult = P;
580       }
581       break;
582
583     case GeomAbs_Sphere:
584       {
585         ProjLib_Sphere P(mySurface->Sphere());
586         Project(P,myCurve);
587         if ( P.IsDone()) { 
588           // on met dans la pseudo-periode ( car Sphere n'est pas
589           // periodique en V !)
590           P.SetInBounds(myCurve->FirstParameter());
591         }
592         myResult = P;
593       }
594       break;
595
596     case GeomAbs_Torus:
597       {
598         ProjLib_Torus P(mySurface->Torus());
599         Project(P,myCurve);
600         myResult = P;
601       }
602       break;
603
604     case GeomAbs_BezierSurface:
605     case GeomAbs_BSplineSurface:
606       {
607         
608         Standard_Boolean IsTrimmed[2] = {Standard_False, Standard_False};
609         Standard_Real f, l, dt;
610         const Standard_Real eps = 0.01;
611         f = myCurve->FirstParameter();
612         l = myCurve->LastParameter();
613         dt = (l-f)*eps;
614
615         Standard_Real U1=0.,U2=0.,V1=0.,V2=0;
616         const Adaptor3d_Surface& S = mySurface->Surface();
617         U1 = S.FirstUParameter();
618         U2 = S.LastUParameter();
619         V1 = S.FirstVParameter();
620         V2 = S.LastVParameter();
621
622         if(IsoIsDeg(S, U1, GeomAbs_IsoU, 0., myTolerance) ) {
623           //Surface has pole at U = Umin
624           gp_Pnt Pole = mySurface->Value(U1, V1);
625           TrimC3d(myCurve, IsTrimmed, dt, Pole);
626         }
627
628         if(IsoIsDeg(S, U2, GeomAbs_IsoU, 0., myTolerance) ) {
629           //Surface has pole at U = Umax
630           gp_Pnt Pole = mySurface->Value(U2, V1);
631           TrimC3d(myCurve, IsTrimmed, dt, Pole);
632         }
633           
634         if(IsoIsDeg(S, V1, GeomAbs_IsoV, 0., myTolerance) ) {
635           //Surface has pole at V = Vmin
636           gp_Pnt Pole = mySurface->Value(U1, V1);
637           TrimC3d(myCurve, IsTrimmed, dt, Pole);
638         }
639
640         if(IsoIsDeg(S, V2, GeomAbs_IsoV, 0., myTolerance) ) {
641           //Surface has pole at V = Vmax
642           gp_Pnt Pole = mySurface->Value(U1, V2);
643           TrimC3d(myCurve, IsTrimmed, dt, Pole);
644         }
645
646         ProjLib_ComputeApproxOnPolarSurface polar(myCurve, 
647                                                   mySurface,
648                                                   myTolerance);
649
650         Handle(Geom2d_BSplineCurve) aRes = polar.BSpline();
651
652         if(IsTrimmed[0] || IsTrimmed[1]) {
653           if(IsTrimmed[0]) {
654             //Add segment before start of curve
655             f = myCurve->FirstParameter();
656             ExtendC2d(aRes, f, -dt, U1, U2, V1, V2);
657           }
658           if(IsTrimmed[1]) {
659             //Add segment after end of curve
660             l = myCurve->LastParameter();
661             ExtendC2d(aRes, l, dt, U1, U2, V1, V2);
662           }
663         }
664           
665         myResult.SetBSpline(aRes);
666         myResult.Done();
667         myResult.SetType(GeomAbs_BSplineCurve);
668       }
669       break;
670
671     default:
672       {
673         Standard_Boolean IsTrimmed[2] = {Standard_False, Standard_False};
674         Standard_Real f, l, dt;
675         const Standard_Real eps = 0.01;
676         
677         if(mySurface->GetType() == GeomAbs_SurfaceOfRevolution) {
678           //Check possible singularity
679
680           gp_Pnt P = mySurface->AxeOfRevolution().Location();
681           gp_Dir N = mySurface->AxeOfRevolution().Direction();
682
683           gp_Lin L(P, N);
684
685           f = myCurve->FirstParameter();
686           l = myCurve->LastParameter();
687           dt = (l-f)*eps;
688
689           P = myCurve->Value(f);
690           if(L.Distance(P) < Precision::Confusion()) {
691             IsTrimmed[0] = Standard_True;
692             f = f+dt;
693             myCurve = myCurve->Trim(f, l, Precision::Confusion());
694           }
695
696           P = myCurve->Value(l);
697           if(L.Distance(P) < Precision::Confusion()) {
698             IsTrimmed[1] = Standard_True;
699             l = l-dt;
700             myCurve = myCurve->Trim(f, l, Precision::Confusion());
701           }
702         }
703
704         ProjLib_CompProjectedCurve Projector(mySurface,myCurve,
705                                              myTolerance,myTolerance);
706         Handle(ProjLib_HCompProjectedCurve) HProjector = 
707           new ProjLib_HCompProjectedCurve();
708         HProjector->Set(Projector);
709
710         // Normalement, dans le cadre de ProjLib, le resultat 
711         // doit etre une et une seule courbe !!!
712         // De plus, cette courbe ne doit pas etre Single point
713         Standard_Integer NbCurves = Projector.NbCurves();
714         Standard_Real Udeb,Ufin;
715         if (NbCurves > 0) {
716           Projector.Bounds(1,Udeb,Ufin);
717         }
718         else {
719           StdFail_NotDone::Raise("ProjLib CompProjectedCurve Not Done");
720         }
721         // Approximons cette courbe algorithmique.
722         Standard_Boolean Only3d = Standard_False;
723         Standard_Boolean Only2d = Standard_True;
724         GeomAbs_Shape Continuity = GeomAbs_C1;
725         Standard_Integer MaxDegree = 14;
726         Standard_Integer MaxSeg    = 16;
727
728         Approx_CurveOnSurface appr(HProjector, mySurface, Udeb, Ufin, 
729                                    myTolerance, 
730                                    Continuity, MaxDegree, MaxSeg, 
731                                    Only3d, Only2d);
732
733         Handle(Geom2d_BSplineCurve) aRes = appr.Curve2d();
734
735         if(IsTrimmed[0] || IsTrimmed[1]) {
736           // Treatment only for surface of revolution
737           Standard_Real u1, u2, v1, v2;
738           u1 = mySurface->FirstUParameter();
739           u2 = mySurface->LastUParameter();
740           v1 = mySurface->FirstVParameter();
741           v2 = mySurface->LastVParameter();
742           
743           if(IsTrimmed[0]) {
744             //Add segment before start of curve
745             ExtendC2d(aRes, f, -dt, u1, u2, v1, v2);
746           }
747           if(IsTrimmed[1]) {
748             //Add segment after end of curve
749             ExtendC2d(aRes, l, dt, u1, u2, v1, v2);
750           }
751         }
752           
753         myResult.SetBSpline(aRes);
754         myResult.Done();
755         myResult.SetType(GeomAbs_BSplineCurve);
756       }
757   }
758   if ( !myResult.IsDone()) {
759     ProjLib_ComputeApprox Comp( myCurve, mySurface, myTolerance);
760     myResult.Done();
761     
762     // set the type
763     if ( SType == GeomAbs_Plane  &&  CType == GeomAbs_BezierCurve) {
764       myResult.SetType(GeomAbs_BezierCurve);
765       myResult.SetBezier(Comp.Bezier()) ;
766     }
767     else {
768       myResult.SetType(GeomAbs_BSplineCurve);
769       myResult.SetBSpline(Comp.BSpline()) ;
770     }
771     // set the periodicity flag
772     if ( SType == GeomAbs_Plane               && 
773          CType == GeomAbs_BSplineCurve        &&
774          myCurve->IsPeriodic()   ) {
775       myResult.SetPeriodic();
776     }
777     myTolerance = Comp.Tolerance();
778   }
779
780   else {
781     // On remet arbitrairement la tol atteinte a une valeur
782     // petite en attendant mieux. dub lbo 11/03/97
783     myTolerance = Min(myTolerance,Precision::Confusion());
784     
785     // Translate the projected curve to keep the first point
786     // In the canonical boundaries of periodic surfaces.
787     if (mySurface->IsUPeriodic()) {
788       // xf
789       Standard_Real aT1, aT2, aU1, aU2, aUPeriod, aUr, aUm, aUmid, dUm, dUr;
790       GeomAbs_CurveType aTypeR;
791       ProjLib_Projector aResult;
792       //
793       aT1=myCurve->FirstParameter();
794       aT2=myCurve->LastParameter();
795       aU1=mySurface->FirstUParameter();
796       aU2=mySurface->LastUParameter();
797       aUPeriod=mySurface->UPeriod();
798       //
799       aTypeR=myResult.GetType();
800       if ((aU2-aU1)<(aUPeriod-myTolerance) && aTypeR == GeomAbs_Line) {
801         aResult=myResult;
802         aResult.UFrame(aT1, aT2, aU1, aUPeriod);
803         //
804         gp_Lin2d &aLr = (gp_Lin2d &) aResult.Line();
805         aUr=aLr.Location().X();
806         gp_Lin2d &aLm = (gp_Lin2d &) myResult.Line();
807         aUm=aLm.Location().X();
808         //
809         aUmid=0.5*(aU2+aU1);
810         dUm=fabs(aUm-aUmid);
811         dUr=fabs(aUr-aUmid);
812         if (dUr<dUm) {
813           myResult=aResult;
814         }
815       }
816       else {
817         myResult.UFrame(aT1, aT2, aU1, aUPeriod);
818       }
819       //
820       /*
821       myResult.UFrame(myCurve->FirstParameter(),
822                       myCurve->LastParameter(),
823                       mySurface->FirstUParameter(),
824                       mySurface->UPeriod());
825       */
826       //xt
827       //  Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 Begin
828       //  Correct the U isoline in periodical surface
829       // to be inside restriction boundaries.
830       if (myResult.GetType() == GeomAbs_Line) {
831         gp_Lin2d &aLine = (gp_Lin2d &) myResult.Line();
832
833         Standard_Real aPeriod = mySurface->UPeriod();
834         Standard_Real aFUPar  = mySurface->FirstUParameter();
835         Standard_Real aLUPar  = mySurface->LastUParameter();
836
837         // Check if the parametric range is lower then the period.
838         if (aLUPar - aFUPar < aPeriod - myTolerance) {
839           Standard_Real aU = aLine.Location().X();
840
841           if (Abs(aU + aPeriod - aFUPar) < myTolerance ||
842               Abs(aU - aPeriod - aFUPar) < myTolerance) {
843             gp_Pnt2d aNewLoc(aFUPar, aLine.Location().Y());
844
845             aLine.SetLocation(aNewLoc);
846           } else if (Abs(aU + aPeriod - aLUPar) < myTolerance ||
847                      Abs(aU - aPeriod - aLUPar) < myTolerance) {
848             gp_Pnt2d aNewLoc(aLUPar, aLine.Location().Y());
849
850             aLine.SetLocation(aNewLoc);
851           }
852         }
853       }
854     }
855 //  Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 End
856
857     if (mySurface->IsVPeriodic()) {
858       myResult.VFrame(myCurve->FirstParameter(),
859                          myCurve->LastParameter(),
860                          mySurface->FirstVParameter(),
861                          mySurface->VPeriod());
862 //  Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 Begin
863 //  Correct the V isoline in a periodical surface
864 // to be inside restriction boundaries.
865       if (myResult.GetType() == GeomAbs_Line) {
866         gp_Lin2d &aLine = (gp_Lin2d &) myResult.Line();
867
868         Standard_Real aPeriod = mySurface->VPeriod();
869         Standard_Real aFVPar  = mySurface->FirstVParameter();
870         Standard_Real aLVPar  = mySurface->LastVParameter();
871
872         // Check if the parametric range is lower then the period.
873         if (aLVPar - aFVPar < aPeriod - myTolerance) {
874           Standard_Real aV = aLine.Location().Y();
875
876           if (Abs(aV + aPeriod - aFVPar) < myTolerance ||
877               Abs(aV - aPeriod - aFVPar) < myTolerance) {
878             gp_Pnt2d aNewLoc(aLine.Location().X(), aFVPar);
879
880             aLine.SetLocation(aNewLoc);
881           } else if (Abs(aV + aPeriod - aLVPar) < myTolerance ||
882                      Abs(aV - aPeriod - aLVPar) < myTolerance) {
883             gp_Pnt2d aNewLoc(aLine.Location().X(), aLVPar);
884
885             aLine.SetLocation(aNewLoc);
886           }
887         }
888       }
889     }
890 //  Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 End
891   } 
892 }
893
894
895 //=======================================================================
896 //function : GetSurface
897 //purpose  : 
898 //=======================================================================
899
900 const Handle(Adaptor3d_HSurface)& ProjLib_ProjectedCurve::GetSurface() const
901 {
902   return mySurface;
903 }
904
905
906 //=======================================================================
907 //function : GetCurve
908 //purpose  : 
909 //=======================================================================
910
911 const Handle(Adaptor3d_HCurve)& ProjLib_ProjectedCurve::GetCurve() const
912 {
913   return myCurve;
914 }
915
916
917 //=======================================================================
918 //function : GetTolerance
919 //purpose  : 
920 //=======================================================================
921
922 Standard_Real ProjLib_ProjectedCurve::GetTolerance() const 
923 {
924   return myTolerance;
925 }
926
927
928 //=======================================================================
929 //function : FirstParameter
930 //purpose  : 
931 //=======================================================================
932
933 Standard_Real ProjLib_ProjectedCurve::FirstParameter() const 
934 {
935   return myCurve->FirstParameter();
936 }
937
938
939 //=======================================================================
940 //function : LastParameter
941 //purpose  : 
942 //=======================================================================
943
944 Standard_Real ProjLib_ProjectedCurve::LastParameter() const 
945 {
946   return myCurve->LastParameter();
947 }
948
949
950 //=======================================================================
951 //function : Continuity
952 //purpose  : 
953 //=======================================================================
954
955 GeomAbs_Shape ProjLib_ProjectedCurve::Continuity() const
956 {
957   Standard_NotImplemented::Raise("");
958   return GeomAbs_C0;
959 }
960
961
962 //=======================================================================
963 //function : NbIntervals
964 //purpose  : 
965 //=======================================================================
966
967 Standard_Integer ProjLib_ProjectedCurve::NbIntervals(const GeomAbs_Shape ) const 
968 {
969   Standard_NotImplemented::Raise("");
970   return 0;
971 }
972
973
974 //=======================================================================
975 //function : Intervals
976 //purpose  : 
977 //=======================================================================
978
979 //void ProjLib_ProjectedCurve::Intervals(TColStd_Array1OfReal&  T,
980 void ProjLib_ProjectedCurve::Intervals(TColStd_Array1OfReal&  ,
981                                        const GeomAbs_Shape ) const 
982 {
983   Standard_NotImplemented::Raise("");
984 }
985
986
987 //=======================================================================
988 //function : IsClosed
989 //purpose  : 
990 //=======================================================================
991
992 Standard_Boolean ProjLib_ProjectedCurve::IsClosed() const
993 {
994   Standard_NotImplemented::Raise("");
995   return Standard_True;
996 }
997
998
999 //=======================================================================
1000 //function : IsPeriodic
1001 //purpose  : 
1002 //=======================================================================
1003
1004 Standard_Boolean ProjLib_ProjectedCurve::IsPeriodic() const
1005 {
1006   return myResult.IsPeriodic();
1007 }
1008
1009
1010 //=======================================================================
1011 //function : Period
1012 //purpose  : 
1013 //=======================================================================
1014
1015 Standard_Real ProjLib_ProjectedCurve::Period() const
1016 {
1017   Standard_NotImplemented::Raise("");
1018   return 0.;
1019 }
1020
1021
1022 //=======================================================================
1023 //function : Value
1024 //purpose  : 
1025 //=======================================================================
1026
1027 gp_Pnt2d ProjLib_ProjectedCurve::Value(const Standard_Real ) const 
1028 {
1029   Standard_NotImplemented::Raise("");
1030   return gp_Pnt2d(0.,0.);
1031 }
1032
1033
1034 //=======================================================================
1035 //function : D0
1036 //purpose  : 
1037 //=======================================================================
1038
1039 void ProjLib_ProjectedCurve::D0(const Standard_Real , gp_Pnt2d& ) const
1040 {
1041   Standard_NotImplemented::Raise("");
1042 }
1043
1044
1045 //=======================================================================
1046 //function : D1
1047 //purpose  : 
1048 //=======================================================================
1049
1050 void ProjLib_ProjectedCurve::D1(const Standard_Real ,
1051                                       gp_Pnt2d&     , 
1052                                       gp_Vec2d&     ) const 
1053 {
1054   Standard_NotImplemented::Raise("");
1055 }
1056
1057
1058 //=======================================================================
1059 //function : D2
1060 //purpose  : 
1061 //=======================================================================
1062
1063 void ProjLib_ProjectedCurve::D2(const Standard_Real , 
1064                                       gp_Pnt2d&     , 
1065                                       gp_Vec2d&     , 
1066                                       gp_Vec2d&     ) const 
1067 {
1068   Standard_NotImplemented::Raise("");
1069 }
1070
1071
1072 //=======================================================================
1073 //function : D3
1074 //purpose  : 
1075 //=======================================================================
1076
1077 void ProjLib_ProjectedCurve::D3(const Standard_Real, 
1078                                       gp_Pnt2d&, 
1079                                       gp_Vec2d&, 
1080                                       gp_Vec2d&, 
1081                                       gp_Vec2d&) const 
1082 {
1083   Standard_NotImplemented::Raise("");
1084 }
1085
1086
1087 //=======================================================================
1088 //function : DN
1089 //purpose  : 
1090 //=======================================================================
1091
1092 gp_Vec2d ProjLib_ProjectedCurve::DN(const Standard_Real, 
1093                                     const Standard_Integer) const 
1094 {
1095   Standard_NotImplemented::Raise("");
1096   return gp_Vec2d(0.,0.);
1097 }
1098
1099
1100 //=======================================================================
1101 //function : Resolution
1102 //purpose  : 
1103 //=======================================================================
1104
1105 Standard_Real ProjLib_ProjectedCurve::Resolution(const Standard_Real) const 
1106 {
1107   Standard_NotImplemented::Raise("");
1108   return 0.;
1109 }
1110     
1111
1112 //=======================================================================
1113 //function : GetType
1114 //purpose  : 
1115 //=======================================================================
1116
1117 GeomAbs_CurveType ProjLib_ProjectedCurve::GetType() const
1118 {
1119   return myResult.GetType();
1120 }
1121
1122
1123 //=======================================================================
1124 //function : Line
1125 //purpose  : 
1126 //=======================================================================
1127
1128 gp_Lin2d ProjLib_ProjectedCurve::Line() const
1129 {
1130   return myResult.Line();
1131 }
1132
1133
1134 //=======================================================================
1135 //function : Circle
1136 //purpose  : 
1137 //=======================================================================
1138
1139 gp_Circ2d ProjLib_ProjectedCurve::Circle() const
1140 {
1141   return myResult.Circle();
1142 }
1143
1144
1145 //=======================================================================
1146 //function : Ellipse
1147 //purpose  : 
1148 //=======================================================================
1149
1150 gp_Elips2d ProjLib_ProjectedCurve::Ellipse() const
1151 {
1152   return myResult.Ellipse();
1153 }
1154
1155
1156 //=======================================================================
1157 //function : Hyperbola
1158 //purpose  : 
1159 //=======================================================================
1160
1161 gp_Hypr2d ProjLib_ProjectedCurve::Hyperbola() const
1162 {
1163   return myResult.Hyperbola();
1164 }
1165
1166
1167 //=======================================================================
1168 //function : Parabola
1169 //purpose  : 
1170 //=======================================================================
1171
1172 gp_Parab2d ProjLib_ProjectedCurve::Parabola() const
1173 {
1174   return myResult.Parabola();
1175 }
1176
1177
1178
1179 //=======================================================================
1180 //function : Degree
1181 //purpose  : 
1182 //=======================================================================
1183
1184 Standard_Integer ProjLib_ProjectedCurve::Degree() const
1185 {
1186   Standard_NoSuchObject_Raise_if 
1187     ( (GetType() != GeomAbs_BSplineCurve) &&
1188       (GetType() != GeomAbs_BezierCurve),
1189      "ProjLib_ProjectedCurve:Degree");
1190   if (GetType() == GeomAbs_BSplineCurve) {
1191     return myResult.BSpline()->Degree();
1192   }
1193   else if (GetType() == GeomAbs_BezierCurve) {
1194     return myResult.Bezier()->Degree();
1195   }
1196
1197   // portage WNT
1198   return 0;
1199 }
1200
1201 //=======================================================================
1202 //function : IsRational
1203 //purpose  : 
1204 //=======================================================================
1205
1206 Standard_Boolean ProjLib_ProjectedCurve::IsRational() const 
1207 {
1208   Standard_NoSuchObject_Raise_if 
1209     ( (GetType() != GeomAbs_BSplineCurve) &&
1210       (GetType() != GeomAbs_BezierCurve),
1211      "ProjLib_ProjectedCurve:IsRational");
1212   if (GetType() == GeomAbs_BSplineCurve) {
1213     return myResult.BSpline()->IsRational();
1214   }
1215   else if (GetType() == GeomAbs_BezierCurve) {
1216     return myResult.Bezier()->IsRational();
1217   }
1218   // portage WNT
1219   return Standard_False;
1220 }
1221
1222 //=======================================================================
1223 //function : NbPoles
1224 //purpose  : 
1225 //=======================================================================
1226
1227 Standard_Integer ProjLib_ProjectedCurve::NbPoles() const
1228 {
1229   Standard_NoSuchObject_Raise_if 
1230     ( (GetType() != GeomAbs_BSplineCurve) &&
1231       (GetType() != GeomAbs_BezierCurve)   
1232      ,"ProjLib_ProjectedCurve:NbPoles"  );
1233   if (GetType() == GeomAbs_BSplineCurve) {
1234     return myResult.BSpline()->NbPoles();
1235   }
1236   else if (GetType() == GeomAbs_BezierCurve) {
1237     return myResult.Bezier()->NbPoles();
1238   }
1239
1240   // portage WNT
1241   return 0;
1242 }
1243
1244 //=======================================================================
1245 //function : NbKnots
1246 //purpose  : 
1247 //=======================================================================
1248
1249 Standard_Integer ProjLib_ProjectedCurve::NbKnots() const 
1250 {
1251   Standard_NoSuchObject_Raise_if ( GetType() != GeomAbs_BSplineCurve, 
1252                                   "ProjLib_ProjectedCurve:NbKnots");
1253   return myResult.BSpline()->NbKnots();
1254 }
1255
1256 //=======================================================================
1257 //function : Bezier
1258 //purpose  : 
1259 //=======================================================================
1260
1261 Handle(Geom2d_BezierCurve) ProjLib_ProjectedCurve::Bezier() const 
1262 {
1263  return myResult.Bezier() ;
1264 }
1265
1266 //=======================================================================
1267 //function : BSpline
1268 //purpose  : 
1269 //=======================================================================
1270
1271 Handle(Geom2d_BSplineCurve) ProjLib_ProjectedCurve::BSpline() const 
1272 {
1273  return myResult.BSpline() ;
1274 }
1275 //=======================================================================
1276 //function : Trim
1277 //purpose  : 
1278 //=======================================================================
1279
1280 Handle(Adaptor2d_HCurve2d) ProjLib_ProjectedCurve::Trim 
1281 //(const Standard_Real First,
1282 // const Standard_Real Last,
1283 // const Standard_Real Tolerance) const 
1284 (const Standard_Real ,
1285  const Standard_Real ,
1286  const Standard_Real ) const 
1287 {
1288   Standard_NotImplemented::Raise("");
1289   return NULL ;
1290 }
1291