// Created on: 1995-07-18 // Created by: Modelistation // Copyright (c) 1995-1999 Matra Datavision // Copyright (c) 1999-2014 OPEN CASCADE SAS // // This file is part of Open CASCADE Technology software library. // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License version 2.1 as published // by the Free Software Foundation, with special exception defined in the file // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT // distribution for complete text of the license and disclaimer of any warranty. // // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const Standard_Real MaxParamVal = 1.0e+10; const Standard_Real aBorderDivisor = 1.0e+4; const Standard_Real HyperbolaLimit = 23.; //ln(MaxParamVal) static Standard_Boolean IsQuadric(const GeomAbs_SurfaceType theSType) { if (theSType == GeomAbs_Plane) return Standard_True; if (theSType == GeomAbs_Cylinder) return Standard_True; if (theSType == GeomAbs_Cone) return Standard_True; if (theSType == GeomAbs_Sphere) return Standard_True; if (theSType == GeomAbs_Torus) return Standard_True; return Standard_False; } static Standard_Boolean IsConic(const GeomAbs_CurveType theCType) { if (theCType == GeomAbs_Line) return Standard_True; if (theCType == GeomAbs_Circle) return Standard_True; if (theCType == GeomAbs_Ellipse) return Standard_True; if (theCType == GeomAbs_Hyperbola) return Standard_True; if (theCType == GeomAbs_Parabola) return Standard_True; return Standard_False; } // restrict maximal parameter on hyperbola to avoid FPE static Standard_Real GetCurvMaxParamVal (const Adaptor3d_Curve& theC) { if (theC.GetType() == GeomAbs_Hyperbola) { return HyperbolaLimit; } if (theC.GetType() == GeomAbs_OffsetCurve) { Handle(Geom_Curve) aBC (theC.OffsetCurve()->BasisCurve()); Handle(Geom_TrimmedCurve) aTC = Handle(Geom_TrimmedCurve)::DownCast (aBC); if (! aTC.IsNull()) { aBC = aTC->BasisCurve(); } if (aBC->IsKind (STANDARD_TYPE(Geom_Hyperbola))) return HyperbolaLimit; } return MaxParamVal; } // restrict maximal parameter on surfaces based on hyperbola to avoid FPE static void GetSurfMaxParamVals (const Adaptor3d_Surface& theS, Standard_Real& theUmax, Standard_Real& theVmax) { theUmax = theVmax = MaxParamVal; if (theS.GetType() == GeomAbs_SurfaceOfExtrusion) { theUmax = GetCurvMaxParamVal (theS.BasisCurve()->Curve()); } else if (theS.GetType() == GeomAbs_SurfaceOfRevolution) { theVmax = GetCurvMaxParamVal (theS.BasisCurve()->Curve()); } else if (theS.GetType() == GeomAbs_OffsetSurface) { GetSurfMaxParamVals (theS.BasisSurface()->Surface(), theUmax, theVmax); } } //======================================================================= //function : Extrema_GenExtCS //purpose : //======================================================================= Extrema_GenExtCS::Extrema_GenExtCS() : myDone(Standard_False), mytmin(0.0), mytsup(0.0), myumin(0.0), myusup(0.0), myvmin(0.0), myvsup(0.0), mytsample(0), myusample(0), myvsample(0), mytol1(0.0), mytol2(0.0), myS(NULL) { } //======================================================================= //function : Extrema_GenExtCS //purpose : //======================================================================= Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, const Adaptor3d_Surface& S, const Standard_Integer NbT, const Standard_Integer NbU, const Standard_Integer NbV, const Standard_Real Tol1, const Standard_Real Tol2) { Initialize(S, NbU, NbV, Tol2); Perform(C, NbT, Tol1); } //======================================================================= //function : Extrema_GenExtCS //purpose : //======================================================================= Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C, const Adaptor3d_Surface& S, const Standard_Integer NbT, const Standard_Integer NbU, const Standard_Integer NbV, const Standard_Real tmin, const Standard_Real tsup, const Standard_Real Umin, const Standard_Real Usup, const Standard_Real Vmin, const Standard_Real Vsup, const Standard_Real Tol1, const Standard_Real Tol2) { Initialize(S, NbU, NbV, Umin,Usup,Vmin,Vsup,Tol2); Perform(C, NbT, tmin, tsup, Tol1); } //======================================================================= //function : Initialize //purpose : //======================================================================= void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, const Standard_Integer NbU, const Standard_Integer NbV, const Standard_Real Tol2) { myumin = S.FirstUParameter(); myusup = S.LastUParameter(); myvmin = S.FirstVParameter(); myvsup = S.LastVParameter(); Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,Tol2); } //======================================================================= //function : Initialize //purpose : //======================================================================= void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, const Standard_Integer NbU, const Standard_Integer NbV, const Standard_Real Umin, const Standard_Real Usup, const Standard_Real Vmin, const Standard_Real Vsup, const Standard_Real Tol2) { myS = (Adaptor3d_SurfacePtr)&S; myusample = NbU; myvsample = NbV; myumin = Umin; myusup = Usup; myvmin = Vmin; myvsup = Vsup; mytol2 = Tol2; Standard_Real umaxpar, vmaxpar; GetSurfMaxParamVals(*myS, umaxpar, vmaxpar); if(Precision::IsInfinite (myusup)) { myusup = umaxpar; } if(Precision::IsInfinite (myumin)) { myumin = -umaxpar; } if(Precision::IsInfinite (myvsup)) { myvsup = vmaxpar; } if(Precision::IsInfinite (myvmin)) { myvmin = -vmaxpar; } Standard_Real du = (myusup - myumin) / aBorderDivisor; Standard_Real dv = (myvsup - myvmin) / aBorderDivisor; const Standard_Real aMinU = myumin + du; const Standard_Real aMinV = myvmin + dv; const Standard_Real aMaxU = myusup - du; const Standard_Real aMaxV = myvsup - dv; const Standard_Real aStepSU = (aMaxU - aMinU) / myusample; const Standard_Real aStepSV = (aMaxV - aMinV) / myvsample; mySurfPnts = new TColgp_HArray2OfPnt (0, myusample, 0, myvsample); Standard_Real aSU = aMinU; for (Standard_Integer aSUI = 0; aSUI <= myusample; aSUI++, aSU += aStepSU) { Standard_Real aSV = aMinV; for (Standard_Integer aSVI = 0; aSVI <= myvsample; aSVI++, aSV += aStepSV) { mySurfPnts->ChangeValue (aSUI, aSVI) = myS->Value (aSU, aSV); } } } //======================================================================= //function : Perform //purpose : //======================================================================= void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C, const Standard_Integer NbT, const Standard_Real Tol1) { mytmin = C.FirstParameter(); mytsup = C.LastParameter(); Perform(C, NbT, mytmin, mytsup,Tol1); } //======================================================================= //function : Perform //purpose : //======================================================================= void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, const Standard_Integer NbT, const Standard_Real tmin, const Standard_Real tsup, const Standard_Real Tol1) { myDone = Standard_False; myF.Initialize(C,*myS); mytmin = tmin; mytsup = tsup; mytol1 = Tol1; mytsample = NbT; // Modif de lvt pour trimer la surface non pas aux infinis mais a +/- 10000 Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin; Standard_Real aCMaxVal = GetCurvMaxParamVal (C); if (Precision::IsInfinite(mytsup)){ mytsup = aCMaxVal; } if (Precision::IsInfinite(mytmin)){ mytmin = -aCMaxVal; } // Standard_Integer aNbVar = 3; GeomAbs_SurfaceType aSType = myS->GetType(); if (IsQuadric(aSType)) { aNbVar = 1; } else { GeomAbs_CurveType aCType = C.GetType(); if (IsConic(aCType)) { aNbVar = 2; } } math_Vector Tol(1, 3), TUV(1, 3), TUVinf(1, 3), TUVsup(1, 3); // Tol(1) = mytol1; Tol(2) = mytol2; Tol(3) = mytol2; // TUVinf(1) = mytmin; TUVinf(2) = trimumin; TUVinf(3) = trimvmin; // TUVsup(1) = mytsup; TUVsup(2) = trimusup; TUVsup(3) = trimvsup; // // Number of particles used in PSO algorithm (particle swarm optimization). const Standard_Integer aNbParticles = 48; // if (aNbVar == 3) { GlobMinGenCS(C, aNbParticles, TUVinf, TUVsup, TUV); } else if (aNbVar == 2) { GlobMinConicS(C, aNbParticles, TUVinf, TUVsup, TUV); } else { GlobMinCQuadric(C, aNbParticles, TUVinf, TUVsup, TUV); } // Find min approximation math_FunctionSetRoot anA(myF, Tol); anA.Perform(myF, TUV, TUVinf, TUVsup); myDone = Standard_True; } //======================================================================= //function : GlobMinGenCS //purpose : //======================================================================= void Extrema_GenExtCS::GlobMinGenCS(const Adaptor3d_Curve& theC, const Standard_Integer theNbParticles, const math_Vector& theTUVinf, const math_Vector& theTUVsup, math_Vector& theTUV) { math_PSOParticlesPool aParticles(theNbParticles, 3); math_Vector aMinTUV(1, 3); aMinTUV = theTUVinf + (theTUVsup - theTUVinf) / aBorderDivisor; math_Vector aMaxTUV(1, 3); aMaxTUV = theTUVsup - (theTUVsup - theTUVinf) / aBorderDivisor; Standard_Real aStepCU = (aMaxTUV(1) - aMinTUV(1)) / mytsample; Standard_Real aStepSU = (aMaxTUV(2) - aMinTUV(2)) / myusample; Standard_Real aStepSV = (aMaxTUV(3) - aMinTUV(3)) / myvsample; // Correct number of curve samples in case of low resolution Standard_Integer aNewCsample = mytsample; Standard_Real aScaleFactor = 5.0; Standard_Real aResolutionCU = aStepCU / theC.Resolution(1.0); Standard_Real aMinResolution = aScaleFactor * Min(aResolutionCU, Min(aStepSU / myS->UResolution(1.0), aStepSV / myS->VResolution(1.0))); if (aMinResolution > Epsilon(1.0)) { if (aResolutionCU > aMinResolution) { const Standard_Integer aMaxNbNodes = 50; aNewCsample = Min(aMaxNbNodes, RealToInt(mytsample * aResolutionCU / aMinResolution)); aStepCU = (aMaxTUV(1) - aMinTUV(1)) / aNewCsample; } } // Pre-compute curve sample points. TColgp_Array1OfPnt aCurvPnts(0, aNewCsample); Standard_Real aCU1 = aMinTUV(1); for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCU1 += aStepCU) aCurvPnts.SetValue(aCUI, theC.Value(aCU1)); PSO_Particle* aParticle = aParticles.GetWorstParticle(); // Select specified number of particles from pre-computed set of samples Standard_Real aSU = aMinTUV(2); for (Standard_Integer aSUI = 0; aSUI <= myusample; aSUI++, aSU += aStepSU) { Standard_Real aSV = aMinTUV(3); for (Standard_Integer aSVI = 0; aSVI <= myvsample; aSVI++, aSV += aStepSV) { Standard_Real aCU2 = aMinTUV(1); for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCU2 += aStepCU) { Standard_Real aSqDist = mySurfPnts->Value(aSUI, aSVI).SquareDistance(aCurvPnts.Value(aCUI)); if (aSqDist < aParticle->Distance) { aParticle->Position[0] = aCU2; aParticle->Position[1] = aSU; aParticle->Position[2] = aSV; aParticle->BestPosition[0] = aCU2; aParticle->BestPosition[1] = aSU; aParticle->BestPosition[2] = aSV; aParticle->Distance = aSqDist; aParticle->BestDistance = aSqDist; aParticle = aParticles.GetWorstParticle(); } } } } math_Vector aStep(1, 3); aStep(1) = aStepCU; aStep(2) = aStepSU; aStep(3) = aStepSV; // Find min approximation Standard_Real aValue; Extrema_GlobOptFuncCS aFunc(&theC, myS); math_PSO aPSO(&aFunc, theTUVinf, theTUVsup, aStep); aPSO.Perform(aParticles, theNbParticles, aValue, theTUV); } //======================================================================= //function : GlobMinConicS //purpose : //======================================================================= void Extrema_GenExtCS::GlobMinConicS(const Adaptor3d_Curve& theC, const Standard_Integer theNbParticles, const math_Vector& theTUVinf, const math_Vector& theTUVsup, math_Vector& theTUV) { Standard_Integer aNbVar = 2; math_Vector anUVinf(1, aNbVar), anUVsup(1, aNbVar), anUV(1, aNbVar); Standard_Integer i; for (i = 1; i <= aNbVar; ++i) { anUVinf(i) = theTUVinf(i + 1); anUVsup(i) = theTUVsup(i + 1); } // math_PSOParticlesPool aParticles(theNbParticles, aNbVar); math_Vector aMinUV(1, aNbVar); aMinUV = anUVinf + (anUVsup - anUVinf) / aBorderDivisor; math_Vector aMaxUV(1, aNbVar); aMaxUV = anUVsup - (anUVsup - anUVinf) / aBorderDivisor; //Increase numbers of UV samples to improve searching global minimum Standard_Integer anAddsample = Max(mytsample / 2, 3); Standard_Integer anUsample = myusample + anAddsample; Standard_Integer aVsample = myvsample + anAddsample; // Standard_Real aStepSU = (aMaxUV(1) - aMinUV(1)) / anUsample; Standard_Real aStepSV = (aMaxUV(2) - aMinUV(2)) / aVsample; // Extrema_GlobOptFuncConicS aFunc(myS, anUVinf(1), anUVsup(1), anUVinf(2), anUVsup(2)); aFunc.LoadConic(&theC, theTUVinf(1), theTUVsup(1)); PSO_Particle* aParticle = aParticles.GetWorstParticle(); // Select specified number of particles from pre-computed set of samples Standard_Real aSU = aMinUV(1); for (Standard_Integer aSUI = 0; aSUI <= anUsample; aSUI++, aSU += aStepSU) { anUV(1) = aSU; Standard_Real aSV = aMinUV(2); for (Standard_Integer aSVI = 0; aSVI <= aVsample; aSVI++, aSV += aStepSV) { anUV(2) = aSV; Standard_Real aSqDist; if (!aFunc.Value(anUV, aSqDist)) { aSqDist = Precision::Infinite(); } if (aSqDist < aParticle->Distance) { aParticle->Position[0] = aSU; aParticle->Position[1] = aSV; aParticle->BestPosition[0] = aSU; aParticle->BestPosition[1] = aSV; aParticle->Distance = aSqDist; aParticle->BestDistance = aSqDist; aParticle = aParticles.GetWorstParticle(); } } } math_Vector aStep(1, aNbVar); aStep(1) = aStepSU; aStep(2) = aStepSV; // Find min approximation Standard_Real aValue; math_PSO aPSO(&aFunc, anUVinf, anUVsup, aStep); aPSO.Perform(aParticles, theNbParticles, aValue, anUV); // Standard_Real aCT = aFunc.ConicParameter(anUV); if (theC.IsPeriodic()) { if (aCT < theTUVinf(1) - Precision::PConfusion() || aCT > theTUVsup(1) + Precision::PConfusion()) { aCT = ElCLib::InPeriod(aCT, theTUVinf(1), theTUVinf(1) + 2. * M_PI); } } theTUV(1) = aCT; theTUV(2) = anUV(1); theTUV(3) = anUV(2); } //======================================================================= //function : GlobMinCQuadric //purpose : //======================================================================= void Extrema_GenExtCS::GlobMinCQuadric(const Adaptor3d_Curve& theC, const Standard_Integer theNbParticles, const math_Vector& theTUVinf, const math_Vector& theTUVsup, math_Vector& theTUV) { Standard_Integer aNbVar = 1; math_Vector aTinf(1, aNbVar), aTsup(1, aNbVar), aT(1, aNbVar); aTinf(1) = theTUVinf(1); aTsup(1) = theTUVsup(1); // math_PSOParticlesPool aParticles(theNbParticles, aNbVar); math_Vector aMinT(1, aNbVar); aMinT = aTinf + (aTsup - aTinf) / aBorderDivisor; math_Vector aMaxT(1, aNbVar); aMaxT = aTsup - (aTsup - aTinf) / aBorderDivisor; // //Increase numbers of curve samples to improve searching global minimum //because dimension of optimisation task is redused const Standard_Integer aMaxNbNodes = 50; Standard_Integer aNewCsample = mytsample; Standard_Integer anAddsample = Max(myusample / 2, 3); aNewCsample += anAddsample; aNewCsample = Min(aNewCsample, aMaxNbNodes); // // Correct number of curve samples in case of low resolution Standard_Real aStepCT = (aMaxT(1) - aMinT(1)) / aNewCsample; Standard_Real aStepSU = (theTUVsup(2) - theTUVinf(2)) / myusample; Standard_Real aStepSV = (theTUVsup(3) - theTUVinf(3)) / myvsample; Standard_Real aScaleFactor = 5.0; Standard_Real aResolutionCU = aStepCT / theC.Resolution(1.0); Standard_Real aMinResolution = aScaleFactor * Min(aResolutionCU, Min(aStepSU / myS->UResolution(1.0), aStepSV / myS->VResolution(1.0))); if (aMinResolution > Epsilon(1.0)) { if (aResolutionCU > aMinResolution) { aNewCsample = Min(aMaxNbNodes, RealToInt(aNewCsample * aResolutionCU / aMinResolution)); aStepCT = (aMaxT(1) - aMinT(1)) / aNewCsample; } } // Extrema_GlobOptFuncCQuadric aFunc(&theC, aTinf(1), aTsup(1)); aFunc.LoadQuad(myS, theTUVinf(2), theTUVsup(2), theTUVinf(3), theTUVsup(3)); PSO_Particle* aParticle = aParticles.GetWorstParticle(); // Select specified number of particles from pre-computed set of samples Standard_Real aCT = aMinT(1); for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCT += aStepCT) { aT(1) = aCT; Standard_Real aSqDist; if (!aFunc.Value(aT, aSqDist)) { aSqDist = Precision::Infinite(); } if (aSqDist < aParticle->Distance) { aParticle->Position[0] = aCT; aParticle->BestPosition[0] = aCT; aParticle->Distance = aSqDist; aParticle->BestDistance = aSqDist; aParticle = aParticles.GetWorstParticle(); } } // math_Vector aStep(1, aNbVar); aStep(1) = aStepCT; // Find min approximation Standard_Real aValue; math_PSO aPSO(&aFunc, aTinf, aTsup, aStep); aPSO.Perform(aParticles, theNbParticles, aValue, aT); // math_Vector anUV(1, 2); aFunc.QuadricParameters(aT, anUV); if (myS->IsUPeriodic()) { if (anUV(1) < theTUVinf(2) - Precision::PConfusion() || anUV(1)> theTUVsup(2) + Precision::PConfusion()) { anUV(1) = ElCLib::InPeriod(anUV(1), theTUVinf(2), theTUVinf(2) + 2. * M_PI); } } // if (myS->IsVPeriodic()) { if (anUV(2) < theTUVinf(3) - Precision::PConfusion() || anUV(2)> theTUVsup(3) + Precision::PConfusion()) { anUV(2) = ElCLib::InPeriod(anUV(2), theTUVinf(3), theTUVinf(3) + 2. * M_PI); } } // theTUV(1) = aT(1); theTUV(2) = anUV(1); theTUV(3) = anUV(2); } //======================================================================= //function : IsDone //purpose : //======================================================================= Standard_Boolean Extrema_GenExtCS::IsDone() const { return myDone; } //======================================================================= //function : NbExt //purpose : //======================================================================= Standard_Integer Extrema_GenExtCS::NbExt() const { if (!IsDone()) { throw StdFail_NotDone(); } return myF.NbExt(); } //======================================================================= //function : Value //purpose : //======================================================================= Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const { if (N < 1 || N > NbExt()) { throw Standard_OutOfRange(); } return myF.SquareDistance(N); } //======================================================================= //function : PointOnCurve //purpose : //======================================================================= const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const { if (N < 1 || N > NbExt()) { throw Standard_OutOfRange(); } return myF.PointOnCurve(N); } //======================================================================= //function : PointOnSurface //purpose : //======================================================================= const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const { if (N < 1 || N > NbExt()) { throw Standard_OutOfRange(); } return myF.PointOnSurface(N); } //======================================================================= //function : BidonSurface //purpose : //======================================================================= Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const { return (Adaptor3d_SurfacePtr)0L; } //======================================================================= //function : BidonCurve //purpose : //======================================================================= Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const { return (Adaptor3d_CurvePtr)0L; }