1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
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
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134
16
17 #include <ElCLib.hxx>
18 #include <ElSLib.hxx>
19 #include <Extrema_ExtElC.hxx>
20 #include <Extrema_ExtElCS.hxx>
21 #include <Extrema_ExtPElC.hxx>
22 #include <Extrema_ExtPElS.hxx>
23 #include <Extrema_POnCurv.hxx>
24 #include <Extrema_POnSurf.hxx>
25 #include <gp_Circ.hxx>
26 #include <gp_Cone.hxx>
27 #include <gp_Cylinder.hxx>
28 #include <gp_Hypr.hxx>
29 #include <gp_Lin.hxx>
30 #include <gp_Pln.hxx>
31 #include <gp_Sphere.hxx>
32 #include <gp_Torus.hxx>
33 #include <gp_Vec.hxx>
34 #include <IntAna_IntConicQuad.hxx>
35 #include <IntAna_Quadric.hxx>
36 #include <IntAna_QuadQuadGeo.hxx>
37 #include <Precision.hxx>
38 #include <Standard_NotImplemented.hxx>
39 #include <Standard_OutOfRange.hxx>
40 #include <StdFail_NotDone.hxx>
41 #include <TColStd_ListOfInteger.hxx>
42
Extrema_ExtElCS()43 Extrema_ExtElCS::Extrema_ExtElCS()
44 {
45 myDone = Standard_False;
46 myIsPar = Standard_False;
47 myNbExt = 0;
48 }
49
50
Extrema_ExtElCS(const gp_Lin & C,const gp_Pln & S)51 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
52 const gp_Pln& S)
53 {
54 Perform(C, S);
55 }
56
57
58
Perform(const gp_Lin & C,const gp_Pln & S)59 void Extrema_ExtElCS::Perform(const gp_Lin& C,
60 const gp_Pln& S)
61 {
62 myDone = Standard_True;
63 myIsPar = Standard_False;
64 myNbExt = 0;
65
66 if (C.Direction().IsNormal(S.Axis().Direction(),
67 Precision::Angular()))
68 {
69 mySqDist = new TColStd_HArray1OfReal(1, 1);
70 mySqDist->SetValue(1, S.SquareDistance(C));
71 myIsPar = Standard_True;
72 myNbExt = 1;
73 }
74 }
75
76
Extrema_ExtElCS(const gp_Lin & C,const gp_Cylinder & S)77 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
78 const gp_Cylinder& S)
79 {
80 Perform(C, S);
81 }
82
83
84
Perform(const gp_Lin & C,const gp_Cylinder & S)85 void Extrema_ExtElCS::Perform(const gp_Lin& C,
86 const gp_Cylinder& S)
87 {
88 myDone = Standard_False;
89 myNbExt = 0;
90 myIsPar = Standard_False;
91
92 gp_Ax3 Pos = S.Position();
93
94 Standard_Boolean isParallel = Standard_False;
95
96 Standard_Real radius = S.Radius();
97 Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
98 if (Extrem.IsParallel())
99 {
100 isParallel = Standard_True;
101 }
102 else
103 {
104 Standard_Integer i, aStartIdx = 0;
105
106 Extrema_POnCurv myPOnC1, myPOnC2;
107 Extrem.Points(1, myPOnC1, myPOnC2);
108 gp_Pnt PonAxis = myPOnC1.Value();
109 gp_Pnt PC = myPOnC2.Value();
110
111 // line intersects the cylinder
112 if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
113 {
114 IntAna_Quadric theQuadric(S);
115 IntAna_IntConicQuad Inters(C, theQuadric);
116 if (Inters.IsDone() && Inters.IsInQuadric())
117 {
118 isParallel = Standard_True;
119 }
120 else if (Inters.IsDone())
121 {
122 myNbExt = Inters.NbPoints();
123 aStartIdx = myNbExt;
124 if (myNbExt > 0)
125 {
126 // Not more than 2 additional points from perpendiculars.
127 mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
128 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
129 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
130 Standard_Real u, v, w;
131 for (i = 1; i <= myNbExt; i++)
132 {
133 mySqDist->SetValue(i, 0.);
134 gp_Pnt P_int = Inters.Point(i);
135 w = Inters.ParamOnConic(i);
136 Extrema_POnCurv PonC(w, P_int);
137 myPoint1->SetValue(i, PonC);
138 ElSLib::CylinderParameters(Pos, radius, P_int, u, v);
139 Extrema_POnSurf PonS(u, v, P_int);
140 myPoint2->SetValue(i, PonS);
141 }
142 }
143 }
144 }
145 else
146 {
147 // line is tangent or outside of the cylinder
148 Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
149 if (ExPS.IsDone())
150 {
151 if (aStartIdx == 0)
152 {
153 myNbExt = ExPS.NbExt();
154 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
155 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
156 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
157 }
158 else
159 myNbExt += ExPS.NbExt();
160
161 for (i = aStartIdx + 1; i <= myNbExt; i++)
162 {
163 myPoint1->SetValue(i, myPOnC2);
164 myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
165 mySqDist->SetValue(i, (myPOnC2.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
166 }
167 }
168 }
169
170 myDone = Standard_True;
171 }
172
173 if (isParallel)
174 {
175 // Line direction is similar to cylinder axis of rotation.
176 // The case is possible when either extrema returned parallel status
177 // or Intersection tool returned infinite number of solutions.
178 // This is possible due to Intersection algorithm uses more precise
179 // characteristics to consider given geometries parallel.
180 // In the latter case there may be several extremas, thus we look for
181 // the one with the lowest distance and use it as a final solution.
182 mySqDist = new TColStd_HArray1OfReal(1, 1);
183 Standard_Real aDist = Extrem.SquareDistance(1);
184 const Standard_Integer aNbExt = Extrem.NbExt();
185 for (Standard_Integer i = 2; i <= aNbExt; i++)
186 {
187 const Standard_Real aD = Extrem.SquareDistance(i);
188 if (aD < aDist)
189 {
190 aDist = aD;
191 }
192 }
193
194 aDist = sqrt(aDist) - radius;
195 mySqDist->SetValue(1, aDist * aDist);
196 myDone = Standard_True;
197 myIsPar = Standard_True;
198 myNbExt = 1;
199 }
200 }
201
202
203
Extrema_ExtElCS(const gp_Lin & C,const gp_Cone & S)204 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
205 const gp_Cone& S)
206 { Perform(C, S);}
207
208
209
210 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
211 // const gp_Cone& S)
Perform(const gp_Lin &,const gp_Cone &)212 void Extrema_ExtElCS::Perform(const gp_Lin& ,
213 const gp_Cone& )
214 {
215 throw Standard_NotImplemented();
216
217 }
218
219
220
Extrema_ExtElCS(const gp_Lin & C,const gp_Sphere & S)221 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
222 const gp_Sphere& S)
223 {
224 Perform(C, S);
225 }
226
227
228
Perform(const gp_Lin & C,const gp_Sphere & S)229 void Extrema_ExtElCS::Perform(const gp_Lin& C,
230 const gp_Sphere& S)
231 {
232 // In case of intersection - return four points:
233 // 2 intersection points and 2 perpendicular.
234 // No intersection - only min and max.
235
236 myDone = Standard_False;
237 myNbExt = 0;
238 myIsPar = Standard_False;
239 Standard_Integer aStartIdx = 0;
240
241 gp_Pnt aCenter = S.Location();
242
243 Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());
244
245 Standard_Integer i;
246 if (Extrem.IsDone() &&
247 Extrem.NbExt() > 0)
248 {
249 Extrema_POnCurv myPOnC1 = Extrem.Point(1);
250 if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
251 {
252 IntAna_IntConicQuad aLinSphere(C, S);
253 if (aLinSphere.IsDone())
254 {
255 myNbExt = aLinSphere.NbPoints();
256 aStartIdx = myNbExt;
257 // Not more than 2 additional points from perpendiculars.
258 mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
259 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
260 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
261
262 for (i = 1; i <= myNbExt; i++)
263 {
264 Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));
265
266 Standard_Real u,v;
267 ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
268 Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));
269
270 myPoint1->SetValue(i, aCPnt);
271 myPoint2->SetValue(i, aSPnt);
272 mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
273 }
274 }
275 }
276
277 Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
278 if (ExPS.IsDone())
279 {
280 if (aStartIdx == 0)
281 {
282 myNbExt = ExPS.NbExt();
283 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
284 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
285 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
286 }
287 else
288 myNbExt += ExPS.NbExt();
289
290 for (i = aStartIdx + 1; i <= myNbExt; i++)
291 {
292 myPoint1->SetValue(i, myPOnC1);
293 myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
294 mySqDist->SetValue(i,(myPOnC1.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
295 }
296 }
297 }
298 myDone = Standard_True;
299 }
300
301
Extrema_ExtElCS(const gp_Lin & C,const gp_Torus & S)302 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
303 const gp_Torus& S)
304 { Perform(C, S);}
305
306
307
308 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
309 // const gp_Torus& S)
Perform(const gp_Lin &,const gp_Torus &)310 void Extrema_ExtElCS::Perform(const gp_Lin& ,
311 const gp_Torus& )
312 {
313 throw Standard_NotImplemented();
314
315 }
316
317
318 // Circle-?
319
Extrema_ExtElCS(const gp_Circ & C,const gp_Pln & S)320 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
321 const gp_Pln& S)
322 {
323 Perform(C, S);
324 }
325
326
327
Perform(const gp_Circ & C,const gp_Pln & S)328 void Extrema_ExtElCS::Perform(const gp_Circ& C,
329 const gp_Pln& S)
330 {
331 myDone = Standard_True;
332 myIsPar = Standard_False;
333 myNbExt = 0;
334
335 gp_Ax2 Pos = C.Position();
336 gp_Dir NCirc = Pos.Direction();
337 gp_Dir NPln = S.Axis().Direction();
338
339 Standard_Boolean isParallel = Standard_False;
340
341 if (NCirc.IsParallel(NPln, Precision::Angular())) {
342 isParallel = Standard_True;
343 }
344 else {
345
346 gp_Dir ExtLine = NCirc ^ NPln;
347 ExtLine = ExtLine ^ NCirc;
348 //
349 gp_Dir XDir = Pos.XDirection();
350 Standard_Real T[2];
351 T[0] = XDir.AngleWithRef(ExtLine, NCirc);
352 if(T[0] < 0.)
353 {
354 //Put in period
355 T[0] += M_PI;
356 }
357 T[1] = T[0] + M_PI;
358 //
359 myNbExt = 2;
360 //Check intersection
361 IntAna_IntConicQuad anInter(C, S,
362 Precision::Angular(),
363 Precision::Confusion());
364
365 if (anInter.IsDone() && anInter.IsInQuadric())
366 {
367 isParallel = Standard_True;
368 }
369 else if (anInter.IsDone())
370 {
371 if(anInter.NbPoints() > 1)
372 {
373 myNbExt += anInter.NbPoints();
374 }
375 }
376
377 if (!isParallel)
378 {
379 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
380 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
381 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
382
383 Standard_Integer i;
384 gp_Pnt PC, PP;
385 Standard_Real U, V;
386 Extrema_POnCurv POnC;
387 Extrema_POnSurf POnS;
388 for (i = 0; i < 2; ++i)
389 {
390 PC = ElCLib::CircleValue(T[i], C.Position(), C.Radius());
391 POnC.SetValues(T[i], PC);
392 myPoint1->SetValue(i + 1, POnC);
393 ElSLib::PlaneParameters(S.Position(), PC, U, V);
394 PP = ElSLib::PlaneValue(U, V, S.Position());
395 POnS.SetParameters(U, V, PP);
396 myPoint2->SetValue(i + 1, POnS);
397 mySqDist->SetValue(i + 1, PC.SquareDistance(PP));
398 }
399 //
400 if (myNbExt > 2)
401 {
402 //Add intersection points
403 for (i = 1; i <= anInter.NbPoints(); ++i)
404 {
405 Standard_Real t = anInter.ParamOnConic(i);
406 PC = ElCLib::CircleValue(t, C.Position(), C.Radius());
407 POnC.SetValues(t, PC);
408 myPoint1->SetValue(i + 2, POnC);
409 ElSLib::PlaneParameters(S.Position(), PC, U, V);
410 PP = ElSLib::PlaneValue(U, V, S.Position());
411 POnS.SetParameters(U, V, PP);
412 myPoint2->SetValue(i + 2, POnS);
413 mySqDist->SetValue(i + 2, PC.SquareDistance(PP));
414 }
415 }
416 }
417 }
418
419 if (isParallel)
420 {
421 mySqDist = new TColStd_HArray1OfReal(1, 1);
422 mySqDist->SetValue(1, S.SquareDistance(C.Location()));
423 myIsPar = Standard_True;
424 myNbExt = 1;
425 }
426 }
427
Extrema_ExtElCS(const gp_Circ & C,const gp_Cylinder & S)428 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
429 const gp_Cylinder& S)
430 {
431 Perform(C, S);
432 }
433
434
435
436 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 Begin
437 // Implementation of the method.
Perform(const gp_Circ & C,const gp_Cylinder & S)438 void Extrema_ExtElCS::Perform(const gp_Circ& C,
439 const gp_Cylinder& S)
440 {
441 myDone = Standard_False;
442 myIsPar = Standard_False;
443 myNbExt = 0;
444
445 // Get an axis line of the cylinder.
446 gp_Lin anAxis(S.Axis());
447
448 // Compute extrema between the circle and the line.
449 Extrema_ExtElC anExtC(anAxis, C, 0.);
450
451 if (!anExtC.IsDone())
452 return;
453
454 Standard_Boolean isParallel = Standard_False;
455
456 if (anExtC.IsParallel()) {
457 isParallel = Standard_True;
458 } else {
459 Standard_Integer aNbExt = anExtC.NbExt();
460 Standard_Integer i;
461 Standard_Integer aCurI = 1;
462 Standard_Real aTolConf = Precision::Confusion();
463 Standard_Real aCylRad = S.Radius();
464
465 // Check whether two objects have intersection points
466 IntAna_Quadric aCylQuad(S);
467 IntAna_IntConicQuad aCircCylInter(C, aCylQuad);
468 Standard_Integer aNbInter = 0;
469 if (aCircCylInter.IsDone() && aCircCylInter.IsInQuadric())
470 {
471 isParallel = Standard_True;
472 }
473 else if (aCircCylInter.IsDone())
474 {
475 aNbInter = aCircCylInter.NbPoints();
476 }
477
478 if (!isParallel)
479 {
480 // Compute the extremas.
481 myNbExt = 2 * aNbExt + aNbInter;
482 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
483 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
484 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
485
486 for (i = 1; i <= aNbExt; i++) {
487 Extrema_POnCurv aPOnAxis;
488 Extrema_POnCurv aPOnCirc;
489 Standard_Real aSqDist = anExtC.SquareDistance(i);
490 Standard_Real aDist = sqrt(aSqDist);
491
492 anExtC.Points(i, aPOnAxis, aPOnCirc);
493
494 if (aSqDist <= (aTolConf * aTolConf)) {
495 myNbExt -= 2;
496 continue;
497 }
498
499 gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
500 Standard_Real aShift[2] = {aDist + aCylRad, aDist - aCylRad};
501 Standard_Integer j;
502
503 for (j = 0; j < 2; j++) {
504 gp_Vec aVec(aDir);
505 gp_Pnt aPntOnCyl;
506
507 aVec.Multiply(aShift[j]);
508 aPntOnCyl = aPOnCirc.Value().Translated(aVec);
509
510 Standard_Real aU;
511 Standard_Real aV;
512
513 ElSLib::Parameters(S, aPntOnCyl, aU, aV);
514
515 Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);
516
517 myPoint1->SetValue(aCurI, aPOnCirc);
518 myPoint2->SetValue(aCurI, aPOnSurf);
519 mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
520 }
521 }
522
523 // Adding intersection points to the list of extremas
524 for (i=1; i<=aNbInter; i++)
525 {
526 Standard_Real aU;
527 Standard_Real aV;
528
529 gp_Pnt aInterPnt = aCircCylInter.Point(i);
530
531 aU = ElCLib::Parameter(C, aInterPnt);
532 Extrema_POnCurv aPOnCirc(aU, aInterPnt);
533
534 ElSLib::Parameters(S, aInterPnt, aU, aV);
535 Extrema_POnSurf aPOnCyl(aU, aV, aInterPnt);
536 myPoint1->SetValue(aCurI, aPOnCirc);
537 myPoint2->SetValue(aCurI, aPOnCyl);
538 mySqDist->SetValue(aCurI++, 0.0);
539 }
540 }
541 }
542
543 myDone = Standard_True;
544
545 if (isParallel)
546 {
547 // The case is possible when either extrema returned parallel status
548 // or Intersection tool returned infinite number of solutions.
549 // This is possible due to Intersection algorithm uses more precise
550 // characteristics to consider given geometries parallel.
551 // In the latter case there may be several extremas, thus we look for
552 // the one with the lowest distance and use it as a final solution.
553
554 myIsPar = Standard_True;
555 myNbExt = 1;
556 mySqDist = new TColStd_HArray1OfReal(1, 1);
557 Standard_Real aDist = anExtC.SquareDistance(1);
558
559 const Standard_Integer aNbExt = anExtC.NbExt();
560 for (Standard_Integer i = 2; i <= aNbExt; i++)
561 {
562 const Standard_Real aD = anExtC.SquareDistance(i);
563 if (aD < aDist)
564 {
565 aDist = aD;
566 }
567 }
568
569 aDist = sqrt(aDist) - S.Radius();
570 mySqDist->SetValue(1, aDist * aDist);
571 }
572 }
573 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 End
574
575
576
Extrema_ExtElCS(const gp_Circ & C,const gp_Cone & S)577 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
578 const gp_Cone& S)
579 { Perform(C, S);}
580
581
582
583 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
584 // const gp_Cone& S)
Perform(const gp_Circ &,const gp_Cone &)585 void Extrema_ExtElCS::Perform(const gp_Circ& ,
586 const gp_Cone& )
587 {
588 throw Standard_NotImplemented();
589
590 }
591
592
593
594 //=======================================================================
595 //function : Extrema_ExtElCS
596 //purpose : Circle/Sphere
597 //=======================================================================
Extrema_ExtElCS(const gp_Circ & C,const gp_Sphere & S)598 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
599 const gp_Sphere& S)
600 {
601 Perform(C, S);
602 }
603
604 //=======================================================================
605 //function : Perform
606 //purpose : Circle/Sphere
607 //=======================================================================
Perform(const gp_Circ & C,const gp_Sphere & S)608 void Extrema_ExtElCS::Perform(const gp_Circ& C,
609 const gp_Sphere& S)
610 {
611 myDone = Standard_False;
612 myIsPar = Standard_False;
613 myNbExt = 0;
614
615 if (gp_Lin(C.Axis()).SquareDistance(S.Location()) < Precision::SquareConfusion())
616 {
617 // Circle and sphere are parallel
618 myIsPar = Standard_True;
619 myDone = Standard_True;
620 myNbExt = 1;
621
622 // Compute distance from circle to the sphere
623 Standard_Real aSqDistLoc = C.Location().SquareDistance(S.Location());
624 Standard_Real aSqDist = aSqDistLoc + C.Radius() * C.Radius();
625 Standard_Real aDist = sqrt(aSqDist) - S.Radius();
626 mySqDist = new TColStd_HArray1OfReal(1, 1);
627 mySqDist->SetValue(1, aDist * aDist);
628 return;
629 }
630
631 // Intersect sphere with circle's plane
632 gp_Pln CPln(C.Location(), C.Axis().Direction());
633 IntAna_QuadQuadGeo anInter(CPln, S);
634 if (!anInter.IsDone())
635 // not done
636 return;
637
638 if (anInter.TypeInter() != IntAna_Circle)
639 {
640 // Intersection is empty or just a point.
641 // The parallel case has already been considered,
642 // thus, here we have to find only one minimal solution
643 myNbExt = 1;
644 myDone = Standard_True;
645
646 mySqDist = new TColStd_HArray1OfReal(1, 1);
647 myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
648 myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
649
650 // Compute parameter on circle
651 const Standard_Real aT = ElCLib::Parameter(C, S.Location());
652 // Compute point on circle
653 gp_Pnt aPOnC = ElCLib::Value(aT, C);
654
655 // Compute parameters on sphere
656 Standard_Real aU, aV;
657 ElSLib::Parameters(S, aPOnC, aU, aV);
658 // Compute point on sphere
659 gp_Pnt aPOnS = ElSLib::Value(aU, aV, S);
660
661 // Save solution
662 myPoint1->SetValue(1, Extrema_POnCurv(aT, aPOnC));
663 myPoint2->SetValue(1, Extrema_POnSurf(aU, aV, aPOnS));
664 mySqDist->SetValue(1, aPOnC.SquareDistance(aPOnS));
665 return;
666 }
667
668 // Here, the intersection is a circle
669
670 // Intersection circle
671 gp_Circ aCInt = anInter.Circle(1);
672
673 // Perform intersection of the input circle with the intersection circle
674 Extrema_ExtElC anExtC(C, aCInt);
675 Standard_Boolean isExtremaCircCircValid = anExtC.IsDone() // Check if intersection is done
676 && !anExtC.IsParallel() // Parallel case has already been considered
677 && anExtC.NbExt() > 0; // Check that some solutions have been found
678 if (!isExtremaCircCircValid)
679 // not done
680 return;
681
682 myDone = Standard_True;
683
684 // Few solutions
685 Standard_Real aNbExt = anExtC.NbExt();
686 // Find the minimal distance
687 Standard_Real aMinSqDist = ::RealLast();
688 for (Standard_Integer i = 1; i <= aNbExt; ++i)
689 {
690 Standard_Real aSqDist = anExtC.SquareDistance(i);
691 if (aSqDist < aMinSqDist)
692 aMinSqDist = aSqDist;
693 }
694
695 // Collect all solutions close to the minimal one
696 TColStd_ListOfInteger aSols;
697 for (Standard_Integer i = 1; i <= aNbExt; ++i)
698 {
699 Standard_Real aDiff = anExtC.SquareDistance(i) - aMinSqDist;
700 if (aDiff < Precision::SquareConfusion())
701 aSols.Append(i);
702 }
703
704 // Save all minimal solutions
705 myNbExt = aSols.Extent();
706
707 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
708 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
709 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
710
711 TColStd_ListIteratorOfListOfInteger it(aSols);
712 for (Standard_Integer iSol = 1; it.More(); it.Next(), ++iSol)
713 {
714 Extrema_POnCurv P1, P2;
715 anExtC.Points(it.Value(), P1, P2);
716
717 // Compute parameters on sphere
718 Standard_Real aU, aV;
719 ElSLib::Parameters(S, P1.Value(), aU, aV);
720 // Compute point on sphere
721 gp_Pnt aPOnS = ElSLib::Value(aU, aV, S);
722
723 // Save solution
724 myPoint1->SetValue(iSol, P1);
725 myPoint2->SetValue(iSol, Extrema_POnSurf(aU, aV, aPOnS));
726 mySqDist->SetValue(iSol, P1.Value().SquareDistance(aPOnS));
727 }
728 }
729
Extrema_ExtElCS(const gp_Circ & C,const gp_Torus & S)730 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
731 const gp_Torus& S)
732 { Perform(C, S);}
733
734
735
736 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
737 // const gp_Torus& S)
Perform(const gp_Circ &,const gp_Torus &)738 void Extrema_ExtElCS::Perform(const gp_Circ& ,
739 const gp_Torus& )
740 {
741 throw Standard_NotImplemented();
742
743 }
744
Extrema_ExtElCS(const gp_Hypr & C,const gp_Pln & S)745 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
746 const gp_Pln& S)
747 {
748 Perform(C, S);
749 }
750
751
752
Perform(const gp_Hypr & C,const gp_Pln & S)753 void Extrema_ExtElCS::Perform(const gp_Hypr& C,
754 const gp_Pln& S)
755 {
756 myDone = Standard_True;
757 myIsPar = Standard_False;
758 myNbExt = 0;
759
760 gp_Ax2 Pos = C.Position();
761 gp_Dir NHypr = Pos.Direction();
762 gp_Dir NPln = S.Axis().Direction();
763
764 if (NHypr.IsParallel(NPln, Precision::Angular())) {
765
766 mySqDist = new TColStd_HArray1OfReal(1, 1);
767 mySqDist->SetValue(1, S.SquareDistance(C.Location()));
768 myIsPar = Standard_True;
769 myNbExt = 1;
770 }
771 else {
772
773 gp_Dir XDir = Pos.XDirection();
774 gp_Dir YDir = Pos.YDirection();
775
776 Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir));
777 Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir));
778
779 if(Abs(B) > Abs(A)) {
780 Standard_Real T = -0.5 * Log((A+B)/(B-A));
781 gp_Pnt Ph = ElCLib::HyperbolaValue(T, Pos, C.MajorRadius(), C.MinorRadius());
782 Extrema_POnCurv PC(T, Ph);
783 myPoint1 = new Extrema_HArray1OfPOnCurv(1,1);
784 myPoint1->SetValue(1, PC);
785
786 mySqDist = new TColStd_HArray1OfReal(1, 1);
787 mySqDist->SetValue(1, S.SquareDistance(Ph));
788
789 Standard_Real U, V;
790 ElSLib::PlaneParameters(S.Position(), Ph, U, V);
791 gp_Pnt Pp = ElSLib::PlaneValue(U, V, S.Position());
792 Extrema_POnSurf PS(U, V, Pp);
793 myPoint2 = new Extrema_HArray1OfPOnSurf(1,1);
794 myPoint2->SetValue(1, PS);
795
796 myNbExt = 1;
797 }
798 }
799 }
800
801
IsDone() const802 Standard_Boolean Extrema_ExtElCS::IsDone() const
803 {
804 return myDone;
805 }
806
807
NbExt() const808 Standard_Integer Extrema_ExtElCS::NbExt() const
809 {
810 if (!IsDone()) throw StdFail_NotDone();
811 return myNbExt;
812 }
813
SquareDistance(const Standard_Integer N) const814 Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
815 {
816 if (N < 1 || N > NbExt())
817 {
818 throw Standard_OutOfRange();
819 }
820
821 return mySqDist->Value(N);
822 }
823
824
Points(const Standard_Integer N,Extrema_POnCurv & P1,Extrema_POnSurf & P2) const825 void Extrema_ExtElCS::Points(const Standard_Integer N,
826 Extrema_POnCurv& P1,
827 Extrema_POnSurf& P2) const
828 {
829 if (N < 1 || N > NbExt())
830 {
831 throw Standard_OutOfRange();
832 }
833
834 P1 = myPoint1->Value(N);
835 P2 = myPoint2->Value(N);
836 }
837
838
IsParallel() const839 Standard_Boolean Extrema_ExtElCS::IsParallel() const
840 {
841 if (!IsDone())
842 {
843 throw StdFail_NotDone();
844 }
845 return myIsPar;
846 }
847