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