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