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 const&) {
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& theDeriv1, Standard_Real* const theDeriv2 = 0) const
131 if (!CheckParameter(theX))
133 return Standard_False;
137 gp_Vec aDC1, aDC2, aDCC1, aDCC2;
141 myCurve1.D1(theX, aP1, aDC1);
142 myCurve2.D1(theX, aP2, aDC2);
146 myCurve1.D2(theX, aP1, aDC1, aDCC1);
147 myCurve2.D2(theX, aP2, aDC2, aDCC2);
150 const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
152 theDeriv1 = -2.0*aVec1.Dot(aVec2);
156 const gp_Vec aVec3(aDCC2 - aDCC1);
157 *theDeriv2 = -2.0*(aVec2.SquareMagnitude() + aVec1.Dot(aVec3));
160 catch(Standard_Failure const&)
162 return Standard_False;
165 return Standard_True;
168 //returns value and gradient
169 virtual Standard_Boolean Values(const math_Vector& theX,
170 Standard_Real& theVal,
171 math_Vector& theGrad)
173 if (!Value(theX, theVal))
175 return Standard_False;
178 if (!Gradient(theX, theGrad)) {
179 return Standard_False;
182 return Standard_True;
185 //returns value, gradient and hessian
186 virtual Standard_Boolean Values(const math_Vector& theX,
187 Standard_Real& theVal,
188 math_Vector& theGrad,
189 math_Matrix& theHessian)
191 if (!Value(theX, theVal))
193 return Standard_False;
196 if (!Derive(theX(1), theGrad(1), &theHessian(1, 1)))
198 return Standard_False;
201 return Standard_True;
204 Standard_Real FirstParameter() const
210 Standard_Real LastParameter() const
216 GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
218 //checks if the function can be computed when its parameter is
220 Standard_Boolean CheckParameter(const Standard_Real theParam) const
222 return ((myFirst <= theParam) && (theParam <= myLast));
225 const Adaptor3d_Curve& myCurve1;
226 const Adaptor3d_Curve& myCurve2;
227 const Standard_Real myFirst;
228 const Standard_Real myLast;
231 //=======================================================================
232 //class : GeomLib_CheckCurveOnSurface_Local
233 //purpose : Created for parallelization possibility only
234 //=======================================================================
235 class GeomLib_CheckCurveOnSurface_Local
238 GeomLib_CheckCurveOnSurface_Local(
239 const Handle(Geom_Curve)& theCurve3D,
240 const Handle(Geom2d_Curve)& theCurve2D,
241 const Handle(Geom_Surface)& theSurface,
242 const TColStd_Array1OfReal& theIntervalsArr,
243 const Standard_Real theEpsilonRange,
244 const Standard_Integer theNbParticles):
245 myCurve3D(theCurve3D),
246 myCurve2D(theCurve2D),
247 mySurface(theSurface),
248 mySubIntervals(theIntervalsArr),
249 myEpsilonRange(theEpsilonRange),
250 myNbParticles(theNbParticles),
251 myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1),
252 myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1)
256 void operator()(const Standard_Integer& theIndex) const
258 //For every sub-interval (which is set by mySubIntervals array) this method
259 //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function.
260 //This optimal value will be put in corresponding (depending on theIndex - the
261 //identificator of the current interval in mySubIntervals array) cell of
262 //myArrOfDist and myArrOfParam arrays.
263 const GeomAdaptor_Curve anAC(myCurve3D);
264 const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D);
265 const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface);
267 const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
269 GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
270 mySubIntervals.Value(theIndex),
271 mySubIntervals.Value(theIndex+1));
273 Standard_Real aMinDist = RealLast(), aPar = 0.0;
274 if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
276 myArrOfDist(theIndex) = RealLast();
277 myArrOfParam(theIndex) = aFunc.FirstParameter();
281 myArrOfDist(theIndex) = aMinDist;
282 myArrOfParam(theIndex) = aPar;
285 //Returns optimal value (inverse of square of maximal distance)
286 void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
288 //This method looks for the minimal value of myArrOfDist.
290 const Standard_Integer aStartInd = myArrOfDist.Lower();
291 theMinimalValue = myArrOfDist(aStartInd);
292 theParameter = myArrOfParam(aStartInd);
293 for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++)
295 if(myArrOfDist(i) < theMinimalValue)
297 theMinimalValue = myArrOfDist(i);
298 theParameter = myArrOfParam(i);
304 GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&);
305 const Handle(Geom_Curve)& myCurve3D;
306 const Handle(Geom2d_Curve)& myCurve2D;
307 const Handle(Geom_Surface)& mySurface;
309 const TColStd_Array1OfReal& mySubIntervals;
310 const Standard_Real myEpsilonRange;
311 const Standard_Integer myNbParticles;
312 mutable NCollection_Array1<Standard_Real> myArrOfDist;
313 mutable NCollection_Array1<Standard_Real> myArrOfParam;
316 //=======================================================================
317 //function : GeomLib_CheckCurveOnSurface
319 //=======================================================================
320 GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
325 myMaxDistance(RealLast()),
327 myTolRange(Precision::PConfusion())
331 //=======================================================================
332 //function : GeomLib_CheckCurveOnSurface
334 //=======================================================================
335 GeomLib_CheckCurveOnSurface::
336 GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve,
337 const Handle(Geom_Surface)& theSurface,
338 const Standard_Real theFirst,
339 const Standard_Real theLast,
340 const Standard_Real theTolRange):
342 mySurface(theSurface),
346 myMaxDistance(RealLast()),
348 myTolRange(theTolRange)
352 //=======================================================================
355 //=======================================================================
356 void GeomLib_CheckCurveOnSurface::Init()
363 myMaxDistance = RealLast();
364 myMaxParameter = 0.0;
365 myTolRange = Precision::PConfusion();
368 //=======================================================================
371 //=======================================================================
372 void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve,
373 const Handle(Geom_Surface)& theSurface,
374 const Standard_Real theFirst,
375 const Standard_Real theLast,
376 const Standard_Real theTolRange)
379 mySurface = theSurface;
383 myMaxDistance = RealLast();
384 myMaxParameter = 0.0;
385 myTolRange = theTolRange;
388 //=======================================================================
391 //=======================================================================
393 //After fixing bug # 26365, this fragment should be deleted
394 //(together the text "#ifdef HAVE_TBB")
396 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
397 const Standard_Boolean)
399 const Standard_Boolean isTheMTDisabled = Standard_True;
401 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
402 const Standard_Boolean isTheMTDisabled)
405 if( myCurve.IsNull() ||
406 mySurface.IsNull() ||
413 if(((myCurve->FirstParameter() - myFirst) > myTolRange) ||
414 ((myCurve->LastParameter() - myLast) < -myTolRange) ||
415 ((thePCurve->FirstParameter() - myFirst) > myTolRange) ||
416 ((thePCurve->LastParameter() - myLast) < -myTolRange))
422 const Standard_Real anEpsilonRange = 1.e-3;
424 Standard_Integer aNbParticles = 3;
426 //Polynomial function with degree n has not more than n-1 maxima and
427 //minima (degree of 1st derivative is equal to n-1 => 1st derivative has
428 //no greater than n-1 roots). Consequently, this function has
429 //maximum n monotonicity intervals. That is a good idea to try to put
430 //at least one particle in every monotonicity interval. Therefore,
431 //number of particles should be equal to n.
433 const Standard_Integer aNbSubIntervals =
434 FillSubIntervals( myCurve, thePCurve,
435 myFirst, myLast, aNbParticles);
446 TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
447 FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
449 GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
450 mySurface, anIntervals, anEpsilonRange, aNbParticles);
452 OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
454 aComp.OptimalValues(myMaxDistance, myMaxParameter);
456 myMaxDistance = sqrt(Abs(myMaxDistance));
458 catch (Standard_Failure const&) {
463 //=======================================================================
464 // Function : FillSubIntervals
465 // purpose : Divides [theFirst, theLast] interval on parts
466 // in order to make searching-algorithm more precisely
467 // (fills theSubIntervals array).
468 // Returns number of subintervals.
469 //=======================================================================
470 Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d,
471 const Handle(Geom2d_Curve)& theCurve2d,
472 const Standard_Real theFirst,
473 const Standard_Real theLast,
474 Standard_Integer &theNbParticles,
475 TColStd_Array1OfReal* const theSubIntervals)
477 const Standard_Integer aMaxKnots = 101;
478 const Standard_Real anArrTempC[2] = {theFirst, theLast};
479 const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
482 Handle(Geom2d_BSplineCurve) aBS2DCurv;
483 Handle(Geom_BSplineCurve) aBS3DCurv;
484 Standard_Boolean isTrimmed3D = Standard_False, isTrimmed2D = Standard_False;
487 if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
489 aBS3DCurv = Handle(Geom_BSplineCurve)::
490 DownCast(Handle(Geom_TrimmedCurve)::
491 DownCast(theCurve3d)->BasisCurve());
492 isTrimmed3D = Standard_True;
496 aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
499 if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
501 aBS2DCurv = Handle(Geom2d_BSplineCurve)::
502 DownCast(Handle(Geom2d_TrimmedCurve)::
503 DownCast(theCurve2d)->BasisCurve());
504 isTrimmed2D = Standard_True;
508 aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
511 Handle(TColStd_HArray1OfReal) anArrKnots3D, anArrKnots2D;
513 if(!aBS3DCurv.IsNull())
515 if(aBS3DCurv->NbKnots() <= aMaxKnots)
517 anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
521 Standard_Integer KnotCount;
526 const TColStd_Array1OfReal& aKnots = aBS3DCurv->Knots();
527 for(i = aBS3DCurv->FirstUKnotIndex(); i <= aBS3DCurv->LastUKnotIndex(); ++i)
529 if(aKnots(i) > theFirst && aKnots(i) < theLast)
538 KnotCount = aBS3DCurv->LastUKnotIndex() - aBS3DCurv->FirstUKnotIndex() + 1;
540 if(KnotCount <= aMaxKnots)
542 anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots());
546 anArrKnots3D = new TColStd_HArray1OfReal(1, aMaxKnots);
547 anArrKnots3D->SetValue(1, theFirst);
548 anArrKnots3D->SetValue(aMaxKnots, theLast);
550 Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
551 Standard_Real t = theFirst + dt;
552 for(i = 2; i < aMaxKnots; ++i, t += dt)
554 anArrKnots3D->SetValue(i, t);
561 anArrKnots3D = new TColStd_HArray1OfReal(anArrTemp);
563 if(!aBS2DCurv.IsNull())
565 if(aBS2DCurv->NbKnots() <= aMaxKnots)
567 anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
571 Standard_Integer KnotCount;
576 const TColStd_Array1OfReal& aKnots = aBS2DCurv->Knots();
577 for(i = aBS2DCurv->FirstUKnotIndex(); i <= aBS2DCurv->LastUKnotIndex(); ++i)
579 if(aKnots(i) > theFirst && aKnots(i) < theLast)
588 KnotCount = aBS2DCurv->LastUKnotIndex() - aBS2DCurv->FirstUKnotIndex() + 1;
590 if(KnotCount <= aMaxKnots)
592 anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots());
596 anArrKnots2D = new TColStd_HArray1OfReal(1, aMaxKnots);
597 anArrKnots2D->SetValue(1, theFirst);
598 anArrKnots2D->SetValue(aMaxKnots, theLast);
600 Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1);
601 Standard_Real t = theFirst + dt;
602 for(i = 2; i < aMaxKnots; ++i, t += dt)
604 anArrKnots2D->SetValue(i, t);
611 anArrKnots2D = new TColStd_HArray1OfReal(anArrTemp);
615 Standard_Integer aNbSubIntervals = 1;
620 const Standard_Integer anIndMax3D = anArrKnots3D->Upper(),
621 anIndMax2D = anArrKnots2D->Upper();
623 Standard_Integer anIndex3D = anArrKnots3D->Lower(),
624 anIndex2D = anArrKnots2D->Lower();
627 theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
629 while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
631 const Standard_Real aVal3D = anArrKnots3D->Value(anIndex3D),
632 aVal2D = anArrKnots2D->Value(anIndex2D);
633 const Standard_Real aDelta = aVal3D - aVal2D;
635 if(aDelta < Precision::PConfusion())
637 if((aVal3D > theFirst) && (aVal3D < theLast))
642 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
647 if(-aDelta < Precision::PConfusion())
654 if((aVal2D > theFirst) && (aVal2D < theLast))
659 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
667 theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
669 if(!aBS3DCurv.IsNull())
671 theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
674 if(!aBS2DCurv.IsNull())
676 theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
679 catch(Standard_Failure const&)
682 std::cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
683 "FillSubIntervals(): Incorrect filling!" << std::endl;
689 return aNbSubIntervals;
692 //=======================================================================
693 //class : PSO_Perform
694 //purpose : Searches minimal distance with math_PSO class
695 //=======================================================================
696 Standard_Boolean PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
697 const math_Vector &theParInf,
698 const math_Vector &theParSup,
699 const Standard_Real theEpsilon,
700 const Standard_Integer theNbParticles,
701 Standard_Real& theBestValue,
702 math_Vector &theOutputParam)
704 const Standard_Real aDeltaParam = theParSup(1) - theParInf(1);
705 if(aDeltaParam < Precision::PConfusion())
706 return Standard_False;
708 math_Vector aStepPar(1, 1);
709 aStepPar(1) = theEpsilon*aDeltaParam;
711 math_PSOParticlesPool aParticles(theNbParticles, 1);
713 //They are used for finding a position of theNbParticles worst places
714 const Standard_Integer aNbControlPoints = 3*theNbParticles;
716 const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1);
717 Standard_Integer aCount = 1;
718 for(Standard_Real aPrm = theParInf(1); aCount <= aNbControlPoints; aCount++,
719 aPrm = (aCount == aNbControlPoints)? theParSup(1) : aPrm+aStep)
721 Standard_Real aVal = RealLast();
722 if(!theFunction.Value(aPrm, aVal))
725 PSO_Particle* aParticle = aParticles.GetWorstParticle();
727 if(aVal > aParticle->BestDistance)
730 aParticle->Position[0] = aPrm;
731 aParticle->BestPosition[0] = aPrm;
732 aParticle->Distance = aVal;
733 aParticle->BestDistance = aVal;
736 math_PSO aPSO(&theFunction, theParInf, theParSup, aStepPar);
737 aPSO.Perform(aParticles, theNbParticles, theBestValue, theOutputParam);
739 return Standard_True;
742 //=======================================================================
743 //class : MinComputing
744 //purpose : Performs computing minimal value
745 //=======================================================================
746 Standard_Boolean MinComputing (
747 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
748 const Standard_Real theEpsilon, //1.0e-3
749 const Standard_Integer theNbParticles,
750 Standard_Real& theBestValue,
751 Standard_Real& theBestParameter)
758 math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1);
759 aParInf(1) = theFunction.FirstParameter();
760 aParSup(1) = theFunction.LastParameter();
761 theBestParameter = aParInf(1);
762 theBestValue = RealLast();
764 if(!PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
765 theBestValue, anOutputParam))
768 std::cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << std::endl;
770 return Standard_False;
773 theBestParameter = anOutputParam(1);
775 //Here, anOutputParam contains parameter, which is near to optimal.
776 //It needs to be more precise. Precision is made by math_NewtonMinimum.
777 math_NewtonMinimum aMinSol(theFunction);
778 aMinSol.Perform(theFunction, anOutputParam);
780 if(aMinSol.IsDone() && (aMinSol.GetStatus() == math_OK))
781 {//math_NewtonMinimum has precised the value. We take it.
782 aMinSol.Location(anOutputParam);
783 theBestParameter = anOutputParam(1);
784 theBestValue = aMinSol.Minimum();
787 {//Use math_PSO again but on smaller range.
788 const Standard_Real aStep = theEpsilon*(aParSup(1) - aParInf(1));
789 aParInf(1) = theBestParameter - 0.5*aStep;
790 aParSup(1) = theBestParameter + 0.5*aStep;
792 Standard_Real aValue = RealLast();
793 if(PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles,
794 aValue, anOutputParam))
796 if(aValue < theBestValue)
798 theBestValue = aValue;
799 theBestParameter = anOutputParam(1);
804 catch(Standard_Failure const&)
807 std::cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << std::endl;
809 return Standard_False;
812 return Standard_True;