86743bcc92a057b19d28efa62cc4aad9bb7e840b
[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-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 // 24-Aug-95 : xab removed C1 and C2 test : appeller  D1 et D2 
18 //             avec discernement !
19 // 19-09-97  : JPI correction derivee seconde
20
21 #include <CSLib_Offset.hxx>
22 #include <Geom_BezierCurve.hxx>
23 #include <Geom_BSplineCurve.hxx>
24 #include <Geom_Circle.hxx>
25 #include <Geom_Curve.hxx>
26 #include <Geom_Ellipse.hxx>
27 #include <Geom_Geometry.hxx>
28 #include <Geom_Hyperbola.hxx>
29 #include <Geom_Line.hxx>
30 #include <Geom_OffsetCurve.hxx>
31 #include <Geom_Parabola.hxx>
32 #include <Geom_TrimmedCurve.hxx>
33 #include <Geom_UndefinedDerivative.hxx>
34 #include <Geom_UndefinedValue.hxx>
35 #include <gp.hxx>
36 #include <gp_Dir.hxx>
37 #include <gp_Pnt.hxx>
38 #include <gp_Trsf.hxx>
39 #include <gp_Vec.hxx>
40 #include <gp_XYZ.hxx>
41 #include <Standard_ConstructionError.hxx>
42 #include <Standard_NoSuchObject.hxx>
43 #include <Standard_NotImplemented.hxx>
44 #include <Standard_RangeError.hxx>
45 #include <Standard_Type.hxx>
46
47 typedef Geom_OffsetCurve         OffsetCurve;
48 typedef Geom_Curve               Curve;
49 typedef gp_Dir  Dir;
50 typedef gp_Pnt  Pnt;
51 typedef gp_Trsf Trsf;
52 typedef gp_Vec  Vec;
53 typedef gp_XYZ  XYZ;
54
55 //ordre de derivation maximum pour la recherche de la premiere 
56 //derivee non nulle
57 static const int maxDerivOrder = 3;
58 static const Standard_Real MinStep   = 1e-7;
59 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
60
61
62 static gp_Vec dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative
63 // Recalculate derivatives in the singular point
64 // Returns true if the direction of derivatives is changed
65 static Standard_Boolean AdjustDerivative(
66     const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative, Standard_Real theU, gp_Vec& theD1,
67     gp_Vec& theD2 = dummyDerivative, gp_Vec& theD3 = dummyDerivative, gp_Vec& theD4 = dummyDerivative);
68
69
70
71 //=======================================================================
72 //function : Copy
73 //purpose  : 
74 //=======================================================================
75
76 Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const {
77
78  Handle(Geom_OffsetCurve) C;
79  C = new OffsetCurve (basisCurve, offsetValue, direction);
80  return C;
81 }
82
83
84
85 //=======================================================================
86 //function : Geom_OffsetCurve
87 //purpose  : Basis curve cannot be an Offset curve or trimmed from
88 //            offset curve.
89 //=======================================================================
90
91 Geom_OffsetCurve::Geom_OffsetCurve (const Handle(Geom_Curve)& theCurve,
92                                     const Standard_Real       theOffset,
93                                     const gp_Dir&             theDir,
94                                     const Standard_Boolean    isTheNotCheckC0)
95  : direction(theDir), offsetValue(theOffset)
96 {
97   SetBasisCurve (theCurve, isTheNotCheckC0);
98 }
99
100
101 //=======================================================================
102 //function : Reverse
103 //purpose  : 
104 //=======================================================================
105
106 void Geom_OffsetCurve::Reverse ()
107
108   basisCurve->Reverse();
109   offsetValue = -offsetValue;
110 }
111
112
113 //=======================================================================
114 //function : ReversedParameter
115 //purpose  : 
116 //=======================================================================
117
118 Standard_Real Geom_OffsetCurve::ReversedParameter( const Standard_Real U) const 
119 {
120   return basisCurve->ReversedParameter( U);
121 }
122
123 //=======================================================================
124 //function : Direction
125 //purpose  : 
126 //=======================================================================
127
128 const gp_Dir& Geom_OffsetCurve::Direction () const               
129   { return direction; }
130
131 //=======================================================================
132 //function : SetDirection
133 //purpose  : 
134 //=======================================================================
135
136 void Geom_OffsetCurve::SetDirection (const Dir& V)     
137   { direction = V; }
138
139 //=======================================================================
140 //function : SetOffsetValue
141 //purpose  : 
142 //=======================================================================
143
144 void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D)   
145   { offsetValue = D; }
146
147
148 //=======================================================================
149 //function : IsPeriodic
150 //purpose  : 
151 //=======================================================================
152
153 Standard_Boolean Geom_OffsetCurve::IsPeriodic () const
154 {
155   return basisCurve->IsPeriodic();
156 }
157
158 //=======================================================================
159 //function : Period
160 //purpose  : 
161 //=======================================================================
162
163 Standard_Real Geom_OffsetCurve::Period () const
164 {
165   return basisCurve->Period();
166 }
167
168 //=======================================================================
169 //function : SetBasisCurve
170 //purpose  : 
171 //=======================================================================
172
173 void Geom_OffsetCurve::SetBasisCurve (const Handle(Geom_Curve)& C,
174                                       const Standard_Boolean isNotCheckC0)
175 {
176   const Standard_Real aUf = C->FirstParameter(),
177                       aUl = C->LastParameter();
178   Handle(Geom_Curve) aCheckingCurve =  Handle(Geom_Curve)::DownCast(C->Copy());
179   Standard_Boolean isTrimmed = Standard_False;
180
181   while(aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ||
182         aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
183   {
184     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
185     {
186       Handle(Geom_TrimmedCurve) aTrimC = 
187                 Handle(Geom_TrimmedCurve)::DownCast(aCheckingCurve);
188       aCheckingCurve = aTrimC->BasisCurve();
189       isTrimmed = Standard_True;
190     }
191
192     if (aCheckingCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
193     {
194       Handle(Geom_OffsetCurve) aOC = 
195             Handle(Geom_OffsetCurve)::DownCast(aCheckingCurve);
196       aCheckingCurve = aOC->BasisCurve();
197       Standard_Real PrevOff = aOC->Offset();
198       gp_Vec V1(aOC->Direction());
199       gp_Vec V2(direction);
200       gp_Vec Vdir(PrevOff*V1 + offsetValue*V2);
201
202       if (offsetValue >= 0.)
203       {
204         offsetValue = Vdir.Magnitude();
205         direction.SetXYZ(Vdir.XYZ());
206       }
207       else
208       {
209         offsetValue = -Vdir.Magnitude();
210         direction.SetXYZ((-Vdir).XYZ());
211       }
212     }
213   }
214   
215   myBasisCurveContinuity = aCheckingCurve->Continuity();
216
217   Standard_Boolean isC0 = !isNotCheckC0 &&
218                           (myBasisCurveContinuity == GeomAbs_C0);
219
220   // Basis curve must be at least C1
221   if (isC0 && aCheckingCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve)))
222   {
223     Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCheckingCurve);
224     if(aBC->IsG1(aUf, aUl, MyAngularToleranceForG1))
225     {
226       //Checking if basis curve has more smooth (C1, G2 and above) is not done.
227       //It can be done in case of need.
228       myBasisCurveContinuity = GeomAbs_G1;
229       isC0 = Standard_False;
230     }
231
232     // Raise exception if still C0
233     if (isC0)
234       Standard_ConstructionError::Raise("Offset on C0 curve");
235   }
236   //
237   if(isTrimmed)
238   {
239     basisCurve = new Geom_TrimmedCurve(aCheckingCurve, aUf, aUl);
240   } 
241   else
242   {
243     basisCurve = aCheckingCurve;
244   }
245 }
246
247
248
249 //=======================================================================
250 //function : BasisCurve
251 //purpose  : 
252 //=======================================================================
253
254 Handle(Geom_Curve) Geom_OffsetCurve::BasisCurve () const 
255
256   return basisCurve;
257 }
258
259
260 //=======================================================================
261 //function : Continuity
262 //purpose  : 
263 //=======================================================================
264
265 GeomAbs_Shape Geom_OffsetCurve::Continuity () const {
266
267   GeomAbs_Shape OffsetShape=GeomAbs_C0;
268   switch (myBasisCurveContinuity) {
269     case GeomAbs_C0 : OffsetShape = GeomAbs_C0;       break;
270     case GeomAbs_C1 : OffsetShape = GeomAbs_C0;       break;
271     case GeomAbs_C2 : OffsetShape = GeomAbs_C1;       break;
272     case GeomAbs_C3 : OffsetShape = GeomAbs_C2;       break;
273     case GeomAbs_CN : OffsetShape = GeomAbs_CN;       break;
274     case GeomAbs_G1 : OffsetShape = GeomAbs_G1;       break;
275     case GeomAbs_G2 : OffsetShape = GeomAbs_G2;       break;
276   }
277   return OffsetShape;
278 }
279
280
281 //=======================================================================
282 //function : D0
283 //purpose  : 
284 //=======================================================================
285
286 void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const 
287 {
288   gp_Pnt PBasis;
289   gp_Vec VBasis;
290   D0(U,P,PBasis,VBasis);
291 }
292
293
294 //=======================================================================
295 //function : D1
296 //purpose  : 
297 //=======================================================================
298
299 void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const 
300 {
301   gp_Pnt PBasis;
302   gp_Vec V1Basis,V2Basis;
303   D1(U,P,PBasis,V1,V1Basis,V2Basis);
304 }
305
306
307
308 //=======================================================================
309 //function : D2
310 //purpose  : 
311 //=======================================================================
312
313 void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const 
314 {
315   gp_Pnt PBasis;
316   gp_Vec V1Basis,V2Basis,V3Basis;
317   D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis);
318 }
319
320
321 //=======================================================================
322 //function : D3
323 //purpose  : 
324 //=======================================================================
325
326 void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& theP, Vec& theV1, Vec& theV2, Vec& theV3) const
327 {
328    // P(u) = p(u) + Offset * Ndir / R
329    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
330
331    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
332
333    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
334    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
335
336    //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
337    //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
338    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
339    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
340
341   Standard_Boolean IsDirectionChange = Standard_False;
342
343   basisCurve->D3 (theU, theP, theV1, theV2, theV3);
344   Vec aV4 = basisCurve->DN (theU, 4);
345   if(theV1.SquareMagnitude() <= gp::Resolution())
346     IsDirectionChange = AdjustDerivative(basisCurve, 4, theU, theV1, theV2, theV3, aV4);
347
348   CSLib_Offset::D3(theP, theV1, theV2, theV3, aV4, direction, offsetValue,
349                    IsDirectionChange, theP, theV1, theV2, theV3);
350 }
351
352
353 //=======================================================================
354 //function : DN
355 //purpose  : 
356 //=======================================================================
357
358 Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const 
359   {
360   Standard_RangeError_Raise_if (N < 1, "Exception: "
361                               "Geom_OffsetCurve::DN(...). N<1.");
362
363   gp_Vec VN, Vtemp;
364   gp_Pnt Ptemp;
365   switch (N)
366     {
367     case 1:
368       D1( U, Ptemp, VN);
369       break;
370     case 2:
371       D2( U, Ptemp, Vtemp, VN);
372       break;
373     case 3:
374       D3( U, Ptemp, Vtemp, Vtemp, VN);
375       break;
376     default:
377       Standard_NotImplemented::Raise("Exception: "
378         "Derivative order is greater than 3. Cannot compute of derivative.");
379   }
380   
381   return VN;
382 }
383
384 //=======================================================================
385 //function : D0
386 //purpose  : 
387 //=======================================================================
388
389 void  Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP,
390                            gp_Pnt& thePbasis, gp_Vec& theV1basis)const 
391 {
392   basisCurve->D1(theU, thePbasis, theV1basis);
393   Standard_Boolean IsDirectionChange = Standard_False;
394   if(theV1basis.SquareMagnitude() <= gp::Resolution())
395     IsDirectionChange = AdjustDerivative(basisCurve, 1, theU, theV1basis);
396
397   CSLib_Offset::D0(thePbasis, theV1basis, direction, offsetValue, IsDirectionChange, theP);
398 }
399
400 //=======================================================================
401 //function : D1
402 //purpose  : 
403 //=======================================================================
404
405 void Geom_OffsetCurve::D1 ( const Standard_Real theU, 
406                             Pnt& theP , Pnt& thePBasis ,
407                             Vec& theV1, Vec& theV1basis, Vec& theV2basis) const {
408
409    // P(u) = p(u) + Offset * Ndir / R
410    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
411
412    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
413
414   basisCurve->D2 (theU, thePBasis, theV1basis, theV2basis);
415
416   Standard_Boolean IsDirectionChange = Standard_False;
417   if(theV1basis.SquareMagnitude() <= gp::Resolution())
418     IsDirectionChange = AdjustDerivative(basisCurve, 2, theU, theV1basis, theV2basis);
419
420   CSLib_Offset::D1(thePBasis, theV1basis, theV2basis, direction, offsetValue, IsDirectionChange, theP, theV1);
421 }
422
423
424 //=======================================================================
425 //function : D2
426 //purpose  : 
427 //=======================================================================
428
429 void Geom_OffsetCurve::D2 (const Standard_Real theU,
430                            Pnt& theP, Pnt& thePBasis,
431                            Vec& theV1, Vec& theV2,
432                            Vec& theV1basis, Vec& theV2basis, Vec& theV3basis) const
433 {
434    // P(u) = p(u) + Offset * Ndir / R
435    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
436
437    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
438
439    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
440    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
441
442   Standard_Boolean IsDirectionChange = Standard_False;
443
444   basisCurve->D3 (theU, thePBasis, theV1basis, theV2basis, theV3basis);
445
446   if(theV1basis.SquareMagnitude() <= gp::Resolution())
447     IsDirectionChange = AdjustDerivative(basisCurve, 3, theU, theV1basis, theV2basis, theV3basis);
448
449   CSLib_Offset::D2(thePBasis, theV1basis, theV2basis, theV3basis, direction, offsetValue,
450                    IsDirectionChange, theP, theV1, theV2);
451 }
452
453
454 //=======================================================================
455 //function : FirstParameter
456 //purpose  : 
457 //=======================================================================
458
459 Standard_Real Geom_OffsetCurve::FirstParameter () const {
460
461    return basisCurve->FirstParameter();
462 }
463
464
465 //=======================================================================
466 //function : LastParameter
467 //purpose  : 
468 //=======================================================================
469
470 Standard_Real Geom_OffsetCurve::LastParameter () const {
471
472    return basisCurve->LastParameter();
473 }
474
475
476 //=======================================================================
477 //function : Offset
478 //purpose  : 
479 //=======================================================================
480
481 Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
482
483 //=======================================================================
484 //function : Value
485 //purpose  : 
486 //=======================================================================
487
488 void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP, 
489                               Pnt& thePbasis,  Vec& theV1basis) const 
490 {
491   if (myBasisCurveContinuity == GeomAbs_C0)
492     Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!");
493
494   basisCurve->D1(theU, thePbasis, theV1basis);
495   D0(theU,theP);
496 }
497
498
499 //=======================================================================
500 //function : IsClosed
501 //purpose  : 
502 //=======================================================================
503
504 Standard_Boolean Geom_OffsetCurve::IsClosed () const 
505
506   gp_Pnt PF,PL;
507   D0(FirstParameter(),PF);
508   D0(LastParameter(),PL);
509   return ( PF.Distance(PL) <= gp::Resolution());
510 }
511
512
513
514 //=======================================================================
515 //function : IsCN
516 //purpose  : 
517 //=======================================================================
518
519 Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const {
520
521    Standard_RangeError_Raise_if (N < 0, " ");
522    return basisCurve->IsCN (N + 1);
523 }
524
525
526 //=======================================================================
527 //function : Transform
528 //purpose  : 
529 //=======================================================================
530
531 void Geom_OffsetCurve::Transform (const Trsf& T) {
532
533   basisCurve->Transform (T);
534   direction.Transform(T);
535   offsetValue *= T.ScaleFactor();
536 }
537
538 //=======================================================================
539 //function : TransformedParameter
540 //purpose  : 
541 //=======================================================================
542
543 Standard_Real Geom_OffsetCurve::TransformedParameter(const Standard_Real U,
544                                                      const gp_Trsf& T) const
545 {
546   return basisCurve->TransformedParameter(U,T);
547 }
548
549 //=======================================================================
550 //function : ParametricTransformation
551 //purpose  : 
552 //=======================================================================
553
554 Standard_Real Geom_OffsetCurve::ParametricTransformation(const gp_Trsf& T)
555 const
556 {
557   return basisCurve->ParametricTransformation(T);
558 }
559
560 //=======================================================================
561 //function : GetBasisCurveContinuity
562 //purpose  : 
563 //=======================================================================
564 GeomAbs_Shape Geom_OffsetCurve::GetBasisCurveContinuity() const
565 {
566   return myBasisCurveContinuity;
567 }
568
569
570 // ============= Auxiliary functions ===================
571 Standard_Boolean AdjustDerivative(const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative,
572                                   Standard_Real theU, gp_Vec& theD1, gp_Vec& theD2,
573                                   gp_Vec& theD3, gp_Vec& theD4)
574 {
575   static const Standard_Real aTol = gp::Resolution();
576
577   Standard_Boolean IsDirectionChange = Standard_False;
578   const Standard_Real anUinfium   = theCurve->FirstParameter();
579   const Standard_Real anUsupremum = theCurve->LastParameter();
580
581   const Standard_Real DivisionFactor = 1.e-3;
582   Standard_Real du;
583   if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
584     du = 0.0;
585   else
586     du = anUsupremum - anUinfium;
587
588   const Standard_Real aDelta = Max(du * DivisionFactor, MinStep);
589
590   //Derivative is approximated by Taylor-series
591   Standard_Integer anIndex = 1; //Derivative order
592   gp_Vec V;
593
594   do
595   {
596     V =  theCurve->DN(theU, ++anIndex);
597   }
598   while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder);
599
600   Standard_Real u;
601
602   if(theU-anUinfium < aDelta)
603     u = theU+aDelta;
604   else
605     u = theU-aDelta;
606
607   gp_Pnt P1, P2;
608   theCurve->D0(Min(theU, u), P1);
609   theCurve->D0(Max(theU, u), P2);
610
611   gp_Vec V1(P1, P2);
612   IsDirectionChange = V.Dot(V1) < 0.0;
613   Standard_Real aSign = IsDirectionChange ? -1.0 : 1.0;
614
615   theD1 = V * aSign;
616   gp_Vec* aDeriv[3] = {&theD2, &theD3, &theD4};
617   for (Standard_Integer i = 1; i < theMaxDerivative; i++)
618     *(aDeriv[i-1]) = theCurve->DN(theU, anIndex + i) * aSign;
619
620   return IsDirectionChange;
621 }