// 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. //Modified by MPS : 21-05-97 PRO 7598 // Le nombre de solutions rendu est mySqDist.Length() et non // mySqDist.Length()/2 // ajout des valeurs absolues dans le test d'orthogonalite de // GetStateNumber() /*----------------------------------------------------------------------------- Fonctions permettant de rechercher une distance extremale entre 2 courbes C1 et C2 (en partant de points approches C1(u0) et C2(v0)). Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par l'algorithme math_FunctionSetRoot. Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont: { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| } { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| } Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1 et F2 sont egales a: { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2 { Dvf1(u,v) = Dv.Du/||Du|| { Duf2(u,v) = -Du.Dv/||Dv|| { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2 ----------------------------------------------------------------------------*/ #include static const Standard_Real MinTol = 1.e-20; static const Standard_Real TolFactor = 1.e-12; static const Standard_Real MinStep = 1e-7; static const Standard_Integer MaxOrder = 3; //============================================================================= Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C) { const Standard_Integer NPoint = 10; Standard_Real aStartParam, anEndParam; if(C==myC1) { aStartParam = myUinfium; anEndParam = myUsupremum; } else if(C==myC2) { aStartParam = myVinfium; anEndParam = myVsupremum; } else { //Warning: No curve for tolerance computing! return MinTol; } const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint; Standard_Integer aNum = 0; Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative //(it is computed with using NPoint point) do { Standard_Real u = aStartParam + aNum*aStep; //parameter for every point if(u > anEndParam) u = anEndParam; Pnt Ptemp; //empty point (is not used below) Vec VDer; // 1st derivative vector Tool1::D1(*((Curve1*)C), u, Ptemp, VDer); Standard_Real vm = VDer.Magnitude(); if(vm > aMax) aMax = vm; } while(++aNum < NPoint+1); return Max(aMax*TolFactor,MinTol); } //============================================================================= Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol) { math_Vector V1(1,2), V2(1,2); V1(1) = 0.0; V2(1) = 0.0; V1(2) = 0.0; V2(2) = 0.0; SubIntervalInitialize(V1, V2); myMaxDerivOrderC1 = 0; myTolC1=MinTol; myMaxDerivOrderC2 = 0; myTolC2=MinTol; } //============================================================================= Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1, const Curve2& C2, const Standard_Real thetol) : myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2), myTol (thetol) { math_Vector V1(1,2), V2(1,2); V1(1) = Tool1::FirstParameter(*((Curve1*)myC1)); V2(1) = Tool1::LastParameter(*((Curve1*)myC1)); V1(2) = Tool2::FirstParameter(*((Curve2*)myC2)); V2(2) = Tool2::LastParameter(*((Curve2*)myC2)); SubIntervalInitialize(V1, V2); switch(Tool1::GetType(*((Curve1*)myC1))) { case GeomAbs_BezierCurve: case GeomAbs_BSplineCurve: case GeomAbs_OffsetCurve: case GeomAbs_OtherCurve: myMaxDerivOrderC1 = MaxOrder; myTolC1 = SearchOfTolerance((Standard_Address)&C1); break; default: myMaxDerivOrderC1 = 0; myTolC1=MinTol; break; } switch(Tool2::GetType(*((Curve2*)myC2))) { case GeomAbs_BezierCurve: case GeomAbs_BSplineCurve: case GeomAbs_OffsetCurve: case GeomAbs_OtherCurve: myMaxDerivOrderC2 = MaxOrder; myTolC2 = SearchOfTolerance((Standard_Address)&C2); break; default: myMaxDerivOrderC2 = 0; myTolC2=MinTol; break; } } //============================================================================= void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C) { Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()") if (theRank == 1) { myC1 = (Standard_Address)&C; switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType()) { case GeomAbs_BezierCurve: case GeomAbs_BSplineCurve: case GeomAbs_OffsetCurve: case GeomAbs_OtherCurve: myMaxDerivOrderC1 = MaxOrder; myTolC1 = SearchOfTolerance((Standard_Address)&C); break; default: myMaxDerivOrderC1 = 0; myTolC1=MinTol; break; } } else if (theRank == 2) { myC2 = (Standard_Address)&C; switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType()) { case GeomAbs_BezierCurve: case GeomAbs_BSplineCurve: case GeomAbs_OffsetCurve: case GeomAbs_OtherCurve: myMaxDerivOrderC2 = MaxOrder; myTolC2 = SearchOfTolerance((Standard_Address)&C); break; default: myMaxDerivOrderC2 = 0; myTolC2=MinTol; break; } } } //============================================================================= Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F) { myU = UV(1); myV = UV(2); Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu); Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv); Vec P1P2 (myP1,myP2); Standard_Real Ndu = myDu.Magnitude(); if(myMaxDerivOrderC1 != 0) { if (Ndu <= myTolC1) { //Derivative is approximated by Taylor-series const Standard_Real DivisionFactor = 1.e-3; Standard_Real du; if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) du = 0.0; else du = myUsupremum-myUinfium; const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); Standard_Integer n = 1; //Derivative order Vec V; Standard_Boolean IsDeriveFound; do { V = Tool1::DN(*((Curve1*)myC1),myU,++n); Ndu = V.Magnitude(); IsDeriveFound = (Ndu > myTolC1); } while(!IsDeriveFound && n < myMaxDerivOrderC1); if(IsDeriveFound) { Standard_Real u; if(myU-myUinfium < aDelta) u = myU+aDelta; else u = myU-aDelta; Pnt P1, P2; Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1); Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2); Vec V1(P1,P2); Standard_Real aDirFactor = V.Dot(V1); if(aDirFactor < 0.0) myDu = -V; else myDu = V; }//if(IsDeriveFound) else { //Derivative is approximated by three points Pnt Ptemp; //(0,0,0)-coordinate Pnt P1, P2, P3; Standard_Boolean IsParameterGrown; if(myU-myUinfium < 2*aDelta) { Tool1::D0(*((Curve1*)myC1),myU,P1); Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2); Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3); IsParameterGrown = Standard_True; } else { Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1); Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2); Tool1::D0(*((Curve1*)myC1),myU,P3); IsParameterGrown = Standard_False; } Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); if(IsParameterGrown) myDu=-3*V1+4*V2-V3; else myDu=V1-4*V2+3*V3; }//else of if(IsDeriveFound) Ndu = myDu.Magnitude(); }//if (Ndu <= myTolC1) condition }//if(myMaxDerivOrder != 0) if (Ndu <= MinTol) { //Warning: 1st derivative of C1 is equal to zero! return Standard_False; } Standard_Real Ndv = myDv.Magnitude(); if(myMaxDerivOrderC2 != 0) { if (Ndv <= myTolC2) { const Standard_Real DivisionFactor = 1.e-3; Standard_Real dv; if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst())) dv = 0.0; else dv = myVsupremum-myVinfium; const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep); //Derivative is approximated by Taylor-series Standard_Integer n = 1; //Derivative order Vec V; Standard_Boolean IsDeriveFound; do { V = Tool2::DN(*((Curve2*)myC2),myV,++n); Ndv = V.Magnitude(); IsDeriveFound = (Ndv > myTolC2); } while(!IsDeriveFound && n < myMaxDerivOrderC2); if(IsDeriveFound) { Standard_Real v; if(myV-myVinfium < aDelta) v = myV+aDelta; else v = myV-aDelta; Pnt P1, P2; Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1); Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2); Vec V1(P1,P2); Standard_Real aDirFactor = V.Dot(V1); if(aDirFactor < 0.0) myDv = -V; else myDv = V; }//if(IsDeriveFound) else { //Derivative is approximated by three points Pnt Ptemp; //(0,0,0)-coordinate Pnt P1, P2, P3; Standard_Boolean IsParameterGrown; if(myV-myVinfium < 2*aDelta) { Tool2::D0(*((Curve2*)myC2),myV,P1); Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2); Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3); IsParameterGrown = Standard_True; } else { Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1); Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2); Tool2::D0(*((Curve2*)myC2),myV,P3); IsParameterGrown = Standard_False; } Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); if(IsParameterGrown) myDv=-3*V1+4*V2-V3; else myDv=V1-4*V2+3*V3; }//else of if(IsDeriveFound) Ndv = myDv.Magnitude(); }//if (Ndv <= myTolC2) }//if(myMaxDerivOrder != 0) if (Ndv <= MinTol) { //1st derivative of C2 is equal to zero! return Standard_False; } F(1) = P1P2.Dot(myDu)/Ndu; F(2) = P1P2.Dot(myDv)/Ndv; return Standard_True; } //============================================================================= Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV, math_Matrix& Df) { math_Vector F(1,2); return Values(UV,F,Df); } //============================================================================= Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV, math_Vector& F, math_Matrix& Df) { myU = UV(1); myV = UV(2); if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv { //Warning: No function value found! return Standard_False; }//if(Value(UV, F) == Standard_False) Vec Du, Dv, Duu, Dvv; Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu); Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv); //Calling of "Value(...)" function change class member values. //After running it is necessary to return to previous values. const Standard_Real myU_old = myU, myV_old = myV; const Pnt myP1_old = myP1, myP2_old = myP2; const Vec myDu_old = myDu, myDv_old = myDv; //Attention: aDelta value must be greater than same value for "Value(...)" // function to avoid of points' collisions. const Standard_Real DivisionFactor = 0.01; Standard_Real du; if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) du = 0.0; else du = myUsupremum-myUinfium; const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep); Standard_Real dv; if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst())) dv = 0.0; else dv = myVsupremum-myVinfium; const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); Vec P1P2 (myP1,myP2); if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) { //Derivative is approximated by three points math_Vector FF1(1,2), FF2(1,2), FF3(1,2); Standard_Real F1, F2, F3; /////////////////////////// Search of DF1_u derivative (begin) /////////////////// if(myU-myUinfium < 2*aDeltaU) { F1=F(1); math_Vector UV2(1,2), UV3(1,2); UV2(1)=myU+aDeltaU; UV2(2)=myV; UV3(1)=myU+2*aDeltaU; UV3(2)=myV; if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F2 = FF2(1); F3 = FF3(1); Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); }//if(myU-myUinfium < 2*aDeltaU) else { F3 = F(1); math_Vector UV2(1,2), UV1(1,2); UV2(1)=myU-aDeltaU; UV2(2)=myV; UV1(1)=myU-2*aDeltaU; UV1(2)=myV; if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F1 = FF1(1); F2 = FF2(1); Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); }//else of if(myU-myUinfium < 2*aDeltaU) condition /////////////////////////// Search of DF1_u derivative (end) /////////////////// //Return to previous values myU = myU_old; myV = myV_old; /////////////////////////// Search of DF1_v derivative (begin) /////////////////// if(myV-myVinfium < 2*aDeltaV) { F1=F(1); math_Vector UV2(1,2), UV3(1,2); UV2(1)=myU; UV2(2)=myV+aDeltaV; UV3(1)=myU; UV3(2)=myV+2*aDeltaV; if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F2 = FF2(1); F3 = FF3(1); Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); }//if(myV-myVinfium < 2*aDeltaV) else { F3 = F(1); math_Vector UV2(1,2), UV1(1,2); UV2(1)=myU; UV2(2)=myV-aDeltaV; UV1(1)=myU; UV1(2)=myV-2*aDeltaV; if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F1 = FF1(1); F2 = FF2(1); Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); }//else of if(myV-myVinfium < 2*aDeltaV) /////////////////////////// Search of DF1_v derivative (end) /////////////////// //Return to previous values myU = myU_old; myV = myV_old; myP1 = myP1_old, myP2 = myP2_old; myDu = myDu_old, myDv = myDv_old; }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) else { const Standard_Real Ndu = myDu.Magnitude(); Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu)); Df(1,2) = myDv.Dot(myDu)/Ndu; }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) { //Derivative is approximated by three points math_Vector FF1(1,2), FF2(1,2), FF3(1,2); Standard_Real F1, F2, F3; /////////////////////////// Search of DF2_v derivative (begin) /////////////////// if(myV-myVinfium < 2*aDeltaV) { F1=F(2); math_Vector UV2(1,2), UV3(1,2); UV2(1)=myU; UV2(2)=myV+aDeltaV; UV3(1)=myU; UV3(2)=myV+2*aDeltaV; if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F2 = FF2(2); F3 = FF3(2); Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); }//if(myV-myVinfium < 2*aDeltaV) else { F3 = F(2); math_Vector UV2(1,2), UV1(1,2); UV2(1)=myU; UV2(2)=myV-aDeltaV; UV1(1)=myU; UV1(2)=myV-2*aDeltaV; if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F1 = FF1(2); F2 = FF2(2); Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); }//else of if(myV-myVinfium < 2*aDeltaV) /////////////////////////// Search of DF2_v derivative (end) /////////////////// //Return to previous values myU = myU_old; myV = myV_old; /////////////////////////// Search of DF2_u derivative (begin) /////////////////// if(myU-myUinfium < 2*aDeltaU) { F1=F(2); math_Vector UV2(1,2), UV3(1,2); UV2(1)=myU+aDeltaU; UV2(2)=myV; UV3(1)=myU+2*aDeltaU; UV3(2)=myV; if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) { //There are many points close to singularity points and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F2 = FF2(2); F3 = FF3(2); Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); }//if(myU-myUinfium < 2*aDelta) else { F3 = F(2); math_Vector UV2(1,2), UV1(1,2); UV2(1)=myU-aDeltaU; UV2(2)=myV; UV1(1)=myU-2*aDeltaU; UV1(2)=myV; if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) { //There are many points close to singularity points //and which have zero-derivative. //Try to decrease aDelta variable's value. return Standard_False; } F1 = FF1(2); F2 = FF2(2); Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); }//else of if(myU-myUinfium < 2*aDeltaU) /////////////////////////// Search of DF2_u derivative (end) /////////////////// //Return to previous values myU = myU_old; myV = myV_old; myP1 = myP1_old; myP2 = myP2_old; myDu = myDu_old; myDv = myDv_old; }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) else { Standard_Real Ndv = myDv.Magnitude(); Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv)); Df(2,1) = -myDu.Dot(myDv)/Ndv; }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) return Standard_True; }//end of function //============================================================================= Standard_Integer Extrema_FuncExtCC::GetStateNumber () { Vec Du (myDu), Dv (myDv); Vec P1P2 (myP1, myP2); Standard_Real mod = Du.Magnitude(); if(mod > myTolC1) { Du /= mod; } mod = Dv.Magnitude(); if(mod > myTolC2) { Dv /= mod; } if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) { mySqDist.Append(myP1.SquareDistance(myP2)); myPoints.Append(POnC(myU, myP1)); myPoints.Append(POnC(myV, myP2)); } return 0; } //============================================================================= void Extrema_FuncExtCC::Points (const Standard_Integer N, POnC& P1, POnC& P2) const { P1 = myPoints.Value(2*N-1); P2 = myPoints.Value(2*N); } //============================================================================= void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound, const math_Vector& theSupBound) { myUinfium = theInfBound(1); myUsupremum = theSupBound(1); myVinfium = theInfBound(2); myVsupremum = theSupBound(2); } //=============================================================================