Warning, /include/opencascade/Extrema_FuncExtPC.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 #include <Standard_TypeMismatch.hxx>
0016 #include <Precision.hxx>
0017 static const Standard_Real TolFactor = 1.e-12;
0018 static const Standard_Real MinTol = 1.e-20;
0019 static const Standard_Real MinStep = 1e-7;
0020
0021 static const Standard_Integer MaxOrder = 3;
0022
0023
0024
0025 /*-----------------------------------------------------------------------------
0026 Fonction permettant de rechercher une distance extremale entre un point P et
0027 une courbe C (en partant d'un point approche C(u0)).
0028 Cette classe herite de math_FunctionWithDerivative et est utilisee par
0029 les algorithmes math_FunctionRoot et math_FunctionRoots.
0030 Si on note D1c et D2c les derivees premiere et seconde:
0031 { F(u) = (C(u)-P).D1c(u)/ ||Du||}
0032 { DF(u) = ||Du|| + (C(u)-P).D2c(u)/||Du|| - F(u)*Duu*Du/||Du||**2
0033
0034
0035 { F(u) = (C(u)-P).D1c(u) }
0036 { DF(u) = D1c(u).D1c(u) + (C(u)-P).D2c(u)
0037 = ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
0038 ----------------------------------------------------------------------------*/
0039
0040 Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
0041 {
0042 const Standard_Integer NPoint = 10;
0043 const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
0044
0045 Standard_Integer aNum = 0;
0046 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
0047 //(it is computed with using NPoint point)
0048
0049 do
0050 {
0051 Standard_Real u = myUinfium + aNum*aStep; //parameter for every point
0052 if(u > myUsupremum)
0053 u = myUsupremum;
0054
0055 Pnt Ptemp; //empty point (is not used below)
0056 Vec VDer; // 1st derivative vector
0057 Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
0058
0059 if(Precision::IsInfinite(VDer.X()) || Precision::IsInfinite(VDer.Y()))
0060 {
0061 continue;
0062 }
0063
0064 Standard_Real vm = VDer.Magnitude();
0065 if(vm > aMax)
0066 aMax = vm;
0067 }
0068 while(++aNum < NPoint+1);
0069
0070 return Max(aMax*TolFactor,MinTol);
0071
0072 }
0073
0074 //=============================================================================
0075
0076 Extrema_FuncExtPC::Extrema_FuncExtPC():
0077 myU(0.),
0078 myD1f(0.)
0079 {
0080 myPinit = Standard_False;
0081 myCinit = Standard_False;
0082 myD1Init = Standard_False;
0083
0084 SubIntervalInitialize(0.0,0.0);
0085 myMaxDerivOrder = 0;
0086 myTol=MinTol;
0087 }
0088
0089 //=============================================================================
0090
0091 Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
0092 const Curve& C): myU(0.), myD1f(0.)
0093 {
0094 myP = P;
0095 myC = (Standard_Address)&C;
0096 myPinit = Standard_True;
0097 myCinit = Standard_True;
0098 myD1Init = Standard_False;
0099
0100 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
0101 Tool::LastParameter(*((Curve*)myC)));
0102
0103 switch(Tool::GetType(*((Curve*)myC)))
0104 {
0105 case GeomAbs_BezierCurve:
0106 case GeomAbs_BSplineCurve:
0107 case GeomAbs_OffsetCurve:
0108 case GeomAbs_OtherCurve:
0109 myMaxDerivOrder = MaxOrder;
0110 myTol = SearchOfTolerance();
0111 break;
0112 default:
0113 myMaxDerivOrder = 0;
0114 myTol=MinTol;
0115 break;
0116 }
0117 }
0118 //=============================================================================
0119
0120 void Extrema_FuncExtPC::Initialize(const Curve& C)
0121 {
0122 myC = (Standard_Address)&C;
0123 myCinit = Standard_True;
0124 myPoint.Clear();
0125 mySqDist.Clear();
0126 myIsMin.Clear();
0127
0128 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
0129 Tool::LastParameter(*((Curve*)myC)));
0130
0131 switch(Tool::GetType(*((Curve*)myC)))
0132 {
0133 case GeomAbs_BezierCurve:
0134 case GeomAbs_BSplineCurve:
0135 case GeomAbs_OffsetCurve:
0136 case GeomAbs_OtherCurve:
0137 myMaxDerivOrder = MaxOrder;
0138 myTol = SearchOfTolerance();
0139 break;
0140 default:
0141 myMaxDerivOrder = 0;
0142 myTol=MinTol;
0143 break;
0144 }
0145 }
0146
0147 //=============================================================================
0148
0149 void Extrema_FuncExtPC::SetPoint(const Pnt& P)
0150 {
0151 myP = P;
0152 myPinit = Standard_True;
0153 myPoint.Clear();
0154 mySqDist.Clear();
0155 myIsMin.Clear();
0156 }
0157
0158 //=============================================================================
0159
0160
0161 Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
0162 {
0163 if (!myPinit || !myCinit)
0164 throw Standard_TypeMismatch("No init");
0165
0166 myU = U;
0167 Vec D1c;
0168 Tool::D1(*((Curve*)myC),myU,myPc,D1c);
0169
0170 if(Precision::IsInfinite(D1c.X()) || Precision::IsInfinite(D1c.Y()))
0171 {
0172 F = Precision::Infinite();
0173 return Standard_False;
0174 }
0175
0176 Standard_Real Ndu = D1c.Magnitude();
0177
0178 if(myMaxDerivOrder != 0)
0179 {
0180 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
0181 {
0182 const Standard_Real DivisionFactor = 1.e-3;
0183 Standard_Real du;
0184 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
0185 du = 0.0;
0186 else
0187 du = myUsupremum-myUinfium;
0188
0189 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
0190 //Derivative is approximated by Taylor-series
0191
0192 Standard_Integer n = 1; //Derivative order
0193 Vec V;
0194 Standard_Boolean IsDeriveFound;
0195
0196 do
0197 {
0198 V = Tool::DN(*((Curve*)myC),myU,++n);
0199 Ndu = V.Magnitude();
0200 IsDeriveFound = (Ndu > myTol);
0201 }
0202 while(!IsDeriveFound && n < myMaxDerivOrder);
0203
0204 if(IsDeriveFound)
0205 {
0206 Standard_Real u;
0207
0208 if(myU-myUinfium < aDelta)
0209 u = myU+aDelta;
0210 else
0211 u = myU-aDelta;
0212
0213 Pnt P1, P2;
0214 Tool::D0(*((Curve*)myC),Min(myU, u),P1);
0215 Tool::D0(*((Curve*)myC),Max(myU, u),P2);
0216
0217 Vec V1(P1,P2);
0218 Standard_Real aDirFactor = V.Dot(V1);
0219
0220 if(aDirFactor < 0.0)
0221 D1c = -V;
0222 else
0223 D1c = V;
0224 }//if(IsDeriveFound)
0225 else
0226 {
0227 //Derivative is approximated by three points
0228
0229 Pnt Ptemp; //(0,0,0)-coordinate
0230 Pnt P1, P2, P3;
0231 Standard_Boolean IsParameterGrown;
0232
0233 if(myU-myUinfium < 2*aDelta)
0234 {
0235 Tool::D0(*((Curve*)myC),myU,P1);
0236 Tool::D0(*((Curve*)myC),myU+aDelta,P2);
0237 Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
0238 IsParameterGrown = Standard_True;
0239 }
0240 else
0241 {
0242 Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
0243 Tool::D0(*((Curve*)myC),myU-aDelta,P2);
0244 Tool::D0(*((Curve*)myC),myU,P3);
0245 IsParameterGrown = Standard_False;
0246 }
0247
0248 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
0249
0250 if(IsParameterGrown)
0251 D1c=-3*V1+4*V2-V3;
0252 else
0253 D1c=V1-4*V2+3*V3;
0254 }
0255 Ndu = D1c.Magnitude();
0256 }//(if (Ndu <= myTol)) condition
0257 }//if(myMaxDerivOrder != 0)
0258
0259 if (Ndu <= MinTol)
0260 {
0261 //Warning: 1st derivative is equal to zero!
0262 return Standard_False;
0263 }
0264
0265 Vec PPc (myP,myPc);
0266 F = PPc.Dot(D1c)/Ndu;
0267 return Standard_True;
0268 }
0269
0270 //=============================================================================
0271
0272 Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_Real& D1f)
0273 {
0274 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
0275 Standard_Real F;
0276 return Values(U,F,D1f); /* on fait appel a Values pour simplifier la
0277 sauvegarde de l'etat. */
0278 }
0279 //=============================================================================
0280
0281 Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U,
0282 Standard_Real& F,
0283 Standard_Real& D1f)
0284 {
0285 if (!myPinit || !myCinit)
0286 throw Standard_TypeMismatch("No init");
0287
0288 Pnt myPc_old = myPc, myP_old = myP;
0289
0290 if(Value(U,F) == Standard_False)
0291 {
0292 //Warning: No function value found!;
0293
0294 myD1Init = Standard_False;
0295 return Standard_False;
0296 }
0297
0298 myU = U;
0299 myPc = myPc_old;
0300 myP = myP_old;
0301
0302 Vec D1c,D2c;
0303 Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
0304
0305 Standard_Real Ndu = D1c.Magnitude();
0306 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
0307 {
0308 //Derivative is approximated by three points
0309
0310 //Attention: aDelta value must be greater than same value for "Value(...)"
0311 // function to avoid of points' collisions.
0312 const Standard_Real DivisionFactor = 0.01;
0313 Standard_Real du;
0314 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
0315 du = 0.0;
0316 else
0317 du = myUsupremum-myUinfium;
0318
0319 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
0320
0321 Standard_Real F1, F2, F3;
0322
0323 if(myU-myUinfium < 2*aDelta)
0324 {
0325 F1=F;
0326 //const Standard_Real U1 = myU;
0327 const Standard_Real U2 = myU + aDelta;
0328 const Standard_Real U3 = myU + aDelta * 2.0;
0329
0330 if(!((Value(U2,F2)) && (Value(U3,F3))))
0331 {
0332 //There are many points close to singularity points and
0333 //which have zero-derivative. Try to decrease aDelta variable's value.
0334
0335 myD1Init = Standard_False;
0336 return Standard_False;
0337 }
0338
0339 //After calling of Value(...) function variable myU will be redeterminated.
0340 //So we must return it previous value.
0341 D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
0342 }
0343 else
0344 {
0345 F3 = F;
0346 const Standard_Real U1 = myU - aDelta * 2.0;
0347 const Standard_Real U2 = myU - aDelta;
0348 //const Standard_Real U3 = myU;
0349
0350 if(!((Value(U2,F2)) && (Value(U1,F1))))
0351 {
0352 //There are many points close to singularity points and
0353 //which have zero-derivative. Try to decrease aDelta variable's value.
0354 myD1Init = Standard_False;
0355 return Standard_False;
0356 }
0357 //After calling of Value(...) function variable myU will be redeterminated.
0358 //So we must return it previous value.
0359 D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
0360 }
0361 myU = U;
0362 myPc = myPc_old;
0363 myP = myP_old;
0364 }
0365 else
0366 {
0367 Vec PPc (myP,myPc);
0368 D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
0369 }
0370
0371 myD1f = D1f;
0372
0373 myD1Init = Standard_True;
0374 return Standard_True;
0375 }
0376 //=============================================================================
0377
0378 Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
0379 {
0380 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
0381 mySqDist.Append(myPc.SquareDistance(myP));
0382
0383 // It is necessary to always compute myD1f.
0384 myD1Init = Standard_True;
0385 Standard_Real FF, DD;
0386 Values(myU, FF, DD);
0387
0388 Standard_Integer IntVal = 0;
0389 if (myD1f > 0.0)
0390 {
0391 IntVal = 1;
0392 }
0393
0394 myIsMin.Append(IntVal);
0395 myPoint.Append(POnC(myU,myPc));
0396 return 0;
0397 }
0398 //=============================================================================
0399
0400 Standard_Integer Extrema_FuncExtPC::NbExt () const { return mySqDist.Length(); }
0401 //=============================================================================
0402
0403 Standard_Real Extrema_FuncExtPC::SquareDistance (const Standard_Integer N) const
0404 {
0405 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
0406 return mySqDist.Value(N);
0407 }
0408 //=============================================================================
0409 Standard_Boolean Extrema_FuncExtPC::IsMin (const Standard_Integer N) const
0410 {
0411 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
0412 return (myIsMin.Value(N) == 1);
0413 }
0414 //=============================================================================
0415 const POnC & Extrema_FuncExtPC::Point (const Standard_Integer N) const
0416 {
0417 if (!myPinit || !myCinit) throw Standard_TypeMismatch();
0418 return myPoint.Value(N);
0419 }
0420 //=============================================================================
0421
0422 void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
0423 {
0424 myUinfium = theUfirst;
0425 myUsupremum = theUlast;
0426 }
