573f438516588f9803963df4937c64b0c5a2b4d8
[occt.git] / src / Geom / Geom_OffsetCurve.cxx
1 // Created on: 1991-06-25
2 // Created by: JCV
3 // Copyright (c) 1991-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 // 24-Aug-95 : xab removed C1 and C2 test : appeller  D1 et D2 
23 //             avec discernement !
24 // 19-09-97  : JPI correction derivee seconde
25
26
27 #include <Geom_OffsetCurve.ixx>
28
29
30 #include <gp.hxx>
31 #include <gp_XYZ.hxx>
32 #include <Geom_Line.hxx>
33 #include <Geom_Circle.hxx>
34 #include <Geom_Ellipse.hxx>
35 #include <Geom_Hyperbola.hxx>
36 #include <Geom_Parabola.hxx>
37 #include <Geom_BezierCurve.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <Geom_BSplineCurve.hxx>
40 #include <Geom_Geometry.hxx>
41
42 #include <Geom_UndefinedDerivative.hxx>
43 #include <Geom_UndefinedValue.hxx>
44 #include <Standard_ConstructionError.hxx>
45 #include <Standard_RangeError.hxx>
46 #include <Standard_NotImplemented.hxx>
47
48 typedef Geom_OffsetCurve         OffsetCurve;
49 typedef Handle(Geom_OffsetCurve) Handle(OffsetCurve);
50 typedef Geom_Curve               Curve;
51 typedef Handle(Geom_Curve)       Handle(Curve);
52 typedef Handle(Geom_Geometry)    Handle(Geometry);
53
54 typedef gp_Dir  Dir;
55 typedef gp_Pnt  Pnt;
56 typedef gp_Trsf Trsf;
57 typedef gp_Vec  Vec;
58 typedef gp_XYZ  XYZ;
59
60
61
62
63
64 static const int MaxDegree = 9;
65    //ordre de derivation maximum pour la recherche de la premiere 
66    //derivee non nulle
67
68
69
70 //=======================================================================
71 //function : Copy
72 //purpose  : 
73 //=======================================================================
74
75 Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
76
77  Handle(OffsetCurve) C;
78  C = new OffsetCurve (basisCurve, offsetValue, direction);
79  return C;
80 }
81
82
83
84 //=======================================================================
85 //function : Geom_OffsetCurve
86 //purpose  : 
87 //=======================================================================
88
89 Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Curve)& C,
90                                     const Standard_Real           Offset, 
91                                     const Dir&           V      )
92  : direction(V), offsetValue(Offset) {
93   
94   if (C->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve)) {
95     Handle(OffsetCurve) OC = Handle(OffsetCurve)::DownCast(C->Copy());
96     if ((OC->BasisCurve())->Continuity() == GeomAbs_C0)  
97       Standard_ConstructionError::Raise();
98
99     basisCurve = Handle(Curve)::DownCast((OC->BasisCurve())->Copy());
100     Standard_Real PrevOff = OC->Offset();
101     gp_Vec V1(OC->Direction());
102     gp_Vec V2(direction);
103     gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
104
105     if (Offset >= 0.) {
106       offsetValue = Vdir.Magnitude();
107       direction.SetXYZ(Vdir.XYZ());
108     } else {
109       offsetValue = -Vdir.Magnitude();
110       direction.SetXYZ((-Vdir).XYZ());
111     }
112   } else {
113     if (C->Continuity() == GeomAbs_C0) Standard_ConstructionError::Raise();
114     basisCurve = Handle(Curve)::DownCast(C->Copy()); // DownCast: 10-03-93
115   }
116 }
117
118
119 //=======================================================================
120 //function : Reverse
121 //purpose  : 
122 //=======================================================================
123
124 void Geom_OffsetCurve::Reverse ()
125
126   basisCurve->Reverse();
127   offsetValue = -offsetValue;
128 }
129
130
131 //=======================================================================
132 //function : ReversedParameter
133 //purpose  : 
134 //=======================================================================
135
136 Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const 
137 {
138   return basisCurve->ReversedParameter( U);
139 }
140
141 //=======================================================================
142 //function : Direction
143 //purpose  : 
144 //=======================================================================
145
146 const gp_Dir& Geom_OffsetCurve::Direction () const               
147   { return direction; }
148
149 //=======================================================================
150 //function : SetDirection
151 //purpose  : 
152 //=======================================================================
153
154 void Geom_OffsetCurve::SetDirection (const Dir& V)     
155   { direction = V; }
156
157 //=======================================================================
158 //function : SetOffsetValue
159 //purpose  : 
160 //=======================================================================
161
162 void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)   
163   { offsetValue = D; }
164
165
166 //=======================================================================
167 //function : IsPeriodic
168 //purpose  : 
169 //=======================================================================
170
171 Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
172 {
173   return basisCurve->IsPeriodic();
174 }
175
176 //=======================================================================
177 //function : Period
178 //purpose  : 
179 //=======================================================================
180
181 Standard_Real Geom_OffsetCurve::Period () const
182 {
183   return basisCurve->Period();
184 }
185
186 //=======================================================================
187 //function : SetBasisCurve
188 //purpose  : 
189 //=======================================================================
190
191 void Geom_OffsetCurve::SetBasisCurve (const Handle(Curve)& C) {
192
193   if (C->Continuity() == GeomAbs_C0) Standard_ConstructionError::Raise();
194   basisCurve = Handle(Curve)::DownCast(C->Copy()); // DownCast: 10-03-93
195 }
196
197
198
199 //=======================================================================
200 //function : BasisCurve
201 //purpose  : 
202 //=======================================================================
203
204 Handle(Curve) Geom_OffsetCurve::BasisCurve () const 
205
206   return basisCurve;
207 }
208
209
210 //=======================================================================
211 //function : Continuity
212 //purpose  : 
213 //=======================================================================
214
215 GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
216
217   GeomAbs_Shape OffsetShape=GeomAbs_C0;
218   switch (basisCurve->Continuity()) {
219     case GeomAbs_C0 : OffsetShape = GeomAbs_C0;       break;
220     case GeomAbs_C1 : OffsetShape = GeomAbs_C0;       break;
221     case GeomAbs_C2 : OffsetShape = GeomAbs_C1;       break;
222     case GeomAbs_C3 : OffsetShape = GeomAbs_C2;       break;
223     case GeomAbs_CN : OffsetShape = GeomAbs_CN;       break;
224     case GeomAbs_G1 : OffsetShape = GeomAbs_G1;       break;
225     case GeomAbs_G2 : OffsetShape = GeomAbs_G2;       break;
226   }
227   return OffsetShape;
228 }
229
230
231 //=======================================================================
232 //function : D0
233 //purpose  : 
234 //=======================================================================
235
236 void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const 
237 {
238   gp_Pnt PBasis;
239   gp_Vec VBasis;
240   D0(U,P,PBasis,VBasis);
241 }
242
243
244 //=======================================================================
245 //function : D1
246 //purpose  : 
247 //=======================================================================
248
249 void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const 
250 {
251   gp_Pnt PBasis;
252   gp_Vec V1Basis,V2Basis;
253   D1(U,P,PBasis,V1,V1Basis,V2Basis);
254 }
255
256
257
258 //=======================================================================
259 //function : D2
260 //purpose  : 
261 //=======================================================================
262
263 void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const 
264 {
265   gp_Pnt PBasis;
266   gp_Vec V1Basis,V2Basis,V3Basis;
267   D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
268 }
269
270
271 //=======================================================================
272 //function : D3
273 //purpose  : 
274 //=======================================================================
275
276 void Geom_OffsetCurve::D3 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2, Vec& V3) 
277 const {
278
279
280    // P(u) = p(u) + Offset * Ndir / R
281    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
282
283    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
284
285    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
286    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
287
288    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
289    //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
290    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
291    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
292
293
294
295   basisCurve->D3 (U, P, V1, V2, V3);
296   Vec V4 = basisCurve->DN (U, 4);
297   Standard_Integer Index = 2;
298   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
299     V1 = basisCurve->DN (U, Index);
300     Index++;
301   }
302   if (Index != 2) {
303     V2 = basisCurve->DN (U, Index);
304     V3 = basisCurve->DN (U, Index + 1);
305     V4 = basisCurve->DN (U, Index + 2);
306   }
307   XYZ OffsetDir = direction.XYZ();
308   XYZ Ndir      = (V1.XYZ()).Crossed (OffsetDir);
309   XYZ DNdir     = (V2.XYZ()).Crossed (OffsetDir);
310   XYZ D2Ndir    = (V3.XYZ()).Crossed (OffsetDir);
311   XYZ D3Ndir    = (V4.XYZ()).Crossed (OffsetDir);
312   Standard_Real R2  = Ndir.SquareModulus();
313   Standard_Real R   = Sqrt (R2);
314   Standard_Real R3  = R2 * R;
315   Standard_Real R4  = R2 * R2;
316   Standard_Real R5  = R3 * R2;
317   Standard_Real R6  = R3 * R3;
318   Standard_Real R7  = R5 * R2;
319   Standard_Real Dr  = Ndir.Dot (DNdir);
320   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
321   Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
322   if (R7 <= gp::Resolution()) {
323     if (R6 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
324     // V3 = P"' (U) :
325     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
326     D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
327     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 -
328                                  15.0*Dr*Dr*Dr/R6 - D3r));
329     D3Ndir.Multiply (offsetValue/R);
330     V3.Add (Vec(D3Ndir));
331     // V2 = P" (U) :
332     Standard_Real R4 = R2 * R2;
333     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
334     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
335     D2Ndir.Multiply (offsetValue / R);
336     V2.Add (Vec(D2Ndir));
337     // V1 = P' (U) :
338     DNdir.Multiply(R);
339     DNdir.Subtract (Ndir.Multiplied (Dr/R));
340     DNdir.Multiply (offsetValue/R2);
341     V1.Add (Vec(DNdir));
342   }
343   else {
344     // V3 = P"' (U) :
345     D3Ndir.Divide (R);
346     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
347     D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
348     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 
349                 15.0*Dr*Dr*Dr/R7 - D3r));
350     D3Ndir.Multiply (offsetValue);
351     V3.Add (Vec(D3Ndir));
352     // V2 = P" (U) :
353     D2Ndir.Divide (R);
354     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
355     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
356     D2Ndir.Multiply (offsetValue);
357     V2.Add (Vec(D2Ndir));
358     // V1 = P' (U) :
359     DNdir.Multiply (offsetValue/R);
360     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
361     V1.Add (Vec(DNdir));
362   }
363   //P (U) :
364   Ndir.Multiply (offsetValue/R);
365   Ndir.Add (P.XYZ());
366   P.SetXYZ (Ndir);
367 }
368
369
370 //=======================================================================
371 //function : DN
372 //purpose  : 
373 //=======================================================================
374
375 Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const {
376
377
378   if (N < 1)                  Standard_RangeError::Raise();
379   XYZ OffsetDir = direction.XYZ();
380   Pnt P;
381   Vec V1, V2, dummy;
382   if (N == 1) {
383     basisCurve->D2 (U, P, V1, V2);
384     Standard_Integer Index = 2;
385     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
386       V1 = basisCurve->DN (U, Index);
387       Index++;
388     }
389     if (Index != 2)  V2 = basisCurve->DN (U, Index);
390     XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
391     XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
392     Standard_Real R2 = Ndir.SquareModulus();
393     Standard_Real R  = Sqrt (R2);
394     Standard_Real R3 = R * R2;
395     Standard_Real Dr = Ndir.Dot (DNdir);
396     if (R3 <= gp::Resolution()) {
397       if (R2 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
398       Ndir.Multiply (Dr/R);
399       DNdir.Multiply(R);
400       DNdir.Subtract (Ndir);
401       DNdir.Multiply (offsetValue/R2);
402       V1.Add (Vec(DNdir));
403     }
404     else {
405       Ndir.Multiply (offsetValue * Dr / R3);
406       DNdir.Multiply (offsetValue/R);
407       DNdir.Subtract (Ndir);        
408       V1.Add (Vec(DNdir));
409     }
410     dummy = V1;
411   }
412   else if (N == 2) {
413     Vec V3;
414     basisCurve->D3 (U, P, V1, V2, V3);
415     Standard_Integer Index = 2;
416     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
417       V1 = basisCurve->DN (U, Index);
418       Index++;
419     }
420     if (Index != 2) {
421       V2 = basisCurve->DN (U, Index);
422       V3 = basisCurve->DN (U, Index + 1);
423     }
424     XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
425     XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
426     XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
427     Standard_Real R2  = Ndir.SquareModulus();
428     Standard_Real R   = Sqrt (R2);
429     Standard_Real R3  = R2 * R;   Standard_Real R4 = R2 * R2;  Standard_Real R5 = R3 * R2;
430     Standard_Real Dr  = Ndir.Dot (DNdir);
431     Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
432     if (R5 <= gp::Resolution()) {
433       if (R4 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
434       Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r/R2));
435       DNdir.Multiply (2.0 * Dr / R2);
436       D2Ndir.Subtract (DNdir);
437       D2Ndir.Subtract (Ndir);
438       D2Ndir.Multiply (offsetValue / R);
439       V2.Add (Vec(D2Ndir));
440     }
441     else {
442       Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r / R2));
443       DNdir.Multiply (2.0 * Dr / R2);
444       D2Ndir.Divide (R);
445       D2Ndir.Subtract (DNdir);
446       D2Ndir.Subtract (Ndir);
447       D2Ndir.Multiply (offsetValue);
448       V2.Add (Vec(D2Ndir));
449     }
450     dummy = V2;
451   }
452   else if (N == 3) {
453     Vec V3;
454     basisCurve->D3 (U, P, V1, V2, V3);
455     Vec V4 = basisCurve->DN (U, 4);
456     Standard_Integer Index = 2;
457     while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
458       V1 = basisCurve->DN (U, Index);
459       Index++;
460     }
461     if (Index != 2) {
462       V2 = basisCurve->DN (U, Index);
463       V3 = basisCurve->DN (U, Index + 1);
464       V4 = basisCurve->DN (U, Index + 2);
465     }
466     XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
467     XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
468     XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
469     XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir);
470     Standard_Real R2  = Ndir.SquareModulus();
471     Standard_Real R   = Sqrt (R2);  Standard_Real R3 = R2 * R;  Standard_Real R4 = R2 * R2;
472     Standard_Real R5  = R3 * R2;  Standard_Real R6 = R3 * R3;  Standard_Real R7 = R5 * R2;
473     Standard_Real Dr  = Ndir.Dot (DNdir);
474     Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
475     Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
476     if (R7 <= gp::Resolution()) {
477       if (R6 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
478       D2Ndir.Multiply (3.0 * Dr / R2);
479       DNdir.Multiply (3.0 * ((D2r/R2) + (Dr*Dr)/R4));
480       Ndir.Multiply (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r);
481       D3Ndir.Subtract (D2Ndir);
482       D3Ndir.Subtract (DNdir);
483       D3Ndir.Add (Ndir);
484       D3Ndir.Multiply (offsetValue/R);
485       V3.Add (Vec(D3Ndir));
486     }
487     else {
488       D2Ndir.Multiply (3.0 * Dr / R3);
489       DNdir.Multiplied (3.0 * ((D2r/R3) + (Dr*Dr/R5)));
490       Ndir.Multiply (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r);
491       D3Ndir.Divide (R);
492       D3Ndir.Subtract (D2Ndir);
493       D3Ndir.Subtract (DNdir);
494       D3Ndir.Add (Ndir);
495       D3Ndir.Multiply (offsetValue);
496       V3.Add (Vec(D3Ndir));
497     }
498     dummy = V3;
499   }
500   else { Standard_NotImplemented::Raise();  }
501
502   return dummy;
503 }
504
505 //=======================================================================
506 //function : D0
507 //purpose  : 
508 //=======================================================================
509
510 void  Geom_OffsetCurve::D0(const Standard_Real U,
511                            gp_Pnt& P,
512                            gp_Pnt& Pbasis,
513                            gp_Vec& V1basis)const 
514 {
515
516   basisCurve->D1 (U, Pbasis, V1basis);
517   Standard_Integer Index = 2;
518   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
519     V1basis = basisCurve->DN (U, Index);
520     Index++;
521   }
522   XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
523   Standard_Real R = Ndir.Modulus();
524   if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
525   Ndir.Multiply (offsetValue/R);
526   Ndir.Add (Pbasis.XYZ());
527   P.SetXYZ(Ndir);
528 }
529
530 //=======================================================================
531 //function : D1
532 //purpose  : 
533 //=======================================================================
534
535 void Geom_OffsetCurve::D1 ( const Standard_Real U, 
536                             Pnt& P , Pnt& PBasis ,
537                             Vec& V1, Vec& V1basis, Vec& V2basis) const {
538
539    // P(u) = p(u) + Offset * Ndir / R
540    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
541
542    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
543
544   basisCurve->D2 (U, PBasis, V1basis, V2basis);
545   V1 = V1basis;
546   Vec V2 = V2basis;
547   Standard_Integer Index = 2;
548   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
549     V1 = basisCurve->DN (U, Index);
550     Index++;
551   }
552   if (Index != 2)  V2 = basisCurve->DN (U, Index);
553   XYZ OffsetDir = direction.XYZ();
554   XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
555   XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
556   Standard_Real R2 = Ndir.SquareModulus();
557   Standard_Real R  = Sqrt (R2);
558   Standard_Real R3 = R * R2;
559   Standard_Real Dr = Ndir.Dot (DNdir);
560   if (R3 <= gp::Resolution()) {
561     //We try another computation but the stability is not very good.
562     if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
563     DNdir.Multiply(R);
564     DNdir.Subtract (Ndir.Multiplied (Dr/R));
565     DNdir.Multiply (offsetValue/R2);
566     V1.Add (Vec(DNdir));
567   }
568   else {
569     // Same computation as IICURV in EUCLID-IS because the stability is
570     // better
571     DNdir.Multiply (offsetValue/R);
572     DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3));        
573     V1.Add (Vec(DNdir));
574   }
575   Ndir.Multiply (offsetValue/R);
576   Ndir.Add (PBasis.XYZ());
577   P.SetXYZ (Ndir);
578 }
579
580
581 //=======================================================================
582 //function : D2
583 //purpose  : 
584 //=======================================================================
585
586 void Geom_OffsetCurve::D2 (const Standard_Real U, 
587                            Pnt& P      , Pnt& PBasis ,
588                            Vec& V1     , Vec& V2     , 
589                            Vec& V1basis, Vec& V2basis, Vec& V3basis) const {
590
591    // P(u) = p(u) + Offset * Ndir / R
592    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
593
594    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
595
596    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
597    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
598
599   basisCurve->D3 (U, PBasis, V1basis, V2basis, V3basis);
600   Standard_Integer Index = 2;
601   V1     = V1basis;
602   V2     = V2basis;
603   Vec V3 = V3basis;
604   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
605     V1 = basisCurve->DN (U, Index);
606     Index++;
607   }
608   if (Index != 2) {
609     V2 = basisCurve->DN (U, Index);
610     V3 = basisCurve->DN (U, Index + 1);
611   }
612   XYZ OffsetDir = direction.XYZ();
613   XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
614   XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
615   XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
616   Standard_Real R2    = Ndir.SquareModulus();
617   Standard_Real R     = Sqrt (R2);
618   Standard_Real R3    = R2 * R;
619   Standard_Real R4    = R2 * R2;
620   Standard_Real R5    = R3 * R2;
621   Standard_Real Dr    = Ndir.Dot (DNdir);
622   Standard_Real D2r   = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
623   if (R5 <= gp::Resolution()) {
624     //We try another computation but the stability is not very good
625     //dixit ISG.
626     if (R4 <= gp::Resolution())  Geom_UndefinedDerivative::Raise();
627     // V2 = P" (U) :
628     Standard_Real R4 = R2 * R2;
629     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
630     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
631     D2Ndir.Multiply (offsetValue / R);
632     V2.Add (Vec(D2Ndir));
633     // V1 = P' (U) :
634     DNdir.Multiply(R);
635     DNdir.Subtract (Ndir.Multiplied (Dr/R));
636     DNdir.Multiply (offsetValue/R2);
637     V1.Add (Vec(DNdir));
638   }
639   else {
640     // Same computation as IICURV in EUCLID-IS because the stability is
641     // better.
642     // V2 = P" (U) :
643     D2Ndir.Multiply (offsetValue/R);
644     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
645     D2Ndir.Add (Ndir.Multiplied (
646                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
647                                      )
648                     );
649     V2.Add (Vec(D2Ndir));
650     // V1 = P' (U) :
651     DNdir.Multiply (offsetValue/R);
652     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
653     V1.Add (Vec(DNdir));
654   }
655   //P (U) :
656   Ndir.Multiply (offsetValue/R);
657   Ndir.Add (PBasis.XYZ());
658   P.SetXYZ (Ndir);
659 }
660
661
662 //=======================================================================
663 //function : FirstParameter
664 //purpose  : 
665 //=======================================================================
666
667 Standard_Real Geom_OffsetCurve::FirstParameter () const {
668
669    return basisCurve->FirstParameter();
670 }
671
672
673 //=======================================================================
674 //function : LastParameter
675 //purpose  : 
676 //=======================================================================
677
678 Standard_Real Geom_OffsetCurve::LastParameter () const {
679
680    return basisCurve->LastParameter();
681 }
682
683
684 //=======================================================================
685 //function : Offset
686 //purpose  : 
687 //=======================================================================
688
689 Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
690
691 //=======================================================================
692 //function : Value
693 //purpose  : 
694 //=======================================================================
695
696 void Geom_OffsetCurve::Value (
697 const Standard_Real U, Pnt& P, Pnt& PBasis,  Vec& V1basis) const {
698
699   if (basisCurve->Continuity() == GeomAbs_C0)  Geom_UndefinedValue::Raise();
700   basisCurve->D1 (U, PBasis, V1basis);
701   Standard_Integer Index = 2;
702   while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
703     V1basis = basisCurve->DN (U, Index);
704     Index++;
705   }
706   XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
707   Standard_Real R = Ndir.Modulus();
708   if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
709   Ndir.Multiply (offsetValue/R);
710   Ndir.Add (PBasis.XYZ());
711   P.SetXYZ (Ndir);
712 }
713
714
715 //=======================================================================
716 //function : IsClosed
717 //purpose  : 
718 //=======================================================================
719
720 Standard_Boolean Geom_OffsetCurve::IsClosed () const 
721
722   gp_Pnt PF,PL;
723   D0(FirstParameter(),PF);
724   D0(LastParameter(),PL);
725   return ( PF.Distance(PL) <= gp::Resolution());
726 }
727
728
729
730 //=======================================================================
731 //function : IsCN
732 //purpose  : 
733 //=======================================================================
734
735 Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
736
737    Standard_RangeError_Raise_if (N < 0, " ");
738    return basisCurve->IsCN (N + 1);
739 }
740
741
742 //=======================================================================
743 //function : Transform
744 //purpose  : 
745 //=======================================================================
746
747 void Geom_OffsetCurve::Transform (const Trsf& T) {
748
749   basisCurve->Transform (T);
750   direction.Transform(T);
751   offsetValue *= T.ScaleFactor();
752 }
753
754 //=======================================================================
755 //function : TransformedParameter
756 //purpose  : 
757 //=======================================================================
758
759 Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
760                                                      const gp_Trsf& T) const
761 {
762   return basisCurve->TransformedParameter(U,T);
763 }
764
765 //=======================================================================
766 //function : ParametricTransformation
767 //purpose  : 
768 //=======================================================================
769
770 Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
771 const
772 {
773   return basisCurve->ParametricTransformation(T);
774 }