0026506: Change class BRepLib_CheckCurveOnSurface
[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
37 class GeomLib_CheckCurveOnSurface_TargetFunc;
38
39 static 
40 Standard_Boolean MinComputing(
41                 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
42                 const Standard_Real theEpsilon, //1.0e-3
43                 const Standard_Integer theNbParticles,
44                 Standard_Real& theBestValue,
45                 Standard_Real& theBestParameter);
46
47 static Standard_Integer FillSubIntervals( const Handle(Geom_Curve)& theCurve3d,
48                                           const Handle(Geom2d_Curve)& theCurve2d,
49                                           const Standard_Real theFirst,
50                                           const Standard_Real theLast,
51                                           Standard_Integer &theNbParticles,
52                                           TColStd_Array1OfReal* const theSubIntervals = 0);
53
54 //=======================================================================
55 //class   : GeomLib_CheckCurveOnSurface_TargetFunc
56 //purpose : Target function (to be minimized)
57 //=======================================================================
58 class GeomLib_CheckCurveOnSurface_TargetFunc :
59   public math_MultipleVarFunctionWithHessian
60 {
61  public:
62   GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
63                                           const Adaptor3d_Curve& theAdCS,
64                                           const Standard_Real theFirst,
65                                           const Standard_Real theLast):
66   myCurve1(theC3D),
67   myCurve2(theAdCS),
68   myFirst(theFirst),
69   myLast(theLast)
70   {
71   }
72   
73   //returns the number of parameters of the function
74   //(the function is one-dimension).
75   virtual Standard_Integer NbVariables() const {
76     return 1;
77   }
78   
79   //returns value of the function when parameters are equal to theX
80   virtual Standard_Boolean Value(const math_Vector& theX,
81                                  Standard_Real& theFVal)
82   {
83     return Value(theX(1), theFVal);
84   }
85
86   //returns value of the one-dimension-function when parameter
87   //is equal to theX
88   Standard_Boolean Value( const  Standard_Real theX,
89                           Standard_Real& theFVal) const
90   {
91     try
92     {
93       OCC_CATCH_SIGNALS
94       if (!CheckParameter(theX))
95         return Standard_False;
96
97       const gp_Pnt  aP1(myCurve1.Value(theX)),
98                     aP2(myCurve2.Value(theX));
99       
100       theFVal = -1.0*aP1.SquareDistance(aP2);
101     }
102     catch(Standard_Failure) {
103       return Standard_False;
104     }
105     //
106     return Standard_True;
107   }
108
109   //see analogical method for abstract owner class math_MultipleVarFunction
110   virtual Standard_Integer GetStateNumber()
111   {
112     return 0;
113   }
114   
115   //returns the gradient of the function when parameters are
116   //equal to theX
117   virtual Standard_Boolean Gradient(const math_Vector& theX,
118                                     math_Vector& theGrad)
119   {
120     return Derive(theX(1), theGrad(1));
121   }
122
123   //returns 1st derivative of the the one-dimension-function when
124   //parameter is equal to theX
125   Standard_Boolean Derive(const Standard_Real theX, Standard_Real& theDeriv) const
126   {
127     try
128     {
129       OCC_CATCH_SIGNALS
130       if (!CheckParameter(theX))
131       {
132         return Standard_False;
133       }
134       //
135       gp_Pnt aP1, aP2;
136       gp_Vec aDC1, aDC2;
137       //
138       myCurve1.D1(theX, aP1, aDC1);
139       myCurve2.D1(theX, aP2, aDC2);
140
141       const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
142       //
143       theDeriv = -2.0*aVec1.Dot(aVec2);
144     }
145     catch(Standard_Failure)
146     {
147       return Standard_False;
148     }
149
150     return Standard_True;
151   }
152   
153   //returns value and gradient   
154   virtual Standard_Boolean Values(const math_Vector& theX,
155                                   Standard_Real& theVal,
156                                   math_Vector& theGrad) 
157   {
158     if (!Value(theX, theVal))
159     {
160       return Standard_False;
161     }
162     //
163     if (!Gradient(theX, theGrad)) {
164       return Standard_False;
165     }
166     //
167     return Standard_True;
168   }
169
170   //returns value, gradient and hessian
171   virtual Standard_Boolean Values(const math_Vector& theX,
172                                   Standard_Real& theVal,
173                                   math_Vector& theGrad,
174                                   math_Matrix& theHessian)
175   {
176     if (!Value(theX, theVal))
177     {
178       return Standard_False;
179     }
180     //
181     if (!Gradient(theX, theGrad))
182     {
183       return Standard_False;
184     }
185     //
186     theHessian(1,1) = theGrad(1);
187     //
188     return Standard_True;
189   }
190   //
191   Standard_Real FirstParameter() const
192   {
193     return myFirst;
194   }
195
196   //
197   Standard_Real LastParameter() const
198   {
199     return myLast;
200   }
201   
202  private:
203   GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
204
205   //checks if the function can be computed when its parameter is
206   //equal to theParam
207    Standard_Boolean CheckParameter(const Standard_Real theParam) const
208    {
209      return ((myFirst <= theParam) && (theParam <= myLast));
210    }
211
212    const Adaptor3d_Curve& myCurve1;
213    const Adaptor3d_Curve& myCurve2;
214    const Standard_Real myFirst;
215    const Standard_Real myLast;
216 };
217
218 //=======================================================================
219 //class   : GeomLib_CheckCurveOnSurface_Local
220 //purpose : Created for parallelization possibility only
221 //=======================================================================
222 class GeomLib_CheckCurveOnSurface_Local
223 {
224 public:
225   GeomLib_CheckCurveOnSurface_Local(
226               const Handle(Geom_Curve)& theCurve3D,
227               const Handle(Geom2d_Curve)& theCurve2D,
228               const Handle(Geom_Surface)& theSurface,
229               const TColStd_Array1OfReal& theIntervalsArr,
230               const Standard_Real theEpsilonRange,
231               const Standard_Integer theNbParticles):
232   myCurve3D(theCurve3D),
233   myCurve2D(theCurve2D),
234   mySurface(theSurface),
235   mySubIntervals(theIntervalsArr),
236   myEpsilonRange(theEpsilonRange),
237   myNbParticles(theNbParticles),
238   myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1),
239   myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1)
240   {
241   }
242   
243   void operator()(const Standard_Integer& theIndex) const
244   {
245     //For every sub-interval (which is set by mySubIntervals array) this method
246     //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function.
247     //This optimal value will be put in corresponding (depending on theIndex - the
248     //identificator of the current interval in mySubIntervals array) cell of 
249     //myArrOfDist and myArrOfParam arrays.
250     const GeomAdaptor_Curve anAC(myCurve3D);
251     const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D);
252     const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface);
253
254     const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
255
256     GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
257                                                   mySubIntervals.Value(theIndex),
258                                                   mySubIntervals.Value(theIndex+1));
259
260     Standard_Real aMinDist = RealLast(), aPar = 0.0;
261     if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
262     {
263       myArrOfDist(theIndex) = RealLast();
264       myArrOfParam(theIndex) = aFunc.FirstParameter();
265       return;
266     }
267
268     myArrOfDist(theIndex) = aMinDist;
269     myArrOfParam(theIndex) = aPar;
270   }
271
272   //Returns optimal value (inverse of square of maximal distance)
273   void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
274   {
275     //This method looks for the minimal value of myArrOfDist.
276
277     const Standard_Integer aStartInd = myArrOfDist.Lower();
278     theMinimalValue = myArrOfDist(aStartInd);
279     theParameter = myArrOfParam(aStartInd);
280     for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++)
281     {
282       if(myArrOfDist(i) < theMinimalValue)
283       {
284         theMinimalValue = myArrOfDist(i);
285         theParameter = myArrOfParam(i);
286       }
287     }
288   }
289
290 private:
291   GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&);
292   const Handle(Geom_Curve)& myCurve3D;
293   const Handle(Geom2d_Curve)& myCurve2D;
294   const Handle(Geom_Surface)& mySurface;
295
296   const TColStd_Array1OfReal& mySubIntervals;
297   const Standard_Real myEpsilonRange;
298   const Standard_Integer myNbParticles;
299   mutable NCollection_Array1<Standard_Real> myArrOfDist;
300   mutable NCollection_Array1<Standard_Real> myArrOfParam;
301 };
302
303 //=======================================================================
304 //function : GeomLib_CheckCurveOnSurface
305 //purpose  : 
306 //=======================================================================
307 GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
308 :
309   myFirst(0.),
310   myLast(0.),
311   myErrorStatus(0),
312   myMaxDistance(RealLast()),
313   myMaxParameter(0.)
314 {
315 }
316
317 //=======================================================================
318 //function : GeomLib_CheckCurveOnSurface
319 //purpose  : 
320 //=======================================================================
321 GeomLib_CheckCurveOnSurface::
322   GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve,
323                               const Handle(Geom_Surface)& theSurface,
324                               const Standard_Real theFirst,
325                               const Standard_Real theLast):
326   myCurve(theCurve),
327   mySurface(theSurface),
328   myFirst(theFirst),
329   myLast(theLast),
330   myErrorStatus(0),
331   myMaxDistance(RealLast()),
332   myMaxParameter(0.)
333 {
334 }
335
336 //=======================================================================
337 //function : Init
338 //purpose  : 
339 //=======================================================================
340 void GeomLib_CheckCurveOnSurface::Init()
341 {
342   myCurve.Nullify();
343   mySurface.Nullify();
344   myFirst = 0.0;
345   myLast = 0.0;
346   myErrorStatus = 0;
347   myMaxDistance = RealLast();
348   myMaxParameter = 0.0;
349 }
350
351 //=======================================================================
352 //function : Init
353 //purpose  : 
354 //=======================================================================
355 void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve,
356                                         const Handle(Geom_Surface)& theSurface,
357                                         const Standard_Real theFirst,
358                                         const Standard_Real theLast)
359 {
360   myCurve = theCurve;
361   mySurface = theSurface;
362   myFirst = theFirst;
363   myLast = theLast;
364   myErrorStatus = 0;
365   myMaxDistance = RealLast();
366   myMaxParameter = 0.0;
367 }
368
369 //=======================================================================
370 //function : Perform
371 //purpose  : 
372 //=======================================================================
373
374 #ifndef HAVE_TBB
375 //After fixing bug # 26365, this fragment should be deleted
376 //(together the text "#ifdef HAVE_TBB")
377
378 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
379                                           const Standard_Boolean)
380 {
381   const Standard_Boolean isTheMTDisabled = Standard_True;
382 #else
383 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
384                                           const Standard_Boolean isTheMTDisabled)
385 {
386 #endif
387   if( myCurve.IsNull() ||
388       mySurface.IsNull() ||
389       thePCurve.IsNull())
390   {
391     myErrorStatus = 1;
392     return;
393   }
394
395   if( (myCurve->FirstParameter() > myFirst) ||
396       (myCurve->LastParameter() < myLast) ||
397       (thePCurve->FirstParameter() > myFirst) ||
398       (thePCurve->LastParameter() < myLast))
399   {
400     myErrorStatus = 2;
401     return;
402   }
403
404   const Standard_Real anEpsilonRange = 1.e-3;
405
406   Standard_Integer aNbParticles = 3;
407
408   //Polynomial function with degree n has not more than n-1 maxima and
409   //minima (degree of 1st derivative is equal to n-1 => 1st derivative has
410   //no greater than n-1 roots). Consequently, this function has
411   //maximum n monotonicity intervals. That is a good idea to try to put
412   //at least one particle in every monotonicity interval. Therefore,
413   //number of particles should be equal to n. 
414
415   const Standard_Integer aNbSubIntervals = 
416                               FillSubIntervals( myCurve, thePCurve,
417                                                 myFirst, myLast, aNbParticles);
418
419   if(!aNbSubIntervals)
420   {
421     myErrorStatus = 3;
422     return;
423   }
424
425   try {
426     OCC_CATCH_SIGNALS
427
428     TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
429     FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
430
431     GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
432                                 mySurface, anIntervals, anEpsilonRange, aNbParticles);
433
434     OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
435
436     aComp.OptimalValues(myMaxDistance, myMaxParameter);
437
438     myMaxDistance = sqrt(Abs(myMaxDistance));
439   }
440   catch (Standard_Failure) {
441     myErrorStatus = 3;
442   }
443 }
444
445 //=======================================================================
446 // Function : FillSubIntervals
447 // purpose : Divides [theFirst, theLast] interval on parts
448 //            in order to make searching-algorithm more precisely
449 //            (fills theSubIntervals array).
450 //            Returns number of subintervals.
451 //=======================================================================
452 Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d,
453                                   const Handle(Geom2d_Curve)& theCurve2d,
454                                   const Standard_Real theFirst,
455                                   const Standard_Real theLast,
456                                   Standard_Integer &theNbParticles,
457                                   TColStd_Array1OfReal* const theSubIntervals)
458 {
459   const Standard_Real anArrTempC[2] = {theFirst, theLast};
460   const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
461
462   theNbParticles = 3;
463   Handle(Geom2d_BSplineCurve) aBS2DCurv;
464   Handle(Geom_BSplineCurve) aBS3DCurv;
465
466   //
467   if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
468   {
469     aBS3DCurv = Handle(Geom_BSplineCurve)::
470                       DownCast(Handle(Geom_TrimmedCurve)::
471                       DownCast(theCurve3d)->BasisCurve());
472   }
473   else
474   {
475     aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
476   }
477
478   if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
479   {
480     aBS2DCurv = Handle(Geom2d_BSplineCurve)::
481                       DownCast(Handle(Geom2d_TrimmedCurve)::
482                       DownCast(theCurve2d)->BasisCurve());
483   }
484   else
485   {
486     aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
487   }
488
489   const TColStd_Array1OfReal &anArrKnots3D = !aBS3DCurv.IsNull() ? 
490                                               aBS3DCurv->Knots() :
491                                               anArrTemp;
492   const TColStd_Array1OfReal &anArrKnots2D = !aBS2DCurv.IsNull() ?
493                                               aBS2DCurv->Knots() :
494                                               anArrTemp;
495
496   Standard_Integer aNbSubIntervals = 1;
497
498   try
499   {
500     OCC_CATCH_SIGNALS
501     const Standard_Integer  anIndMax3D = anArrKnots3D.Upper(),
502                             anIndMax2D = anArrKnots2D.Upper();
503
504     Standard_Integer  anIndex3D = anArrKnots3D.Lower(),
505                       anIndex2D = anArrKnots2D.Lower();
506
507     if(theSubIntervals)
508       theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
509
510     while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
511     {
512       const Standard_Real aVal3D = anArrKnots3D.Value(anIndex3D),
513                           aVal2D = anArrKnots2D.Value(anIndex2D);
514       const Standard_Real aDelta = aVal3D - aVal2D;
515
516       if(aDelta < Precision::PConfusion())
517       {//aVal3D <= aVal2D
518         if((aVal3D > theFirst) && (aVal3D < theLast))
519         {
520           aNbSubIntervals++;
521         
522           if(theSubIntervals)
523             theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
524         }
525
526         anIndex3D++;
527
528         if(-aDelta < Precision::PConfusion())
529         {//aVal3D == aVal2D
530           anIndex2D++;
531         }
532       }
533       else
534       {//aVal2D < aVal3D
535         if((aVal2D > theFirst) && (aVal2D < theLast))
536         {
537           aNbSubIntervals++;
538           
539           if(theSubIntervals)
540             theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
541         }
542
543         anIndex2D++;
544       }
545     }
546
547     if(theSubIntervals)
548       theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
549
550     if(!aBS3DCurv.IsNull())
551     {
552       theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
553     }
554
555     if(!aBS2DCurv.IsNull())
556     {
557       theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
558     }
559   }
560   catch(Standard_Failure)
561   {
562 #ifdef OCCT_DEBUG
563     cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
564             "FillSubIntervals(): Incorrect filling!" << endl;
565 #endif
566
567     aNbSubIntervals = 0;
568   }
569
570   return aNbSubIntervals;
571 }
572
573 //=======================================================================
574 //class   : MinComputing
575 //purpose : Performs computing minimal value
576 //=======================================================================
577 Standard_Boolean MinComputing (
578                 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
579                 const Standard_Real theEpsilon, //1.0e-3
580                 const Standard_Integer theNbParticles,
581                 Standard_Real& theBestValue,
582                 Standard_Real& theBestParameter)
583 {
584   try
585   {
586     OCC_CATCH_SIGNALS
587
588     //They are used for finding a position of theNbParticles worst places
589     const Standard_Integer aNbControlPoints = 3*theNbParticles;
590     //
591     math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1), aStepPar(1,1);
592     aParInf(1) = theFunction.FirstParameter();
593     aParSup(1) = theFunction.LastParameter();
594     theBestParameter = aParInf(1);
595     theBestValue = RealLast();
596
597     const Standard_Real aDeltaParam = aParSup(1) - aParInf(1);
598     if(aDeltaParam < Precision::PConfusion())
599         return Standard_False;
600
601     aStepPar(1) = theEpsilon*aDeltaParam;
602
603     math_PSOParticlesPool aParticles(theNbParticles, 1);
604
605     const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1);
606     Standard_Integer aCount = 1;
607     for(Standard_Real aPrm = aParInf(1); aCount <= aNbControlPoints; aCount++,
608                       aPrm = (aCount == aNbControlPoints)? aParSup(1) : aPrm+aStep)
609     {
610       Standard_Real aVal = RealLast();
611       theFunction.Value(aPrm, aVal);
612
613       PSO_Particle* aParticle = aParticles.GetWorstParticle();
614
615       if(aVal > aParticle->BestDistance)
616         continue;
617
618       aParticle->Position[0] = aPrm;
619       aParticle->BestPosition[0] = aPrm;
620       aParticle->Distance     = aVal;
621       aParticle->BestDistance = aVal;
622     }
623
624     math_PSO aPSO(&theFunction, aParInf, aParSup, aStepPar);
625     aPSO.Perform(aParticles, theNbParticles, theBestValue, anOutputParam);
626
627     //Here, anOutputParam contains parameter, which is near to optimal.
628     //It needs to be more precise. Precision is made by math_NewtonMinimum.
629     math_NewtonMinimum anA(theFunction);
630     anA.Perform(theFunction, anOutputParam);
631
632     if(!anA.IsDone())
633     {
634 #ifdef OCCT_DEBUG
635       cout << "BRepLib_CheckCurveOnSurface::Compute(): No solution found!" << endl;
636 #endif
637       return Standard_False;
638     }
639
640     anA.Location(anOutputParam);
641     theBestParameter =  anOutputParam(1);
642     theBestValue = anA.Minimum();
643   }
644   catch(Standard_Failure)
645   {
646 #ifdef OCCT_DEBUG
647     cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
648 #endif
649     return Standard_False;
650   }
651
652   return Standard_True;
653 }