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