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