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>
37 class GeomLib_CheckCurveOnSurface_TargetFunc;
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);
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);
54 //=======================================================================
55 //class : GeomLib_CheckCurveOnSurface_TargetFunc
56 //purpose : Target function (to be minimized)
57 //=======================================================================
58 class GeomLib_CheckCurveOnSurface_TargetFunc :
59 public math_MultipleVarFunctionWithHessian
62 GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
63 const Adaptor3d_Curve& theAdCS,
64 const Standard_Real theFirst,
65 const Standard_Real theLast):
73 //returns the number of parameters of the function
74 //(the function is one-dimension).
75 virtual Standard_Integer NbVariables() const {
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)
83 return Value(theX(1), theFVal);
86 //returns value of the one-dimension-function when parameter
88 Standard_Boolean Value( const Standard_Real theX,
89 Standard_Real& theFVal) const
94 if (!CheckParameter(theX))
95 return Standard_False;
97 const gp_Pnt aP1(myCurve1.Value(theX)),
98 aP2(myCurve2.Value(theX));
100 theFVal = -1.0*aP1.SquareDistance(aP2);
102 catch(Standard_Failure) {
103 return Standard_False;
106 return Standard_True;
109 //see analogical method for abstract owner class math_MultipleVarFunction
110 virtual Standard_Integer GetStateNumber()
115 //returns the gradient of the function when parameters are
117 virtual Standard_Boolean Gradient(const math_Vector& theX,
118 math_Vector& theGrad)
120 return Derive(theX(1), theGrad(1));
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
130 if (!CheckParameter(theX))
132 return Standard_False;
138 myCurve1.D1(theX, aP1, aDC1);
139 myCurve2.D1(theX, aP2, aDC2);
141 const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
143 theDeriv = -2.0*aVec1.Dot(aVec2);
145 catch(Standard_Failure)
147 return Standard_False;
150 return Standard_True;
153 //returns value and gradient
154 virtual Standard_Boolean Values(const math_Vector& theX,
155 Standard_Real& theVal,
156 math_Vector& theGrad)
158 if (!Value(theX, theVal))
160 return Standard_False;
163 if (!Gradient(theX, theGrad)) {
164 return Standard_False;
167 return Standard_True;
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)
176 if (!Value(theX, theVal))
178 return Standard_False;
181 if (!Gradient(theX, theGrad))
183 return Standard_False;
186 theHessian(1,1) = theGrad(1);
188 return Standard_True;
191 Standard_Real FirstParameter() const
197 Standard_Real LastParameter() const
203 GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
205 //checks if the function can be computed when its parameter is
207 Standard_Boolean CheckParameter(const Standard_Real theParam) const
209 return ((myFirst <= theParam) && (theParam <= myLast));
212 const Adaptor3d_Curve& myCurve1;
213 const Adaptor3d_Curve& myCurve2;
214 const Standard_Real myFirst;
215 const Standard_Real myLast;
218 //=======================================================================
219 //class : GeomLib_CheckCurveOnSurface_Local
220 //purpose : Created for parallelization possibility only
221 //=======================================================================
222 class GeomLib_CheckCurveOnSurface_Local
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)
243 void operator()(const Standard_Integer& theIndex) const
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);
254 const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
256 GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
257 mySubIntervals.Value(theIndex),
258 mySubIntervals.Value(theIndex+1));
260 Standard_Real aMinDist = RealLast(), aPar = 0.0;
261 if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
263 myArrOfDist(theIndex) = RealLast();
264 myArrOfParam(theIndex) = aFunc.FirstParameter();
268 myArrOfDist(theIndex) = aMinDist;
269 myArrOfParam(theIndex) = aPar;
272 //Returns optimal value (inverse of square of maximal distance)
273 void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
275 //This method looks for the minimal value of myArrOfDist.
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++)
282 if(myArrOfDist(i) < theMinimalValue)
284 theMinimalValue = myArrOfDist(i);
285 theParameter = myArrOfParam(i);
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;
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;
303 //=======================================================================
304 //function : GeomLib_CheckCurveOnSurface
306 //=======================================================================
307 GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
312 myMaxDistance(RealLast()),
317 //=======================================================================
318 //function : GeomLib_CheckCurveOnSurface
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):
327 mySurface(theSurface),
331 myMaxDistance(RealLast()),
336 //=======================================================================
339 //=======================================================================
340 void GeomLib_CheckCurveOnSurface::Init()
347 myMaxDistance = RealLast();
348 myMaxParameter = 0.0;
351 //=======================================================================
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)
361 mySurface = theSurface;
365 myMaxDistance = RealLast();
366 myMaxParameter = 0.0;
369 //=======================================================================
372 //=======================================================================
375 //After fixing bug # 26365, this fragment should be deleted
376 //(together the text "#ifdef HAVE_TBB")
378 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
379 const Standard_Boolean)
381 const Standard_Boolean isTheMTDisabled = Standard_True;
383 void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
384 const Standard_Boolean isTheMTDisabled)
387 if( myCurve.IsNull() ||
388 mySurface.IsNull() ||
395 if( (myCurve->FirstParameter() > myFirst) ||
396 (myCurve->LastParameter() < myLast) ||
397 (thePCurve->FirstParameter() > myFirst) ||
398 (thePCurve->LastParameter() < myLast))
404 const Standard_Real anEpsilonRange = 1.e-3;
406 Standard_Integer aNbParticles = 3;
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.
415 const Standard_Integer aNbSubIntervals =
416 FillSubIntervals( myCurve, thePCurve,
417 myFirst, myLast, aNbParticles);
428 TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
429 FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
431 GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
432 mySurface, anIntervals, anEpsilonRange, aNbParticles);
434 OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
436 aComp.OptimalValues(myMaxDistance, myMaxParameter);
438 myMaxDistance = sqrt(Abs(myMaxDistance));
440 catch (Standard_Failure) {
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)
459 const Standard_Real anArrTempC[2] = {theFirst, theLast};
460 const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
463 Handle(Geom2d_BSplineCurve) aBS2DCurv;
464 Handle(Geom_BSplineCurve) aBS3DCurv;
467 if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
469 aBS3DCurv = Handle(Geom_BSplineCurve)::
470 DownCast(Handle(Geom_TrimmedCurve)::
471 DownCast(theCurve3d)->BasisCurve());
475 aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
478 if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
480 aBS2DCurv = Handle(Geom2d_BSplineCurve)::
481 DownCast(Handle(Geom2d_TrimmedCurve)::
482 DownCast(theCurve2d)->BasisCurve());
486 aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
489 const TColStd_Array1OfReal &anArrKnots3D = !aBS3DCurv.IsNull() ?
492 const TColStd_Array1OfReal &anArrKnots2D = !aBS2DCurv.IsNull() ?
496 Standard_Integer aNbSubIntervals = 1;
501 const Standard_Integer anIndMax3D = anArrKnots3D.Upper(),
502 anIndMax2D = anArrKnots2D.Upper();
504 Standard_Integer anIndex3D = anArrKnots3D.Lower(),
505 anIndex2D = anArrKnots2D.Lower();
508 theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
510 while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
512 const Standard_Real aVal3D = anArrKnots3D.Value(anIndex3D),
513 aVal2D = anArrKnots2D.Value(anIndex2D);
514 const Standard_Real aDelta = aVal3D - aVal2D;
516 if(aDelta < Precision::PConfusion())
518 if((aVal3D > theFirst) && (aVal3D < theLast))
523 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
528 if(-aDelta < Precision::PConfusion())
535 if((aVal2D > theFirst) && (aVal2D < theLast))
540 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
548 theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
550 if(!aBS3DCurv.IsNull())
552 theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
555 if(!aBS2DCurv.IsNull())
557 theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
560 catch(Standard_Failure)
563 cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
564 "FillSubIntervals(): Incorrect filling!" << endl;
570 return aNbSubIntervals;
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)
588 //They are used for finding a position of theNbParticles worst places
589 const Standard_Integer aNbControlPoints = 3*theNbParticles;
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();
597 const Standard_Real aDeltaParam = aParSup(1) - aParInf(1);
598 if(aDeltaParam < Precision::PConfusion())
599 return Standard_False;
601 aStepPar(1) = theEpsilon*aDeltaParam;
603 math_PSOParticlesPool aParticles(theNbParticles, 1);
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)
610 Standard_Real aVal = RealLast();
611 theFunction.Value(aPrm, aVal);
613 PSO_Particle* aParticle = aParticles.GetWorstParticle();
615 if(aVal > aParticle->BestDistance)
618 aParticle->Position[0] = aPrm;
619 aParticle->BestPosition[0] = aPrm;
620 aParticle->Distance = aVal;
621 aParticle->BestDistance = aVal;
624 math_PSO aPSO(&theFunction, aParInf, aParSup, aStepPar);
625 aPSO.Perform(aParticles, theNbParticles, theBestValue, anOutputParam);
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);
635 cout << "BRepLib_CheckCurveOnSurface::Compute(): No solution found!" << endl;
637 return Standard_False;
640 anA.Location(anOutputParam);
641 theBestParameter = anOutputParam(1);
642 theBestValue = anA.Minimum();
644 catch(Standard_Failure)
647 cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
649 return Standard_False;
652 return Standard_True;