1 // Copyright (c) 2015-... OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 #include <CSLib_Offset.hxx>
16 #include <gp_Dir2d.hxx>
19 #include <Geom_UndefinedValue.hxx>
20 #include <Geom_UndefinedDerivative.hxx>
21 #include <Geom2d_UndefinedValue.hxx>
22 #include <Geom2d_UndefinedDerivative.hxx>
24 #include <Standard_NullValue.hxx>
27 // ========== Offset values for 2D curves ==========
29 void CSLib_Offset::D0(const gp_Pnt2d& theBasePoint,
30 const gp_Vec2d& theBaseDeriv,
31 Standard_Real theOffset,
32 Standard_Boolean , // unused
33 gp_Pnt2d& theResPoint)
35 if (theBaseDeriv.SquareMagnitude() <= gp::Resolution())
36 Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector "
37 "because tangent vector has zero-magnitude!");
39 gp_Dir2d aNormal(theBaseDeriv.Y(), -theBaseDeriv.X());
40 theResPoint.SetCoord(theBasePoint.X() + aNormal.X() * theOffset,
41 theBasePoint.Y() + aNormal.Y() * theOffset);
44 void CSLib_Offset::D1(const gp_Pnt2d& theBasePoint,
45 const gp_Vec2d& theBaseD1,
46 const gp_Vec2d& theBaseD2,
47 Standard_Real theOffset,
48 Standard_Boolean , // unused
49 gp_Pnt2d& theResPoint,
50 gp_Vec2d& theResDeriv)
52 // P(u) = p(u) + Offset * Ndir / R
53 // with R = || p' ^ Z|| and Ndir = P' ^ Z
55 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
57 gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
58 gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
59 Standard_Real R2 = Ndir.SquareModulus();
60 Standard_Real R = Sqrt (R2);
61 Standard_Real R3 = R * R2;
62 Standard_Real Dr = Ndir.Dot(DNdir);
63 if (R3 <= gp::Resolution())
65 if (R2 <= gp::Resolution())
66 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
67 //We try another computation but the stability is not very good.
69 DNdir.Subtract(Ndir.Multiplied(Dr / R));
70 DNdir.Multiply(theOffset / R2);
74 // Same computation as IICURV in EUCLID-IS because the stability is better
75 DNdir.Multiply(theOffset / R);
76 DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3));
80 D0(theBasePoint, theBaseD1, theOffset, Standard_False, theResPoint);
82 theResDeriv = theBaseD1.Added(gp_Vec2d(DNdir));
85 void CSLib_Offset::D2(const gp_Pnt2d& theBasePoint,
86 const gp_Vec2d& theBaseD1,
87 const gp_Vec2d& theBaseD2,
88 const gp_Vec2d& theBaseD3,
89 Standard_Real theOffset,
90 Standard_Boolean theIsDirectionChange,
91 gp_Pnt2d& theResPoint,
95 // P(u) = p(u) + Offset * Ndir / R
96 // with R = || p' ^ Z|| and Ndir = P' ^ Z
98 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
100 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
101 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
103 gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
104 gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
105 gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X());
106 Standard_Real R2 = Ndir.SquareModulus();
107 Standard_Real R = Sqrt(R2);
108 Standard_Real R3 = R2 * R;
109 Standard_Real R4 = R2 * R2;
110 Standard_Real R5 = R3 * R2;
111 Standard_Real Dr = Ndir.Dot(DNdir);
112 Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot (DNdir);
113 if (R5 <= gp::Resolution())
115 if (R4 <= gp::Resolution())
116 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
117 //We try another computation but the stability is not very good dixit ISG.
119 D2Ndir.Subtract(DNdir.Multiplied (2.0 * Dr / R2));
120 D2Ndir.Add(Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
121 D2Ndir.Multiply(theOffset / R);
125 DNdir.Subtract(Ndir.Multiplied(Dr / R));
126 DNdir.Multiply(theOffset / R2);
130 // Same computation as IICURV in EUCLID-IS because the stability is better.
132 D2Ndir.Multiply(theOffset / R);
133 D2Ndir.Subtract(DNdir.Multiplied (2.0 * theOffset * Dr / R3));
134 D2Ndir.Add (Ndir.Multiplied(theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
137 DNdir.Multiply(theOffset / R);
138 DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3));
142 D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint);
144 theResD1 = theBaseD1.Added(gp_Vec2d(DNdir));
146 if (theIsDirectionChange)
147 theResD2 = -theBaseD2;
149 theResD2 = theBaseD2;
150 theResD2.Add(gp_Vec2d(D2Ndir));
153 void CSLib_Offset::D3(const gp_Pnt2d& theBasePoint,
154 const gp_Vec2d& theBaseD1,
155 const gp_Vec2d& theBaseD2,
156 const gp_Vec2d& theBaseD3,
157 const gp_Vec2d& theBaseD4,
158 Standard_Real theOffset,
159 Standard_Boolean theIsDirectionChange,
160 gp_Pnt2d& theResPoint,
165 // P(u) = p(u) + Offset * Ndir / R
166 // with R = || p' ^ Z|| and Ndir = P' ^ Z
168 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
170 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
171 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
173 // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2 ) * D2Ndir -
174 // (3.0 * D2r / R2) * DNdir) + (3.0 * Dr * Dr / R4) * DNdir -
175 // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
176 // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
178 gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
179 gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
180 gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X());
181 gp_XY D3Ndir(theBaseD4.Y(), -theBaseD4.X());
182 Standard_Real R2 = Ndir.SquareModulus();
183 Standard_Real R = Sqrt (R2);
184 Standard_Real R3 = R2 * R;
185 Standard_Real R4 = R2 * R2;
186 Standard_Real R5 = R3 * R2;
187 Standard_Real R6 = R3 * R3;
188 Standard_Real R7 = R5 * R2;
189 Standard_Real Dr = Ndir.Dot(DNdir);
190 Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot (DNdir);
191 Standard_Real D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
193 if (R7 <= gp::Resolution())
195 if (R6 <= gp::Resolution())
196 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
197 //We try another computation but the stability is not very good dixit ISG.
199 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * theOffset * Dr / R2));
201 (DNdir.Multiplied ((3.0 * theOffset) * ((D2r/R2) + (Dr*Dr)/R4))));
202 D3Ndir.Add (Ndir.Multiplied (
203 (theOffset * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))));
204 D3Ndir.Multiply (theOffset/R);
206 Standard_Real R4 = R2 * R2;
207 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
208 D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
209 D2Ndir.Multiply (theOffset / R);
212 DNdir.Subtract (Ndir.Multiplied (Dr/R));
213 DNdir.Multiply (theOffset/R2);
217 // Same computation as IICURV in EUCLID-IS because the stability is better.
219 D3Ndir.Multiply (theOffset/R);
220 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * theOffset * Dr / R3));
221 D3Ndir.Subtract (DNdir.Multiplied (
222 ((3.0 * theOffset) * ((D2r/R3) + (Dr*Dr)/R5))) );
223 D3Ndir.Add (Ndir.Multiplied (
224 (theOffset * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))));
226 D2Ndir.Multiply (theOffset/R);
227 D2Ndir.Subtract (DNdir.Multiplied (2.0 * theOffset * Dr / R3));
228 D2Ndir.Subtract (Ndir.Multiplied (
229 theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
231 DNdir.Multiply (theOffset/R);
232 DNdir.Subtract (Ndir.Multiplied (theOffset*Dr/R3));
236 D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint);
238 theResD1 = theBaseD1.Added(gp_Vec2d(DNdir));
240 theResD2 = theBaseD2.Added(gp_Vec2d(D2Ndir));
242 if (theIsDirectionChange)
243 theResD3 = -theBaseD3;
245 theResD3 = theBaseD3;
246 theResD3.Add(gp_Vec2d(D2Ndir));
250 // ========== Offset values for 3D curves ==========
252 void CSLib_Offset::D0(const gp_Pnt& theBasePoint,
253 const gp_Vec& theBaseDeriv,
254 const gp_Dir& theOffsetDirection,
255 Standard_Real theOffsetValue,
256 Standard_Boolean , // unused
259 gp_XYZ Ndir = (theBaseDeriv.XYZ()).Crossed(theOffsetDirection.XYZ());
260 Standard_Real R = Ndir.Modulus();
261 if (R <= gp::Resolution())
262 Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector "
263 "because tangent vector has zero-magnitude!");
265 Ndir.Multiply(theOffsetValue / R);
266 Ndir.Add(theBasePoint.XYZ());
267 theResPoint.SetXYZ(Ndir);
270 void CSLib_Offset::D1(const gp_Pnt& theBasePoint,
271 const gp_Vec& theBaseD1,
272 const gp_Vec& theBaseD2,
273 const gp_Dir& theOffsetDirection,
274 Standard_Real theOffsetValue,
275 Standard_Boolean , // unused
279 // P(u) = p(u) + Offset * Ndir / R
280 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
282 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
284 gp_XYZ Ndir = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
285 gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
286 Standard_Real R2 = Ndir.SquareModulus();
287 Standard_Real R = Sqrt (R2);
288 Standard_Real R3 = R * R2;
289 Standard_Real Dr = Ndir.Dot (DNdir);
290 if (R3 <= gp::Resolution()) {
291 if (R2 <= gp::Resolution())
292 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
293 //We try another computation but the stability is not very good.
295 DNdir.Subtract(Ndir.Multiplied(Dr / R));
296 DNdir.Multiply(theOffsetValue / R2);
299 // Same computation as IICURV in EUCLID-IS because the stability is
301 DNdir.Multiply(theOffsetValue / R);
302 DNdir.Subtract(Ndir.Multiplied(theOffsetValue * Dr / R3));
306 D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, Standard_False, theResPoint);
308 theResDeriv = theBaseD1.Added(gp_Vec(DNdir));
311 void CSLib_Offset::D2(const gp_Pnt& theBasePoint,
312 const gp_Vec& theBaseD1,
313 const gp_Vec& theBaseD2,
314 const gp_Vec& theBaseD3,
315 const gp_Dir& theOffsetDirection,
316 Standard_Real theOffsetValue,
317 Standard_Boolean theIsDirectionChange,
322 // P(u) = p(u) + Offset * Ndir / R
323 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
325 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
327 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
328 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
330 gp_XYZ Ndir = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
331 gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
332 gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ());
333 Standard_Real R2 = Ndir.SquareModulus();
334 Standard_Real R = Sqrt (R2);
335 Standard_Real R3 = R2 * R;
336 Standard_Real R4 = R2 * R2;
337 Standard_Real R5 = R3 * R2;
338 Standard_Real Dr = Ndir.Dot (DNdir);
339 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
341 if (R5 <= gp::Resolution()) {
342 if (R4 <= gp::Resolution())
343 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
344 //We try another computation but the stability is not very good
347 Standard_Real R4 = R2 * R2;
348 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
349 D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
350 D2Ndir.Multiply (theOffsetValue / R);
354 DNdir.Subtract (Ndir.Multiplied (Dr/R));
355 DNdir.Multiply (theOffsetValue/R2);
358 // Same computation as IICURV in EUCLID-IS because the stability is
361 D2Ndir.Multiply (theOffsetValue/R);
362 D2Ndir.Subtract (DNdir.Multiplied (2.0 * theOffsetValue * Dr / R3));
363 D2Ndir.Add (Ndir.Multiplied (theOffsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
366 DNdir.Multiply (theOffsetValue/R);
367 DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3));
371 D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint);
373 theResD1 = theBaseD1.Added(gp_Vec(DNdir));
375 if (theIsDirectionChange)
376 theResD2 = -theBaseD2;
378 theResD2 = theBaseD2;
379 theResD2.Add(gp_Vec(D2Ndir));
382 void CSLib_Offset::D3(const gp_Pnt& theBasePoint,
383 const gp_Vec& theBaseD1,
384 const gp_Vec& theBaseD2,
385 const gp_Vec& theBaseD3,
386 const gp_Vec& theBaseD4,
387 const gp_Dir& theOffsetDirection,
388 Standard_Real theOffsetValue,
389 Standard_Boolean theIsDirectionChange,
395 // P(u) = p(u) + Offset * Ndir / R
396 // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
398 // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R))
400 // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
401 // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
403 //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
404 // (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
405 // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
406 // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
408 gp_XYZ Ndir = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
409 gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
410 gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ());
411 gp_XYZ D3Ndir = (theBaseD4.XYZ()).Crossed(theOffsetDirection.XYZ());
412 Standard_Real R2 = Ndir.SquareModulus();
413 Standard_Real R = Sqrt (R2);
414 Standard_Real R3 = R2 * R;
415 Standard_Real R4 = R2 * R2;
416 Standard_Real R5 = R3 * R2;
417 Standard_Real R6 = R3 * R3;
418 Standard_Real R7 = R5 * R2;
419 Standard_Real Dr = Ndir.Dot (DNdir);
420 Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
421 Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
422 if (R7 <= gp::Resolution()) {
423 if (R6 <= gp::Resolution())
424 Standard_NullValue::Raise("CSLib_Offset: Null derivative");
426 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
427 D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
428 D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r));
429 D3Ndir.Multiply (theOffsetValue/R);
431 Standard_Real R4 = R2 * R2;
432 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
433 D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
434 D2Ndir.Multiply (theOffsetValue / R);
437 DNdir.Subtract (Ndir.Multiplied (Dr/R));
438 DNdir.Multiply (theOffsetValue/R2);
443 D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
444 D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
445 D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r));
446 D3Ndir.Multiply (theOffsetValue);
449 D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
450 D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
451 D2Ndir.Multiply (theOffsetValue);
453 DNdir.Multiply (theOffsetValue/R);
454 DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3));
458 D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint);
460 theResD1 = theBaseD1.Added(gp_Vec(DNdir));
462 theResD2 = theBaseD2.Added(gp_Vec(D2Ndir));
464 if (theIsDirectionChange)
465 theResD3 = -theBaseD3;
467 theResD3 = theBaseD3;
468 theResD3.Add(gp_Vec(D2Ndir));