1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134
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>
31 #include <gp_Sphere.hxx>
32 #include <gp_Torus.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>
43 Extrema_ExtElCS::Extrema_ExtElCS()
45 myDone = Standard_False;
46 myIsPar = Standard_False;
51 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
59 void Extrema_ExtElCS::Perform(const gp_Lin& C,
62 myDone = Standard_True;
63 myIsPar = Standard_False;
66 if (C.Direction().IsNormal(S.Axis().Direction(),
67 Precision::Angular()))
69 mySqDist = new TColStd_HArray1OfReal(1, 1);
70 mySqDist->SetValue(1, S.SquareDistance(C));
71 myIsPar = Standard_True;
77 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
85 void Extrema_ExtElCS::Perform(const gp_Lin& C,
88 myDone = Standard_False;
90 myIsPar = Standard_False;
92 gp_Ax3 Pos = S.Position();
94 Standard_Boolean isParallel = Standard_False;
96 Standard_Real radius = S.Radius();
97 Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
98 if (Extrem.IsParallel())
100 isParallel = Standard_True;
104 Standard_Integer i, aStartIdx = 0;
106 Extrema_POnCurv myPOnC1, myPOnC2;
107 Extrem.Points(1, myPOnC1, myPOnC2);
108 gp_Pnt PonAxis = myPOnC1.Value();
109 gp_Pnt PC = myPOnC2.Value();
111 // line intersects the cylinder
112 if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
114 IntAna_Quadric theQuadric(S);
115 IntAna_IntConicQuad Inters(C, theQuadric);
116 if (Inters.IsDone() && Inters.IsInQuadric())
118 isParallel = Standard_True;
120 else if (Inters.IsDone())
122 myNbExt = Inters.NbPoints();
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++)
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);
147 // line is tangent or outside of the cylinder
148 Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
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);
159 myNbExt += ExPS.NbExt();
161 for (i = aStartIdx + 1; i <= myNbExt; i++)
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()));
170 myDone = Standard_True;
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++)
187 const Standard_Real aD = Extrem.SquareDistance(i);
194 aDist = sqrt(aDist) - radius;
195 mySqDist->SetValue(1, aDist * aDist);
196 myDone = Standard_True;
197 myIsPar = Standard_True;
204 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
210 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
212 void Extrema_ExtElCS::Perform(const gp_Lin& ,
215 throw Standard_NotImplemented();
221 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
229 void Extrema_ExtElCS::Perform(const gp_Lin& C,
232 // In case of intersection - return four points:
233 // 2 intersection points and 2 perpendicular.
234 // No intersection - only min and max.
236 myDone = Standard_False;
238 myIsPar = Standard_False;
239 Standard_Integer aStartIdx = 0;
241 gp_Pnt aCenter = S.Location();
243 Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());
246 if (Extrem.IsDone() &&
249 Extrema_POnCurv myPOnC1 = Extrem.Point(1);
250 if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
252 IntAna_IntConicQuad aLinSphere(C, S);
253 if (aLinSphere.IsDone())
255 myNbExt = aLinSphere.NbPoints();
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);
262 for (i = 1; i <= myNbExt; i++)
264 Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));
267 ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
268 Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));
270 myPoint1->SetValue(i, aCPnt);
271 myPoint2->SetValue(i, aSPnt);
272 mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
277 Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
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);
288 myNbExt += ExPS.NbExt();
290 for (i = aStartIdx + 1; i <= myNbExt; i++)
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()));
298 myDone = Standard_True;
302 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
308 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
309 // const gp_Torus& S)
310 void Extrema_ExtElCS::Perform(const gp_Lin& ,
313 throw Standard_NotImplemented();
320 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
328 void Extrema_ExtElCS::Perform(const gp_Circ& C,
331 myDone = Standard_True;
332 myIsPar = Standard_False;
335 gp_Ax2 Pos = C.Position();
336 gp_Dir NCirc = Pos.Direction();
337 gp_Dir NPln = S.Axis().Direction();
339 Standard_Boolean isParallel = Standard_False;
341 if (NCirc.IsParallel(NPln, Precision::Angular())) {
342 isParallel = Standard_True;
346 gp_Dir ExtLine = NCirc ^ NPln;
347 ExtLine = ExtLine ^ NCirc;
349 gp_Dir XDir = Pos.XDirection();
351 T[0] = XDir.AngleWithRef(ExtLine, NCirc);
361 IntAna_IntConicQuad anInter(C, S,
362 Precision::Angular(),
363 Precision::Confusion());
365 if (anInter.IsDone() && anInter.IsInQuadric())
367 isParallel = Standard_True;
369 else if (anInter.IsDone())
371 if(anInter.NbPoints() > 1)
373 myNbExt += anInter.NbPoints();
379 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
380 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
381 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
386 Extrema_POnCurv POnC;
387 Extrema_POnSurf POnS;
388 for (i = 0; i < 2; ++i)
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));
402 //Add intersection points
403 for (i = 1; i <= anInter.NbPoints(); ++i)
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));
421 mySqDist = new TColStd_HArray1OfReal(1, 1);
422 mySqDist->SetValue(1, S.SquareDistance(C.Location()));
423 myIsPar = Standard_True;
428 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
429 const gp_Cylinder& S)
436 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 Begin
437 // Implementation of the method.
438 void Extrema_ExtElCS::Perform(const gp_Circ& C,
439 const gp_Cylinder& S)
441 myDone = Standard_False;
442 myIsPar = Standard_False;
445 // Get an axis line of the cylinder.
446 gp_Lin anAxis(S.Axis());
448 // Compute extrema between the circle and the line.
449 Extrema_ExtElC anExtC(anAxis, C, 0.);
451 if (!anExtC.IsDone())
454 Standard_Boolean isParallel = Standard_False;
456 if (anExtC.IsParallel()) {
457 isParallel = Standard_True;
459 Standard_Integer aNbExt = anExtC.NbExt();
461 Standard_Integer aCurI = 1;
462 Standard_Real aTolConf = Precision::Confusion();
463 Standard_Real aCylRad = S.Radius();
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())
471 isParallel = Standard_True;
473 else if (aCircCylInter.IsDone())
475 aNbInter = aCircCylInter.NbPoints();
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);
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);
492 anExtC.Points(i, aPOnAxis, aPOnCirc);
494 if (aSqDist <= (aTolConf * aTolConf)) {
499 gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
500 Standard_Real aShift[2] = {aDist + aCylRad, aDist - aCylRad};
503 for (j = 0; j < 2; j++) {
507 aVec.Multiply(aShift[j]);
508 aPntOnCyl = aPOnCirc.Value().Translated(aVec);
513 ElSLib::Parameters(S, aPntOnCyl, aU, aV);
515 Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);
517 myPoint1->SetValue(aCurI, aPOnCirc);
518 myPoint2->SetValue(aCurI, aPOnSurf);
519 mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
523 // Adding intersection points to the list of extremas
524 for (i=1; i<=aNbInter; i++)
529 gp_Pnt aInterPnt = aCircCylInter.Point(i);
531 aU = ElCLib::Parameter(C, aInterPnt);
532 Extrema_POnCurv aPOnCirc(aU, aInterPnt);
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);
543 myDone = Standard_True;
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.
554 myIsPar = Standard_True;
556 mySqDist = new TColStd_HArray1OfReal(1, 1);
557 Standard_Real aDist = anExtC.SquareDistance(1);
559 const Standard_Integer aNbExt = anExtC.NbExt();
560 for (Standard_Integer i = 2; i <= aNbExt; i++)
562 const Standard_Real aD = anExtC.SquareDistance(i);
569 aDist = sqrt(aDist) - S.Radius();
570 mySqDist->SetValue(1, aDist * aDist);
573 // Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 End
577 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
583 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
585 void Extrema_ExtElCS::Perform(const gp_Circ& ,
588 throw Standard_NotImplemented();
594 //=======================================================================
595 //function : Extrema_ExtElCS
596 //purpose : Circle/Sphere
597 //=======================================================================
598 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
604 //=======================================================================
606 //purpose : Circle/Sphere
607 //=======================================================================
608 void Extrema_ExtElCS::Perform(const gp_Circ& C,
611 myDone = Standard_False;
612 myIsPar = Standard_False;
615 if (gp_Lin(C.Axis()).SquareDistance(S.Location()) < Precision::SquareConfusion())
617 // Circle and sphere are parallel
618 myIsPar = Standard_True;
619 myDone = Standard_True;
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);
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())
638 if (anInter.TypeInter() != IntAna_Circle)
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
644 myDone = Standard_True;
646 mySqDist = new TColStd_HArray1OfReal(1, 1);
647 myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
648 myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
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);
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);
662 myPoint1->SetValue(1, Extrema_POnCurv(aT, aPOnC));
663 myPoint2->SetValue(1, Extrema_POnSurf(aU, aV, aPOnS));
664 mySqDist->SetValue(1, aPOnC.SquareDistance(aPOnS));
668 // Here, the intersection is a circle
670 // Intersection circle
671 gp_Circ aCInt = anInter.Circle(1);
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)
682 myDone = Standard_True;
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)
690 Standard_Real aSqDist = anExtC.SquareDistance(i);
691 if (aSqDist < aMinSqDist)
692 aMinSqDist = aSqDist;
695 // Collect all solutions close to the minimal one
696 TColStd_ListOfInteger aSols;
697 for (Standard_Integer i = 1; i <= aNbExt; ++i)
699 Standard_Real aDiff = anExtC.SquareDistance(i) - aMinSqDist;
700 if (aDiff < Precision::SquareConfusion())
704 // Save all minimal solutions
705 myNbExt = aSols.Extent();
707 mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
708 myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
709 myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
711 TColStd_ListIteratorOfListOfInteger it(aSols);
712 for (Standard_Integer iSol = 1; it.More(); it.Next(), ++iSol)
714 Extrema_POnCurv P1, P2;
715 anExtC.Points(it.Value(), P1, P2);
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);
724 myPoint1->SetValue(iSol, P1);
725 myPoint2->SetValue(iSol, Extrema_POnSurf(aU, aV, aPOnS));
726 mySqDist->SetValue(iSol, P1.Value().SquareDistance(aPOnS));
730 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
736 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
737 // const gp_Torus& S)
738 void Extrema_ExtElCS::Perform(const gp_Circ& ,
741 throw Standard_NotImplemented();
745 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
753 void Extrema_ExtElCS::Perform(const gp_Hypr& C,
756 myDone = Standard_True;
757 myIsPar = Standard_False;
760 gp_Ax2 Pos = C.Position();
761 gp_Dir NHypr = Pos.Direction();
762 gp_Dir NPln = S.Axis().Direction();
764 if (NHypr.IsParallel(NPln, Precision::Angular())) {
766 mySqDist = new TColStd_HArray1OfReal(1, 1);
767 mySqDist->SetValue(1, S.SquareDistance(C.Location()));
768 myIsPar = Standard_True;
773 gp_Dir XDir = Pos.XDirection();
774 gp_Dir YDir = Pos.YDirection();
776 Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir));
777 Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir));
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);
786 mySqDist = new TColStd_HArray1OfReal(1, 1);
787 mySqDist->SetValue(1, S.SquareDistance(Ph));
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);
802 Standard_Boolean Extrema_ExtElCS::IsDone() const
808 Standard_Integer Extrema_ExtElCS::NbExt() const
810 if (!IsDone()) throw StdFail_NotDone();
814 Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
816 if (N < 1 || N > NbExt())
818 throw Standard_OutOfRange();
821 return mySqDist->Value(N);
825 void Extrema_ExtElCS::Points(const Standard_Integer N,
827 Extrema_POnSurf& P2) const
829 if (N < 1 || N > NbExt())
831 throw Standard_OutOfRange();
834 P1 = myPoint1->Value(N);
835 P2 = myPoint2->Value(N);
839 Standard_Boolean Extrema_ExtElCS::IsParallel() const
843 throw StdFail_NotDone();