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 #include <Standard_TypeMismatch.hxx>
16 #include <Precision.hxx>
17 static const Standard_Real TolFactor = 1.e-12;
18 static const Standard_Real MinTol = 1.e-20;
19 static const Standard_Real MinStep = 1e-7;
21 static const Standard_Integer MaxOrder = 3;
25 /*-----------------------------------------------------------------------------
26 Fonction permettant de rechercher une distance extremale entre un point P et
27 une courbe C (en partant d'un point approche C(u0)).
28 Cette classe herite de math_FunctionWithDerivative et est utilisee par
29 les algorithmes math_FunctionRoot et math_FunctionRoots.
30 Si on note D1c et D2c les derivees premiere et seconde:
31 { F(u) = (C(u)-P).D1c(u)/ ||Du||}
32 { DF(u) = ||Du|| + (C(u)-P).D2c(u)/||Du|| - F(u)*Duu*Du/||Du||**2
35 { F(u) = (C(u)-P).D1c(u) }
36 { DF(u) = D1c(u).D1c(u) + (C(u)-P).D2c(u)
37 = ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
38 ----------------------------------------------------------------------------*/
40 Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
42 const Standard_Integer NPoint = 10;
43 const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
45 Standard_Integer aNum = 0;
46 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
47 //(it is computed with using NPoint point)
51 Standard_Real u = myUinfium + aNum*aStep; //parameter for every point
55 Pnt Ptemp; //empty point (is not used below)
56 Vec VDer; // 1st derivative vector
57 Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
59 if(Precision::IsInfinite(VDer.X()) || Precision::IsInfinite(VDer.Y()))
64 Standard_Real vm = VDer.Magnitude();
68 while(++aNum < NPoint+1);
70 return Max(aMax*TolFactor,MinTol);
74 //=============================================================================
76 Extrema_FuncExtPC::Extrema_FuncExtPC():
80 myPinit = Standard_False;
81 myCinit = Standard_False;
82 myD1Init = Standard_False;
84 SubIntervalInitialize(0.0,0.0);
89 //=============================================================================
91 Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
92 const Curve& C): myU(0.), myD1f(0.)
95 myC = (Standard_Address)&C;
96 myPinit = Standard_True;
97 myCinit = Standard_True;
98 myD1Init = Standard_False;
100 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
101 Tool::LastParameter(*((Curve*)myC)));
103 switch(Tool::GetType(*((Curve*)myC)))
105 case GeomAbs_BezierCurve:
106 case GeomAbs_BSplineCurve:
107 case GeomAbs_OffsetCurve:
108 case GeomAbs_OtherCurve:
109 myMaxDerivOrder = MaxOrder;
110 myTol = SearchOfTolerance();
118 //=============================================================================
120 void Extrema_FuncExtPC::Initialize(const Curve& C)
122 myC = (Standard_Address)&C;
123 myCinit = Standard_True;
128 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
129 Tool::LastParameter(*((Curve*)myC)));
131 switch(Tool::GetType(*((Curve*)myC)))
133 case GeomAbs_BezierCurve:
134 case GeomAbs_BSplineCurve:
135 case GeomAbs_OffsetCurve:
136 case GeomAbs_OtherCurve:
137 myMaxDerivOrder = MaxOrder;
138 myTol = SearchOfTolerance();
147 //=============================================================================
149 void Extrema_FuncExtPC::SetPoint(const Pnt& P)
152 myPinit = Standard_True;
158 //=============================================================================
161 Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
163 if (!myPinit || !myCinit)
164 throw Standard_TypeMismatch("No init");
168 Tool::D1(*((Curve*)myC),myU,myPc,D1c);
170 if(Precision::IsInfinite(D1c.X()) || Precision::IsInfinite(D1c.Y()))
172 F = Precision::Infinite();
173 return Standard_False;
176 Standard_Real Ndu = D1c.Magnitude();
178 if(myMaxDerivOrder != 0)
180 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
182 const Standard_Real DivisionFactor = 1.e-3;
184 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
187 du = myUsupremum-myUinfium;
189 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
190 //Derivative is approximated by Taylor-series
192 Standard_Integer n = 1; //Derivative order
194 Standard_Boolean IsDeriveFound;
198 V = Tool::DN(*((Curve*)myC),myU,++n);
200 IsDeriveFound = (Ndu > myTol);
202 while(!IsDeriveFound && n < myMaxDerivOrder);
208 if(myU-myUinfium < aDelta)
214 Tool::D0(*((Curve*)myC),Min(myU, u),P1);
215 Tool::D0(*((Curve*)myC),Max(myU, u),P2);
218 Standard_Real aDirFactor = V.Dot(V1);
227 //Derivative is approximated by three points
229 Pnt Ptemp; //(0,0,0)-coordinate
231 Standard_Boolean IsParameterGrown;
233 if(myU-myUinfium < 2*aDelta)
235 Tool::D0(*((Curve*)myC),myU,P1);
236 Tool::D0(*((Curve*)myC),myU+aDelta,P2);
237 Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
238 IsParameterGrown = Standard_True;
242 Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
243 Tool::D0(*((Curve*)myC),myU-aDelta,P2);
244 Tool::D0(*((Curve*)myC),myU,P3);
245 IsParameterGrown = Standard_False;
248 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
255 Ndu = D1c.Magnitude();
256 }//(if (Ndu <= myTol)) condition
257 }//if(myMaxDerivOrder != 0)
261 //Warning: 1st derivative is equal to zero!
262 return Standard_False;
266 F = PPc.Dot(D1c)/Ndu;
267 return Standard_True;
270 //=============================================================================
272 Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_Real& D1f)
274 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
276 return Values(U,F,D1f); /* on fait appel a Values pour simplifier la
277 sauvegarde de l'etat. */
279 //=============================================================================
281 Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U,
285 if (!myPinit || !myCinit)
286 throw Standard_TypeMismatch("No init");
288 Pnt myPc_old = myPc, myP_old = myP;
290 if(Value(U,F) == Standard_False)
292 //Warning: No function value found!;
294 myD1Init = Standard_False;
295 return Standard_False;
303 Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
305 Standard_Real Ndu = D1c.Magnitude();
306 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
308 //Derivative is approximated by three points
310 //Attention: aDelta value must be greater than same value for "Value(...)"
311 // function to avoid of points' collisions.
312 const Standard_Real DivisionFactor = 0.01;
314 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
317 du = myUsupremum-myUinfium;
319 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
321 Standard_Real F1, F2, F3;
323 if(myU-myUinfium < 2*aDelta)
326 //const Standard_Real U1 = myU;
327 const Standard_Real U2 = myU + aDelta;
328 const Standard_Real U3 = myU + aDelta * 2.0;
330 if(!((Value(U2,F2)) && (Value(U3,F3))))
332 //There are many points close to singularity points and
333 //which have zero-derivative. Try to decrease aDelta variable's value.
335 myD1Init = Standard_False;
336 return Standard_False;
339 //After calling of Value(...) function variable myU will be redeterminated.
340 //So we must return it previous value.
341 D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
346 const Standard_Real U1 = myU - aDelta * 2.0;
347 const Standard_Real U2 = myU - aDelta;
348 //const Standard_Real U3 = myU;
350 if(!((Value(U2,F2)) && (Value(U1,F1))))
352 //There are many points close to singularity points and
353 //which have zero-derivative. Try to decrease aDelta variable's value.
354 myD1Init = Standard_False;
355 return Standard_False;
357 //After calling of Value(...) function variable myU will be redeterminated.
358 //So we must return it previous value.
359 D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
368 D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
373 myD1Init = Standard_True;
374 return Standard_True;
376 //=============================================================================
378 Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
380 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
381 mySqDist.Append(myPc.SquareDistance(myP));
383 // It is necessary to always compute myD1f.
384 myD1Init = Standard_True;
385 Standard_Real FF, DD;
388 Standard_Integer IntVal = 0;
394 myIsMin.Append(IntVal);
395 myPoint.Append(POnC(myU,myPc));
398 //=============================================================================
400 Standard_Integer Extrema_FuncExtPC::NbExt () const { return mySqDist.Length(); }
401 //=============================================================================
403 Standard_Real Extrema_FuncExtPC::SquareDistance (const Standard_Integer N) const
405 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
406 return mySqDist.Value(N);
408 //=============================================================================
409 Standard_Boolean Extrema_FuncExtPC::IsMin (const Standard_Integer N) const
411 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
412 return (myIsMin.Value(N) == 1);
414 //=============================================================================
415 const POnC & Extrema_FuncExtPC::Point (const Standard_Integer N) const
417 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
418 return myPoint.Value(N);
420 //=============================================================================
422 void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
424 myUinfium = theUfirst;
425 myUsupremum = theUlast;