e8e186ddbe49ac89c2550c47ad07b1e32b550ab0
[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& 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 }