src/Extrema/Extrema_FuncExtPC.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   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
   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.
   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 //
  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.
  18
  19 #include <Standard_TypeMismatch.hxx>
  20 #include <Precision.hxx>
  21 static const Standard_Real TolFactor = 1.e-12;
  22 static const Standard_Real MinTol    = 1.e-20;
  23 static const Standard_Real MinStep   = 1e-7;
  24
  25 static const Standard_Integer MaxOrder = 3;
  26
  27
  28
  29 /*-----------------------------------------------------------------------------
  30 Fonction permettant de rechercher une distance extremale entre un point P et
  31 une courbe C (en partant d'un point approche C(u0)).
  32 Cette classe herite de math_FunctionWithDerivative et est utilisee par
  33 les algorithmes math_FunctionRoot et math_FunctionRoots.
  34 Si 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
  37
  38
  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) }
  42 ----------------------------------------------------------------------------*/
  43
  44 Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
  45 {
  46   const Standard_Integer NPoint = 10;
  47   const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
  48
  49   Standard_Integer aNum = 0;
  50   Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
  51   //(it is computed with using NPoint point)
  52
  53   do
  54   {
  55     Standard_Real u = myUinfium + aNum*aStep;    //parameter for every point
  56     if(u > myUsupremum)
  57       u = myUsupremum;
  58
  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;
  65   }
  66   while(++aNum < NPoint+1);
  67
  68   return Max(aMax*TolFactor,MinTol);
  69
  70 }
  71
  72 //=============================================================================
  73
  74 Extrema_FuncExtPC::Extrema_FuncExtPC():
  75 myU(0.),
  76 myD1f(0.)
  77 {
  78   myPinit = Standard_False;
  79   myCinit = Standard_False;
  80   myD1Init = Standard_False;
  81
  82   SubIntervalInitialize(0.0,0.0);
  83   myMaxDerivOrder = 0;
  84   myTol=MinTol;
  85 }
  86
  87 //=============================================================================
  88
  89 Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
  90                                       const Curve& C): myU(0.), myD1f(0.)
  91 {
  92   myP = P;
  93   myC = (Standard_Address)&C;
  94   myPinit = Standard_True;
  95   myCinit = Standard_True;
  96   myD1Init = Standard_False;
  97
  98   SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
  99     Tool::LastParameter(*((Curve*)myC)));
 100
 101   switch(Tool::GetType(*((Curve*)myC)))
 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;
 113   }
 114 }
 115 //=============================================================================
 116
 117 void Extrema_FuncExtPC::Initialize(const Curve& C)
 118 {
 119   myC = (Standard_Address)&C;
 120   myCinit = Standard_True;
 121   myPoint.Clear();
 122   mySqDist.Clear();
 123   myIsMin.Clear();
 124
 125   SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
 126     Tool::LastParameter(*((Curve*)myC)));
 127
 128   switch(Tool::GetType(*((Curve*)myC)))
 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;
 140   }
 141 }
 142
 143 //=============================================================================
 144
 145 void 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
 156
 157 Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
 158 {
 159   if (!myPinit || !myCinit)
 160     Standard_TypeMismatch::Raise("No init");
 161
 162   myU = U;
 163   Vec D1c;
 164   Tool::D1(*((Curve*)myC),myU,myPc,D1c);
 165   Standard_Real Ndu = D1c.Magnitude();
 166
 167   if(myMaxDerivOrder != 0)
 168   {
 169     if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
 170     {
 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;
 177
 178       const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 179       //Derivative is approximated by Taylor-series
 180
 181       Standard_Integer n = 1; //Derivative order
 182       Vec V;
 183       Standard_Boolean IsDeriveFound;
 184
 185       do
 186       {
 187         V = Tool::DN(*((Curve*)myC),myU,++n);
 188         Ndu = V.Magnitude();
 189         IsDeriveFound = (Ndu > myTol);
 190       }
 191       while(!IsDeriveFound && n < myMaxDerivOrder);
 192
 193       if(IsDeriveFound)
 194       {
 195         Standard_Real u;
 196
 197         if(myU-myUinfium < aDelta)
 198           u = myU+aDelta;
 199         else
 200           u = myU-aDelta;
 201
 202         Pnt P1, P2;
 203         Tool::D0(*((Curve*)myC),Min(myU, u),P1);
 204         Tool::D0(*((Curve*)myC),Max(myU, u),P2);
 205
 206         Vec V1(P1,P2);
 207         Standard_Real aDirFactor = V.Dot(V1);
 208
 209         if(aDirFactor < 0.0)
 210           D1c = -V;
 211         else
 212           D1c = V;
 213       }//if(IsDeriveFound)
 214       else
 215       {
 216         //Derivative is approximated by three points
 217
 218         Pnt Ptemp; //(0,0,0)-coordinate
 219         Pnt P1, P2, P3;
 220         Standard_Boolean IsParameterGrown;
 221
 222         if(myU-myUinfium < 2*aDelta)
 223         {
 224           Tool::D0(*((Curve*)myC),myU,P1);
 225           Tool::D0(*((Curve*)myC),myU+aDelta,P2);
 226           Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
 227           IsParameterGrown = Standard_True;
 228         }
 229         else
 230         {
 231           Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
 232           Tool::D0(*((Curve*)myC),myU-aDelta,P2);
 233           Tool::D0(*((Curve*)myC),myU,P3);
 234           IsParameterGrown = Standard_False;
 235         }
 236
 237         Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
 238
 239         if(IsParameterGrown)
 240           D1c=-3*V1+4*V2-V3;
 241         else
 242           D1c=V1-4*V2+3*V3;
 243       }
 244       Ndu = D1c.Magnitude();
 245     }//(if (Ndu <= myTol)) condition
 246   }//if(myMaxDerivOrder != 0)
 247
 248   if (Ndu <= MinTol)
 249   {
 250     //Warning: 1st derivative is equal to zero!
 251     return Standard_False;
 252   }
 253
 254   Vec PPc (myP,myPc);
 255   F = PPc.Dot(D1c)/Ndu;
 256   return Standard_True;
 257 }
 258
 259 //=============================================================================
 260
 261 Standard_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
 266                            sauvegarde de l'etat. */
 267 }
 268 //=============================================================================
 269
 270 Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U,
 271                                             Standard_Real& F,
 272                                             Standard_Real& D1f)
 273 {
 274   if (!myPinit || !myCinit)
 275     Standard_TypeMismatch::Raise("No init");
 276
 277   Pnt myPc_old = myPc, myP_old = myP;
 278
 279   if(Value(U,F) == Standard_False)
 280   {
 281     //Warning: No function value found!;
 282
 283     myD1Init = Standard_False;
 284     return Standard_False;
 285   }
 286
 287   myU = U;
 288   myPc = myPc_old;
 289   myP = myP_old;
 290
 291   Vec D1c,D2c;
 292   Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
 293
 294   Standard_Real Ndu = D1c.Magnitude();
 295   if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
 296   {
 297     //Derivative is approximated by three points
 298
 299     //Attention:  aDelta value must be greater than same value for "Value(...)"
 300     //            function to avoid of points' collisions.
 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;
 307
 308     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 309
 310     Standard_Real F1, F2, F3;
 311
 312     if(myU-myUinfium < 2*aDelta)
 313     {
 314       F1=F;
 315       //const Standard_Real U1 = myU;
 316       const Standard_Real U2 = myU + aDelta;
 317       const Standard_Real U3 = myU + aDelta * 2.0;
 318
 319       if(!((Value(U2,F2)) && (Value(U3,F3))))
 320       {
 321         //There are many points close to singularity points and
 322         //which have zero-derivative. Try to decrease aDelta variable's value.
 323
 324         myD1Init = Standard_False;
 325         return Standard_False;
 326       }
 327
 328       //After calling of Value(...) function variable myU will be redeterminated.
 329       //So we must return it previous value.
 330       D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
 331     }
 332     else
 333     {
 334       F3 = F;
 335       const Standard_Real U1 = myU - aDelta * 2.0;
 336       const Standard_Real U2 = myU - aDelta;
 337       //const Standard_Real U3 = myU;
 338
 339       if(!((Value(U2,F2)) && (Value(U1,F1))))
 340       {
 341         //There are many points close to singularity points and
 342         //which have zero-derivative. Try to decrease aDelta variable's value.
 343         myD1Init = Standard_False;
 344         return Standard_False;
 345       }
 346       //After calling of Value(...) function variable myU will be redeterminated.
 347       //So we must return it previous value.
 348       D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
 349     }
 350     myU = U;
 351     myPc = myPc_old;
 352     myP = myP_old;
 353   }
 354   else
 355   {
 356     Vec PPc (myP,myPc);
 357     D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
 358   }
 359
 360   myD1f = D1f;
 361
 362   myD1Init = Standard_True;
 363   return Standard_True;
 364 }
 365 //=============================================================================
 366
 367 Standard_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
 388 Standard_Integer Extrema_FuncExtPC::NbExt () const { return mySqDist.Length(); }
 389 //=============================================================================
 390
 391 Standard_Real Extrema_FuncExtPC::SquareDistance (const Standard_Integer N) const
 392 {
 393   if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
 394   return mySqDist.Value(N);
 395 }
 396 //=============================================================================
 397 Standard_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 //=============================================================================
 403 const POnC & Extrema_FuncExtPC::Point (const Standard_Integer N) const
 404 {
 405   if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
 406   return myPoint.Value(N);
 407 }
 408 //=============================================================================
 409
 410 void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
 411 {
 412   myUinfium = theUfirst;
 413   myUsupremum = theUlast;
 414 }