src/GeomLib/GeomLib_CheckCurveOnSurface.cxx

   1 // Created by: Nikolai BUKHALOV
   2 // Copyright (c) 2015 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   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.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <Adaptor2d_HCurve2d.hxx>
  16 #include <Adaptor3d_Curve.hxx>
  17 #include <Adaptor3d_CurveOnSurface.hxx>
  18 #include <Adaptor3d_HSurface.hxx>
  19 #include <Geom_BSplineCurve.hxx>
  20 #include <Geom_TrimmedCurve.hxx>
  21 #include <Geom2d_BSplineCurve.hxx>
  22 #include <Geom2d_TrimmedCurve.hxx>
  23 #include <Geom2dAdaptor_GHCurve.hxx>
  24 #include <GeomAdaptor_Curve.hxx>
  25 #include <GeomAdaptor_HSurface.hxx>
  26 #include <GeomLib_CheckCurveOnSurface.hxx>
  27 #include <gp_Pnt.hxx>
  28 #include <math_Matrix.hxx>
  29 #include <math_MultipleVarFunctionWithHessian.hxx>
  30 #include <math_NewtonMinimum.hxx>
  31 #include <math_PSO.hxx>
  32 #include <math_PSOParticlesPool.hxx>
  33 #include <OSD_Parallel.hxx>
  34 #include <Standard_ErrorHandler.hxx>
  35 #include <TColStd_Array1OfReal.hxx>
  36 #include <TColStd_HArray1OfReal.hxx>
  37
  38 class GeomLib_CheckCurveOnSurface_TargetFunc;
  39
  40 static
  41 Standard_Boolean MinComputing(
  42                 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
  43                 const Standard_Real theEpsilon, //1.0e-3
  44                 const Standard_Integer theNbParticles,
  45                 Standard_Real& theBestValue,
  46                 Standard_Real& theBestParameter);
  47
  48 static Standard_Integer FillSubIntervals( const Handle(Geom_Curve)& theCurve3d,
  49                                           const Handle(Geom2d_Curve)& theCurve2d,
  50                                           const Standard_Real theFirst,
  51                                           const Standard_Real theLast,
  52                                           Standard_Integer &theNbParticles,
  53                                           TColStd_Array1OfReal* const theSubIntervals = 0);
  54
  55 //=======================================================================
  56 //class   : GeomLib_CheckCurveOnSurface_TargetFunc
  57 //purpose : Target function (to be minimized)
  58 //=======================================================================
  59 class GeomLib_CheckCurveOnSurface_TargetFunc :
  60   public math_MultipleVarFunctionWithHessian
  61 {
  62  public:
  63   GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
  64                                           const Adaptor3d_Curve& theAdCS,
  65                                           const Standard_Real theFirst,
  66                                           const Standard_Real theLast):
  67   myCurve1(theC3D),
  68   myCurve2(theAdCS),
  69   myFirst(theFirst),
  70   myLast(theLast)
  71   {
  72   }
  73
  74   //returns the number of parameters of the function
  75   //(the function is one-dimension).
  76   virtual Standard_Integer NbVariables() const {
  77     return 1;
  78   }
  79
  80   //returns value of the function when parameters are equal to theX
  81   virtual Standard_Boolean Value(const math_Vector& theX,
  82                                  Standard_Real& theFVal)
  83   {
  84     return Value(theX(1), theFVal);
  85   }
  86
  87   //returns value of the one-dimension-function when parameter
  88   //is equal to theX
  89   Standard_Boolean Value( const  Standard_Real theX,
  90                           Standard_Real& theFVal) const
  91   {
  92     try
  93     {
  94       OCC_CATCH_SIGNALS
  95       if (!CheckParameter(theX))
  96         return Standard_False;
  97
  98       const gp_Pnt  aP1(myCurve1.Value(theX)),
  99                     aP2(myCurve2.Value(theX));
 100
 101       theFVal = -1.0*aP1.SquareDistance(aP2);
 102     }
 103     catch(Standard_Failure const&) {
 104       return Standard_False;
 105     }
 106     //
 107     return Standard_True;
 108   }
 109
 110   //see analogical method for abstract owner class math_MultipleVarFunction
 111   virtual Standard_Integer GetStateNumber()
 112   {
 113     return 0;
 114   }
 115
 116   //returns the gradient of the function when parameters are
 117   //equal to theX
 118   virtual Standard_Boolean Gradient(const math_Vector& theX,
 119                                     math_Vector& theGrad)
 120   {
 121     return Derive(theX(1), theGrad(1));
 122   }
 123
 124   //returns 1st derivative of the the one-dimension-function when
 125   //parameter is equal to theX
 126   Standard_Boolean Derive(const Standard_Real theX, Standard_Real& theDeriv1, Standard_Real* const theDeriv2 = 0) const
 127   {
 128     try
 129     {
 130       OCC_CATCH_SIGNALS
 131       if (!CheckParameter(theX))
 132       {
 133         return Standard_False;
 134       }
 135       //
 136       gp_Pnt aP1, aP2;
 137       gp_Vec aDC1, aDC2, aDCC1, aDCC2;
 138       //
 139       if (!theDeriv2)
 140       {
 141         myCurve1.D1(theX, aP1, aDC1);
 142         myCurve2.D1(theX, aP2, aDC2);
 143       }
 144       else
 145       {
 146         myCurve1.D2(theX, aP1, aDC1, aDCC1);
 147         myCurve2.D2(theX, aP2, aDC2, aDCC2);
 148       }
 149
 150       const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
 151       //
 152       theDeriv1 = -2.0*aVec1.Dot(aVec2);
 153
 154       if (theDeriv2)
 155       {
 156         const gp_Vec aVec3(aDCC2 - aDCC1);
 157         *theDeriv2 = -2.0*(aVec2.SquareMagnitude() + aVec1.Dot(aVec3));
 158       }
 159     }
 160     catch(Standard_Failure const&)
 161     {
 162       return Standard_False;
 163     }
 164
 165     return Standard_True;
 166   }
 167
 168   //returns value and gradient
 169   virtual Standard_Boolean Values(const math_Vector& theX,
 170                                   Standard_Real& theVal,
 171                                   math_Vector& theGrad)
 172   {
 173     if (!Value(theX, theVal))
 174     {
 175       return Standard_False;
 176     }
 177     //
 178     if (!Gradient(theX, theGrad)) {
 179       return Standard_False;
 180     }
 181     //
 182     return Standard_True;
 183   }
 184
 185   //returns value, gradient and hessian
 186   virtual Standard_Boolean Values(const math_Vector& theX,
 187                                   Standard_Real& theVal,
 188                                   math_Vector& theGrad,
 189                                   math_Matrix& theHessian)
 190   {
 191     if (!Value(theX, theVal))
 192     {
 193       return Standard_False;
 194     }
 195     //
 196     if (!Derive(theX(1), theGrad(1), &theHessian(1, 1)))
 197     {
 198       return Standard_False;
 199     }
 200     //
 201     return Standard_True;
 202   }
 203   //
 204   Standard_Real FirstParameter() const
 205   {
 206     return myFirst;
 207   }
 208
 209   //
 210   Standard_Real LastParameter() const
 211   {
 212     return myLast;
 213   }
 214
 215  private:
 216   GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
 217
 218   //checks if the function can be computed when its parameter is
 219   //equal to theParam
 220    Standard_Boolean CheckParameter(const Standard_Real theParam) const
 221    {
 222      return ((myFirst <= theParam) && (theParam <= myLast));
 223    }
 224
 225    const Adaptor3d_Curve& myCurve1;
 226    const Adaptor3d_Curve& myCurve2;
 227    const Standard_Real myFirst;
 228    const Standard_Real myLast;
 229 };
 230
 231 //=======================================================================
 232 //class   : GeomLib_CheckCurveOnSurface_Local
 233 //purpose : Created for parallelization possibility only
 234 //=======================================================================
 235 class GeomLib_CheckCurveOnSurface_Local
 236 {
 237 public:
 238   GeomLib_CheckCurveOnSurface_Local(
 239               const Handle(Geom_Curve)& theCurve3D,
 240               const Handle(Geom2d_Curve)& theCurve2D,
 241               const Handle(Geom_Surface)& theSurface,
 242               const TColStd_Array1OfReal& theIntervalsArr,
 243               const Standard_Real theEpsilonRange,
 244               const Standard_Integer theNbParticles):
 245   myCurve3D(theCurve3D),
 246   myCurve2D(theCurve2D),
 247   mySurface(theSurface),
 248   mySubIntervals(theIntervalsArr),
 249   myEpsilonRange(theEpsilonRange),
 250   myNbParticles(theNbParticles),
 251   myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1),
 252   myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1)
 253   {
 254   }
 255
 256   void operator()(const Standard_Integer& theIndex) const
 257   {
 258     //For every sub-interval (which is set by mySubIntervals array) this method
 259     //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function.
 260     //This optimal value will be put in corresponding (depending on theIndex - the
 261     //identificator of the current interval in mySubIntervals array) cell of
 262     //myArrOfDist and myArrOfParam arrays.
 263     const GeomAdaptor_Curve anAC(myCurve3D);
 264     const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D);
 265     const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface);
 266
 267     const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
 268
 269     GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
 270                                                   mySubIntervals.Value(theIndex),
 271                                                   mySubIntervals.Value(theIndex+1));
 272
 273     Standard_Real aMinDist = RealLast(), aPar = 0.0;
 274     if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
 275     {
 276       myArrOfDist(theIndex) = RealLast();
 277       myArrOfParam(theIndex) = aFunc.FirstParameter();
 278       return;
 279     }
 280
 281     myArrOfDist(theIndex) = aMinDist;
 282     myArrOfParam(theIndex) = aPar;
 283   }
 284
 285   //Returns optimal value (inverse of square of maximal distance)
 286   void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
 287   {
 288     //This method looks for the minimal value of myArrOfDist.
 289
 290     const Standard_Integer aStartInd = myArrOfDist.Lower();
 291     theMinimalValue = myArrOfDist(aStartInd);
 292     theParameter = myArrOfParam(aStartInd);
 293     for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++)
 294     {
 295       if(myArrOfDist(i) < theMinimalValue)
 296       {
 297         theMinimalValue = myArrOfDist(i);
 298         theParameter = myArrOfParam(i);
 299       }
 300     }
 301   }
 302
 303 private:
 304   GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&);
 305   const Handle(Geom_Curve)& myCurve3D;
 306   const Handle(Geom2d_Curve)& myCurve2D;
 307   const Handle(Geom_Surface)& mySurface;
 308
 309   const TColStd_Array1OfReal& mySubIntervals;
 310   const Standard_Real myEpsilonRange;
 311   const Standard_Integer myNbParticles;
 312   mutable NCollection_Array1<Standard_Real> myArrOfDist;
 313   mutable NCollection_Array1<Standard_Real> myArrOfParam;
 314 };
 315
 316 //=======================================================================
 317 //function : GeomLib_CheckCurveOnSurface
 318 //purpose  :
 319 //=======================================================================
 320 GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
 321 :
 322   myFirst(0.),
 323   myLast(0.),
 324   myErrorStatus(0),
 325   myMaxDistance(RealLast()),
 326   myMaxParameter(0.),
 327   myTolRange(Precision::PConfusion())
 328 {
 329 }
 330
 331 //=======================================================================
 332 //function : GeomLib_CheckCurveOnSurface
 333 //purpose  :
 334 //=======================================================================
 335 GeomLib_CheckCurveOnSurface::
 336   GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve,
 337                               const Handle(Geom_Surface)& theSurface,
 338                               const Standard_Real theFirst,
 339                               const Standard_Real theLast,
 340                               const Standard_Real theTolRange):
 341   myCurve(theCurve),
 342   mySurface(theSurface),
 343   myFirst(theFirst),
 344   myLast(theLast),
 345   myErrorStatus(0),
 346   myMaxDistance(RealLast()),
 347   myMaxParameter(0.),
 348   myTolRange(theTolRange)
 349 {
 350 }
 351
 352 //=======================================================================
 353 //function : Init
 354 //purpose  :
 355 //=======================================================================
 356 void GeomLib_CheckCurveOnSurface::Init()
 357 {
 358   myCurve.Nullify();
 359   mySurface.Nullify();
 360   myFirst = 0.0;
 361   myLast = 0.0;
 362   myErrorStatus = 0;
 363   myMaxDistance = RealLast();
 364   myMaxParameter = 0.0;
 365   myTolRange = Precision::PConfusion();
 366 }
 367
 368 //=======================================================================
 369 //function : Init
 370 //purpose  :
 371 //=======================================================================
 372 void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve,
 373                                         const Handle(Geom_Surface)& theSurface,
 374                                         const Standard_Real theFirst,
 375                                         const Standard_Real theLast,
 376                                         const Standard_Real theTolRange)
 377 {
 378   myCurve = theCurve;
 379   mySurface = theSurface;
 380   myFirst = theFirst;
 381   myLast = theLast;
 382   myErrorStatus = 0;
 383   myMaxDistance = RealLast();
 384   myMaxParameter = 0.0;
 385   myTolRange = theTolRange;
 386 }
 387
 388 //=======================================================================
 389 //function : Perform
 390 //purpose  :
 391 //=======================================================================
 392 #ifndef HAVE_TBB
 393 //After fixing bug # 26365, this fragment should be deleted
 394 //(together the text "#ifdef HAVE_TBB")
 395
 396 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
 397                                           const Standard_Boolean)
 398 {
 399   const Standard_Boolean isTheMTDisabled = Standard_True;
 400 #else
 401 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
 402                                           const Standard_Boolean isTheMTDisabled)
 403 {
 404 #endif
 405   if( myCurve.IsNull() ||
 406       mySurface.IsNull() ||
 407       thePCurve.IsNull())
 408   {
 409     myErrorStatus = 1;
 410     return;
 411   }
 412
 413   if(((myCurve->FirstParameter() - myFirst) > myTolRange) ||
 414      ((myCurve->LastParameter() - myLast) < -myTolRange) ||
 415      ((thePCurve->FirstParameter() - myFirst) > myTolRange) ||
 416      ((thePCurve->LastParameter() - myLast) < -myTolRange))
 417   {
 418     myErrorStatus = 2;
 419     return;
 420   }
 421
 422   const Standard_Real anEpsilonRange = 1.e-3;
 423
 424   Standard_Integer aNbParticles = 3;
 425
 426   //Polynomial function with degree n has not more than n-1 maxima and
 427   //minima (degree of 1st derivative is equal to n-1 => 1st derivative has
 428   //no greater than n-1 roots). Consequently, this function has
 429   //maximum n monotonicity intervals. That is a good idea to try to put
 430   //at least one particle in every monotonicity interval. Therefore,
 431   //number of particles should be equal to n.
 432
 433   const Standard_Integer aNbSubIntervals =
 434                               FillSubIntervals( myCurve, thePCurve,
 435                                                 myFirst, myLast, aNbParticles);
 436
 437   if(!aNbSubIntervals)
 438   {
 439     myErrorStatus = 3;
 440     return;
 441   }
 442
 443   try {
 444     OCC_CATCH_SIGNALS
 445
 446     TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
 447     FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
 448
 449     GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
 450                                 mySurface, anIntervals, anEpsilonRange, aNbParticles);
 451
 452     OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
 453
 454     aComp.OptimalValues(myMaxDistance, myMaxParameter);
 455
 456     myMaxDistance = sqrt(Abs(myMaxDistance));
 457   }
 458   catch (Standard_Failure const&) {
 459     myErrorStatus = 3;
 460   }
 461 }
 462
 463 //=======================================================================
 464 // Function : FillSubIntervals
 465 // purpose : Divides [theFirst, theLast] interval on parts
 466 //            in order to make searching-algorithm more precisely
 467 //            (fills theSubIntervals array).
 468 //            Returns number of subintervals.
 469 //=======================================================================
 470 Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d,
 471                                   const Handle(Geom2d_Curve)& theCurve2d,
 472                                   const Standard_Real theFirst,
 473                                   const Standard_Real theLast,
 474                                   Standard_Integer &theNbParticles,
 475                                   TColStd_Array1OfReal* const theSubIntervals)
 476 {
 477   const Standard_Integer aMaxKnots = 101;
 478   const Standard_Real anArrTempC[2] = {theFirst, theLast};
 479   const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
 480
 481   theNbParticles = 3;
 482   Handle(Geom2d_BSplineCurve) aBS2DCurv;
 483   Handle(Geom_BSplineCurve) aBS3DCurv;
 484   Standard_Boolean isTrimmed3D = Standard_False, isTrimmed2D = Standard_False;
 485
 486   //
 487   if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
 488   {
 489     aBS3DCurv = Handle(Geom_BSplineCurve)::
 490                       DownCast(Handle(Geom_TrimmedCurve)::
 491                       DownCast(theCurve3d)->BasisCurve());
 492     isTrimmed3D = Standard_True;
 493   }
 494   else
 495   {
 496     aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
 497   }
 498
 499   if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
 500   {
 501     aBS2DCurv = Handle(Geom2d_BSplineCurve)::
 502                       DownCast(Handle(Geom2d_TrimmedCurve)::
 503                       DownCast(theCurve2d)->BasisCurve());
 504     isTrimmed2D = Standard_True;
 505   }
 506   else
 507   {
 508     aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
 509   }
 510
 511   Handle(TColStd_HArray1OfReal) anArrKnots3D,  anArrKnots2D;
 512
 513   if(!aBS3DCurv.IsNull())
 514   {
 515     if(aBS3DCurv->NbKnots() <= aMaxKnots)
 516     {
 517       anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
 518     }
 519     else
 520     {
 521       Standard_Integer KnotCount;
 522       if(isTrimmed3D)
 523       {
 524         Standard_Integer i;
 525         KnotCount = 0;
 526         const TColStd_Array1OfReal& aKnots = aBS3DCurv->Knots();
 527         for(i = aBS3DCurv->FirstUKnotIndex(); i <= aBS3DCurv->LastUKnotIndex(); ++i)
 528         {
 529           if(aKnots(i) > theFirst && aKnots(i) < theLast)
 530           {
 531             ++KnotCount;
 532           }
 533         }
 534         KnotCount += 2;
 535       }
 536       else
 537       {
 538         KnotCount = aBS3DCurv->LastUKnotIndex() - aBS3DCurv->FirstUKnotIndex() + 1;
 539       }
 540       if(KnotCount <= aMaxKnots)
 541       {
 542         anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
 543       }
 544       else
 545       {
 546         anArrKnots3D = new TColStd_HArray1OfReal(1, aMaxKnots);
 547         anArrKnots3D->SetValue(1, theFirst);
 548         anArrKnots3D->SetValue(aMaxKnots, theLast);
 549         Standard_Integer i;
 550         Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
 551         Standard_Real t = theFirst + dt;
 552         for(i = 2; i < aMaxKnots; ++i, t += dt)
 553         {
 554           anArrKnots3D->SetValue(i, t);
 555         }
 556       }
 557     }
 558   }
 559   else
 560   {
 561     anArrKnots3D = new TColStd_HArray1OfReal(anArrTemp);
 562   }
 563   if(!aBS2DCurv.IsNull())
 564   {
 565     if(aBS2DCurv->NbKnots() <= aMaxKnots)
 566     {
 567       anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
 568     }
 569     else
 570     {
 571       Standard_Integer KnotCount;
 572       if(isTrimmed2D)
 573       {
 574         Standard_Integer i;
 575         KnotCount = 0;
 576         const TColStd_Array1OfReal& aKnots = aBS2DCurv->Knots();
 577         for(i = aBS2DCurv->FirstUKnotIndex(); i <= aBS2DCurv->LastUKnotIndex(); ++i)
 578         {
 579           if(aKnots(i) > theFirst && aKnots(i) < theLast)
 580           {
 581             ++KnotCount;
 582           }
 583         }
 584         KnotCount += 2;
 585       }
 586       else
 587       {
 588         KnotCount = aBS2DCurv->LastUKnotIndex() - aBS2DCurv->FirstUKnotIndex() + 1;
 589       }
 590       if(KnotCount <= aMaxKnots)
 591       {
 592         anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
 593       }
 594       else
 595       {
 596         anArrKnots2D = new TColStd_HArray1OfReal(1, aMaxKnots);
 597         anArrKnots2D->SetValue(1, theFirst);
 598         anArrKnots2D->SetValue(aMaxKnots, theLast);
 599         Standard_Integer i;
 600         Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
 601         Standard_Real t = theFirst + dt;
 602         for(i = 2; i < aMaxKnots; ++i, t += dt)
 603         {
 604           anArrKnots2D->SetValue(i, t);
 605         }
 606       }
 607     }
 608   }
 609   else
 610   {
 611     anArrKnots2D = new TColStd_HArray1OfReal(anArrTemp);
 612   }
 613
 614
 615   Standard_Integer aNbSubIntervals = 1;
 616
 617   try
 618   {
 619     OCC_CATCH_SIGNALS
 620     const Standard_Integer  anIndMax3D = anArrKnots3D->Upper(),
 621                             anIndMax2D = anArrKnots2D->Upper();
 622
 623     Standard_Integer  anIndex3D = anArrKnots3D->Lower(),
 624                       anIndex2D = anArrKnots2D->Lower();
 625
 626     if(theSubIntervals)
 627       theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
 628
 629     while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
 630     {
 631       const Standard_Real aVal3D = anArrKnots3D->Value(anIndex3D),
 632                           aVal2D = anArrKnots2D->Value(anIndex2D);
 633       const Standard_Real aDelta = aVal3D - aVal2D;
 634
 635       if(aDelta < Precision::PConfusion())
 636       {//aVal3D <= aVal2D
 637         if((aVal3D > theFirst) && (aVal3D < theLast))
 638         {
 639           aNbSubIntervals++;
 640
 641           if(theSubIntervals)
 642             theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
 643         }
 644
 645         anIndex3D++;
 646
 647         if(-aDelta < Precision::PConfusion())
 648         {//aVal3D == aVal2D
 649           anIndex2D++;
 650         }
 651       }
 652       else
 653       {//aVal2D < aVal3D
 654         if((aVal2D > theFirst) && (aVal2D < theLast))
 655         {
 656           aNbSubIntervals++;
 657
 658           if(theSubIntervals)
 659             theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
 660         }
 661
 662         anIndex2D++;
 663       }
 664     }
 665
 666     if(theSubIntervals)
 667       theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
 668
 669     if(!aBS3DCurv.IsNull())
 670     {
 671       theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
 672     }
 673
 674     if(!aBS2DCurv.IsNull())
 675     {
 676       theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
 677     }
 678   }
 679   catch(Standard_Failure const&)
 680   {
 681 #ifdef OCCT_DEBUG
 682     std::cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
 683             "FillSubIntervals(): Incorrect filling!" << std::endl;
 684 #endif
 685
 686     aNbSubIntervals = 0;
 687   }
 688
 689   return aNbSubIntervals;
 690 }
 691
 692 //=======================================================================
 693 //class   : PSO_Perform
 694 //purpose : Searches minimal distance with math_PSO class
 695 //=======================================================================
 696 Standard_Boolean PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
 697                              const math_Vector &theParInf,
 698                              const math_Vector &theParSup,
 699                              const Standard_Real theEpsilon,
 700                              const Standard_Integer theNbParticles,
 701                              Standard_Real& theBestValue,
 702                              math_Vector &theOutputParam)
 703 {
 704   const Standard_Real aDeltaParam = theParSup(1) - theParInf(1);
 705   if(aDeltaParam < Precision::PConfusion())
 706     return Standard_False;
 707
 708   math_Vector aStepPar(1, 1);
 709   aStepPar(1) = theEpsilon*aDeltaParam;
 710
 711   math_PSOParticlesPool aParticles(theNbParticles, 1);
 712
 713   //They are used for finding a position of theNbParticles worst places
 714   const Standard_Integer aNbControlPoints = 3*theNbParticles;
 715
 716   const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1);
 717   Standard_Integer aCount = 1;
 718   for(Standard_Real aPrm = theParInf(1); aCount <= aNbControlPoints; aCount++,
 719     aPrm = (aCount == aNbControlPoints)? theParSup(1) : aPrm+aStep)
 720   {
 721     Standard_Real aVal = RealLast();
 722     if(!theFunction.Value(aPrm, aVal))
 723       continue;
 724
 725     PSO_Particle* aParticle = aParticles.GetWorstParticle();
 726
 727     if(aVal > aParticle->BestDistance)
 728       continue;
 729
 730     aParticle->Position[0] = aPrm;
 731     aParticle->BestPosition[0] = aPrm;
 732     aParticle->Distance     = aVal;
 733     aParticle->BestDistance = aVal;
 734   }
 735
 736   math_PSO aPSO(&theFunction, theParInf, theParSup, aStepPar);
 737   aPSO.Perform(aParticles, theNbParticles, theBestValue, theOutputParam);
 738
 739   return Standard_True;
 740 }
 741
 742 //=======================================================================
 743 //class   : MinComputing
 744 //purpose : Performs computing minimal value
 745 //=======================================================================
 746 Standard_Boolean MinComputing (
 747                 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
 748                 const Standard_Real theEpsilon, //1.0e-3
 749                 const Standard_Integer theNbParticles,
 750                 Standard_Real& theBestValue,
 751                 Standard_Real& theBestParameter)
 752 {
 753   try
 754   {
 755     OCC_CATCH_SIGNALS
 756
 757     //
 758     math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1);
 759     aParInf(1) = theFunction.FirstParameter();
 760     aParSup(1) = theFunction.LastParameter();
 761     theBestParameter = aParInf(1);
 762     theBestValue = RealLast();
 763
 764     if(!PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
 765                     theBestValue, anOutputParam))
 766     {
 767 #ifdef OCCT_DEBUG
 768       std::cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << std::endl;
 769 #endif
 770       return Standard_False;
 771     }
 772
 773     theBestParameter = anOutputParam(1);
 774
 775     //Here, anOutputParam contains parameter, which is near to optimal.
 776     //It needs to be more precise. Precision is made by math_NewtonMinimum.
 777     math_NewtonMinimum aMinSol(theFunction);
 778     aMinSol.Perform(theFunction, anOutputParam);
 779
 780     if(aMinSol.IsDone() && (aMinSol.GetStatus() == math_OK))
 781     {//math_NewtonMinimum has precised the value. We take it.
 782       aMinSol.Location(anOutputParam);
 783       theBestParameter =  anOutputParam(1);
 784       theBestValue = aMinSol.Minimum();
 785     }
 786     else
 787     {//Use math_PSO again but on smaller range.
 788       const Standard_Real aStep = theEpsilon*(aParSup(1) - aParInf(1));
 789       aParInf(1) = theBestParameter - 0.5*aStep;
 790       aParSup(1) = theBestParameter + 0.5*aStep;
 791
 792       Standard_Real aValue = RealLast();
 793       if(PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
 794                      aValue, anOutputParam))
 795       {
 796         if(aValue < theBestValue)
 797         {
 798           theBestValue = aValue;
 799           theBestParameter = anOutputParam(1);
 800         }
 801       }
 802     }
 803   }
 804   catch(Standard_Failure const&)
 805   {
 806 #ifdef OCCT_DEBUG
 807     std::cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << std::endl;
 808 #endif
 809     return Standard_False;
 810   }
 811
 812   return Standard_True;
 813 }