1 // Created by: Nikolai BUKHALOV
2 // Copyright (c) 2015 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
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>
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>
38 class GeomLib_CheckCurveOnSurface_TargetFunc;
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);
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);
55 //=======================================================================
56 //class : GeomLib_CheckCurveOnSurface_TargetFunc
57 //purpose : Target function (to be minimized)
58 //=======================================================================
59 class GeomLib_CheckCurveOnSurface_TargetFunc :
60 public math_MultipleVarFunctionWithHessian
63 GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
64 const Adaptor3d_Curve& theAdCS,
65 const Standard_Real theFirst,
66 const Standard_Real theLast):
74 //returns the number of parameters of the function
75 //(the function is one-dimension).
76 virtual Standard_Integer NbVariables() const {
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)
84 return Value(theX(1), theFVal);
87 //returns value of the one-dimension-function when parameter
89 Standard_Boolean Value( const Standard_Real theX,
90 Standard_Real& theFVal) const
95 if (!CheckParameter(theX))
96 return Standard_False;
98 const gp_Pnt aP1(myCurve1.Value(theX)),
99 aP2(myCurve2.Value(theX));
101 theFVal = -1.0*aP1.SquareDistance(aP2);
103 catch(Standard_Failure) {
104 return Standard_False;
107 return Standard_True;
110 //see analogical method for abstract owner class math_MultipleVarFunction
111 virtual Standard_Integer GetStateNumber()
116 //returns the gradient of the function when parameters are
118 virtual Standard_Boolean Gradient(const math_Vector& theX,
119 math_Vector& theGrad)
121 return Derive(theX(1), theGrad(1));
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
131 if (!CheckParameter(theX))
133 return Standard_False;
139 myCurve1.D1(theX, aP1, aDC1);
140 myCurve2.D1(theX, aP2, aDC2);
142 const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
144 theDeriv = -2.0*aVec1.Dot(aVec2);
146 catch(Standard_Failure)
148 return Standard_False;
151 return Standard_True;
154 //returns value and gradient
155 virtual Standard_Boolean Values(const math_Vector& theX,
156 Standard_Real& theVal,
157 math_Vector& theGrad)
159 if (!Value(theX, theVal))
161 return Standard_False;
164 if (!Gradient(theX, theGrad)) {
165 return Standard_False;
168 return Standard_True;
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)
177 if (!Value(theX, theVal))
179 return Standard_False;
182 if (!Gradient(theX, theGrad))
184 return Standard_False;
187 theHessian(1,1) = theGrad(1);
189 return Standard_True;
192 Standard_Real FirstParameter() const
198 Standard_Real LastParameter() const
204 GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
206 //checks if the function can be computed when its parameter is
208 Standard_Boolean CheckParameter(const Standard_Real theParam) const
210 return ((myFirst <= theParam) && (theParam <= myLast));
213 const Adaptor3d_Curve& myCurve1;
214 const Adaptor3d_Curve& myCurve2;
215 const Standard_Real myFirst;
216 const Standard_Real myLast;
219 //=======================================================================
220 //class : GeomLib_CheckCurveOnSurface_Local
221 //purpose : Created for parallelization possibility only
222 //=======================================================================
223 class GeomLib_CheckCurveOnSurface_Local
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)
244 void operator()(const Standard_Integer& theIndex) const
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);
255 const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
257 GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
258 mySubIntervals.Value(theIndex),
259 mySubIntervals.Value(theIndex+1));
261 Standard_Real aMinDist = RealLast(), aPar = 0.0;
262 if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
264 myArrOfDist(theIndex) = RealLast();
265 myArrOfParam(theIndex) = aFunc.FirstParameter();
269 myArrOfDist(theIndex) = aMinDist;
270 myArrOfParam(theIndex) = aPar;
273 //Returns optimal value (inverse of square of maximal distance)
274 void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
276 //This method looks for the minimal value of myArrOfDist.
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++)
283 if(myArrOfDist(i) < theMinimalValue)
285 theMinimalValue = myArrOfDist(i);
286 theParameter = myArrOfParam(i);
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;
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;
304 //=======================================================================
305 //function : GeomLib_CheckCurveOnSurface
307 //=======================================================================
308 GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
313 myMaxDistance(RealLast()),
315 myTolRange(Precision::PConfusion())
319 //=======================================================================
320 //function : GeomLib_CheckCurveOnSurface
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):
330 mySurface(theSurface),
334 myMaxDistance(RealLast()),
336 myTolRange(theTolRange)
340 //=======================================================================
343 //=======================================================================
344 void GeomLib_CheckCurveOnSurface::Init()
351 myMaxDistance = RealLast();
352 myMaxParameter = 0.0;
353 myTolRange = Precision::PConfusion();
356 //=======================================================================
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)
367 mySurface = theSurface;
371 myMaxDistance = RealLast();
372 myMaxParameter = 0.0;
373 myTolRange = theTolRange;
376 //=======================================================================
379 //=======================================================================
381 //After fixing bug # 26365, this fragment should be deleted
382 //(together the text "#ifdef HAVE_TBB")
384 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
385 const Standard_Boolean)
387 const Standard_Boolean isTheMTDisabled = Standard_True;
389 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
390 const Standard_Boolean isTheMTDisabled)
393 if( myCurve.IsNull() ||
394 mySurface.IsNull() ||
401 if(((myCurve->FirstParameter() - myFirst) > myTolRange) ||
402 ((myCurve->LastParameter() - myLast) < -myTolRange) ||
403 ((thePCurve->FirstParameter() - myFirst) > myTolRange) ||
404 ((thePCurve->LastParameter() - myLast) < -myTolRange))
410 const Standard_Real anEpsilonRange = 1.e-3;
412 Standard_Integer aNbParticles = 3;
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.
421 const Standard_Integer aNbSubIntervals =
422 FillSubIntervals( myCurve, thePCurve,
423 myFirst, myLast, aNbParticles);
434 TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
435 FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
437 GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
438 mySurface, anIntervals, anEpsilonRange, aNbParticles);
440 OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
442 aComp.OptimalValues(myMaxDistance, myMaxParameter);
444 myMaxDistance = sqrt(Abs(myMaxDistance));
446 catch (Standard_Failure) {
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)
465 const Standard_Integer aMaxKnots = 101;
466 const Standard_Real anArrTempC[2] = {theFirst, theLast};
467 const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
470 Handle(Geom2d_BSplineCurve) aBS2DCurv;
471 Handle(Geom_BSplineCurve) aBS3DCurv;
472 Standard_Boolean isTrimmed3D = Standard_False, isTrimmed2D = Standard_False;
475 if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
477 aBS3DCurv = Handle(Geom_BSplineCurve)::
478 DownCast(Handle(Geom_TrimmedCurve)::
479 DownCast(theCurve3d)->BasisCurve());
480 isTrimmed3D = Standard_True;
484 aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
487 if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
489 aBS2DCurv = Handle(Geom2d_BSplineCurve)::
490 DownCast(Handle(Geom2d_TrimmedCurve)::
491 DownCast(theCurve2d)->BasisCurve());
492 isTrimmed2D = Standard_True;
496 aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
499 Handle(TColStd_HArray1OfReal) anArrKnots3D, anArrKnots2D;
501 if(!aBS3DCurv.IsNull())
503 if(aBS3DCurv->NbKnots() <= aMaxKnots)
505 anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
509 Standard_Integer KnotCount;
514 const TColStd_Array1OfReal& aKnots = aBS3DCurv->Knots();
515 for(i = aBS3DCurv->FirstUKnotIndex(); i <= aBS3DCurv->LastUKnotIndex(); ++i)
517 if(aKnots(i) > theFirst && aKnots(i) < theLast)
526 KnotCount = aBS3DCurv->LastUKnotIndex() - aBS3DCurv->FirstUKnotIndex() + 1;
528 if(KnotCount <= aMaxKnots)
530 anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
534 anArrKnots3D = new TColStd_HArray1OfReal(1, aMaxKnots);
535 anArrKnots3D->SetValue(1, theFirst);
536 anArrKnots3D->SetValue(aMaxKnots, theLast);
538 Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
539 Standard_Real t = theFirst + dt;
540 for(i = 2; i < aMaxKnots; ++i, t += dt)
542 anArrKnots3D->SetValue(i, t);
549 anArrKnots3D = new TColStd_HArray1OfReal(anArrTemp);
551 if(!aBS2DCurv.IsNull())
553 if(aBS2DCurv->NbKnots() <= aMaxKnots)
555 anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
559 Standard_Integer KnotCount;
564 const TColStd_Array1OfReal& aKnots = aBS2DCurv->Knots();
565 for(i = aBS2DCurv->FirstUKnotIndex(); i <= aBS2DCurv->LastUKnotIndex(); ++i)
567 if(aKnots(i) > theFirst && aKnots(i) < theLast)
576 KnotCount = aBS2DCurv->LastUKnotIndex() - aBS2DCurv->FirstUKnotIndex() + 1;
578 if(KnotCount <= aMaxKnots)
580 anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
584 anArrKnots2D = new TColStd_HArray1OfReal(1, aMaxKnots);
585 anArrKnots2D->SetValue(1, theFirst);
586 anArrKnots2D->SetValue(aMaxKnots, theLast);
588 Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
589 Standard_Real t = theFirst + dt;
590 for(i = 2; i < aMaxKnots; ++i, t += dt)
592 anArrKnots2D->SetValue(i, t);
599 anArrKnots2D = new TColStd_HArray1OfReal(anArrTemp);
603 Standard_Integer aNbSubIntervals = 1;
608 const Standard_Integer anIndMax3D = anArrKnots3D->Upper(),
609 anIndMax2D = anArrKnots2D->Upper();
611 Standard_Integer anIndex3D = anArrKnots3D->Lower(),
612 anIndex2D = anArrKnots2D->Lower();
615 theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
617 while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
619 const Standard_Real aVal3D = anArrKnots3D->Value(anIndex3D),
620 aVal2D = anArrKnots2D->Value(anIndex2D);
621 const Standard_Real aDelta = aVal3D - aVal2D;
623 if(aDelta < Precision::PConfusion())
625 if((aVal3D > theFirst) && (aVal3D < theLast))
630 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
635 if(-aDelta < Precision::PConfusion())
642 if((aVal2D > theFirst) && (aVal2D < theLast))
647 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
655 theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
657 if(!aBS3DCurv.IsNull())
659 theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
662 if(!aBS2DCurv.IsNull())
664 theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
667 catch(Standard_Failure)
670 cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
671 "FillSubIntervals(): Incorrect filling!" << endl;
677 return aNbSubIntervals;
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)
692 const Standard_Real aDeltaParam = theParSup(1) - theParInf(1);
693 if(aDeltaParam < Precision::PConfusion())
694 return Standard_False;
696 math_Vector aStepPar(1, 1);
697 aStepPar(1) = theEpsilon*aDeltaParam;
699 math_PSOParticlesPool aParticles(theNbParticles, 1);
701 //They are used for finding a position of theNbParticles worst places
702 const Standard_Integer aNbControlPoints = 3*theNbParticles;
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)
709 Standard_Real aVal = RealLast();
710 if(!theFunction.Value(aPrm, aVal))
713 PSO_Particle* aParticle = aParticles.GetWorstParticle();
715 if(aVal > aParticle->BestDistance)
718 aParticle->Position[0] = aPrm;
719 aParticle->BestPosition[0] = aPrm;
720 aParticle->Distance = aVal;
721 aParticle->BestDistance = aVal;
724 math_PSO aPSO(&theFunction, theParInf, theParSup, aStepPar);
725 aPSO.Perform(aParticles, theNbParticles, theBestValue, theOutputParam);
727 return Standard_True;
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)
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();
752 if(!PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
753 theBestValue, anOutputParam))
756 cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << endl;
758 return Standard_False;
761 theBestParameter = anOutputParam(1);
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);
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();
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;
780 Standard_Real aValue = RealLast();
781 if(PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
782 aValue, anOutputParam))
784 if(aValue < theBestValue)
786 theBestValue = aValue;
787 theBestParameter = anOutputParam(1);
792 catch(Standard_Failure)
795 cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
797 return Standard_False;
800 return Standard_True;