0029611: Misprint in the formula of Hessian computation in the file GeomLib_CheckCurv...
[occt.git] / 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& 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)
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) {
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)
680   {
681 #ifdef OCCT_DEBUG
682     cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
683             "FillSubIntervals(): Incorrect filling!" << 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       cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << 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)
805   {
806 #ifdef OCCT_DEBUG
807     cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
808 #endif
809     return Standard_False;
810   }
811
812   return Standard_True;
813 }