CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
3//
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
8//
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.
11//
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.
18
7fd59977 19#include <Standard_TypeMismatch.hxx>
408ecc39 20#include <Precision.hxx>
32ca7a51 21static const Standard_Real TolFactor = 1.e-12;
22static const Standard_Real MinTol = 1.e-20;
23static const Standard_Real MinStep = 1e-7;
24
25static const Standard_Integer MaxOrder = 3;
26
27
7fd59977 28
29/*-----------------------------------------------------------------------------
28cec2ba 30Fonction permettant de rechercher une distance extremale entre un point P et
7fd59977 31une courbe C (en partant d'un point approche C(u0)).
28cec2ba 32Cette classe herite de math_FunctionWithDerivative et est utilisee par
7fd59977 33les algorithmes math_FunctionRoot et math_FunctionRoots.
28cec2ba 34Si on note D1c et D2c les derivees premiere et seconde:
35{ F(u) = (C(u)-P).D1c(u)/ ||Du||}
36{ DF(u) = ||Du|| + (C(u)-P).D2c(u)/||Du|| - F(u)*Duu*Du/||Du||**2
7fd59977 37
38
28cec2ba 39{ F(u) = (C(u)-P).D1c(u) }
40{ DF(u) = D1c(u).D1c(u) + (C(u)-P).D2c(u)
41= ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
7fd59977 42----------------------------------------------------------------------------*/
43
32ca7a51 44Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
28cec2ba 45{
32ca7a51 46 const Standard_Integer NPoint = 10;
47 const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
28cec2ba 48
32ca7a51 49 Standard_Integer aNum = 0;
50 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
28cec2ba 51 //(it is computed with using NPoint point)
52
32ca7a51 53 do
28cec2ba 54 {
32ca7a51 55 Standard_Real u = myUinfium + aNum*aStep; //parameter for every point
56 if(u > myUsupremum)
57 u = myUsupremum;
28cec2ba 58
32ca7a51 59 Pnt Ptemp; //empty point (is not used below)
60 Vec VDer; // 1st derivative vector
61 Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
62 Standard_Real vm = VDer.Magnitude();
63 if(vm > aMax)
64 aMax = vm;
28cec2ba 65 }
32ca7a51 66 while(++aNum < NPoint+1);
28cec2ba 67
32ca7a51 68 return Max(aMax*TolFactor,MinTol);
28cec2ba 69
70}
32ca7a51 71
72//=============================================================================
7fd59977 73
74Extrema_FuncExtPC::Extrema_FuncExtPC():
75myU(0.),
76myD1f(0.)
77{
78 myPinit = Standard_False;
79 myCinit = Standard_False;
80 myD1Init = Standard_False;
28cec2ba 81
32ca7a51 82 SubIntervalInitialize(0.0,0.0);
83 myMaxDerivOrder = 0;
84 myTol=MinTol;
7fd59977 85}
86
87//=============================================================================
88
89Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
28cec2ba 90 const Curve& C): myU(0.), myD1f(0.)
91{
7fd59977 92 myP = P;
94 myPinit = Standard_True;
95 myCinit = Standard_True;
96 myD1Init = Standard_False;
28cec2ba 97
32ca7a51 98 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
28cec2ba 99 Tool::LastParameter(*((Curve*)myC)));
100
32ca7a51 101 switch(Tool::GetType(*((Curve*)myC)))
28cec2ba 102 {
103 case GeomAbs_BezierCurve:
104 case GeomAbs_BSplineCurve:
105 case GeomAbs_OtherCurve:
106 myMaxDerivOrder = MaxOrder;
107 myTol = SearchOfTolerance();
108 break;
109 default:
110 myMaxDerivOrder = 0;
111 myTol=MinTol;
112 break;
32ca7a51 113 }
28cec2ba 114}
7fd59977 115//=============================================================================
116
117void Extrema_FuncExtPC::Initialize(const Curve& C)
28cec2ba 118{
120 myCinit = Standard_True;
121 myPoint.Clear();
122 mySqDist.Clear();
123 myIsMin.Clear();
28cec2ba 124
32ca7a51 125 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
28cec2ba 126 Tool::LastParameter(*((Curve*)myC)));
127
32ca7a51 128 switch(Tool::GetType(*((Curve*)myC)))
28cec2ba 129 {
130 case GeomAbs_BezierCurve:
131 case GeomAbs_BSplineCurve:
132 case GeomAbs_OtherCurve:
133 myMaxDerivOrder = MaxOrder;
134 myTol = SearchOfTolerance();
135 break;
136 default:
137 myMaxDerivOrder = 0;
138 myTol=MinTol;
139 break;
32ca7a51 140 }
28cec2ba 141}
7fd59977 142
143//=============================================================================
144
145void Extrema_FuncExtPC::SetPoint(const Pnt& P)
146{
147 myP = P;
148 myPinit = Standard_True;
149 myPoint.Clear();
150 mySqDist.Clear();
151 myIsMin.Clear();
152}
153
154//=============================================================================
155
32ca7a51 156
7fd59977 157Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
158{
32ca7a51 159 if (!myPinit || !myCinit)
160 Standard_TypeMismatch::Raise("No init");
28cec2ba 161
7fd59977 162 myU = U;
163 Vec D1c;
164 Tool::D1(*((Curve*)myC),myU,myPc,D1c);
165 Standard_Real Ndu = D1c.Magnitude();
32ca7a51 166
167 if(myMaxDerivOrder != 0)
28cec2ba 168 {
32ca7a51 169 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
28cec2ba 170 {
32ca7a51 171 const Standard_Real DivisionFactor = 1.e-3;
172 Standard_Real du;
173 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
174 du = 0.0;
175 else
176 du = myUsupremum-myUinfium;
28cec2ba 177
32ca7a51 178 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
28cec2ba 179 //Derivative is approximated by Taylor-series
180
32ca7a51 181 Standard_Integer n = 1; //Derivative order
182 Vec V;
183 Standard_Boolean IsDeriveFound;
28cec2ba 184
32ca7a51 185 do
28cec2ba 186 {
32ca7a51 187 V = Tool::DN(*((Curve*)myC),myU,++n);
188 Ndu = V.Magnitude();
189 IsDeriveFound = (Ndu > myTol);
28cec2ba 190 }
32ca7a51 191 while(!IsDeriveFound && n < myMaxDerivOrder);
28cec2ba 192
32ca7a51 193 if(IsDeriveFound)
28cec2ba 194 {
32ca7a51 195 Standard_Real u;
28cec2ba 196
199 else
28cec2ba 201
32ca7a51 202 Pnt P1, P2;
203 Tool::D0(*((Curve*)myC),Min(myU, u),P1);
204 Tool::D0(*((Curve*)myC),Max(myU, u),P2);
28cec2ba 205
32ca7a51 206 Vec V1(P1,P2);
28cec2ba 208
210 D1c = -V;
211 else
212 D1c = V;
28cec2ba 213 }//if(IsDeriveFound)
32ca7a51 214 else
28cec2ba 215 {
216 //Derivative is approximated by three points
32ca7a51 217
218 Pnt Ptemp; //(0,0,0)-coordinate
219 Pnt P1, P2, P3;
220 Standard_Boolean IsParameterGrown;
28cec2ba 221
28cec2ba 223 {
32ca7a51 224 Tool::D0(*((Curve*)myC),myU,P1);
227 IsParameterGrown = Standard_True;
28cec2ba 228 }
32ca7a51 229 else
28cec2ba 230 {
233 Tool::D0(*((Curve*)myC),myU,P3);
234 IsParameterGrown = Standard_False;
28cec2ba 235 }
236
32ca7a51 237 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
28cec2ba 238
32ca7a51 239 if(IsParameterGrown)
240 D1c=-3*V1+4*V2-V3;
241 else
242 D1c=V1-4*V2+3*V3;
28cec2ba 243 }
244 Ndu = D1c.Magnitude();
245 }//(if (Ndu <= myTol)) condition
246 }//if(myMaxDerivOrder != 0)
32ca7a51 247
248 if (Ndu <= MinTol)
28cec2ba 249 {
250 //Warning: 1st derivative is equal to zero!
32ca7a51 251 return Standard_False;
28cec2ba 252 }
253
7fd59977 254 Vec PPc (myP,myPc);
255 F = PPc.Dot(D1c)/Ndu;
256 return Standard_True;
257}
258
259//=============================================================================
260
261Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_Real& D1f)
262{
263 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
264 Standard_Real F;
265 return Values(U,F,D1f); /* on fait appel a Values pour simplifier la
28cec2ba 266 sauvegarde de l'etat. */
7fd59977 267}
268//=============================================================================
269
28cec2ba 270Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U,
271 Standard_Real& F,
272 Standard_Real& D1f)
273{
32ca7a51 274 if (!myPinit || !myCinit)
275 Standard_TypeMismatch::Raise("No init");
28cec2ba 276
32ca7a51 277 Pnt myPc_old = myPc, myP_old = myP;
278
279 if(Value(U,F) == Standard_False)
28cec2ba 280 {
281 //Warning: No function value found!;
32ca7a51 282
283 myD1Init = Standard_False;
284 return Standard_False;
28cec2ba 285 }
286
7fd59977 287 myU = U;
32ca7a51 288 myPc = myPc_old;
289 myP = myP_old;
28cec2ba 290
7fd59977 291 Vec D1c,D2c;
292 Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
293
294 Standard_Real Ndu = D1c.Magnitude();
32ca7a51 295 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
28cec2ba 296 {
297 //Derivative is approximated by three points
32ca7a51 298
28cec2ba 299 //Attention: aDelta value must be greater than same value for "Value(...)"
300 // function to avoid of points' collisions.
32ca7a51 301 const Standard_Real DivisionFactor = 0.01;
302 Standard_Real du;
303 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
304 du = 0.0;
305 else
306 du = myUsupremum-myUinfium;
28cec2ba 307
32ca7a51 308 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
28cec2ba 309
32ca7a51 310 Standard_Real F1, F2, F3;
28cec2ba 311
28cec2ba 313 {
32ca7a51 314 F1=F;
cc9d78db 315 //const Standard_Real U1 = myU;
316 const Standard_Real U2 = myU + aDelta;
317 const Standard_Real U3 = myU + aDelta * 2.0;
28cec2ba 318
32ca7a51 319 if(!((Value(U2,F2)) && (Value(U3,F3))))
28cec2ba 320 {
321 //There are many points close to singularity points and
322 //which have zero-derivative. Try to decrease aDelta variable's value.
323
32ca7a51 324 myD1Init = Standard_False;
325 return Standard_False;
32ca7a51 326 }
28cec2ba 327
328 //After calling of Value(...) function variable myU will be redeterminated.
329 //So we must return it previous value.
331 }
32ca7a51 332 else
28cec2ba 333 {
32ca7a51 334 F3 = F;
cc9d78db 335 const Standard_Real U1 = myU - aDelta * 2.0;
336 const Standard_Real U2 = myU - aDelta;
337 //const Standard_Real U3 = myU;
338
32ca7a51 339 if(!((Value(U2,F2)) && (Value(U1,F1))))
28cec2ba 340 {
341 //There are many points close to singularity points and
342 //which have zero-derivative. Try to decrease aDelta variable's value.
32ca7a51 343 myD1Init = Standard_False;
344 return Standard_False;
32ca7a51 345 }
28cec2ba 346 //After calling of Value(...) function variable myU will be redeterminated.
347 //So we must return it previous value.
32ca7a51 349 }
28cec2ba 350 myU = U;
351 myPc = myPc_old;
352 myP = myP_old;
353 }
32ca7a51 354 else
28cec2ba 355 {
32ca7a51 356 Vec PPc (myP,myPc);
357 D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
28cec2ba 358 }
7fd59977 359
360 myD1f = D1f;
32ca7a51 361
7fd59977 362 myD1Init = Standard_True;
363 return Standard_True;
28cec2ba 364}
7fd59977 365//=============================================================================
366
367Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
368{
369 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
370 mySqDist.Append(myPc.SquareDistance(myP));
371 Standard_Integer IntVal;
372 if (!myD1Init) {
373 myD1Init = Standard_True;
374 Standard_Real FF, DD;
375 Values(myU, FF, DD);
376 }
377 if (!myD1Init) IntVal = 0;
378 else {
379 if (myD1f > 0.) { IntVal = 1; }
380 else { IntVal = 0; }
381 }
382 myIsMin.Append(IntVal);
383 myPoint.Append(POnC(myU,myPc));
384 return 0;
385}
386//=============================================================================
387
388Standard_Integer Extrema_FuncExtPC::NbExt () const { return mySqDist.Length(); }
389//=============================================================================
390
391Standard_Real Extrema_FuncExtPC::SquareDistance (const Standard_Integer N) const
392{
393 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
394 return mySqDist.Value(N);
395}
396//=============================================================================
397Standard_Boolean Extrema_FuncExtPC::IsMin (const Standard_Integer N) const
398{
399 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
400 return (myIsMin.Value(N) == 1);
401}
402//=============================================================================
5d99f2c8 403const POnC & Extrema_FuncExtPC::Point (const Standard_Integer N) const
7fd59977 404{
405 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
406 return myPoint.Value(N);
407}
408//=============================================================================
409
32ca7a51 410void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
28cec2ba 411{
32ca7a51 412 myUinfium = theUfirst;
413 myUsupremum = theUlast;
28cec2ba 414}