projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031642: Visualization - crash in Graphic3d_Structure::SetVisual() on redisplaying...
[occt.git] / src / GccAna / GccAna_Circ2d3Tan.cxx
CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
973c2be1 2// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 3//
973c2be1 4// This file is part of Open CASCADE Technology software library.
b311480e 5//
d5f74e42 6// This library is free software; you can redistribute it and/or modify it under
7// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 8// by the Free Software Foundation, with special exception defined in the file
9// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10// distribution for complete text of the license and disclaimer of any warranty.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
7fd59977 14
42cf5bc1 15
7fd59977 16#include <ElCLib.hxx>
42cf5bc1 17#include <GccAna_Circ2d3Tan.hxx>
7fd59977 18#include <GccAna_Circ2dBisec.hxx>
42cf5bc1 19#include <GccEnt_BadQualifier.hxx>
20#include <GccEnt_QualifiedCirc.hxx>
21#include <GccEnt_QualifiedLin.hxx>
7fd59977 22#include <GccInt_BCirc.hxx>
7fd59977 23#include <GccInt_BElips.hxx>
24#include <GccInt_BHyper.hxx>
42cf5bc1 25#include <GccInt_BLine.hxx>
26#include <GccInt_IType.hxx>
27#include <gp_Circ2d.hxx>
28#include <gp_Dir2d.hxx>
29#include <gp_Lin2d.hxx>
30#include <gp_Pnt2d.hxx>
31#include <IntAna2d_AnaIntersection.hxx>
32#include <IntAna2d_Conic.hxx>
33#include <IntAna2d_IntPoint.hxx>
34#include <math_DirectPolynomialRoots.hxx>
7fd59977 35#include <Standard_OutOfRange.hxx>
42cf5bc1 36#include <StdFail_NotDone.hxx>
7fd59977 37#include <TColStd_Array1OfReal.hxx>
7fd59977 38
39//=========================================================================
0d969553 40// Creation of a circle tangent to three circles. +
7fd59977 41//=========================================================================
7fd59977 42GccAna_Circ2d3Tan::
43 GccAna_Circ2d3Tan (const GccEnt_QualifiedCirc& Qualified1,
44 const GccEnt_QualifiedCirc& Qualified2,
45 const GccEnt_QualifiedCirc& Qualified3,
46 const Standard_Real Tolerance ):
47
48//=========================================================================
0d969553 49// Initialization of fields. +
7fd59977 50//=========================================================================
51
52 cirsol(1,16) ,
53 qualifier1(1,16),
54 qualifier2(1,16),
55 qualifier3(1,16),
56 TheSame1(1,16) ,
57 TheSame2(1,16) ,
58 TheSame3(1,16) ,
59 pnttg1sol(1,16),
60 pnttg2sol(1,16),
61 pnttg3sol(1,16),
62 par1sol(1,16) ,
63 par2sol(1,16) ,
64 par3sol(1,16) ,
65 pararg1(1,16) ,
66 pararg2(1,16) ,
67 pararg3(1,16)
68{
69
70 gp_Dir2d dirx(1.0,0.0);
71 Standard_Real Tol = Abs(Tolerance);
72 WellDone = Standard_False;
73 NbrSol = 0;
74 if (!(Qualified1.IsEnclosed() || Qualified1.IsEnclosing() ||
75 Qualified1.IsOutside() || Qualified1.IsUnqualified()) ||
76 !(Qualified2.IsEnclosed() || Qualified2.IsEnclosing() ||
77 Qualified2.IsOutside() || Qualified2.IsUnqualified()) ||
78 !(Qualified3.IsEnclosed() || Qualified3.IsEnclosing() ||
79 Qualified3.IsOutside() || Qualified3.IsUnqualified())) {
9775fa61 80 throw GccEnt_BadQualifier();
7fd59977 81 return;
82 }
83
84//=========================================================================
0d969553 85// Processing. +
7fd59977 86//=========================================================================
87
88 gp_Circ2d Cir1 = Qualified1.Qualified();
89 gp_Circ2d Cir2 = Qualified2.Qualified();
90 gp_Circ2d Cir3 = Qualified3.Qualified();
91 Standard_Real R1 = Cir1.Radius();
92 Standard_Real R2 = Cir2.Radius();
93 Standard_Real R3 = Cir3.Radius();
94 gp_Pnt2d center1(Cir1.Location());
95 gp_Pnt2d center2(Cir2.Location());
96 gp_Pnt2d center3(Cir3.Location());
97
98 Standard_Real X1 = center1.X();
99 Standard_Real X2 = center2.X();
100 Standard_Real X3 = center3.X();
101
102 Standard_Real Y1 = center1.Y();
103 Standard_Real Y2 = center2.Y();
104 Standard_Real Y3 = center3.Y();
105
106 gp_XY dir2 = center1.XY() - center2.XY();
107 gp_XY dir3 = center1.XY() - center3.XY();
108
109//////////
110 if ((Abs(R1 - R2) <= Tolerance && center1.IsEqual(center2, Tolerance)) ||
111 (Abs(R1 - R3) <= Tolerance && center1.IsEqual(center3, Tolerance)) ||
112 (Abs(R2 - R3) <= Tolerance && center2.IsEqual(center3, Tolerance)))
113 return;
114 else {
115 if (Abs(dir2^dir3) <= Tolerance) {
116 Standard_Real Dist1 = center1.Distance(center2);
117 Standard_Real Dist2 = center1.Distance(center3);
118 Standard_Real Dist3 = center2.Distance(center3);
119 if (Abs(Abs(R1 - R2) - Dist1) <= Tolerance) {
120 if (Abs(Abs(R1 - R3) - Dist2) <= Tolerance) {
121 if (Abs(Abs(R2 - R3) - Dist3) <= Tolerance)
122 return;
123 } else if (Abs(R1 + R3 - Dist2) <= Tolerance) {
124 if (Abs(R2 + R3 - Dist3) <= Tolerance)
125 return;
126 }
127 } else if (Abs(R1 + R2 - Dist1) <= Tolerance) {
128 if (Abs(Abs(R1 - R3) - Dist2) <= Tolerance &&
129 Abs(R2 + R3 - Dist3) <= Tolerance) {
130 } else {
131 if (Abs(Abs(R2 - R3) - Dist3) <= Tolerance &&
132 Abs(R1 + R3 - Dist2) <= Tolerance)
133 return;
134 }
135 }
136 }
137 }
138/////////
139 TColStd_Array1OfReal A2(1, 8), B2(1, 8), C2(1, 8), D2(1, 8), E2(1, 8), F2(1, 8);
140 TColStd_Array1OfReal A3(1, 8), B3(1, 8), C3(1, 8), D3(1, 8), E3(1, 8), F3(1, 8);
141 TColStd_Array1OfReal Beta2(1, 8), Gamma2(1, 8), Delta2(1, 8);
142 TColStd_Array1OfReal Beta3(1, 8), Gamma3(1, 8), Delta3(1, 8);
143 Standard_Real a2, b2, c2, d2, e2, f2;
144 Standard_Real a3, b3, c3, d3, e3, f3;
145 Standard_Real A, B, C, D, E;
146 Standard_Boolean IsSame;
147 Standard_Boolean IsTouch;
148 Standard_Integer FirstIndex;
149
150 Standard_Integer i, j, k, l;
151 TColStd_Array1OfReal xSol(1, 64);
152 TColStd_Array1OfReal ySol(1, 64);
153 TColStd_Array1OfReal rSol(1, 16);
154 TColStd_Array1OfInteger FirstSol(1, 9);
155 TColStd_Array1OfReal xSol1(1, 32);
156 TColStd_Array1OfReal ySol1(1, 32);
157 TColStd_Array1OfReal rSol1(1, 32);
158 TColStd_Array1OfInteger FirstSol1(1, 9);
159 Standard_Real x, y, r;
160 Standard_Real m, n, t, s, v;
161 Standard_Real p, q;
162 Standard_Real Epsilon;
163
164 Standard_Integer CurSol;
165
166//*********************************************************************************************
167//*********************************************************************************************
168
169// Actually we have to find solutions of eight systems of equations:
170// _ _
171// | (X - X1)2 + (Y - Y1)2 = (R - R1)2 | (X - X1)2 + (Y - Y1)2 = (R + R1)2
172// 1) < (X - X2)2 + (Y - Y2)2 = (R - R2)2 2) < (X - X2)2 + (Y - Y2)2 = (R - R2)2
173// \_(X - X3)2 + (Y - Y3)2 = (R - R3)2 \_(X - X3)2 + (Y - Y3)2 = (R - R3)2
174// _ _
175// | (X - X1)2 + (Y - Y1)2 = (R - R1)2 | (X - X1)2 + (Y - Y1)2 = (R - R1)2
176// 3) < (X - X2)2 + (Y - Y2)2 = (R + R2)2 4) < (X - X2)2 + (Y - Y2)2 = (R - R2)2
177// \_(X - X3)2 + (Y - Y3)2 = (R - R3)2 \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
178// _ _
179// | (X - X1)2 + (Y - Y1)2 = (R + R1)2 | (X - X1)2 + (Y - Y1)2 = (R + R1)2
180// 5) < (X - X2)2 + (Y - Y2)2 = (R + R2)2 6) < (X - X2)2 + (Y - Y2)2 = (R - R2)2
181// \_(X - X3)2 + (Y - Y3)2 = (R - R3)2 \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
182// _ _
183// | (X - X1)2 + (Y - Y1)2 = (R - R1)2 | (X - X1)2 + (Y - Y1)2 = (R + R1)2
184// 7) < (X - X2)2 + (Y - Y2)2 = (R + R2)2 8) < (X - X2)2 + (Y - Y2)2 = (R + R2)2
185// \_(X - X3)2 + (Y - Y3)2 = (R + R3)2 \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
186
187// each equation (X - Xi)2 + (Y - Yi)2 = (R +- Ri)2 means that the circle (X,Y,R) is tangent
188// to the circle (Xi,Yi,Ri).
189
190// The number of each system is very important. Further index i shows the numer of the system
191
192// Further Beta, Gamma and Delta are coefficients of the equation:
193// R +- Ri = Beta*X + Gamma*Y + Delta where i=2 or i=3
194
195//*********************************************************************************************
196//*********************************************************************************************
197
198// Verification do two circles touch each other or not
199// if at least one circle touches other one IsTouch become Standard_Standard_True
200
201 if (Abs((X1 - X2)*(X1 - X2) + (Y1 - Y2)*(Y1 - Y2) - (R1 - R2)*(R1 - R2)) <= Tolerance ||
202 Abs((X1 - X2)*(X1 - X2) + (Y1 - Y2)*(Y1 - Y2) - (R1 + R2)*(R1 + R2)) <= Tolerance ||
203 Abs((X1 - X3)*(X1 - X3) + (Y1 - Y3)*(Y1 - Y3) - (R1 - R3)*(R1 - R3)) <= Tolerance ||
204 Abs((X1 - X3)*(X1 - X3) + (Y1 - Y3)*(Y1 - Y3) - (R1 + R3)*(R1 + R3)) <= Tolerance ||
205 Abs((X2 - X3)*(X2 - X3) + (Y2 - Y3)*(Y2 - Y3) - (R2 - R3)*(R2 - R3)) <= Tolerance ||
206 Abs((X2 - X3)*(X2 - X3) + (Y2 - Y3)*(Y2 - Y3) - (R2 + R3)*(R2 + R3)) <= Tolerance)
207 IsTouch = Standard_True;
208 else
209 IsTouch = Standard_False;
210
211// First step:
212// We are searching for Beta, Gamma and Delta coefficients
213// and also coefficients of the system of second order equations:
214// _
215// | a2*x*x +2*b2*x*y + c2*y*y +2*d2*x + 2*e2*y + f2 = 0
216// <
217// \_ a3*x*x +2*b3*x*y + c3*y*y +2*d3*x + 2*e3*y + f3 = 0 ,
218
219// obtained by exclusion of R from source systems.
220
221 for (i = 1; i <= 8; i++) {
222
223// _
224// | (X - X1)2 + (Y - Y1)2 = (R +- R1)2
225// <
226// \_(X - X2)2 + (Y - Y2)2 = (R +- R2)2
227
228 if (i == 1 || i == 4 || i == 5 || i == 8) {
229 if (Abs(R1 - R2) > Tolerance) {
230 Beta2(i) = (X1 - X2)/(R1 - R2);
231 Gamma2(i) = (Y1 - Y2)/(R1 - R2);
232 Delta2(i) = (X2*X2 - X1*X1 + Y2*Y2 - Y1*Y1 + (R1 - R2)*(R1 - R2))/(2*(R1 - R2));
233 }
234 } else {
235 Beta2(i) = (X1 - X2)/(R1 + R2);
236 Gamma2(i) = (Y1 - Y2)/(R1 + R2);
237 Delta2(i) = (X2*X2 - X1*X1 + Y2*Y2 - Y1*Y1 + (R1 + R2)*(R1 + R2))/(2*(R1 + R2));
238 }
239 if ((i == 1 || i == 4 || i == 5 || i == 8) &&
240 (Abs(R1 - R2) <= Tolerance)) {
241// If R1 = R2
242 A2(i) = 0.;
243 B2(i) = 0.;
244 C2(i) = 0.;
245 D2(i) = X2 - X1;
246 E2(i) = Y2 - Y1;
247 F2(i) = X1*X1 - X2*X2 + Y1*Y1 - Y2*Y2;
248 } else {
249 A2(i) = Beta2(i)*Beta2(i) - 1.;
250 B2(i) = Beta2(i)*Gamma2(i);
251 C2(i) = Gamma2(i)*Gamma2(i) - 1.;
252 D2(i) = Beta2(i)*Delta2(i) + X2;
253 E2(i) = Gamma2(i)*Delta2(i) + Y2;
254 F2(i) = Delta2(i)*Delta2(i) - X2*X2 - Y2*Y2;
255 }
256
257// _
258// | (X - X1)2 + (Y - Y1)2 = (R +- R1)2
259// <
260// \_(X - X3)2 + (Y - Y3)2 = (R +- R3)2
261
262 if (i == 1 || i == 3 || i == 6 || i == 8) {
263 if (Abs(R1 - R3) > Tolerance) {
264 Beta3(i) = (X1 - X3)/(R1 - R3);
265 Gamma3(i) = (Y1 - Y3)/(R1 - R3);
266 Delta3(i) = (X3*X3 - X1*X1 + Y3*Y3 - Y1*Y1 + (R1 - R3)*(R1 - R3))/(2*(R1 - R3));
267 }
268 } else {
269 Beta3(i) = (X1 - X3)/(R1 + R3);
270 Gamma3(i) = (Y1 - Y3)/(R1 + R3);
271 Delta3(i) = (X3*X3 - X1*X1 + Y3*Y3 - Y1*Y1 + (R1 + R3)*(R1 + R3))/(2*(R1 + R3));
272 }
273 if ((i == 1 || i == 3 || i == 6 || i == 8) &&
274 (Abs(R1 - R3) <= Tolerance)) {
275 A3(i) = 0.;
276 B3(i) = 0.;
277 C3(i) = 0.;
278 D3(i) = X3 - X1;
279 E3(i) = Y3 - Y1;
280 F3(i) = X1*X1 - X3*X3 + Y1*Y1 - Y3*Y3;
281 } else {
282 A3(i) = Beta3(i)*Beta3(i) - 1.;
283 B3(i) = Beta3(i)*Gamma3(i);
284 C3(i) = Gamma3(i)*Gamma3(i) - 1.;
285 D3(i) = Beta3(i)*Delta3(i) + X3;
286 E3(i) = Gamma3(i)*Delta3(i) + Y3;
287 F3(i) = Delta3(i)*Delta3(i) - X3*X3 - Y3*Y3;
288 }
289 }
290
291// Second step:
292// We are searching for the couple (X,Y) as a solution of the system:
293// _
294// | a2*x*x +2*b2*x*y + c2*y*y +2*d2*x + 2*e2*y + f2 = 0
295// <
296// \_ a3*x*x +2*b3*x*y + c3*y*y +2*d3*x + 2*e3*y + f3 = 0
297
298 CurSol = 1;
299 for (i = 1; i <= 8; i++) {
300 a2 = A2(i); a3 = A3(i);
301 b2 = B2(i); b3 = B3(i);
302 c2 = C2(i); c3 = C3(i);
303 d2 = D2(i); d3 = D3(i);
304 e2 = E2(i); e3 = E3(i);
305 f2 = F2(i); f3 = F3(i);
306
307 FirstSol(i) = CurSol;
308
309// In some cases we know that some systems have no solution in any case due to qualifiers
310 if (((i == 2 || i == 5 || i == 6 || i == 8) &&
311 (Qualified1.IsEnclosed() || Qualified1.IsEnclosing())) ||
312 ((i == 1 || i == 3 || i == 4 || i == 7) && Qualified1.IsOutside()))
313 continue;
314
315 if (((i == 3 || i == 5 || i == 7 || i == 8) &&
316 (Qualified2.IsEnclosed() || Qualified2.IsEnclosing())) ||
317 ((i == 1 || i == 2 || i == 4 || i == 6) && Qualified2.IsOutside()))
318 continue;
319
320 if (((i == 4 || i == 6 || i == 7 || i == 8) &&
321 (Qualified3.IsEnclosed() || Qualified3.IsEnclosing())) ||
322 ((i == 1 || i == 2 || i == 3 || i == 5) && Qualified3.IsOutside()))
323 continue;
324
325// Check is Cir1 a solution of this system or not
326// In that case equations are equal to each other
327 if (Abs(a2 - a3) <= Tolerance && Abs(b2 - b3) <= Tolerance && Abs(c2 - c3) <= Tolerance &&
328 Abs(d2 - d3) <= Tolerance && Abs(e2 - e3) <= Tolerance && Abs(f2 - f3) <= Tolerance) {
329 xSol(CurSol) = X1;
330 ySol(CurSol) = Y1;
331 CurSol++;
332 continue;
333 }
334// 1) a2 = 0
335 if (Abs(a2) <= Tolerance) {
336
337// 1.1) b2y + d2 = 0
338// Searching for solution of the equation Ay2 + By + C = 0
339 A = c2; B = 2.*e2; C = f2;
340 math_DirectPolynomialRoots yRoots(A, B, C);
341 if (yRoots.IsDone() && !yRoots.InfiniteRoots())
342 for (k = 1; k <= yRoots.NbSolutions(); k++) {
343// for each y solution:
344 y = yRoots.Value(k);
345// Searching for solution of the equation Ax2 + Bx + C = 0
346 if (!(k == 2 && Abs(y - yRoots.Value(1)) <= 10*Tolerance) &&
347 Abs(b2*y + d2) <= b2*Tolerance) {
348 A = a3; B = 2*(b3*y + d3); C = c3*(y*y) + 2*e3*y + f3;
349 math_DirectPolynomialRoots xRoots(A, B, C);
350 if (xRoots.IsDone() && !xRoots.InfiniteRoots())
351 for (j = 1; j <= xRoots.NbSolutions(); j++) {
352 x = xRoots.Value(j);
353 if (!(j == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
354 xSol(CurSol) = x;
355 ySol(CurSol) = y;
356 CurSol++;
357 }
358 }
359 }
360 }
361
362// 1.2) b2y + d2 != 0
363 A = a3*c2*c2 - 4*b2*(b3*c2 - b2*c3);
364 B = 4*a3*c2*e2 - 4*b3*(c2*d2 + 2*b2*e2) + 4*b2*(2*c3*d2 - c2*d3 + 2*b2*e3);
365 C = 2*a3*(c2*f2 + 2*e2*e2) - 4*b3*(b2*f2 + 2*e2*d2) + 4*c3*d2*d2 - 4*d3*(c2*d2 + 2*b2*e2)
366 + 16*b2*e3*d2 + 4*b2*b2*f3;
367 D = 4*a3*e2*f2 - 4*b3*d2*f2 - 4*d3*(b2*f2 + 2*d2*e2) + 8*d2*d2*e3 + 8*b2*d2*f3;
368 E = a3*f2*f2 - 4*d2*d3*f2 + 4*d2*d2*f3;
369
370// Searching for solution of the equation Ay4 + By3 + Cy2 + Dy + E = 0
371// Special case: one circle touches other
372 if (IsTouch) {
373// Derivation of the equation Ay4 + By3 + Cy2 + Dy + E
374 math_DirectPolynomialRoots yRoots1(4*A, 3*B, 2*C, D);
375 if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
376 for (k = 1; k <= yRoots1.NbSolutions(); k++) {
377 y = yRoots1.Value(k);
378// Check if this value is already catched
379 IsSame = Standard_False;
380 for (l = 1; l < k; l++)
381 if (Abs(y - yRoots1.Value(l)) <= 10*Tolerance) IsSame = Standard_True;
382
383 Epsilon = (Abs((Abs((Abs(4*A*y) + Abs(3*B))*y) + Abs(2*C))*y) + Abs(D));
384 if (Abs((((A*y + B)*y + C)*y + D)*y + E) <= Epsilon*Tolerance) {
385 if (!IsSame && Abs(b2*y + d2) > b2*Tolerance) {
386 x = -(c2*(y*y) + 2*e2*y + f2)/(2*(b2*y + d2));
387 xSol(CurSol) = x;
388 ySol(CurSol) = y;
389 CurSol++;
390 }
391 }
392 }
393 }
394
395 math_DirectPolynomialRoots yRoots1(A, B, C, D, E);
396 if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
397 for (k = 1; k <= yRoots1.NbSolutions(); k++) {
398 y = yRoots1.Value(k);
399// Check if this value is already catched
400 IsSame = Standard_False;
401 FirstIndex = (i == 1) ? 1 : FirstSol(i);
402 for (l = FirstIndex; l < CurSol; l++)
403 if (Abs(y - ySol(l)) <= 10*Tolerance) IsSame = Standard_True;
404
405 if (!IsSame && Abs(b2*y + d2) > b2*Tolerance) {
406 x = -(c2*(y*y) + 2*e2*y + f2)/(2*(b2*y + d2));
407 xSol(CurSol) = x;
408 ySol(CurSol) = y;
409 CurSol++;
410 }
411 }
412 } else {
413// 2) a2 != 0
414// Coefficients of the equation (sy + v)Sqrt(p2 - q) = (my2 + ny + t)
415 m = 2*a3*b2*b2/(a2*a2) - 2*b2*b3/a2 - a3*c2/a2 + c3;
416 n = 4*a3*b2*d2/(a2*a2) - 2*b3*d2/a2 - 2*b2*d3/a2 - 2*a3*e2/a2 + 2*e3;
417 t = 2*a3*d2*d2/(a2*a2) - 2*d2*d3/a2 - a3*f2/a2 + f3;
418 s = 2*b3 - 2*a3*b2/a2;
419 v = 2*d3 - 2*d2*a3/a2;
420
421//------------------------------------------
422// If s = v = 0
423 if (Abs(s) <= Tolerance && Abs(v) <= Tolerance) {
424 math_DirectPolynomialRoots yRoots(m, n, t);
425 if (yRoots.IsDone() && !yRoots.InfiniteRoots())
426 for (k = 1; k <= yRoots.NbSolutions(); k++) {
427// for each y solution:
428 y = yRoots.Value(k);
429
430 p = -(b2*y + d2)/a2;
431 q = (c2*(y*y) + 2*e2*y + f2)/a2;
432 Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
433 if (!(k == 2 && Abs(y - yRoots.Value(1)) <= 10*Tolerance) &&
434 p*p - q >= -Epsilon*Tolerance) {
435 A = a2;
436 B = 2*(b2*y + d2);
437 C = c2*y*y + 2*e2*y + f2;
438 math_DirectPolynomialRoots xRoots(A, B, C);
439 if (xRoots.IsDone() && !xRoots.InfiniteRoots())
440 for (l = 1; l <= xRoots.NbSolutions(); l++) {
441// for each x solution:
442 x = xRoots.Value(l);
443
444 if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
445 xSol(CurSol) = x;
446 ySol(CurSol) = y;
447 CurSol++;
448 }
449 }
450 }
451 }
452 } else {
453//------------------------------------------
454// If (s*y + v) != 0
455
456 A = s*s*(b2*b2 - a2*c2) - m*m*a2*a2;
457 B = 2*s*v*(b2*b2 - a2*c2) + 2*s*s*(b2*d2 - a2*e2) - 2*m*n*a2*a2;
458 C = v*v*(b2*b2 - a2*c2) + 4*s*v*(b2*d2 - a2*e2) + s*s*(d2*d2 - a2*f2) - a2*a2*(2*m*t + n*n);
459 D = 2*v*v*(b2*d2 - a2*e2) + 2*s*v*(d2*d2 - a2*f2) - 2*n*t*a2*a2;
460 E = v*v*(d2*d2 - a2*f2) - t*t*a2*a2;
461
462// Searching for solution of the equation Ay4 + By3 + Cy2 + Dy + E = 0
463// Special case: one circle touches other
464 if (IsTouch) {
465// Derivation of the equation Ay4 + By3 + Cy2 + Dy + E
466 math_DirectPolynomialRoots yRoots1(4*A, 3*B, 2*C, D);
467 if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
468 for (k = 1; k <= yRoots1.NbSolutions(); k++) {
469 y = yRoots1.Value(k);
470
471 p = -(b2*y + d2)/a2;
472 q = (c2*(y*y) + 2*e2*y + f2)/a2;
473
474// Check if this value is already catched
475 IsSame = Standard_False;
476 FirstIndex = (i == 1) ? 1 : FirstSol(i);
477 for (l = FirstIndex; l < CurSol; l++)
478 if (Abs(y - ySol(l)) <= 10*Tolerance) IsSame = Standard_True;
479
480 Epsilon = (Abs((Abs((Abs(4*A*y) + Abs(3*B))*y) + Abs(2*C))*y) + Abs(D));
481 if (Abs((((A*y + B)*y + C)*y + D)*y + E) <= Epsilon*Tolerance) {
482
483 Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
484 if (!IsSame && p*p - q >= -Epsilon*Tolerance) {
485 A = a2;
486 B = 2*(b2*y + d2);
487 C = c2*y*y + 2*e2*y + f2;
488 math_DirectPolynomialRoots xRoots(A, B, C);
489 if (xRoots.IsDone() && !xRoots.InfiniteRoots())
490 for (l = 1; l <= xRoots.NbSolutions(); l++) {
491// for each x solution:
492 x = xRoots.Value(l);
493
494 if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
495 xSol(CurSol) = x;
496 ySol(CurSol) = y;
497 CurSol++;
498 }
499 }
500 }
501 }
502 }
503 }
504
505 math_DirectPolynomialRoots yRoots(A, B, C, D, E);
506 if (yRoots.IsDone() && !yRoots.InfiniteRoots())
507 for (k = 1; k <= yRoots.NbSolutions(); k++) {
508// for each y solution:
509 y = yRoots.Value(k);
510
511 p = -(b2*y + d2)/a2;
512 q = (c2*(y*y) + 2*e2*y + f2)/a2;
513
514// Check if this value is already catched
515 IsSame = Standard_False;
516 for (l = 1; l < k; l++)
517 if (Abs(y - yRoots.Value(l)) <= 10*Tolerance) IsSame = Standard_True;
518
519 Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
520 if (!IsSame && p*p - q >= -Epsilon*Tolerance) {
521 A = a2;
522 B = 2*(b2*y + d2);
523 C = c2*y*y + 2*e2*y + f2;
524 math_DirectPolynomialRoots xRoots(A, B, C);
525 if (xRoots.IsDone() && !xRoots.InfiniteRoots())
526 for (l = 1; l <= xRoots.NbSolutions(); l++) {
527// for each x solution:
528 x = xRoots.Value(l);
529
530 if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
531 xSol(CurSol) = x;
532 ySol(CurSol) = y;
533 CurSol++;
534 }
535 }
536 }
537 }
538 }
539 }
540 }
541 FirstSol(9) = CurSol;
542
543// Third step:
544// Check of couples (X,Y) and searching for R. R must be great than 0
545 CurSol = 1;
546 for (i = 1; i <= 8; i++) {
547 FirstSol1(i) = CurSol;
548 for (j = FirstSol(i); j < FirstSol(i + 1); j++) {
549 x = xSol(j);
550 y = ySol(j);
551// in some cases when R1 = R2 :
552 if ((i == 1 || i == 4 || i == 5 || i == 8) && (Abs(R1 - R2) <= Tolerance)) {
553 if (i == 1 || i == 4) {
554 r = R1 + Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
555 Epsilon = 10*(2*Abs(r - R2) + Abs(x - X2) + Abs(y - Y2));
556 if (Abs((r - R2)*(r - R2) - (x - X2)*(x - X2) - (y - Y2)*(y - Y2)) <=
557 Epsilon*Tolerance) {
558 xSol1(CurSol) = x;
559 ySol1(CurSol) = y;
560 rSol1(CurSol) = r;
561 CurSol++;
562 }
563 r = R1 - Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
564 Epsilon = 10*(2*Abs(r - R2) + Abs(x - X2) + Abs(y - Y2));
565 if ((r > Tolerance) &&
566 (Abs((r - R2)*(r - R2) - (x - X2)*(x - X2) - (y - Y2)*(y - Y2)) <=
567 Epsilon*Tolerance)) {
568 xSol1(CurSol) = x;
569 ySol1(CurSol) = y;
570 rSol1(CurSol) = r;
571 CurSol++;
572 }
573 } else {
574// i == 5 || i == 8
575 r = - R1 + Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
576 if (r > Tolerance) {
577 xSol1(CurSol) = x;
578 ySol1(CurSol) = y;
579 rSol1(CurSol) = r;
580 CurSol++;
581 }
582 }
583 } else {
584// Other cases
585 if (i == 1 || i == 4) {
586 r = R2 + Beta2(i)*x + Gamma2(i)*y + Delta2(i);
587 if (r > Tolerance) {
588 xSol1(CurSol) = x;
589 ySol1(CurSol) = y;
590 rSol1(CurSol) = r;
591 CurSol++;
592 }
593 }
594 if (i == 5 || i == 8) {
595 r = -R2 - Beta2(i)*x - Gamma2(i)*y - Delta2(i);
596 if (r > Tolerance) {
597 xSol1(CurSol) = x;
598 ySol1(CurSol) = y;
599 rSol1(CurSol) = r;
600 CurSol++;
601 }
602 }
603 if (i == 3 || i == 7) {
604 r = - R2 + Beta2(i)*x + Gamma2(i)*y + Delta2(i);
605 if (r > Tolerance) {
606 xSol1(CurSol) = x;
607 ySol1(CurSol) = y;
608 rSol1(CurSol) = r;
609 CurSol++;
610 }
611 }
612 if (i == 2 || i == 6) {
613 r = R2 - Beta2(i)*x - Gamma2(i)*y - Delta2(i);
614 if (r > Tolerance) {
615 xSol1(CurSol) = x;
616 ySol1(CurSol) = y;
617 rSol1(CurSol) = r;
618 CurSol++;
619 }
620 }
621 }
622 }
623 }
624 FirstSol1(9) = CurSol;
625// Fourth step
626// Check of triplets (X,Y,R).
627 CurSol = 1;
628 for (i = 1; i <= 8; i++) {
629 FirstSol(i) = CurSol;
630 for (j = FirstSol1(i); j < FirstSol1(i + 1); j++) {
631 x = xSol1(j);
632 y = ySol1(j);
633 r = rSol1(j);
634// in some cases when R1 = R3 :
635 if ((i == 1 || i == 3 || i == 6 || i == 8) && Abs(R1 - R3) <= Tolerance) {
636 if (i == 1 || i == 3) {
637 Epsilon = 10*(2*Abs(r - R3) + Abs(x - X3) + Abs(y - Y3));
638 if (Abs((r - R3)*(r - R3) - (x - X3)*(x - X3) - (y - Y3)*(y - Y3)) <=
639 Epsilon*Tolerance) {
640 xSol(CurSol) = x;
641 ySol(CurSol) = y;
642 rSol(CurSol) = r;
643 CurSol++;
644 }
645 } else {
646// i == 6 || i == 8
647 Epsilon = 10*(2*(r + R3) + Abs(x - X3) + Abs(y - Y3));
648 if (Abs((r + R3)*(r + R3) - (x - X3)*(x - X3) - (y - Y3)*(y - Y3)) <=
649 Epsilon*Tolerance) {
650 xSol(CurSol) = x;
651 ySol(CurSol) = y;
652 rSol(CurSol) = r;
653 CurSol++;
654 }
655 }
656 } else {
657// Other cases
658 Epsilon = 10*(Abs(Beta3(i)) + Abs(Gamma3(i)) + 1.);
659 if (i == 1 || i == 3)
660 if (Abs(R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i) - r) <= Epsilon*Tolerance) {
661 xSol(CurSol) = x;
662 ySol(CurSol) = y;
663 rSol(CurSol) = r;
664 CurSol++;
665 }
666 if (i == 6 || i == 8)
667 if (Abs(R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i) + r) <= Epsilon*Tolerance) {
668 xSol(CurSol) = x;
669 ySol(CurSol) = y;
670 rSol(CurSol) = r;
671 CurSol++;
672 }
673 if (i == 4 || i == 7)
674 if (Abs(Beta3(i)*x + Gamma3(i)*y + Delta3(i) - r - R3) <= Epsilon*Tolerance) {
675 xSol(CurSol) = x;
676 ySol(CurSol) = y;
677 rSol(CurSol) = r;
678 CurSol++;
679 }
680 if (i == 2 || i == 5)
681 if (Abs(r - R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i)) <= Epsilon*Tolerance) {
682 xSol(CurSol) = x;
683 ySol(CurSol) = y;
684 rSol(CurSol) = r;
685 CurSol++;
686 }
687 }
688 }
689 }
690 FirstSol(9) = CurSol;
691
692// Fifth step:
693// We have found all solutions. We have to calculate some parameters for each one.
694 for (i = 1 ; i <= 8; i++) {
695 for (j = FirstSol(i); j < FirstSol(i + 1); j++) {
696
697 if ((Qualified1.IsEnclosed() && rSol(j) > R1) ||
698 (Qualified1.IsEnclosing() && rSol(j) < R1))
699 continue;
700 if ((Qualified2.IsEnclosed() && rSol(j) > R2) ||
701 (Qualified2.IsEnclosing() && rSol(j) < R2))
702 continue;
703 if ((Qualified3.IsEnclosed() && rSol(j) > R3) ||
704 (Qualified3.IsEnclosing() && rSol(j) < R3))
705 continue;
706
707 NbrSol++;
708
709// RLE, avoid out of range
710 if (NbrSol > cirsol.Upper()) NbrSol = cirsol.Upper();
711
712 gp_Pnt2d Center = gp_Pnt2d(xSol(j), ySol(j));
713
714 cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),rSol(j));
715
716// ==========================================================
717 Standard_Real distcc1 = Center.Distance(center1);
718 if (!Qualified1.IsUnqualified())
719 qualifier1(NbrSol) = Qualified1.Qualifier();
720 else if (Abs(distcc1 + rSol(j) - R1) <= Tol)
721 qualifier1(NbrSol) = GccEnt_enclosed;
722 else if (Abs(distcc1 - R1 - rSol(j)) <= Tol)
723 qualifier1(NbrSol) = GccEnt_outside;
724 else qualifier1(NbrSol) = GccEnt_enclosing;
725
726 Standard_Real distcc2 = Center.Distance(center1);
727 if (!Qualified2.IsUnqualified())
728 qualifier2(NbrSol) = Qualified2.Qualifier();
729 else if (Abs(distcc2 + rSol(j) - R2) <= Tol)
730 qualifier2(NbrSol) = GccEnt_enclosed;
731 else if (Abs(distcc2 - R2 - rSol(j)) <= Tol)
732 qualifier2(NbrSol) = GccEnt_outside;
733 else qualifier2(NbrSol) = GccEnt_enclosing;
734
735 Standard_Real distcc3 = Center.Distance(center1);
736 if (!Qualified3.IsUnqualified())
737 qualifier3(NbrSol) = Qualified3.Qualifier();
738 else if (Abs(distcc3 + rSol(j) - R3) <= Tol)
739 qualifier3(NbrSol) = GccEnt_enclosed;
740 else if (Abs(distcc3 - R3 - rSol(j)) <= Tol)
741 qualifier3(NbrSol) = GccEnt_outside;
742 else qualifier3(NbrSol) = GccEnt_enclosing;
743
744// ==========================================================
745
746 if (Center.Distance(Cir1.Location()) <= Tolerance)
747 TheSame1(NbrSol) = 1;
748 else {
749 TheSame1(NbrSol) = 0;
750 gp_Dir2d dc;
751 if ((i == 2 || i == 5 || i == 6 || i == 8) || rSol(j) > R1)
752 dc.SetXY(Cir1.Location().XY() - Center.XY());
753 else
754 dc.SetXY(Center.XY() - Cir1.Location().XY());
755 pnttg1sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
756 par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
757 pnttg1sol(NbrSol));
758 pararg1(NbrSol)=ElCLib::Parameter(Cir1,pnttg1sol(NbrSol));
759 }
760
761 if (Center.Distance(Cir2.Location()) <= Tolerance)
762 TheSame2(NbrSol) = 1;
763 else {
764 TheSame2(NbrSol) = 0;
765 gp_Dir2d dc;
766 if ((i == 3 || i == 5 || i == 7 || i == 8) || rSol(j) > R2)
767 dc.SetXY(Cir2.Location().XY() - Center.XY());
768 else
769 dc.SetXY(Center.XY() - Cir2.Location().XY());
770 pnttg2sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
771 par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
772 pnttg2sol(NbrSol));
773 pararg2(NbrSol)=ElCLib::Parameter(Cir2,pnttg2sol(NbrSol));
774 }
775
776 if (Center.Distance(Cir3.Location()) <= Tolerance)
777 TheSame3(NbrSol) = 1;
778 else {
779 TheSame3(NbrSol) = 0;
780 gp_Dir2d dc;
781 if ((i == 4 || i == 6 || i == 7 || i == 8) || rSol(j) > R3)
782 dc.SetXY(Cir3.Location().XY() - Center.XY());
783 else
784 dc.SetXY(Center.XY() - Cir3.Location().XY());
785 pnttg3sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
786 par3sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
787 pnttg3sol(NbrSol));
788 pararg3(NbrSol)=ElCLib::Parameter(Cir3,pnttg3sol(NbrSol));
789 }
790 }
791 }
792 WellDone = Standard_True;
793}
794
795//=========================================================================
796
797Standard_Boolean GccAna_Circ2d3Tan::
798 IsDone () const {
799 return WellDone;
800 }
801
802Standard_Integer GccAna_Circ2d3Tan::
803 NbSolutions () const {
804 return NbrSol;
805 }
806
807gp_Circ2d GccAna_Circ2d3Tan::
808 ThisSolution (const Standard_Integer Index) const
809{
810 if (!WellDone)
9775fa61 811 throw StdFail_NotDone();
7fd59977 812
813 if (Index <= 0 ||Index > NbrSol)
9775fa61 814 throw Standard_OutOfRange();
7fd59977 815
816 return cirsol(Index);
817}
818
819void GccAna_Circ2d3Tan::
820 WhichQualifier(const Standard_Integer Index ,
821 GccEnt_Position& Qualif1 ,
822 GccEnt_Position& Qualif2 ,
823 GccEnt_Position& Qualif3 ) const
824{
9775fa61 825 if (!WellDone) { throw StdFail_NotDone(); }
826 else if (Index <= 0 ||Index > NbrSol) { throw Standard_OutOfRange(); }
7fd59977 827 else {
828 Qualif1 = qualifier1(Index);
829 Qualif2 = qualifier2(Index);
830 Qualif3 = qualifier3(Index);
831 }
832}
833
834void GccAna_Circ2d3Tan::
835 Tangency1 (const Standard_Integer Index,
836 Standard_Real& ParSol,
837 Standard_Real& ParArg,
838 gp_Pnt2d& PntSol) const {
839 if (!WellDone) {
9775fa61 840 throw StdFail_NotDone();
7fd59977 841 }
842 else if (Index <= 0 ||Index > NbrSol) {
9775fa61 843 throw Standard_OutOfRange();
7fd59977 844 }
845 else {
846 if (TheSame1(Index) == 0) {
847 ParSol = par1sol(Index);
848 ParArg = pararg1(Index);
849 PntSol = gp_Pnt2d(pnttg1sol(Index));
850 }
9775fa61 851 else { throw StdFail_NotDone(); }
7fd59977 852 }
853 }
854
855void GccAna_Circ2d3Tan::
856 Tangency2 (const Standard_Integer Index,
857 Standard_Real& ParSol,
858 Standard_Real& ParArg,
859 gp_Pnt2d& PntSol) const{
860 if (!WellDone) {
9775fa61 861 throw StdFail_NotDone();
7fd59977 862 }
863 else if (Index <= 0 ||Index > NbrSol) {
9775fa61 864 throw Standard_OutOfRange();
7fd59977 865 }
866 else {
867 if (TheSame2(Index) == 0) {
868 ParSol = par2sol(Index);
869 ParArg = pararg2(Index);
870 PntSol = gp_Pnt2d(pnttg2sol(Index));
871 }
9775fa61 872 else { throw StdFail_NotDone(); }
7fd59977 873 }
874 }
875
876void GccAna_Circ2d3Tan::
877 Tangency3 (const Standard_Integer Index,
878 Standard_Real& ParSol,
879 Standard_Real& ParArg,
880 gp_Pnt2d& PntSol) const{
881 if (!WellDone) {
9775fa61 882 throw StdFail_NotDone();
7fd59977 883 }
884 else if (Index <= 0 ||Index > NbrSol) {
9775fa61 885 throw Standard_OutOfRange();
7fd59977 886 }
887 else {
888 if (TheSame3(Index) == 0) {
889 ParSol = par3sol(Index);
890 ParArg = pararg3(Index);
891 PntSol = gp_Pnt2d(pnttg3sol(Index));
892 }
9775fa61 893 else { throw StdFail_NotDone(); }
7fd59977 894 }
895 }
896
897Standard_Boolean GccAna_Circ2d3Tan::
898 IsTheSame1 (const Standard_Integer Index) const
899{
900 if (!WellDone)
9775fa61 901 throw StdFail_NotDone();
7fd59977 902
903 if (Index <= 0 ||Index > NbrSol)
9775fa61 904 throw Standard_OutOfRange();
7fd59977 905
906 if (TheSame1(Index) == 0)
907 return Standard_False;
908
909 return Standard_True;
910}
911
912Standard_Boolean GccAna_Circ2d3Tan::
913 IsTheSame2 (const Standard_Integer Index) const
914{
915 if (!WellDone)
9775fa61 916 throw StdFail_NotDone();
7fd59977 917
918 if (Index <= 0 ||Index > NbrSol)
9775fa61 919 throw Standard_OutOfRange();
7fd59977 920
921 if (TheSame2(Index) == 0)
922 return Standard_False;
923
924 return Standard_True;
925}
926
927Standard_Boolean GccAna_Circ2d3Tan::
928 IsTheSame3 (const Standard_Integer Index) const
929{
930 if (!WellDone)
9775fa61 931 throw StdFail_NotDone();
7fd59977 932
933 if (Index <= 0 ||Index > NbrSol)
9775fa61 934 throw Standard_OutOfRange();
7fd59977 935
936 if (TheSame3(Index) == 0)
937 return Standard_False;
938
939 return Standard_True;
940}
941
942
Open CASCADE Technology repository