projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0029712: Extrema algorithm raises exception
[occt.git] / src / Extrema / Extrema_ExtElC2d.cxx
CommitLineData
b311480e 1// Created on: 1994-01-04
2// Created by: Christophe MARION
3// Copyright (c) 1994-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
7fd59977 17
7fd59977 18#include <ElCLib.hxx>
42cf5bc1 19#include <Extrema_ExtElC2d.hxx>
20#include <Extrema_ExtPElC2d.hxx>
21#include <Extrema_POnCurv2d.hxx>
22#include <gp_Circ2d.hxx>
23#include <gp_Elips2d.hxx>
24#include <gp_Hypr2d.hxx>
25#include <gp_Lin2d.hxx>
26#include <gp_Parab2d.hxx>
7fd59977 27#include <math_DirectPolynomialRoots.hxx>
42cf5bc1 28#include <math_TrigonometricFunctionRoots.hxx>
7fd59977 29#include <Precision.hxx>
42cf5bc1 30#include <Standard_NotImplemented.hxx>
31#include <Standard_OutOfRange.hxx>
32#include <StdFail_InfiniteSolutions.hxx>
33#include <StdFail_NotDone.hxx>
7fd59977 34
638ad7f3 35//=======================================================================
36//function : Extrema_ExtElC2d
37//purpose :
38//=======================================================================
39Extrema_ExtElC2d::Extrema_ExtElC2d()
40{
41 myDone = Standard_False;
42 myIsPar = Standard_False;
43 myNbExt = 0;
44
45 for (Standard_Integer i = 0; i < 8; i++)
46 {
47 mySqDist[i] = RealLast();
48 }
49}
7fd59977 50
15a954de 51//=======================================================================
52//function : Extrema_ExtElC2d
53//purpose :
54//=======================================================================
55Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Lin2d& C1,
56 const gp_Lin2d& C2,
57 const Standard_Real)
7fd59977 58/*-----------------------------------------------------------------------------
0d969553
Y		 59Function:
60 Find min distance between 2 straight lines.
61
62Method:
63 Let D1 and D2 be 2 directions of straight lines C1 and C2.
64 2 cases are considered:
65 1- if Angle(D1,D2) < AngTol, the straight lines are parallel.
66 The distance is the distance between any point of C1 and straight line C2.
67 2- if Angle(D1,D2) > AngTol:
68 Let P = C1(u1) and P =C2(u2) the point intersection:
7fd59977 69
70-----------------------------------------------------------------------------*/
71{
72 myDone = Standard_False;
73 myIsPar = Standard_False;
74 myNbExt = 0;
75
15a954de 76 gp_Vec2d D1(C1.Direction());
77 gp_Vec2d D2(C2.Direction());
78 if (D1.IsParallel(D2, Precision::Angular()))
79 {
7fd59977 80 myIsPar = Standard_True;
81 mySqDist[0] = C2.SquareDistance(C1.Location());
638ad7f3 82 myNbExt = 1;
7fd59977 83 }
15a954de 84 else
85 {
86 // Vector from P1 to P2 (P2 - P1).
87 gp_Vec2d aP1P2(C1.Location(), C2.Location());
88
89 // Solve linear system using Cramer's rule:
90 // D1.X * t1 + D2.X * (-t2) = P2.X - P1.X
91 // D1.Y * t1 + D2.Y * (-t2) = P2.Y - P1.Y
92
93 // There is no division by zero since lines are not parallel.
94 Standard_Real aDelim = 1 / (D1^D2);
95
96 Standard_Real aParam1 = (aP1P2 ^ D2) * aDelim;
97 Standard_Real aParam2 = -(D1 ^ aP1P2) * aDelim; // -1.0 coefficient before t2.
98
99 gp_Pnt2d P1 = ElCLib::Value(aParam1, C1);
100 gp_Pnt2d P2 = ElCLib::Value(aParam2, C2);
101
102 mySqDist[myNbExt] = 0.0;
103 myPoint[myNbExt][0] = Extrema_POnCurv2d(aParam1,P1);
104 myPoint[myNbExt][1] = Extrema_POnCurv2d(aParam2,P2);
105 myNbExt = 1;
7fd59977 106 }
15a954de 107
7fd59977 108 myDone = Standard_True;
109}
110//=============================================================================
111
112Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Lin2d& C1,
113 const gp_Circ2d& C2,
114 const Standard_Real)
115/*-----------------------------------------------------------------------------
0d969553
Y		 116Function:
117 Find extreme distances between straight line C1 and circle C2.
118
119Method:
120 Let P1=C1(u1) and P2=C2(u2) be two solution points
121 D the direction of straight line C1
122 T the tangent at point P2;
123 Then, ( P1P2.D = 0. (1)
7fd59977 124 ( P1P2.T = 0. (2)
125-----------------------------------------------------------------------------*/
126{
127 myIsPar = Standard_False;
128 myDone = Standard_False;
129 myNbExt = 0;
130
0d969553 131// Calculate T1 in the reference of the circle ...
7fd59977 132 gp_Dir2d D = C1.Direction();
133 gp_Dir2d x2, y2;
134 x2 = C2.XAxis().Direction();
135 y2 = C2.YAxis().Direction();
136
137 Standard_Real Dx = D.Dot(x2);
138 Standard_Real Dy = D.Dot(y2);
139 Standard_Real U1, teta[2];
140 gp_Pnt2d O1=C1.Location();
7fd59977 141 gp_Pnt2d P1, P2;
7fd59977 142
143 if (Abs(Dy) <= RealEpsilon()) {
c6541a0c 144 teta[0] = M_PI/2.0;
7fd59977 145 }
146 else teta[0] = ATan(-Dx/Dy);
c6541a0c
D		 147 teta[1] = teta[0]+ M_PI;
148 if (teta[0] < 0.0) teta[0] = teta[0] + 2.0*M_PI;
7fd59977 149
150 P2 = ElCLib::Value(teta[0], C2);
151 U1 = (gp_Vec2d(O1, P2)).Dot(D);
152 P1 = ElCLib::Value(U1, C1);
153 mySqDist[myNbExt] = P1.SquareDistance(P2);
154 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
155 myPoint[myNbExt][1] = Extrema_POnCurv2d(teta[0],P2);
156 myNbExt++;
157
158 P2 = ElCLib::Value(teta[1], C2);
159 U1 = (gp_Vec2d(O1, P2)).Dot(D);
160 P1 = ElCLib::Value(U1, C1);
161 mySqDist[myNbExt] = P1.SquareDistance(P2);
162 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
163 myPoint[myNbExt][1] = Extrema_POnCurv2d(teta[1],P2);
164 myNbExt++;
165 myDone = Standard_True;
166}
167
168
169// =============================================================================
170Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Lin2d& C1,
171 const gp_Elips2d& C2)
172{
173 myDone = Standard_True;
174 myIsPar = Standard_False;
175 myDone = Standard_False;
176 myNbExt = 0;
177
0d969553 178// Calculate T1 in the reference of the ellipse ...
7fd59977 179 gp_Dir2d D = C1.Direction();
180 gp_Dir2d x2, y2;
181 x2 = C2.XAxis().Direction();
182 y2 = C2.YAxis().Direction();
183
184 Standard_Real Dx = D.Dot(x2);
185 Standard_Real Dy = D.Dot(y2);
186 Standard_Real U1, teta[2], r1 = C2.MajorRadius(), r2 = C2.MinorRadius();
7fd59977 187 gp_Pnt2d O1=C1.Location(), P1, P2;
7fd59977 188
189 if (Abs(Dy) <= RealEpsilon()) {
c6541a0c 190 teta[0] = M_PI/2.0;
7fd59977 191 }
192 else teta[0] = ATan(-Dx*r2/(Dy*r1));
193
c6541a0c
D		 194 teta[1] = teta[0] + M_PI;
195 if (teta[0] < 0.0) teta[0] += 2.0*M_PI;
7fd59977 196 P2 = ElCLib::Value(teta[0], C2);
197 U1 = (gp_Vec2d(O1, P2)).Dot(D);
198 P1 = ElCLib::Value(U1, C1);
199 mySqDist[myNbExt] = P1.SquareDistance(P2);
200 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
201 myPoint[myNbExt][1] = Extrema_POnCurv2d(teta[0],P2);
202 myNbExt++;
203
204
205 P2 = ElCLib::Value(teta[1], C2);
206 U1 = (gp_Vec2d(O1, P2)).Dot(D);
207 P1 = ElCLib::Value(U1, C1);
208 mySqDist[myNbExt] = P1.SquareDistance(P2);
209 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
210 myPoint[myNbExt][1] = Extrema_POnCurv2d(teta[1],P2);
211 myNbExt++;
212 myDone = Standard_True;
213}
214
215
216
217//=============================================================================
218
219Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Lin2d& C1, const gp_Hypr2d& C2)
220{
221 myIsPar = Standard_False;
222 myDone = Standard_False;
223 myNbExt = 0;
224
0d969553 225// Calculate T1 in the reference of the parabole ...
7fd59977 226 gp_Dir2d D = C1.Direction();
227 gp_Dir2d x2, y2;
228 x2 = C2.XAxis().Direction();
229 y2 = C2.YAxis().Direction();
230 Standard_Real Dx = D.Dot(x2);
231 Standard_Real Dy = D.Dot(y2);
232
233 Standard_Real U1, v2, U2=0, R = C2.MajorRadius(), r = C2.MinorRadius();
234 gp_Pnt2d P1, P2;
235 if (Abs(Dy) < RealEpsilon()) { return;}
236 if (Abs(R - r*Dx/Dy) < RealEpsilon()) return;
237
238 v2 = (R + r*Dx/Dy)/(R - r*Dx/Dy);
239 if (v2 > 0.0) U2 = Log(Sqrt(v2));
240 P2 = ElCLib::Value(U2, C2);
241
242 U1 = (gp_Vec2d(C1.Location(), P2)).Dot(D);
243 P1 = ElCLib::Value(U1, C1);
244 mySqDist[myNbExt] = P1.SquareDistance(P2);
245 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
246 myPoint[myNbExt][1] = Extrema_POnCurv2d(U2,P2);
247 myNbExt++;
248 myDone = Standard_True;
249}
250
251
252
253//============================================================================
254
255Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Lin2d& C1, const gp_Parab2d& C2)
256{
257 myIsPar = Standard_False;
258 myDone = Standard_False;
259 myNbExt = 0;
260
0d969553 261// Calculate T1 in the reference of the parabole ...
7fd59977 262 gp_Dir2d D = C1.Direction();
263 gp_Dir2d x2, y2;
264 x2 = C2.MirrorAxis().Direction();
265 y2 = C2.Axis().YAxis().Direction();
266 Standard_Real Dx = D.Dot(x2);
267 Standard_Real Dy = D.Dot(y2);
268
269 Standard_Real U1, U2, P = C2.Parameter();
270 gp_Pnt2d P1, P2;
271 if (Abs(Dy) < RealEpsilon()) { return; }
272 U2 = Dx*P/Dy;
273 P2 = ElCLib::Value(U2, C2);
274
275 U1 = (gp_Vec2d(C1.Location(), P2)).Dot(D);
276 P1 = ElCLib::Value(U1, C1);
277 mySqDist[myNbExt] = P1.SquareDistance(P2);
278 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1,P1);
279 myPoint[myNbExt][1] = Extrema_POnCurv2d(U2,P2);
280 myNbExt++;
281 myDone = Standard_True;
282}
283
284
285
286//============================================================================
287
288Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Circ2d& C1, const gp_Circ2d& C2)
289{
290 myIsPar = Standard_False;
291 myDone = Standard_False;
292 myNbExt = 0;
293 myDone = Standard_True;
294
295 gp_Pnt2d O1 = C1.Location();
296 gp_Pnt2d O2 = C2.Location();
297
298 gp_Vec2d DO1O2 (O1, O2);
638ad7f3 299 const Standard_Real aSqDCenters = DO1O2.SquareMagnitude();
300 if (aSqDCenters < Precision::SquareConfusion()) {
7fd59977 301 myIsPar = Standard_True;
638ad7f3 302 myNbExt = 1;
303 myDone = Standard_True;
304 const Standard_Real aDR = C1.Radius() - C2.Radius();
305 mySqDist[0] = aDR*aDR;
306 return;
7fd59977 307 }
308
309 Standard_Integer NoSol, kk;
310 Standard_Real U1, U2;
311 Standard_Real r1 = C1.Radius(), r2 = C2.Radius();
312 Standard_Real Usol2[2], Usol1[2];
313 gp_Pnt2d P1[2], P2[2];
638ad7f3 314 gp_Vec2d O1O2(DO1O2/Sqrt(aSqDCenters));
7fd59977 315
316 P1[0] = O1.Translated(r1*O1O2);
317 Usol1[0] = ElCLib::Parameter(C1, P1[0]);
318 P1[1] = O1.Translated(-r1*O1O2);
319 Usol1[1] = ElCLib::Parameter(C1, P1[1]);
320
321 P2[0] = O2.Translated(r2*O1O2);
322 Usol2[0] = ElCLib::Parameter(C2, P2[0]);
323 P2[1] = O2.Translated(-r2*O1O2);
324 Usol2[1] = ElCLib::Parameter(C2, P2[1]);
325
326 for (NoSol = 0; NoSol <= 1; NoSol++) {
327 U1 = Usol1[NoSol];
328 for (kk = 0; kk <= 1; kk++) {
329 U2 = Usol2[kk];
330 mySqDist[myNbExt] = P2[kk].SquareDistance(P1[NoSol]);
331 myPoint[myNbExt][0] = Extrema_POnCurv2d(U1, P1[NoSol]);
332 myPoint[myNbExt][1] = Extrema_POnCurv2d(U2, P2[kk]);
333 myNbExt++;
334 }
335 }
336}
337//===========================================================================
338
339Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Circ2d& C1, const gp_Elips2d& C2)
340{
341 myIsPar = Standard_False;
342 myDone = Standard_False;
343 myNbExt = 0;
344
345 Standard_Integer i, j;
346
347 Extrema_ExtPElC2d ExtElip(C1.Location(), C2,
c6541a0c 348 Precision::Confusion(), 0.0, 2.0*M_PI);
7fd59977 349
350 if (ExtElip.IsDone()) {
351 for (i = 1; i <= ExtElip.NbExt(); i++) {
352 Extrema_ExtPElC2d ExtCirc(ExtElip.Point(i).Value(), C1,
c6541a0c 353 Precision::Confusion(), 0.0, 2.0*M_PI);
7fd59977 354 if (ExtCirc.IsDone()) {
355 for (j = 1; j <= ExtCirc.NbExt(); j++) {
356 mySqDist[myNbExt] = ExtCirc.SquareDistance(j);
357 myPoint[myNbExt][0] = ExtCirc.Point(j);
358 myPoint[myNbExt][1] = ExtElip.Point(i);
359 myNbExt++;
360 }
361 }
362 myDone = Standard_True;
363 }
364 }
365}
366//============================================================================
367
368Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Circ2d& C1, const gp_Hypr2d& C2)
369{
370 myIsPar = Standard_False;
371 myDone = Standard_False;
372 myNbExt = 0;
373
374 Standard_Integer i, j;
375
376 Extrema_ExtPElC2d ExtHyp(C1.Location(), C2, Precision::Confusion(),
377 RealFirst(), RealLast());
378
379 if (ExtHyp.IsDone()) {
380 for (i = 1; i <= ExtHyp.NbExt(); i++) {
381 Extrema_ExtPElC2d ExtCirc(ExtHyp.Point(i).Value(), C1,
c6541a0c 382 Precision::Confusion(), 0.0, 2.0*M_PI);
7fd59977 383 if (ExtCirc.IsDone()) {
384 for (j = 1; j <= ExtCirc.NbExt(); j++) {
385 mySqDist[myNbExt] = ExtCirc.SquareDistance(j);
386 myPoint[myNbExt][0] = ExtCirc.Point(j);
387 myPoint[myNbExt][1] = ExtHyp.Point(i);
388 myNbExt++;
389 }
390 }
391 myDone = Standard_True;
392 }
393 }
394}
395//============================================================================
396
397Extrema_ExtElC2d::Extrema_ExtElC2d (const gp_Circ2d& C1, const gp_Parab2d& C2)
398{
399 myIsPar = Standard_False;
400 myDone = Standard_False;
401 myNbExt = 0;
402
403 Standard_Integer i, j;
404
405 Extrema_ExtPElC2d ExtParab(C1.Location(), C2, Precision::Confusion(),
406 RealFirst(), RealLast());
407
408 if (ExtParab.IsDone()) {
409 for (i = 1; i <= ExtParab.NbExt(); i++) {
410 Extrema_ExtPElC2d ExtCirc(ExtParab.Point(i).Value(),
c6541a0c 411 C1, Precision::Confusion(), 0.0, 2.0*M_PI);
7fd59977 412 if (ExtCirc.IsDone()) {
413 for (j = 1; j <= ExtCirc.NbExt(); j++) {
414 mySqDist[myNbExt] = ExtCirc.SquareDistance(j);
415 myPoint[myNbExt][0] = ExtCirc.Point(j);
416 myPoint[myNbExt][1] = ExtParab.Point(i);
417 myNbExt++;
418 }
419 }
420 myDone = Standard_True;
421 }
422 }
423}
424//============================================================================
425
7fd59977 426Standard_Boolean Extrema_ExtElC2d::IsDone () const { return myDone; }
427//============================================================================
428
429Standard_Boolean Extrema_ExtElC2d::IsParallel () const
430{
9775fa61 431 if (!IsDone()) { throw StdFail_NotDone(); }
7fd59977 432 return myIsPar;
433}
434//============================================================================
435
436Standard_Integer Extrema_ExtElC2d::NbExt () const
437{
638ad7f3 438 if (!IsDone())
439 {
440 throw StdFail_NotDone();
441 }
442
7fd59977 443 return myNbExt;
444}
445//============================================================================
446
447Standard_Real Extrema_ExtElC2d::SquareDistance (const Standard_Integer N) const
448{
638ad7f3 449 if (N < 1 || N > NbExt())
450 {
451 throw Standard_OutOfRange();
7fd59977 452 }
638ad7f3 453
454 return mySqDist[N - 1];
7fd59977 455}
456//============================================================================
457
458void Extrema_ExtElC2d::Points (const Standard_Integer N,
459 Extrema_POnCurv2d& P1,
460 Extrema_POnCurv2d& P2) const
461{
638ad7f3 462 if (IsParallel())
463 {
464 throw StdFail_InfiniteSolutions();
465 }
466
9775fa61 467 if (N < 1 || N > NbExt()) { throw Standard_OutOfRange(); }
7fd59977 468 P1 = myPoint[N-1][0];
469 P2 = myPoint[N-1][1];
470}
471//============================================================================
Open CASCADE Technology repository