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