5f9ea8d453d883ce4e187854092e2bfaa9f1c021
[occt.git] / src / GeomAdaptor / GeomAdaptor_Curve.cxx
1 // Created on: 1993-04-29
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 // 20/02/97 : PMN -> Positionement local sur BSpline (PRO6902)
18 // 10/07/97 : PMN -> Pas de calcul de resolution dans Nb(Intervals)(PRO9248)
19 // 20/10/97 : RBV -> traitement des offset curves
20
21 #define No_Standard_RangeError
22 #define No_Standard_OutOfRange
23
24 #include <GeomAdaptor_Curve.ixx>
25
26 #include <GeomAdaptor_HCurve.hxx>
27 #include <Adaptor3d_HCurve.hxx>
28 #include <BSplCLib.hxx>
29 #include <BSplCLib_Cache.hxx>
30 #include <GeomAbs_Shape.hxx>
31 #include <TColgp_Array1OfPnt.hxx>
32 #include <TColStd_Array1OfReal.hxx>
33 #include <TColStd_Array1OfInteger.hxx>
34 #include <TColStd_HArray1OfInteger.hxx>
35 #include <Precision.hxx>
36 #include <Geom_TrimmedCurve.hxx>
37 #include <Geom_Circle.hxx>
38 #include <Geom_Line.hxx>
39 #include <Geom_TrimmedCurve.hxx>
40 #include <Geom_BezierCurve.hxx>
41 #include <Geom_BSplineCurve.hxx>
42 #include <Geom_Ellipse.hxx>
43 #include <Geom_Parabola.hxx>
44 #include <Geom_Hyperbola.hxx>
45 //#include <GeomConvert_BSplineCurveKnotSplitting.hxx>
46
47 #include <Standard_OutOfRange.hxx>
48 #include <Standard_NoSuchObject.hxx>
49 #include <Standard_NullObject.hxx>
50 #include <Standard_NotImplemented.hxx>
51 #include <Geom_OffsetCurve.hxx>
52 #include <CSLib_Offset.hxx>
53
54 #define myBspl (*((Handle(Geom_BSplineCurve)*)&myCurve))
55 #define PosTol Precision::PConfusion()/2
56
57 static const int maxDerivOrder = 3;
58 static const Standard_Real MinStep   = 1e-7;
59
60 static gp_Vec dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative
61 // Recalculate derivatives in the singular point
62 // Returns true if the direction of derivatives is changed
63 static Standard_Boolean AdjustDerivative(
64     const Handle(Adaptor3d_HCurve)& theAdaptor, Standard_Integer theMaxDerivative, Standard_Real theU, gp_Vec& theD1,
65     gp_Vec& theD2 = dummyDerivative, gp_Vec& theD3 = dummyDerivative, gp_Vec& theD4 = dummyDerivative);
66
67
68 //=======================================================================
69 //function : LocalContinuity
70 //purpose  : Computes the Continuity of a BSplineCurve 
71 //           between the parameters U1 and U2
72 //           The continuity is C(d-m) 
73 //             with   d = degree, 
74 //                    m = max multiplicity of the Knots between U1 and U2
75 //=======================================================================
76
77 GeomAbs_Shape GeomAdaptor_Curve::LocalContinuity(const Standard_Real U1, 
78                                                  const Standard_Real U2) 
79      const 
80 {
81   Standard_NoSuchObject_Raise_if(myTypeCurve!=GeomAbs_BSplineCurve," ");
82   Standard_Integer Nb = myBspl->NbKnots();
83   Standard_Integer Index1 = 0;
84   Standard_Integer Index2 = 0;
85   Standard_Real newFirst, newLast;
86   TColStd_Array1OfReal    TK(1,Nb);
87   TColStd_Array1OfInteger TM(1,Nb);
88   myBspl->Knots(TK);
89   myBspl->Multiplicities(TM);
90   BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,U1,myBspl->IsPeriodic(),
91                             1,Nb,Index1,newFirst);
92   BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,U2,myBspl->IsPeriodic(),
93                             1,Nb,Index2,newLast);
94   if ( Abs(newFirst-TK(Index1+1))<Precision::PConfusion()) { 
95     if (Index1 < Nb) Index1++;
96   }
97   if ( Abs(newLast-TK(Index2))<Precision::PConfusion())
98     Index2--;
99   Standard_Integer MultMax;
100   // attention aux courbes peridiques.
101   if ( (myBspl->IsPeriodic()) && (Index1 == Nb) )
102     Index1 = 1;
103
104   if ( Index2 - Index1 <= 0) {
105     MultMax = 100;  // CN entre 2 Noeuds consecutifs
106   }
107   else {
108     MultMax = TM(Index1+1);
109     for(Standard_Integer i = Index1+1;i<=Index2;i++) {
110       if ( TM(i)>MultMax) MultMax=TM(i);
111     }
112     MultMax = myBspl->Degree() - MultMax;
113   }
114   if ( MultMax <= 0) {
115     return GeomAbs_C0;
116   }
117   else if ( MultMax == 1) {
118     return GeomAbs_C1;
119   } 
120   else if ( MultMax == 2) {
121     return GeomAbs_C2;
122   }
123   else if ( MultMax == 3) {
124     return GeomAbs_C3;
125   }
126   else { 
127     return GeomAbs_CN;
128   }
129 }
130
131
132 //=======================================================================
133 //function : Load
134 //purpose  : 
135 //=======================================================================
136
137 void GeomAdaptor_Curve::load(const Handle(Geom_Curve)& C,
138                              const Standard_Real UFirst,
139                              const Standard_Real ULast)
140 {
141   myFirst = UFirst;
142   myLast  = ULast;
143
144   if ( myCurve != C) {
145     myCurve = C;
146     
147     const Handle(Standard_Type)& TheType = C->DynamicType();
148     if ( TheType == STANDARD_TYPE(Geom_TrimmedCurve)) {
149       Load((*((Handle(Geom_TrimmedCurve)*)&C))->BasisCurve(),UFirst,ULast);
150     }
151     else if ( TheType ==  STANDARD_TYPE(Geom_Circle)) {
152       myTypeCurve = GeomAbs_Circle;
153     }
154     else if ( TheType ==STANDARD_TYPE(Geom_Line)) {
155       myTypeCurve = GeomAbs_Line;
156     }
157     else if ( TheType == STANDARD_TYPE(Geom_Ellipse)) {
158       myTypeCurve = GeomAbs_Ellipse;
159     }
160     else if ( TheType == STANDARD_TYPE(Geom_Parabola)) {
161       myTypeCurve = GeomAbs_Parabola;
162     }
163     else if ( TheType == STANDARD_TYPE(Geom_Hyperbola)) {
164       myTypeCurve = GeomAbs_Hyperbola;
165     }
166     else if ( TheType == STANDARD_TYPE(Geom_BezierCurve)) {
167       myTypeCurve = GeomAbs_BezierCurve;
168     }
169     else if ( TheType == STANDARD_TYPE(Geom_BSplineCurve)) {
170       myTypeCurve = GeomAbs_BSplineCurve;
171       // Create cache for B-spline
172       myCurveCache = new BSplCLib_Cache(myBspl->Degree(), myBspl->IsPeriodic(), 
173         myBspl->KnotSequence(), myBspl->Poles(), myBspl->Weights());
174     }
175     else if ( TheType == STANDARD_TYPE(Geom_OffsetCurve)) {
176       myTypeCurve = GeomAbs_OtherCurve;
177       // Create nested adaptor for base curve
178       Handle(Geom_Curve) aBase = Handle(Geom_OffsetCurve)::DownCast(myCurve)->BasisCurve();
179       myOffsetBaseCurveAdaptor = new GeomAdaptor_HCurve(aBase);
180     }
181     else {
182       myTypeCurve = GeomAbs_OtherCurve;
183     }
184   }
185 }  
186
187 //    --
188 //    --     Global methods - Apply to the whole curve.
189 //    --     
190
191 //=======================================================================
192 //function : Continuity
193 //purpose  : 
194 //=======================================================================
195
196 GeomAbs_Shape GeomAdaptor_Curve::Continuity() const 
197 {
198   if (myTypeCurve == GeomAbs_BSplineCurve)
199     return LocalContinuity(myFirst, myLast);
200
201   if (myCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
202   {
203     const GeomAbs_Shape S =
204       (*((Handle(Geom_OffsetCurve)*)&myCurve))->GetBasisCurveContinuity();
205     switch(S)
206     {
207       case GeomAbs_CN: return GeomAbs_CN;
208       case GeomAbs_C3: return GeomAbs_C2;
209       case GeomAbs_C2: return GeomAbs_C1;
210       case GeomAbs_C1: return GeomAbs_C0; 
211       case GeomAbs_G1: return GeomAbs_G1;
212       case GeomAbs_G2: return GeomAbs_G2;
213       default:
214         Standard_NoSuchObject::Raise("GeomAdaptor_Curve::Continuity");
215     }
216   }
217   else if (myTypeCurve == GeomAbs_OtherCurve) {
218     Standard_NoSuchObject::Raise("GeomAdaptor_Curve::Contunuity");
219   }
220
221   return GeomAbs_CN;
222 }
223
224 //=======================================================================
225 //function : NbIntervals
226 //purpose  : 
227 //=======================================================================
228
229 Standard_Integer GeomAdaptor_Curve::NbIntervals(const GeomAbs_Shape S) const
230 {
231   Standard_Integer myNbIntervals = 1;
232   Standard_Integer NbSplit;
233   if (myTypeCurve == GeomAbs_BSplineCurve) {
234     Standard_Integer FirstIndex = myBspl->FirstUKnotIndex();
235     Standard_Integer LastIndex  = myBspl->LastUKnotIndex();
236     TColStd_Array1OfInteger Inter (1, LastIndex-FirstIndex+1);
237     if ( S > Continuity()) {
238       Standard_Integer Cont;
239       switch ( S) {
240       case GeomAbs_G1:
241       case GeomAbs_G2:
242         Standard_DomainError::Raise("GeomAdaptor_Curve::NbIntervals");
243         break;
244       case GeomAbs_C0:
245         myNbIntervals = 1;
246         break;
247       case GeomAbs_C1:
248       case GeomAbs_C2:
249       case GeomAbs_C3: 
250       case GeomAbs_CN: 
251         {
252           if      ( S == GeomAbs_C1) Cont = 1;
253           else if ( S == GeomAbs_C2) Cont = 2;
254           else if ( S == GeomAbs_C3) Cont = 3;
255           else                       Cont = myBspl->Degree();
256           Standard_Integer FirstIndex = myBspl->FirstUKnotIndex();
257           Standard_Integer LastIndex  = myBspl->LastUKnotIndex();
258           Standard_Integer Degree = myBspl->Degree();
259           Standard_Integer NbKnots = myBspl->NbKnots();
260           TColStd_Array1OfInteger Mults (1, NbKnots);
261           myBspl->Multiplicities (Mults);
262           NbSplit = 1;
263           Standard_Integer Index   = FirstIndex;
264           Inter (NbSplit) = Index;
265           Index++;
266           NbSplit++;
267           while (Index < LastIndex) 
268             {
269               if (Degree - Mults (Index) < Cont) 
270                 {
271                   Inter (NbSplit) = Index;
272                   NbSplit++;
273                 }
274               Index++;
275             }
276           Inter (NbSplit) = Index;
277
278           Standard_Integer NbInt = NbSplit-1;
279           
280           Standard_Integer Nb = myBspl->NbKnots();
281           Standard_Integer Index1 = 0;
282           Standard_Integer Index2 = 0;
283           Standard_Real newFirst, newLast;
284           TColStd_Array1OfReal    TK(1,Nb);
285           TColStd_Array1OfInteger TM(1,Nb);
286           myBspl->Knots(TK);
287           myBspl->Multiplicities(TM);
288           BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,myFirst,
289                                     myBspl->IsPeriodic(),
290                                     1,Nb,Index1,newFirst);
291           BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,myLast,
292                                     myBspl->IsPeriodic(),
293                                     1,Nb,Index2,newLast);
294
295           // On decale eventuellement les indices  
296           // On utilise une "petite" tolerance, la resolution ne doit 
297           // servir que pour les tres longue courbes....(PRO9248)
298           Standard_Real Eps = Min(Resolution(Precision::Confusion()),
299                                   Precision::PConfusion()); 
300           if ( Abs(newFirst-TK(Index1+1))< Eps) Index1++;
301           if ( newLast-TK(Index2)> Eps) Index2++;
302           
303           myNbIntervals = 1;
304           for ( Standard_Integer i=1; i<=NbInt; i++)
305             if (Inter(i)>Index1 && Inter(i)<Index2) myNbIntervals++;
306         }
307         break;
308       }
309     }
310   }
311   
312   else if (myCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))){
313     GeomAbs_Shape BaseS=GeomAbs_C0;
314     switch(S){
315     case GeomAbs_G1:
316     case GeomAbs_G2:
317       Standard_DomainError::Raise("GeomAdaptor_Curve::NbIntervals");
318       break;
319     case GeomAbs_C0: BaseS = GeomAbs_C1; break;
320     case GeomAbs_C1: BaseS = GeomAbs_C2; break;
321     case GeomAbs_C2: BaseS = GeomAbs_C3; break;
322     default: BaseS = GeomAbs_CN;
323     }
324     GeomAdaptor_Curve C
325       ((*((Handle(Geom_OffsetCurve)*)&myCurve))->BasisCurve());
326     // akm 05/04/02 (OCC278)  If our curve is trimmed we must recalculate 
327     //                    the number of intervals obtained from the basis to
328     //              vvv   reflect parameter bounds
329     Standard_Integer iNbBasisInt = C.NbIntervals(BaseS), iInt;
330     if (iNbBasisInt>1)
331     {
332       TColStd_Array1OfReal rdfInter(1,1+iNbBasisInt);
333       C.Intervals(rdfInter,BaseS);
334       for (iInt=1; iInt<=iNbBasisInt; iInt++)
335         if (rdfInter(iInt)>myFirst && rdfInter(iInt)<myLast)
336           myNbIntervals++;
337     }
338     // akm 05/04/02 ^^^
339   }
340   return myNbIntervals;
341 }
342
343 //=======================================================================
344 //function : Intervals
345 //purpose  : 
346 //=======================================================================
347
348 void GeomAdaptor_Curve::Intervals(TColStd_Array1OfReal& T,
349                                   const GeomAbs_Shape S   ) const
350 {
351   Standard_Integer myNbIntervals = 1;
352   Standard_Integer NbSplit;
353   Standard_Real FirstParam = myFirst, LastParam = myLast;
354
355   if (myTypeCurve == GeomAbs_BSplineCurve) 
356     {
357       Standard_Integer FirstIndex = myBspl->FirstUKnotIndex();
358       Standard_Integer LastIndex  = myBspl->LastUKnotIndex();
359       TColStd_Array1OfInteger Inter (1, LastIndex-FirstIndex+1);
360       
361       if ( S > Continuity()) {
362         Standard_Integer Cont;
363         switch ( S) {
364         case GeomAbs_G1:
365         case GeomAbs_G2:
366           Standard_DomainError::Raise("Geom2dAdaptor_Curve::NbIntervals");
367           break;
368         case GeomAbs_C0:
369           myNbIntervals = 1;
370           break;
371         case GeomAbs_C1:
372         case GeomAbs_C2:
373         case GeomAbs_C3: 
374         case GeomAbs_CN: 
375           {
376             if      ( S == GeomAbs_C1) Cont = 1;
377             else if ( S == GeomAbs_C2) Cont = 2;
378             else if ( S == GeomAbs_C3) Cont = 3;
379             else                       Cont = myBspl->Degree();
380             Standard_Integer FirstIndex = myBspl->FirstUKnotIndex();
381             Standard_Integer LastIndex  = myBspl->LastUKnotIndex();
382             Standard_Integer Degree = myBspl->Degree();
383             Standard_Integer NbKnots = myBspl->NbKnots();
384             TColStd_Array1OfInteger Mults (1, NbKnots);
385             myBspl->Multiplicities (Mults);
386             NbSplit = 1;
387             Standard_Integer Index   = FirstIndex;
388             Inter (NbSplit) = Index;
389             Index++;
390             NbSplit++;
391             while (Index < LastIndex) 
392               {
393                 if (Degree - Mults (Index) < Cont) 
394                   {
395                     Inter (NbSplit) = Index;
396                     NbSplit++;
397                   }
398                 Index++;
399               }
400             Inter (NbSplit) = Index;
401             Standard_Integer NbInt = NbSplit-1;
402             //        GeomConvert_BSplineCurveKnotSplitting Convector(myBspl, Cont);
403             //        Standard_Integer NbInt = Convector.NbSplits()-1;
404             //        TColStd_Array1OfInteger Inter(1,NbInt+1);
405             //        Convector.Splitting( Inter);
406             
407             Standard_Integer Nb = myBspl->NbKnots();
408             Standard_Integer Index1 = 0;
409             Standard_Integer Index2 = 0;
410             Standard_Real newFirst, newLast;
411             TColStd_Array1OfReal    TK(1,Nb);
412             TColStd_Array1OfInteger TM(1,Nb);
413             myBspl->Knots(TK);
414             myBspl->Multiplicities(TM);
415             BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,myFirst,
416                                       myBspl->IsPeriodic(),
417                                       1,Nb,Index1,newFirst);
418             BSplCLib::LocateParameter(myBspl->Degree(),TK,TM,myLast,
419                                       myBspl->IsPeriodic(),
420                                       1,Nb,Index2,newLast);
421             FirstParam = newFirst;
422             LastParam = newLast;
423             // On decale eventuellement les indices  
424             // On utilise une "petite" tolerance, la resolution ne doit 
425             // servir que pour les tres longue courbes....(PRO9248)
426             Standard_Real Eps = Min(Resolution(Precision::Confusion()),
427                                     Precision::PConfusion()); 
428             if ( Abs(newFirst-TK(Index1+1))< Eps) Index1++;
429             if ( newLast-TK(Index2)> Eps) Index2++;
430             
431             Inter( 1) = Index1;
432             myNbIntervals = 1;
433             for ( Standard_Integer i=1; i<=NbInt; i++) {
434               if (Inter(i) > Index1 && Inter(i)<Index2 ) {
435                 myNbIntervals++;
436                 Inter(myNbIntervals) = Inter(i);
437               }
438             }
439             Inter(myNbIntervals+1) = Index2;
440             
441             for (Standard_Integer I=1;I<=myNbIntervals+1;I++) {
442               T(I) = TK(Inter(I));
443             }
444           }
445           break;
446         }
447       }
448     }
449
450   else if (myCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))){
451     GeomAbs_Shape BaseS=GeomAbs_C0;
452     switch(S){
453     case GeomAbs_G1:
454     case GeomAbs_G2:
455       Standard_DomainError::Raise("GeomAdaptor_Curve::NbIntervals");
456       break;
457     case GeomAbs_C0: BaseS = GeomAbs_C1; break;
458     case GeomAbs_C1: BaseS = GeomAbs_C2; break;
459     case GeomAbs_C2: BaseS = GeomAbs_C3; break;
460     default: BaseS = GeomAbs_CN;
461     }
462     GeomAdaptor_Curve C
463       ((*((Handle(Geom_OffsetCurve)*)&myCurve))->BasisCurve());
464     // akm 05/04/02 (OCC278)  If our curve is trimmed we must recalculate 
465     //                    the array of intervals obtained from the basis to
466     //              vvv   reflect parameter bounds
467     Standard_Integer iNbBasisInt = C.NbIntervals(BaseS), iInt;
468     if (iNbBasisInt>1)
469     {
470       TColStd_Array1OfReal rdfInter(1,1+iNbBasisInt);
471       C.Intervals(rdfInter,BaseS);
472       for (iInt=1; iInt<=iNbBasisInt; iInt++)
473         if (rdfInter(iInt)>myFirst && rdfInter(iInt)<myLast)
474           T(++myNbIntervals)=rdfInter(iInt);
475     }
476     // old - myNbIntervals = C.NbIntervals(BaseS);
477     // old - C.Intervals(T, BaseS);
478     // akm 05/04/02 ^^^
479   }
480   
481   T( T.Lower() ) = FirstParam;
482   T( T.Lower() + myNbIntervals ) = LastParam;
483 }
484
485 //=======================================================================
486 //function : Trim
487 //purpose  : 
488 //=======================================================================
489
490 Handle(Adaptor3d_HCurve) GeomAdaptor_Curve::Trim(const Standard_Real First,
491                                                  const Standard_Real Last,
492                                                  const Standard_Real /*Tol*/) const
493 {
494   return Handle(GeomAdaptor_HCurve)(new GeomAdaptor_HCurve(myCurve,First,Last));
495 }
496
497
498 //=======================================================================
499 //function : IsClosed
500 //purpose  : 
501 //=======================================================================
502
503 Standard_Boolean GeomAdaptor_Curve::IsClosed() const
504 {
505   if (!Precision::IsPositiveInfinite(myLast) &&
506       !Precision::IsNegativeInfinite(myFirst))
507   {
508     const gp_Pnt Pd = Value(myFirst);
509     const gp_Pnt Pf = Value(myLast);
510     return (Pd.Distance(Pf) <= Precision::Confusion());
511   }
512   return Standard_False;
513 }
514
515 //=======================================================================
516 //function : IsPeriodic
517 //purpose  : 
518 //=======================================================================
519
520 Standard_Boolean GeomAdaptor_Curve::IsPeriodic() const 
521 {
522   return (myCurve->IsPeriodic()? IsClosed() : Standard_False);
523 }
524
525 //=======================================================================
526 //function : Period
527 //purpose  : 
528 //=======================================================================
529
530 Standard_Real GeomAdaptor_Curve::Period() const 
531 {
532   return myCurve->LastParameter() - myCurve->FirstParameter();
533 }
534
535 //=======================================================================
536 //function : RebuildCache
537 //purpose  : 
538 //=======================================================================
539 void GeomAdaptor_Curve::RebuildCache(const Standard_Real theParameter) const
540 {
541   myCurveCache->BuildCache(theParameter, myBspl->Degree(), 
542     myBspl->IsPeriodic(), myBspl->KnotSequence(), 
543     myBspl->Poles(), myBspl->Weights());
544 }
545
546 //=======================================================================
547 //function : Value
548 //purpose  : 
549 //=======================================================================
550
551 gp_Pnt GeomAdaptor_Curve::Value(const Standard_Real U) const
552 {
553   if (myTypeCurve == GeomAbs_BSplineCurve)
554     return ValueBSpline(U);
555   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
556     return ValueOffset(U);
557   return myCurve->Value(U);
558 }
559
560 //=======================================================================
561 //function : ValueBSpline
562 //purpose  : 
563 //=======================================================================
564 gp_Pnt GeomAdaptor_Curve::ValueBSpline(const Standard_Real theU) const
565 {
566   if (theU == myFirst || theU == myLast)
567   {
568     Standard_Integer Ideb = 0, Ifin = 0;
569     if (theU == myFirst) {
570       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
571       if (Ideb<1) Ideb=1;
572       if (Ideb>=Ifin) Ifin = Ideb+1;
573     }
574     if (theU == myLast) {
575       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
576       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
577       if (Ideb>=Ifin) Ideb = Ifin-1;
578     }
579     return myBspl->LocalValue(theU, Ideb, Ifin);
580   }
581   else if (!myCurveCache.IsNull()) // use cached B-spline data
582   {
583     if (!myCurveCache->IsCacheValid(theU))
584       RebuildCache(theU);
585     gp_Pnt aRes;
586     myCurveCache->D0(theU, aRes);
587     return aRes;
588   }
589   return myCurve->Value(theU);
590 }
591
592 //=======================================================================
593 //function : ValueOffset
594 //purpose  : 
595 //=======================================================================
596 gp_Pnt GeomAdaptor_Curve::ValueOffset(const Standard_Real theU) const
597 {
598   gp_Pnt aP;
599   gp_Vec aV;
600   myOffsetBaseCurveAdaptor->D1(theU, aP, aV);
601   Standard_Boolean IsDirectionChange = Standard_False;
602   if(aV.SquareMagnitude() <= gp::Resolution())
603     IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 1, theU, aV);
604
605   Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve);
606   Standard_Real anOffsetVal = anOffC->Offset();
607   const gp_Dir& anOffsetDir = anOffC->Direction();
608
609   CSLib_Offset::D0(aP, aV, anOffsetDir, anOffsetVal, IsDirectionChange, aP);
610   return aP;
611 }
612
613 //=======================================================================
614 //function : D0
615 //purpose  : 
616 //=======================================================================
617
618 void GeomAdaptor_Curve::D0(const Standard_Real U, gp_Pnt& P) const
619 {
620   if (myTypeCurve == GeomAbs_BSplineCurve)
621     D0BSpline(U, P);
622   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
623     D0Offset(U, P);
624   else
625     myCurve->D0(U, P);
626 }
627
628 //=======================================================================
629 //function : D0BSpline
630 //purpose  : 
631 //=======================================================================
632 void GeomAdaptor_Curve::D0BSpline(const Standard_Real theU, gp_Pnt& theP) const
633 {
634   if (theU == myFirst || theU == myLast) 
635   {
636     Standard_Integer Ideb = 0, Ifin = 0;
637     if (theU == myFirst) {
638       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
639       if (Ideb<1) Ideb=1;
640       if (Ideb>=Ifin) Ifin = Ideb+1;
641     }
642     if (theU == myLast) {
643       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
644       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
645       if (Ideb>=Ifin) Ideb = Ifin-1;
646     }
647     myBspl->LocalD0(theU, Ideb, Ifin, theP);
648     return;
649   }
650   else if (!myCurveCache.IsNull()) // use cached B-spline data
651   {
652     if (!myCurveCache->IsCacheValid(theU))
653       RebuildCache(theU);
654     myCurveCache->D0(theU, theP);
655     return;
656   }
657   myCurve->D0(theU, theP);
658 }
659
660 //=======================================================================
661 //function : D0Offset
662 //purpose  : 
663 //=======================================================================
664 void GeomAdaptor_Curve::D0Offset(const Standard_Real theU, gp_Pnt& theP) const
665 {
666   theP = ValueOffset(theU);
667 }
668
669 //=======================================================================
670 //function : D1
671 //purpose  : 
672 //=======================================================================
673
674 void GeomAdaptor_Curve::D1(const Standard_Real U, gp_Pnt& P, gp_Vec& V) const 
675 {
676   if (myTypeCurve == GeomAbs_BSplineCurve)
677     D1BSpline(U, P, V);
678   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
679     D1Offset(U, P, V);
680   else
681     myCurve->D1(U, P, V);
682 }
683
684 //=======================================================================
685 //function : D1BSpline
686 //purpose  : 
687 //=======================================================================
688 void GeomAdaptor_Curve::D1BSpline(const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const
689 {
690   if (theU == myFirst || theU == myLast) 
691   {
692     Standard_Integer Ideb = 0, Ifin = 0;
693     if (theU == myFirst) {
694       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
695       if (Ideb<1) Ideb=1;
696       if (Ideb>=Ifin) Ifin = Ideb+1;
697     }
698     if (theU == myLast) {
699       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
700       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
701       if (Ideb>=Ifin) Ideb = Ifin-1;
702     }
703     myBspl->LocalD1(theU, Ideb, Ifin, theP, theV);
704     return;
705   }
706   else if (!myCurveCache.IsNull()) // use cached B-spline data
707   {
708     if (!myCurveCache->IsCacheValid(theU))
709       RebuildCache(theU);
710     myCurveCache->D1(theU, theP, theV);
711     return;
712   }
713   myCurve->D1(theU, theP, theV);
714 }
715
716 //=======================================================================
717 //function : D1Offset
718 //purpose  : 
719 //=======================================================================
720 void GeomAdaptor_Curve::D1Offset(const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const
721 {
722   gp_Vec aV2;
723   myOffsetBaseCurveAdaptor->D2 (theU, theP, theV, aV2);
724
725   Standard_Boolean IsDirectionChange = Standard_False;
726   if(theV.SquareMagnitude() <= gp::Resolution())
727     IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 2, theU, theV, aV2);
728
729   Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve);
730   Standard_Real anOffsetVal = anOffC->Offset();
731   const gp_Dir& anOffsetDir = anOffC->Direction();
732   CSLib_Offset::D1(theP, theV, aV2, anOffsetDir, anOffsetVal, IsDirectionChange, theP, theV);
733 }
734
735
736 //=======================================================================
737 //function : D2
738 //purpose  : 
739 //=======================================================================
740
741 void GeomAdaptor_Curve::D2(const Standard_Real U, 
742                            gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const 
743 {
744   if (myTypeCurve == GeomAbs_BSplineCurve)
745     D2BSpline(U, P, V1, V2);
746   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
747     D2Offset(U, P, V1, V2);
748   else
749     myCurve->D2(U, P, V1, V2);
750 }
751
752 //=======================================================================
753 //function : D2BSpline
754 //purpose  : 
755 //=======================================================================
756 void GeomAdaptor_Curve::D2BSpline(const Standard_Real theU, gp_Pnt& theP,
757                                   gp_Vec& theV1, gp_Vec& theV2) const
758 {
759   if (theU == myFirst || theU == myLast)
760   {
761     Standard_Integer Ideb = 0, Ifin = 0;
762     if (theU == myFirst) {
763       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
764       if (Ideb<1) Ideb=1;
765       if (Ideb>=Ifin) Ifin = Ideb+1;
766     }
767     if (theU == myLast) {
768       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
769       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
770       if (Ideb>=Ifin) Ideb = Ifin-1;
771     }
772     myBspl->LocalD2(theU, Ideb, Ifin, theP, theV1, theV2);
773     return;
774   }
775   else if (!myCurveCache.IsNull()) // use cached B-spline data
776   {
777     if (!myCurveCache->IsCacheValid(theU))
778       RebuildCache(theU);
779     myCurveCache->D2(theU, theP, theV1, theV2);
780     return;
781   }
782   myCurve->D2(theU, theP, theV1, theV2);
783 }
784
785 //=======================================================================
786 //function : D2Offset
787 //purpose  : 
788 //=======================================================================
789 void GeomAdaptor_Curve::D2Offset(const Standard_Real theU, gp_Pnt& theP,
790                                   gp_Vec& theV1, gp_Vec& theV2) const
791 {
792   gp_Vec V3;
793   myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, V3);
794
795   Standard_Boolean IsDirectionChange = Standard_False;
796   if(theV1.SquareMagnitude() <= gp::Resolution())
797     IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 3, theU, theV1, theV2, V3);
798
799   Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve);
800   Standard_Real anOffsetVal = anOffC->Offset();
801   const gp_Dir& anOffsetDir = anOffC->Direction();
802   CSLib_Offset::D2(theP, theV1, theV2, V3, anOffsetDir, anOffsetVal, IsDirectionChange, theP, theV1, theV2);
803 }
804
805 //=======================================================================
806 //function : D3
807 //purpose  : 
808 //=======================================================================
809
810 void GeomAdaptor_Curve::D3(const Standard_Real U, 
811                            gp_Pnt& P, gp_Vec& V1, 
812                            gp_Vec& V2, gp_Vec& V3) const 
813 {
814   if (myTypeCurve == GeomAbs_BSplineCurve)
815     D3BSpline(U, P, V1, V2, V3);
816   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
817     D3Offset(U, P, V1, V2, V3);
818   else
819     myCurve->D3(U, P, V1, V2, V3);
820 }
821
822 //=======================================================================
823 //function : D3BSpline
824 //purpose  : 
825 //=======================================================================
826 void GeomAdaptor_Curve::D3BSpline(const Standard_Real theU,
827                                   gp_Pnt& theP, gp_Vec& theV1,
828                                   gp_Vec& theV2, gp_Vec& theV3) const
829 {
830   if (theU == myFirst || theU == myLast)
831   {
832     Standard_Integer Ideb = 0, Ifin = 0;
833     if (theU == myFirst) {
834       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
835       if (Ideb<1) Ideb=1;
836       if (Ideb>=Ifin) Ifin = Ideb+1;
837     }
838     if (theU == myLast) {
839       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
840       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
841       if (Ideb>=Ifin) Ideb = Ifin-1;
842     }
843     myBspl->LocalD3(theU, Ideb, Ifin, theP, theV1, theV2, theV3);
844     return;
845   }
846   else if (!myCurveCache.IsNull()) // use cached B-spline data
847   {
848     if (!myCurveCache->IsCacheValid(theU))
849       RebuildCache(theU);
850     myCurveCache->D3(theU, theP, theV1, theV2, theV3);
851     return;
852   }
853   myCurve->D3(theU, theP, theV1, theV2, theV3);
854 }
855
856 //=======================================================================
857 //function : D3Offset
858 //purpose  : 
859 //=======================================================================
860 void GeomAdaptor_Curve::D3Offset(const Standard_Real theU,
861                                  gp_Pnt& theP, gp_Vec& theV1,
862                                  gp_Vec& theV2, gp_Vec& theV3) const
863 {
864   myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, theV3);
865   gp_Vec V4 = myOffsetBaseCurveAdaptor->DN(theU, 4);
866
867   Standard_Boolean IsDirectionChange = Standard_False;
868   if(theV1.SquareMagnitude() <= gp::Resolution())
869     IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 4, theU, theV1, theV2, theV3, V4);
870
871   Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve);
872   Standard_Real anOffsetVal = anOffC->Offset();
873   const gp_Dir& anOffsetDir = anOffC->Direction();
874   CSLib_Offset::D3(theP, theV1, theV2, theV3, V4, anOffsetDir, anOffsetVal, IsDirectionChange,
875                    theP, theV1, theV2, theV3);
876 }
877
878 //=======================================================================
879 //function : DN
880 //purpose  : 
881 //=======================================================================
882
883 gp_Vec GeomAdaptor_Curve::DN(const Standard_Real    U, 
884                              const Standard_Integer N) const 
885 {
886   if (myTypeCurve == GeomAbs_BSplineCurve)
887     return DNBSpline(U, N);
888   else if (myCurve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
889     return DNOffset(U, N);
890
891   return myCurve->DN(U, N);
892 }
893
894 gp_Vec GeomAdaptor_Curve::DNBSpline(const Standard_Real    U,
895                                     const Standard_Integer N) const
896 {
897   if ((U==myFirst || U==myLast))
898   {
899     Standard_Integer Ideb = 0, Ifin = 0;
900     if (U==myFirst) {
901       myBspl->LocateU(myFirst, PosTol, Ideb, Ifin);
902       if (Ideb<1) Ideb=1;
903       if (Ideb>=Ifin) Ifin = Ideb+1;
904     }
905     if (U==myLast) {
906       myBspl->LocateU(myLast, PosTol, Ideb, Ifin);
907       if (Ifin>myBspl->NbKnots()) Ifin = myBspl->NbKnots();
908       if (Ideb>=Ifin) Ideb = Ifin-1;
909     } 
910     return myBspl->LocalDN( U, Ideb, Ifin, N);
911   }
912   return myCurve->DN( U, N);
913 }
914
915 gp_Vec GeomAdaptor_Curve::DNOffset(const Standard_Real    U,
916                                    const Standard_Integer N) const
917 {
918   gp_Pnt aPnt;
919   gp_Vec aVec, aVN;
920
921   switch (N)
922   {
923   case 1:
924     D1Offset(U, aPnt, aVN);
925     break;
926   case 2:
927     D2Offset(U, aPnt, aVec, aVN);
928     break;
929   case 3:
930     D3Offset(U, aPnt, aVec, aVec, aVN);
931     break;
932   default:
933     aVN = myCurve->DN(U, N);
934   }
935   return aVN;
936 }
937
938 //=======================================================================
939 //function : Resolution
940 //purpose  : 
941 //=======================================================================
942
943 Standard_Real GeomAdaptor_Curve::Resolution(const Standard_Real R3D) const
944 {
945   switch ( myTypeCurve) {
946   case GeomAbs_Line :
947     return R3D;
948   case GeomAbs_Circle: {
949     Standard_Real R = (*((Handle(Geom_Circle)*)&myCurve))->Circ().Radius();
950     if ( R > R3D/2. )
951       return 2*ASin(R3D/(2*R));
952     else
953       return 2*M_PI;
954   }
955   case GeomAbs_Ellipse: {
956     return R3D / (*((Handle(Geom_Ellipse)*)&myCurve))->MajorRadius();
957   }
958   case GeomAbs_BezierCurve: {
959     Standard_Real res;
960     (*((Handle(Geom_BezierCurve)*)&myCurve))->Resolution(R3D,res);
961     return res;
962   }
963   case GeomAbs_BSplineCurve: {
964     Standard_Real res;
965     (*((Handle(Geom_BSplineCurve)*)&myCurve))->Resolution(R3D,res);
966     return res;
967   }
968   default:
969     return Precision::Parametric(R3D);
970   }  
971 }
972
973
974 //    --
975 //    --     The following methods must  be called when GetType returned
976 //    --     the corresponding type.
977 //    --     
978
979 //=======================================================================
980 //function : Line
981 //purpose  : 
982 //=======================================================================
983
984 gp_Lin GeomAdaptor_Curve::Line() const 
985 {
986   Standard_NoSuchObject_Raise_if(myTypeCurve != GeomAbs_Line, "");
987   return (*((Handle(Geom_Line)*)&myCurve))->Lin();  
988 }
989
990 //=======================================================================
991 //function : Circle
992 //purpose  : 
993 //=======================================================================
994
995 gp_Circ  GeomAdaptor_Curve::Circle() const 
996 {
997   Standard_NoSuchObject_Raise_if(myTypeCurve != GeomAbs_Circle, "");
998   return (*((Handle(Geom_Circle)*)&myCurve))->Circ();
999 }
1000
1001 //=======================================================================
1002 //function : Ellipse
1003 //purpose  : 
1004 //=======================================================================
1005
1006 gp_Elips GeomAdaptor_Curve::Ellipse() const 
1007 {
1008   Standard_NoSuchObject_Raise_if(myTypeCurve != GeomAbs_Ellipse, "");
1009   return (*((Handle(Geom_Ellipse)*)&myCurve))->Elips();
1010 }
1011
1012 //=======================================================================
1013 //function : Hyperbola
1014 //purpose  : 
1015 //=======================================================================
1016
1017 gp_Hypr GeomAdaptor_Curve::Hyperbola() const 
1018 {
1019   Standard_NoSuchObject_Raise_if(myTypeCurve != GeomAbs_Hyperbola, "");
1020   return (*((Handle(Geom_Hyperbola)*)&myCurve))->Hypr();  
1021 }
1022
1023 //=======================================================================
1024 //function : Parabola
1025 //purpose  : 
1026 //=======================================================================
1027
1028 gp_Parab GeomAdaptor_Curve::Parabola() const 
1029 {
1030   Standard_NoSuchObject_Raise_if(myTypeCurve != GeomAbs_Parabola, "");
1031   return (*((Handle(Geom_Parabola)*)&myCurve))->Parab();
1032 }
1033
1034 //=======================================================================
1035 //function : Degree
1036 //purpose  : 
1037 //=======================================================================
1038
1039 Standard_Integer GeomAdaptor_Curve::Degree() const
1040 {
1041   if (myTypeCurve == GeomAbs_BezierCurve)
1042     return (*((Handle(Geom_BezierCurve)*)&myCurve))->Degree();
1043   else if (myTypeCurve == GeomAbs_BSplineCurve)
1044     return (*((Handle(Geom_BSplineCurve)*)&myCurve))->Degree();
1045   else
1046     Standard_NoSuchObject::Raise();
1047   // portage WNT 
1048   return 0;
1049 }
1050
1051 //=======================================================================
1052 //function : IsRational
1053 //purpose  : 
1054 //=======================================================================
1055
1056 Standard_Boolean GeomAdaptor_Curve::IsRational() const {
1057   switch( myTypeCurve) {
1058   case GeomAbs_BSplineCurve:
1059     return (*((Handle(Geom_BSplineCurve)*)&myCurve))->IsRational();
1060   case GeomAbs_BezierCurve:
1061     return (*((Handle(Geom_BezierCurve)*)&myCurve))->IsRational();
1062   default:
1063     return Standard_False;
1064   }
1065 }
1066
1067 //=======================================================================
1068 //function : NbPoles
1069 //purpose  : 
1070 //=======================================================================
1071
1072 Standard_Integer GeomAdaptor_Curve::NbPoles() const
1073 {
1074   if (myTypeCurve == GeomAbs_BezierCurve)
1075     return (*((Handle(Geom_BezierCurve)*)&myCurve))->NbPoles();
1076   else if (myTypeCurve == GeomAbs_BSplineCurve)
1077     return (*((Handle(Geom_BSplineCurve)*)&myCurve))->NbPoles();
1078   else
1079     Standard_NoSuchObject::Raise();
1080   // portage WNT
1081   return 0;
1082 }
1083
1084 //=======================================================================
1085 //function : NbKnots
1086 //purpose  : 
1087 //=======================================================================
1088
1089 Standard_Integer GeomAdaptor_Curve::NbKnots() const
1090 {
1091   if ( myTypeCurve != GeomAbs_BSplineCurve)
1092     Standard_NoSuchObject::Raise("GeomAdaptor_Curve::NbKnots");
1093   return (*((Handle(Geom_BSplineCurve)*)&myCurve))->NbKnots();
1094 }
1095
1096 //=======================================================================
1097 //function : Bezier
1098 //purpose  : 
1099 //=======================================================================
1100
1101 Handle(Geom_BezierCurve) GeomAdaptor_Curve::Bezier() const 
1102 {
1103  if ( myTypeCurve != GeomAbs_BezierCurve)
1104     Standard_NoSuchObject::Raise("GeomAdaptor_Curve::Bezier");
1105   return *((Handle(Geom_BezierCurve)*)&myCurve);
1106 }
1107
1108 //=======================================================================
1109 //function : BSpline
1110 //purpose  : 
1111 //=======================================================================
1112
1113 Handle(Geom_BSplineCurve) GeomAdaptor_Curve::BSpline() const 
1114 {
1115  if ( myTypeCurve != GeomAbs_BSplineCurve)
1116     Standard_NoSuchObject::Raise("GeomAdaptor_Curve::BSpline");
1117
1118   return *((Handle(Geom_BSplineCurve)*)&myCurve);
1119 }
1120
1121
1122 // ============= Auxiliary functions ===================
1123 Standard_Boolean AdjustDerivative(const Handle(Adaptor3d_HCurve)& theAdaptor, Standard_Integer theMaxDerivative,
1124                                   Standard_Real theU, gp_Vec& theD1, gp_Vec& theD2,
1125                                   gp_Vec& theD3, gp_Vec& theD4)
1126 {
1127   static const Standard_Real aTol = gp::Resolution();
1128
1129   Standard_Boolean IsDirectionChange = Standard_False;
1130   const Standard_Real anUinfium   = theAdaptor->FirstParameter();
1131   const Standard_Real anUsupremum = theAdaptor->LastParameter();
1132
1133   const Standard_Real DivisionFactor = 1.e-3;
1134   Standard_Real du;
1135   if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
1136     du = 0.0;
1137   else
1138     du = anUsupremum - anUinfium;
1139
1140   const Standard_Real aDelta = Max(du * DivisionFactor, MinStep);
1141
1142   //Derivative is approximated by Taylor-series
1143   Standard_Integer anIndex = 1; //Derivative order
1144   gp_Vec V;
1145
1146   do
1147   {
1148     V =  theAdaptor->DN(theU, ++anIndex);
1149   }
1150   while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder);
1151
1152   Standard_Real u;
1153
1154   if(theU-anUinfium < aDelta)
1155     u = theU+aDelta;
1156   else
1157     u = theU-aDelta;
1158
1159   gp_Pnt P1, P2;
1160   theAdaptor->D0(Min(theU, u), P1);
1161   theAdaptor->D0(Max(theU, u), P2);
1162
1163   gp_Vec V1(P1, P2);
1164   IsDirectionChange = V.Dot(V1) < 0.0;
1165   Standard_Real aSign = IsDirectionChange ? -1.0 : 1.0;
1166
1167   theD1 = V * aSign;
1168   gp_Vec* aDeriv[3] = {&theD2, &theD3, &theD4};
1169   for (Standard_Integer i = 1; i < theMaxDerivative; i++)
1170     *(aDeriv[i-1]) = theAdaptor->DN(theU, anIndex + i) * aSign;
1171
1172   return IsDirectionChange;
1173 }