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