1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
19 #include <Extrema_ExtPElC.ixx>
20 #include <StdFail_NotDone.hxx>
21 #include <math_DirectPolynomialRoots.hxx>
22 #include <math_TrigonometricFunctionRoots.hxx>
24 #include <Standard_OutOfRange.hxx>
25 #include <Standard_NotImplemented.hxx>
26 #include <Precision.hxx>
29 //=============================================================================
31 Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; }
32 //=============================================================================
34 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
36 const Standard_Real Tol,
37 const Standard_Real Uinf,
38 const Standard_Real Usup)
40 Perform(P, L, Tol, Uinf, Usup);
43 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
45 const Standard_Real Tol,
46 const Standard_Real Uinf,
47 const Standard_Real Usup)
49 myDone = Standard_False;
51 gp_Vec V1 = gp_Vec(L.Direction());
52 gp_Pnt OR = L.Location();
54 Standard_Real Mydist = V1.Dot(V);
55 if ((Mydist >= Uinf-Tol) &&
56 (Mydist <= Usup+Tol)){
58 gp_Pnt MyP = OR.Translated(Mydist*V1);
59 Extrema_POnCurv MyPOnCurve(Mydist, MyP);
60 mySqDist[0] = P.SquareDistance(MyP);
61 myPoint[0] = MyPOnCurve;
62 myIsMin[0] = Standard_True;
64 myDone = Standard_True;
71 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
73 const Standard_Real Tol,
74 const Standard_Real Uinf,
75 const Standard_Real Usup)
77 Perform(P, C, Tol, Uinf, Usup);
80 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
82 const Standard_Real Tol,
83 const Standard_Real Uinf,
84 const Standard_Real Usup)
85 /*-----------------------------------------------------------------------------
87 Find values of parameter u such as:
88 - dist(P,C(u)) pass by an extrema,
93 1- Projection of point P in the plane of the circle,
94 2- Calculation of u solutions in [0.,2.*M_PI]:
95 Let Pp, the projected point and
96 O, the center of the circle;
98 - if Pp is mixed with 0, there is an infinite number of solutions;
99 IsDone() renvoie Standard_False.
100 - otherwise, 2 points are solutions for the complete circle:
101 . Us1 = angle(OPp,OX) corresponds to the minimum,
102 . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
103 ( Us1 - M_PI otherwise;
104 Us2 corresponds to the maximum.
105 3- Calculate the extrema in [Uinf,Usup].
106 -----------------------------------------------------------------------------*/
108 myDone = Standard_False;
111 // 1- Projection of the point P in the plane of circle -> Pp ...
113 gp_Pnt O = C.Location();
114 gp_Vec Axe (C.Axis().Direction());
115 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
116 gp_Pnt Pp = P.Translated(Trsl);
118 // 2- Calculate u solutions in [0.,2.*PI] ...
121 if (OPp.Magnitude() < Tol) { return; }
122 Standard_Real Usol[2];
123 Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
124 Usol[1] = Usol[0] + M_PI;
126 Standard_Real myuinf = Uinf;
127 //modified by NIZNHY-PKV Fri Apr 20 15:03:28 2001 f
128 //Standard_Real TolU = Tol*C.Radius();
129 Standard_Real TolU, aR;
131 TolU=Precision::Infinite();
132 if (aR > gp::Resolution()) {
135 //modified by NIZNHY-PKV Fri Apr 20 15:03:32 2001 t
136 ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
137 ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
138 if (((Usol[0]-2*M_PI-Uinf) < TolU) && ((Usol[0]-2*M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
139 if (((Usol[1]-2*M_PI-Uinf) < TolU) && ((Usol[1]-2*M_PI-Uinf) > -TolU)) Usol[1] = Uinf;
142 // 3- Calculate extrema in [Umin,Umax] ...
146 for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
148 if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
149 Cu = ElCLib::Value(Us,C);
150 mySqDist[myNbExt] = Cu.SquareDistance(P);
151 myIsMin[myNbExt] = (NoSol == 0);
152 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
156 myDone = Standard_True;
158 //=============================================================================
160 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
162 const Standard_Real Tol,
163 const Standard_Real Uinf,
164 const Standard_Real Usup)
166 Perform(P, C, Tol, Uinf, Usup);
171 void Extrema_ExtPElC::Perform (const gp_Pnt& P,
173 const Standard_Real Tol,
174 const Standard_Real Uinf,
175 const Standard_Real Usup)
176 /*-----------------------------------------------------------------------------
178 Find values of parameter u so that:
179 - dist(P,C(u)) passes by an extremum,
184 1- Projection of point P in the plane of the ellipse,
185 2- Calculation of the solutions:
186 Let Pp, the projected point; find values u so that:
187 (C(u)-Pp).C'(u) = 0. (1)
188 Let Cos = cos(u) and Sin = sin(u),
189 C(u) = (A*Cos,B*Sin) and Pp = (X,Y);
190 Then, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
191 (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
192 Use algorithm math_TrigonometricFunctionRoots to solve this equation.
193 -----------------------------------------------------------------------------*/
195 myDone = Standard_False;
198 // 1- Projection of point P in the plane of the ellipse -> Pp ...
200 gp_Pnt O = C.Location();
201 gp_Vec Axe (C.Axis().Direction());
202 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
203 gp_Pnt Pp = P.Translated(Trsl);
205 // 2- Calculation of solutions ...
207 Standard_Integer NoSol, NbSol;
208 Standard_Real A = C.MajorRadius();
209 Standard_Real B = C.MinorRadius();
211 Standard_Real OPpMagn = OPp.Magnitude();
212 if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
213 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
214 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
215 // Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
217 Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
218 if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
220 // math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
221 math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
223 if (!Sol.IsDone()) { return; }
226 NbSol = Sol.NbSolutions();
227 for (NoSol = 1; NoSol <= NbSol; NoSol++) {
228 Us = Sol.Value(NoSol);
229 Cu = ElCLib::Value(Us,C);
230 mySqDist[myNbExt] = Cu.SquareDistance(P);
231 myIsMin[myNbExt] = (NoSol == 1);
232 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
235 myDone = Standard_True;
237 //=============================================================================
239 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
241 const Standard_Real Tol,
242 const Standard_Real Uinf,
243 const Standard_Real Usup)
245 Perform(P, C, Tol, Uinf, Usup);
249 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
251 const Standard_Real Tol,
252 const Standard_Real Uinf,
253 const Standard_Real Usup)
254 /*-----------------------------------------------------------------------------
256 Find values of parameter u so that:
257 - dist(P,C(u)) passes by an extremum,
262 1- Projection of point P in the plane of the hyperbola,
263 2- Calculation of solutions:
264 Let Pp, le point projete; on recherche les valeurs u telles que:
265 (C(u)-Pp).C'(u) = 0. (1)
266 Let R and r be the radiuses of the hyperbola,
267 Chu = Cosh(u) and Shu = Sinh(u),
268 C(u) = (R*Chu,r*Shu) and Pp = (X,Y);
269 Then, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
270 (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
272 Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
273 (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
275 ((X*R-Y*r)/2.) * v**(-1) = 0.
276 (2)* v**2 <=> ((R**2+r**2)/4.) * v**4 -
277 ((X*R+Y*r)/2.) * v**3 +
279 ((R**2+r**2)/4.) = 0.
280 Use algorithm math_DirectPolynomialRoots to solve this equation by v.
281 -----------------------------------------------------------------------------*/
283 myDone = Standard_False;
286 // 1- Projection of point P in the plane of hyperbola -> Pp ...
288 gp_Pnt O = C.Location();
289 gp_Vec Axe (C.Axis().Direction());
290 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
291 gp_Pnt Pp = P.Translated(Trsl);
293 // 2- Calculation of solutions ...
295 Standard_Real Tol2 = Tol * Tol;
296 Standard_Real R = C.MajorRadius();
297 Standard_Real r = C.MinorRadius();
299 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
300 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
302 Standard_Real C1 = (R*R+r*r)/4.;
303 math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
304 if (!Sol.IsDone()) { return; }
306 Standard_Real Us, Vs;
307 Standard_Integer NbSol = Sol.NbSolutions();
308 Standard_Boolean DejaEnr;
309 Standard_Integer NoExt;
311 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
312 Vs = Sol.Value(NoSol);
315 if ((Us >= Uinf) && (Us <= Usup)) {
316 Cu = ElCLib::Value(Us,C);
317 DejaEnr = Standard_False;
318 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
319 if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
320 DejaEnr = Standard_True;
326 mySqDist[myNbExt] = Cu.SquareDistance(P);
327 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
328 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
331 } // if ((Us >= Uinf) && (Us <= Usup))
333 } // for (Standard_Integer NoSol = 1; ...
334 myDone = Standard_True;
336 //=============================================================================
338 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
340 const Standard_Real Tol,
341 const Standard_Real Uinf,
342 const Standard_Real Usup)
344 Perform(P, C, Tol, Uinf, Usup);
348 void Extrema_ExtPElC::Perform(const gp_Pnt& P,
350 // const Standard_Real Tol,
351 const Standard_Real ,
352 const Standard_Real Uinf,
353 const Standard_Real Usup)
354 /*-----------------------------------------------------------------------------
356 Find values of parameter u so that:
357 - dist(P,C(u)) pass by an extremum,
362 1- Projection of point P in the plane of the parabola,
363 2- Calculation of solutions:
364 Let Pp, the projected point; find values u so that:
365 (C(u)-Pp).C'(u) = 0. (1)
366 Let F the focus of the parabola,
367 C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
368 Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
369 (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
370 Use algorithm math_DirectPolynomialRoots to solve this equation by U.
371 -----------------------------------------------------------------------------*/
373 myDone = Standard_False;
376 // 1- Projection of point P in the plane of the parabola -> Pp ...
378 gp_Pnt O = C.Location();
379 gp_Vec Axe (C.Axis().Direction());
380 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
381 gp_Pnt Pp = P.Translated(Trsl);
383 // 2- Calculation of solutions ...
385 Standard_Real F = C.Focal();
387 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
388 // Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
389 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
390 math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
391 if (!Sol.IsDone()) { return; }
394 Standard_Integer NbSol = Sol.NbSolutions();
395 Standard_Boolean DejaEnr;
396 Standard_Integer NoExt;
398 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
399 Us = Sol.Value(NoSol);
400 if ((Us >= Uinf) && (Us <= Usup)) {
401 Cu = ElCLib::Value(Us,C);
402 DejaEnr = Standard_False;
403 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
404 if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
405 DejaEnr = Standard_True;
411 mySqDist[myNbExt] = Cu.SquareDistance(P);
412 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
413 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
416 } // if ((Us >= Uinf) && (Us <= Usup))
417 } // for (Standard_Integer NoSol = 1; ...
418 myDone = Standard_True;
420 //=============================================================================
422 Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
423 //=============================================================================
425 Standard_Integer Extrema_ExtPElC::NbExt () const
427 if (!IsDone()) { StdFail_NotDone::Raise(); }
430 //=============================================================================
432 Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
434 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
435 return mySqDist[N-1];
437 //=============================================================================
439 Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
441 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
444 //=============================================================================
446 Extrema_POnCurv Extrema_ExtPElC::Point (const Standard_Integer N) const
448 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
451 //=============================================================================