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.
17 #include <Extrema_ExtPElC.hxx>
18 #include <Extrema_POnCurv.hxx>
19 #include <gp_Circ.hxx>
20 #include <gp_Elips.hxx>
21 #include <gp_Hypr.hxx>
23 #include <gp_Parab.hxx>
25 #include <math_DirectPolynomialRoots.hxx>
26 #include <math_TrigonometricFunctionRoots.hxx>
27 #include <Precision.hxx>
28 #include <Standard_NotImplemented.hxx>
29 #include <Standard_OutOfRange.hxx>
30 #include <StdFail_NotDone.hxx>
32 //=======================================================================
33 //function : Extrema_ExtPElC
35 //=======================================================================
36 Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; }
38 //=======================================================================
39 //function : Extrema_ExtPElC
41 //=======================================================================
42 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
44 const Standard_Real Tol,
45 const Standard_Real Uinf,
46 const Standard_Real Usup)
48 Perform(P, L, Tol, Uinf, Usup);
50 //=======================================================================
53 //=======================================================================
54 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
56 const Standard_Real Tol,
57 const Standard_Real Uinf,
58 const Standard_Real Usup)
60 myDone = Standard_False;
62 gp_Vec V1 (L.Direction());
63 gp_Pnt OR = L.Location();
65 Standard_Real Mydist = V1.Dot(V);
66 if ((Mydist >= Uinf-Tol) &&
67 (Mydist <= Usup+Tol)){
69 gp_Pnt MyP = OR.Translated(Mydist*V1);
70 Extrema_POnCurv MyPOnCurve(Mydist, MyP);
71 mySqDist[0] = P.SquareDistance(MyP);
72 myPoint[0] = MyPOnCurve;
73 myIsMin[0] = Standard_True;
75 myDone = Standard_True;
82 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
84 const Standard_Real Tol,
85 const Standard_Real Uinf,
86 const Standard_Real Usup)
88 Perform(P, C, Tol, Uinf, Usup);
91 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
93 const Standard_Real Tol,
94 const Standard_Real Uinf,
95 const Standard_Real Usup)
96 /*-----------------------------------------------------------------------------
98 Find values of parameter u such as:
99 - dist(P,C(u)) pass by an extrema,
104 1- Projection of point P in the plane of the circle,
105 2- Calculation of u solutions in [0.,2.*M_PI]:
106 Let Pp, the projected point and
107 O, the center of the circle;
109 - if Pp is mixed with 0, there is an infinite number of solutions;
110 IsDone() renvoie Standard_False.
111 - otherwise, 2 points are solutions for the complete circle:
112 . Us1 = angle(OPp,OX) corresponds to the minimum,
113 . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
114 ( Us1 - M_PI otherwise;
115 Us2 corresponds to the maximum.
116 3- Calculate the extrema in [Uinf,Usup].
117 -----------------------------------------------------------------------------*/
119 myDone = Standard_False;
122 // 1- Projection of the point P in the plane of circle -> Pp ...
124 gp_Pnt O = C.Location();
125 gp_Vec Axe (C.Axis().Direction());
126 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
127 gp_Pnt Pp = P.Translated(Trsl);
129 // 2- Calculate u solutions in [0.,2.*PI] ...
132 if (OPp.Magnitude() < Tol) { return; }
133 Standard_Real Usol[2];
134 Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
136 const Standard_Real aAngTol = Precision::Angular();
137 if ( Usol[0] + M_PI < aAngTol )
139 else if ( Usol[0] - M_PI > -aAngTol )
142 Usol[1] = Usol[0] + M_PI;
144 Standard_Real myuinf = Uinf;
145 //Standard_Real TolU = Tol*C.Radius();
146 Standard_Real TolU, aR;
148 TolU=Precision::Infinite();
149 if (aR > gp::Resolution()) {
153 ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
154 ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
155 if (((Usol[0]-2*M_PI-Uinf) < TolU) && ((Usol[0]-2*M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
156 if (((Usol[1]-2*M_PI-Uinf) < TolU) && ((Usol[1]-2*M_PI-Uinf) > -TolU)) Usol[1] = Uinf;
159 // 3- Calculate extrema in [Umin,Umax] ...
163 for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
165 if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
166 Cu = ElCLib::Value(Us,C);
167 mySqDist[myNbExt] = Cu.SquareDistance(P);
168 myIsMin[myNbExt] = (NoSol == 0);
169 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
173 myDone = Standard_True;
175 //=============================================================================
177 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
179 const Standard_Real Tol,
180 const Standard_Real Uinf,
181 const Standard_Real Usup)
183 Perform(P, C, Tol, Uinf, Usup);
188 void Extrema_ExtPElC::Perform (const gp_Pnt& P,
190 const Standard_Real Tol,
191 const Standard_Real Uinf,
192 const Standard_Real Usup)
193 /*-----------------------------------------------------------------------------
195 Find values of parameter u so that:
196 - dist(P,C(u)) passes by an extremum,
201 1- Projection of point P in the plane of the ellipse,
202 2- Calculation of the solutions:
203 Let Pp, the projected point; find values u so that:
204 (C(u)-Pp).C'(u) = 0. (1)
205 Let Cos = cos(u) and Sin = sin(u),
206 C(u) = (A*Cos,B*Sin) and Pp = (X,Y);
207 Then, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
208 (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
209 Use algorithm math_TrigonometricFunctionRoots to solve this equation.
210 -----------------------------------------------------------------------------*/
212 myDone = Standard_False;
215 // 1- Projection of point P in the plane of the ellipse -> Pp ...
217 gp_Pnt O = C.Location();
218 gp_Vec Axe (C.Axis().Direction());
219 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
220 gp_Pnt Pp = P.Translated(Trsl);
222 // 2- Calculation of solutions ...
224 Standard_Integer NoSol, NbSol;
225 Standard_Real A = C.MajorRadius();
226 Standard_Real B = C.MinorRadius();
228 Standard_Real OPpMagn = OPp.Magnitude();
229 if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
230 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
231 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
232 // Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
234 Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
235 if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
237 // math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
238 math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
240 if (!Sol.IsDone()) { return; }
243 NbSol = Sol.NbSolutions();
244 for (NoSol = 1; NoSol <= NbSol; NoSol++) {
245 Us = Sol.Value(NoSol);
246 Cu = ElCLib::Value(Us,C);
247 mySqDist[myNbExt] = Cu.SquareDistance(P);
248 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
249 Cu = ElCLib::Value(Us + 0.1, C);
250 myIsMin[myNbExt] = mySqDist[myNbExt] < Cu.SquareDistance(P);
253 myDone = Standard_True;
255 //=============================================================================
257 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
259 const Standard_Real Tol,
260 const Standard_Real Uinf,
261 const Standard_Real Usup)
263 Perform(P, C, Tol, Uinf, Usup);
267 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
269 const Standard_Real Tol,
270 const Standard_Real Uinf,
271 const Standard_Real Usup)
272 /*-----------------------------------------------------------------------------
274 Find values of parameter u so that:
275 - dist(P,C(u)) passes by an extremum,
280 1- Projection of point P in the plane of the hyperbola,
281 2- Calculation of solutions:
282 Let Pp, le point projete; on recherche les valeurs u telles que:
283 (C(u)-Pp).C'(u) = 0. (1)
284 Let R and r be the radiuses of the hyperbola,
285 Chu = Cosh(u) and Shu = Sinh(u),
286 C(u) = (R*Chu,r*Shu) and Pp = (X,Y);
287 Then, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
288 (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
290 Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
291 (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
293 ((X*R-Y*r)/2.) * v**(-1) = 0.
294 (2)* v**2 <=> ((R**2+r**2)/4.) * v**4 -
295 ((X*R+Y*r)/2.) * v**3 +
297 ((R**2+r**2)/4.) = 0.
298 Use algorithm math_DirectPolynomialRoots to solve this equation by v.
299 -----------------------------------------------------------------------------*/
301 myDone = Standard_False;
304 // 1- Projection of point P in the plane of hyperbola -> Pp ...
306 gp_Pnt O = C.Location();
307 gp_Vec Axe (C.Axis().Direction());
308 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
309 gp_Pnt Pp = P.Translated(Trsl);
311 // 2- Calculation of solutions ...
313 Standard_Real Tol2 = Tol * Tol;
314 Standard_Real R = C.MajorRadius();
315 Standard_Real r = C.MinorRadius();
317 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
318 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
320 Standard_Real C1 = (R*R+r*r)/4.;
321 math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
322 if (!Sol.IsDone()) { return; }
324 Standard_Real Us, Vs;
325 Standard_Integer NbSol = Sol.NbSolutions();
326 Standard_Boolean DejaEnr;
327 Standard_Integer NoExt;
329 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
330 Vs = Sol.Value(NoSol);
333 if ((Us >= Uinf) && (Us <= Usup)) {
334 Cu = ElCLib::Value(Us,C);
335 DejaEnr = Standard_False;
336 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
337 if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
338 DejaEnr = Standard_True;
344 mySqDist[myNbExt] = Cu.SquareDistance(P);
345 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
346 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
349 } // if ((Us >= Uinf) && (Us <= Usup))
351 } // for (Standard_Integer NoSol = 1; ...
352 myDone = Standard_True;
354 //=============================================================================
356 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
358 const Standard_Real Tol,
359 const Standard_Real Uinf,
360 const Standard_Real Usup)
362 Perform(P, C, Tol, Uinf, Usup);
366 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
368 // const Standard_Real Tol,
369 const Standard_Real ,
370 const Standard_Real Uinf,
371 const Standard_Real Usup)
372 /*-----------------------------------------------------------------------------
374 Find values of parameter u so that:
375 - dist(P,C(u)) pass by an extremum,
380 1- Projection of point P in the plane of the parabola,
381 2- Calculation of solutions:
382 Let Pp, the projected point; find values u so that:
383 (C(u)-Pp).C'(u) = 0. (1)
384 Let F the focus of the parabola,
385 C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
386 Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
387 (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
388 Use algorithm math_DirectPolynomialRoots to solve this equation by U.
389 -----------------------------------------------------------------------------*/
391 myDone = Standard_False;
394 // 1- Projection of point P in the plane of the parabola -> Pp ...
396 gp_Pnt O = C.Location();
397 gp_Vec Axe (C.Axis().Direction());
398 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
399 gp_Pnt Pp = P.Translated(Trsl);
401 // 2- Calculation of solutions ...
403 Standard_Real F = C.Focal();
405 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
406 // Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
407 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
408 math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
409 if (!Sol.IsDone()) { return; }
412 Standard_Integer NbSol = Sol.NbSolutions();
413 Standard_Boolean DejaEnr;
414 Standard_Integer NoExt;
416 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
417 Us = Sol.Value(NoSol);
418 if ((Us >= Uinf) && (Us <= Usup)) {
419 Cu = ElCLib::Value(Us,C);
420 DejaEnr = Standard_False;
421 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
422 if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
423 DejaEnr = Standard_True;
429 mySqDist[myNbExt] = Cu.SquareDistance(P);
430 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
431 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
434 } // if ((Us >= Uinf) && (Us <= Usup))
435 } // for (Standard_Integer NoSol = 1; ...
436 myDone = Standard_True;
438 //=============================================================================
440 Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
441 //=============================================================================
443 Standard_Integer Extrema_ExtPElC::NbExt () const
445 if (!IsDone()) { throw StdFail_NotDone(); }
448 //=============================================================================
450 Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
452 if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
453 return mySqDist[N-1];
455 //=============================================================================
457 Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
459 if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
462 //=============================================================================
464 const Extrema_POnCurv& Extrema_ExtPElC::Point (const Standard_Integer N) const
466 if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
469 //=============================================================================