1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 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 //Modified by MPS : 21-05-97 PRO 7598
16 // Le nombre de solutions rendu est mySqDist.Length() et non
17 // mySqDist.Length()/2
18 // ajout des valeurs absolues dans le test d'orthogonalite de
21 /*-----------------------------------------------------------------------------
22 Fonctions permettant de rechercher une distance extremale entre 2 courbes
23 C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
24 Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
25 l'algorithme math_FunctionSetRoot.
26 Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
27 { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
28 { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
29 Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
31 { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
32 { Dvf1(u,v) = Dv.Du/||Du||
33 { Duf2(u,v) = -Du.Dv/||Dv||
34 { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
36 ----------------------------------------------------------------------------*/
38 #include <Precision.hxx>
41 static const Standard_Real MinTol = 1.e-20;
42 static const Standard_Real TolFactor = 1.e-12;
43 static const Standard_Real MinStep = 1e-7;
45 static const Standard_Integer MaxOrder = 3;
49 //=============================================================================
50 Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
52 const Standard_Integer NPoint = 10;
53 Standard_Real aStartParam, anEndParam;
57 aStartParam = myUinfium;
58 anEndParam = myUsupremum;
62 aStartParam = myVinfium;
63 anEndParam = myVsupremum;
67 //Warning: No curve for tolerance computing!
71 const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
73 Standard_Integer aNum = 0;
74 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
75 //(it is computed with using NPoint point)
79 Standard_Real u = aStartParam + aNum*aStep; //parameter for every point
83 Pnt Ptemp; //empty point (is not used below)
84 Vec VDer; // 1st derivative vector
85 Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
86 Standard_Real vm = VDer.Magnitude();
90 while(++aNum < NPoint+1);
92 return Max(aMax*TolFactor,MinTol);
95 //=============================================================================
96 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
98 math_Vector V1(1,2), V2(1,2);
103 SubIntervalInitialize(V1, V2);
104 myMaxDerivOrderC1 = 0;
106 myMaxDerivOrderC2 = 0;
110 //=============================================================================
111 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
113 const Standard_Real thetol) :
114 myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
117 math_Vector V1(1,2), V2(1,2);
118 V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
119 V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
120 V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
121 V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
122 SubIntervalInitialize(V1, V2);
124 switch(Tool1::GetType(*((Curve1*)myC1)))
126 case GeomAbs_BezierCurve:
127 case GeomAbs_BSplineCurve:
128 case GeomAbs_OffsetCurve:
129 case GeomAbs_OtherCurve:
130 myMaxDerivOrderC1 = MaxOrder;
131 myTolC1 = SearchOfTolerance((Standard_Address)&C1);
134 myMaxDerivOrderC1 = 0;
139 switch(Tool2::GetType(*((Curve2*)myC2)))
141 case GeomAbs_BezierCurve:
142 case GeomAbs_BSplineCurve:
143 case GeomAbs_OffsetCurve:
144 case GeomAbs_OtherCurve:
145 myMaxDerivOrderC2 = MaxOrder;
146 myTolC2 = SearchOfTolerance((Standard_Address)&C2);
149 myMaxDerivOrderC2 = 0;
155 //=============================================================================
156 void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
158 Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
162 myC1 = (Standard_Address)&C;
163 switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
165 case GeomAbs_BezierCurve:
166 case GeomAbs_BSplineCurve:
167 case GeomAbs_OffsetCurve:
168 case GeomAbs_OtherCurve:
169 myMaxDerivOrderC1 = MaxOrder;
170 myTolC1 = SearchOfTolerance((Standard_Address)&C);
173 myMaxDerivOrderC1 = 0;
178 else if (theRank == 2)
180 myC2 = (Standard_Address)&C;
181 switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
183 case GeomAbs_BezierCurve:
184 case GeomAbs_BSplineCurve:
185 case GeomAbs_OffsetCurve:
186 case GeomAbs_OtherCurve:
187 myMaxDerivOrderC2 = MaxOrder;
188 myTolC2 = SearchOfTolerance((Standard_Address)&C);
191 myMaxDerivOrderC2 = 0;
198 //=============================================================================
199 Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
203 Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
204 Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
206 Vec P1P2 (myP1,myP2);
208 Standard_Real Ndu = myDu.Magnitude();
211 if(myMaxDerivOrderC1 != 0)
215 //Derivative is approximated by Taylor-series
216 const Standard_Real DivisionFactor = 1.e-3;
218 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
221 du = myUsupremum-myUinfium;
223 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
225 Standard_Integer n = 1; //Derivative order
227 Standard_Boolean IsDeriveFound;
231 V = Tool1::DN(*((Curve1*)myC1),myU,++n);
233 IsDeriveFound = (Ndu > myTolC1);
235 while(!IsDeriveFound && n < myMaxDerivOrderC1);
241 if(myU-myUinfium < aDelta)
247 Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
248 Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
251 Standard_Real aDirFactor = V.Dot(V1);
260 //Derivative is approximated by three points
262 Pnt Ptemp; //(0,0,0)-coordinate
264 Standard_Boolean IsParameterGrown;
266 if(myU-myUinfium < 2*aDelta)
268 Tool1::D0(*((Curve1*)myC1),myU,P1);
269 Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
270 Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
271 IsParameterGrown = Standard_True;
275 Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
276 Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
277 Tool1::D0(*((Curve1*)myC1),myU,P3);
278 IsParameterGrown = Standard_False;
281 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
287 }//else of if(IsDeriveFound)
288 Ndu = myDu.Magnitude();
289 }//if (Ndu <= myTolC1) condition
290 }//if(myMaxDerivOrder != 0)
294 //Warning: 1st derivative of C1 is equal to zero!
295 return Standard_False;
298 Standard_Real Ndv = myDv.Magnitude();
300 if(myMaxDerivOrderC2 != 0)
304 const Standard_Real DivisionFactor = 1.e-3;
306 if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
309 dv = myVsupremum-myVinfium;
311 const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
313 //Derivative is approximated by Taylor-series
315 Standard_Integer n = 1; //Derivative order
317 Standard_Boolean IsDeriveFound;
321 V = Tool2::DN(*((Curve2*)myC2),myV,++n);
323 IsDeriveFound = (Ndv > myTolC2);
325 while(!IsDeriveFound && n < myMaxDerivOrderC2);
331 if(myV-myVinfium < aDelta)
337 Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
338 Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
341 Standard_Real aDirFactor = V.Dot(V1);
350 //Derivative is approximated by three points
352 Pnt Ptemp; //(0,0,0)-coordinate
354 Standard_Boolean IsParameterGrown;
356 if(myV-myVinfium < 2*aDelta)
358 Tool2::D0(*((Curve2*)myC2),myV,P1);
359 Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
360 Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
361 IsParameterGrown = Standard_True;
365 Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
366 Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
367 Tool2::D0(*((Curve2*)myC2),myV,P3);
368 IsParameterGrown = Standard_False;
371 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
377 }//else of if(IsDeriveFound)
379 Ndv = myDv.Magnitude();
380 }//if (Ndv <= myTolC2)
381 }//if(myMaxDerivOrder != 0)
385 //1st derivative of C2 is equal to zero!
386 return Standard_False;
389 F(1) = P1P2.Dot(myDu)/Ndu;
390 F(2) = P1P2.Dot(myDv)/Ndv;
391 return Standard_True;
393 //=============================================================================
395 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
399 return Values(UV,F,Df);
401 //=============================================================================
403 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV,
410 if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
412 //Warning: No function value found!
413 return Standard_False;
414 }//if(Value(UV, F) == Standard_False)
416 Vec Du, Dv, Duu, Dvv;
417 Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
418 Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
420 //Calling of "Value(...)" function change class member values.
421 //After running it is necessary to return to previous values.
422 const Standard_Real myU_old = myU, myV_old = myV;
423 const Pnt myP1_old = myP1, myP2_old = myP2;
424 const Vec myDu_old = myDu, myDv_old = myDv;
427 //Attention: aDelta value must be greater than same value for "Value(...)"
428 // function to avoid of points' collisions.
430 const Standard_Real DivisionFactor = 0.01;
433 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
436 du = myUsupremum-myUinfium;
438 const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
441 if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
444 dv = myVsupremum-myVinfium;
446 const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
448 Vec P1P2 (myP1,myP2);
450 if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
452 //Derivative is approximated by three points
454 math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
455 Standard_Real F1, F2, F3;
457 /////////////////////////// Search of DF1_u derivative (begin) ///////////////////
458 if(myU-myUinfium < 2*aDeltaU)
461 math_Vector UV2(1,2), UV3(1,2);
464 UV3(1)=myU+2*aDeltaU;
466 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
468 //There are many points close to singularity points and which have zero-derivative.
469 //Try to decrease aDelta variable's value.
470 return Standard_False;
476 Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
477 }//if(myU-myUinfium < 2*aDeltaU)
481 math_Vector UV2(1,2), UV1(1,2);
484 UV1(1)=myU-2*aDeltaU;
487 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
489 //There are many points close to singularity points and which have zero-derivative.
490 //Try to decrease aDelta variable's value.
491 return Standard_False;
497 Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
498 }//else of if(myU-myUinfium < 2*aDeltaU) condition
499 /////////////////////////// Search of DF1_u derivative (end) ///////////////////
501 //Return to previous values
505 /////////////////////////// Search of DF1_v derivative (begin) ///////////////////
506 if(myV-myVinfium < 2*aDeltaV)
509 math_Vector UV2(1,2), UV3(1,2);
513 UV3(2)=myV+2*aDeltaV;
515 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
517 //There are many points close to singularity points and which have zero-derivative.
518 //Try to decrease aDelta variable's value.
519 return Standard_False;
524 Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
525 }//if(myV-myVinfium < 2*aDeltaV)
529 math_Vector UV2(1,2), UV1(1,2);
533 UV1(2)=myV-2*aDeltaV;
534 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
536 //There are many points close to singularity points and which have zero-derivative.
537 //Try to decrease aDelta variable's value.
538 return Standard_False;
544 Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
545 }//else of if(myV-myVinfium < 2*aDeltaV)
546 /////////////////////////// Search of DF1_v derivative (end) ///////////////////
548 //Return to previous values
551 myP1 = myP1_old, myP2 = myP2_old;
552 myDu = myDu_old, myDv = myDv_old;
553 }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
556 const Standard_Real Ndu = myDu.Magnitude();
557 Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
558 Df(1,2) = myDv.Dot(myDu)/Ndu;
559 }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
561 if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
563 //Derivative is approximated by three points
565 math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
566 Standard_Real F1, F2, F3;
568 /////////////////////////// Search of DF2_v derivative (begin) ///////////////////
569 if(myV-myVinfium < 2*aDeltaV)
572 math_Vector UV2(1,2), UV3(1,2);
576 UV3(2)=myV+2*aDeltaV;
578 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
580 //There are many points close to singularity points and which have zero-derivative.
581 //Try to decrease aDelta variable's value.
582 return Standard_False;
588 Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
590 }//if(myV-myVinfium < 2*aDeltaV)
594 math_Vector UV2(1,2), UV1(1,2);
598 UV1(2)=myV-2*aDeltaV;
600 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
602 //There are many points close to singularity points and which have zero-derivative.
603 //Try to decrease aDelta variable's value.
604 return Standard_False;
610 Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
611 }//else of if(myV-myVinfium < 2*aDeltaV)
612 /////////////////////////// Search of DF2_v derivative (end) ///////////////////
614 //Return to previous values
618 /////////////////////////// Search of DF2_u derivative (begin) ///////////////////
619 if(myU-myUinfium < 2*aDeltaU)
622 math_Vector UV2(1,2), UV3(1,2);
625 UV3(1)=myU+2*aDeltaU;
627 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
629 //There are many points close to singularity points and which have zero-derivative.
630 //Try to decrease aDelta variable's value.
631 return Standard_False;
637 Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
638 }//if(myU-myUinfium < 2*aDelta)
642 math_Vector UV2(1,2), UV1(1,2);
645 UV1(1)=myU-2*aDeltaU;
648 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
650 //There are many points close to singularity points
651 //and which have zero-derivative.
652 //Try to decrease aDelta variable's value.
653 return Standard_False;
659 Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
660 }//else of if(myU-myUinfium < 2*aDeltaU)
661 /////////////////////////// Search of DF2_u derivative (end) ///////////////////
663 //Return to previous values
670 }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
673 Standard_Real Ndv = myDv.Magnitude();
674 Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
675 Df(2,1) = -myDu.Dot(myDv)/Ndv;
676 }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
678 return Standard_True;
681 //=============================================================================
683 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
685 Vec Du (myDu), Dv (myDv);
686 Vec P1P2 (myP1, myP2);
688 Standard_Real mod = Du.Magnitude();
692 mod = Dv.Magnitude();
697 if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) {
698 mySqDist.Append(myP1.SquareDistance(myP2));
699 myPoints.Append(POnC(myU, myP1));
700 myPoints.Append(POnC(myV, myP2));
704 //=============================================================================
706 void Extrema_FuncExtCC::Points (const Standard_Integer N,
710 P1 = myPoints.Value(2*N-1);
711 P2 = myPoints.Value(2*N);
713 //=============================================================================
715 void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
716 const math_Vector& theSupBound)
718 myUinfium = theInfBound(1);
719 myUsupremum = theSupBound(1);
720 myVinfium = theInfBound(2);
721 myVsupremum = theSupBound(2);
723 //=============================================================================