1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
19 //Modified by MPS : 21-05-97 PRO 7598
20 // Le nombre de solutions rendu est mySqDist.Length() et non
21 // mySqDist.Length()/2
22 // ajout des valeurs absolues dans le test d'orthogonalite de
25 /*-----------------------------------------------------------------------------
26 Fonctions permettant de rechercher une distance extremale entre 2 courbes
27 C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
28 Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
29 l'algorithme math_FunctionSetRoot.
30 Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
31 { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
32 { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
33 Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
35 { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
36 { Dvf1(u,v) = Dv.Du/||Du||
37 { Duf2(u,v) = -Du.Dv/||Dv||
38 { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
40 ----------------------------------------------------------------------------*/
42 #include <Precision.hxx>
45 static const Standard_Real MinTol = 1.e-20;
46 static const Standard_Real TolFactor = 1.e-12;
47 static const Standard_Real MinStep = 1e-7;
49 static const Standard_Integer MaxOrder = 3;
53 //=============================================================================
54 Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
56 const Standard_Integer NPoint = 10;
57 Standard_Real aStartParam, anEndParam;
61 aStartParam = myUinfium;
62 anEndParam = myUsupremum;
66 aStartParam = myVinfium;
67 anEndParam = myVsupremum;
72 cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl;
73 cout << "Warning: No curve for tolerance computing!---"<<endl;
78 const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
80 Standard_Integer aNum = 0;
81 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
82 //(it is computed with using NPoint point)
86 Standard_Real u = aStartParam + aNum*aStep; //parameter for every point
90 Pnt Ptemp; //empty point (is not used below)
91 Vec VDer; // 1st derivative vector
92 Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
93 Standard_Real vm = VDer.Magnitude();
97 while(++aNum < NPoint+1);
99 return Max(aMax*TolFactor,MinTol);
102 //=============================================================================
103 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
105 math_Vector V1(1,2), V2(1,2);
110 SubIntervalInitialize(V1, V2);
111 myMaxDerivOrderC1 = 0;
113 myMaxDerivOrderC2 = 0;
117 //=============================================================================
118 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
120 const Standard_Real thetol) :
121 myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
124 math_Vector V1(1,2), V2(1,2);
125 V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
126 V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
127 V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
128 V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
129 SubIntervalInitialize(V1, V2);
131 switch(Tool1::GetType(*((Curve1*)myC1)))
133 case GeomAbs_BezierCurve:
134 case GeomAbs_BSplineCurve:
135 case GeomAbs_OtherCurve:
136 myMaxDerivOrderC1 = MaxOrder;
137 myTolC1 = SearchOfTolerance((Standard_Address)&C1);
140 myMaxDerivOrderC1 = 0;
145 switch(Tool2::GetType(*((Curve2*)myC2)))
147 case GeomAbs_BezierCurve:
148 case GeomAbs_BSplineCurve:
149 case GeomAbs_OtherCurve:
150 myMaxDerivOrderC2 = MaxOrder;
151 myTolC2 = SearchOfTolerance((Standard_Address)&C2);
154 myMaxDerivOrderC2 = 0;
160 //=============================================================================
161 void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
163 Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
167 myC1 = (Standard_Address)&C;
168 switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
170 case GeomAbs_BezierCurve:
171 case GeomAbs_BSplineCurve:
172 case GeomAbs_OtherCurve:
173 myMaxDerivOrderC1 = MaxOrder;
174 myTolC1 = SearchOfTolerance((Standard_Address)&C);
177 myMaxDerivOrderC1 = 0;
182 else if (theRank == 2)
184 myC2 = (Standard_Address)&C;
185 switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
187 case GeomAbs_BezierCurve:
188 case GeomAbs_BSplineCurve:
189 case GeomAbs_OtherCurve:
190 myMaxDerivOrderC2 = MaxOrder;
191 myTolC2 = SearchOfTolerance((Standard_Address)&C);
194 myMaxDerivOrderC2 = 0;
201 //=============================================================================
202 Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
206 Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
207 Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
209 Vec P1P2 (myP1,myP2);
211 Standard_Real Ndu = myDu.Magnitude();
214 if(myMaxDerivOrderC1 != 0)
218 //Derivative is approximated by Taylor-series
219 const Standard_Real DivisionFactor = 1.e-3;
221 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
224 du = myUsupremum-myUinfium;
226 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
228 Standard_Integer n = 1; //Derivative order
230 Standard_Boolean IsDeriveFound;
234 V = Tool1::DN(*((Curve1*)myC1),myU,++n);
236 IsDeriveFound = (Ndu > myTolC1);
238 while(!IsDeriveFound && n < myMaxDerivOrderC1);
244 if(myU-myUinfium < aDelta)
250 Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
251 Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
254 Standard_Real aDirFactor = V.Dot(V1);
263 //Derivative is approximated by three points
265 Pnt Ptemp; //(0,0,0)-coordinate
267 Standard_Boolean IsParameterGrown;
269 if(myU-myUinfium < 2*aDelta)
271 Tool1::D0(*((Curve1*)myC1),myU,P1);
272 Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
273 Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
274 IsParameterGrown = Standard_True;
278 Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
279 Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
280 Tool1::D0(*((Curve1*)myC1),myU,P3);
281 IsParameterGrown = Standard_False;
284 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
290 }//else of if(IsDeriveFound)
291 Ndu = myDu.Magnitude();
292 }//if (Ndu <= myTolC1) condition
293 }//if(myMaxDerivOrder != 0)
298 cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
299 cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl;
301 return Standard_False;
305 Standard_Real Ndv = myDv.Magnitude();
307 if(myMaxDerivOrderC2 != 0)
311 const Standard_Real DivisionFactor = 1.e-3;
313 if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
316 dv = myVsupremum-myVinfium;
318 const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
320 //Derivative is approximated by Taylor-series
322 Standard_Integer n = 1; //Derivative order
324 Standard_Boolean IsDeriveFound;
328 V = Tool2::DN(*((Curve2*)myC2),myV,++n);
330 IsDeriveFound = (Ndv > myTolC2);
332 while(!IsDeriveFound && n < myMaxDerivOrderC2);
338 if(myV-myVinfium < aDelta)
344 Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
345 Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
348 Standard_Real aDirFactor = V.Dot(V1);
357 //Derivative is approximated by three points
359 Pnt Ptemp; //(0,0,0)-coordinate
361 Standard_Boolean IsParameterGrown;
363 if(myV-myVinfium < 2*aDelta)
365 Tool2::D0(*((Curve2*)myC2),myV,P1);
366 Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
367 Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
368 IsParameterGrown = Standard_True;
372 Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
373 Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
374 Tool2::D0(*((Curve2*)myC2),myV,P3);
375 IsParameterGrown = Standard_False;
378 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
384 }//else of if(IsDeriveFound)
386 Ndv = myDv.Magnitude();
387 }//if (Ndv <= myTolC2)
388 }//if(myMaxDerivOrder != 0)
393 cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
394 cout << "1st derivative of C2 is equal to zero!---"<<endl;
396 return Standard_False;
399 F(1) = P1P2.Dot(myDu)/Ndu;
400 F(2) = P1P2.Dot(myDv)/Ndv;
401 return Standard_True;
403 //=============================================================================
405 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
409 return Values(UV,F,Df);
411 //=============================================================================
413 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV,
420 if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
423 cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl;
424 cout << "Warning: No function value found!---"<<endl;
426 return Standard_False;
427 }//if(Value(UV, F) == Standard_False)
429 Vec Du, Dv, Duu, Dvv;
430 Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
431 Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
433 //Calling of "Value(...)" function change class member values.
434 //After running it is necessary to return to previous values.
435 const Standard_Real myU_old = myU, myV_old = myV;
436 const Pnt myP1_old = myP1, myP2_old = myP2;
437 const Vec myDu_old = myDu, myDv_old = myDv;
440 //Attention: aDelta value must be greater than same value for "Value(...)"
441 // function to avoid of points' collisions.
443 const Standard_Real DivisionFactor = 0.01;
446 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
449 du = myUsupremum-myUinfium;
451 const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
454 if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
457 dv = myVsupremum-myVinfium;
459 const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
461 Vec P1P2 (myP1,myP2);
463 if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
465 //Derivative is approximated by three points
467 math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
468 Standard_Real F1, F2, F3;
470 /////////////////////////// Search of DF1_u derivative (begin) ///////////////////
471 if(myU-myUinfium < 2*aDeltaU)
474 math_Vector UV2(1,2), UV3(1,2);
477 UV3(1)=myU+2*aDeltaU;
479 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
482 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
483 cout << "There are many points close to singularity points "
484 "and which have zero-derivative." << endl;
485 cout << "Try to decrease aDelta variable's value. ---" << endl;
487 return Standard_False;
493 Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
494 }//if(myU-myUinfium < 2*aDeltaU)
498 math_Vector UV2(1,2), UV1(1,2);
501 UV1(1)=myU-2*aDeltaU;
504 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
507 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
508 cout << "There are many points close to singularity points "
509 "and which have zero-derivative." << endl;
510 cout << "Try to decrease aDelta variable's value. ---" << endl;
512 return Standard_False;
518 Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
519 }//else of if(myU-myUinfium < 2*aDeltaU) condition
520 /////////////////////////// Search of DF1_u derivative (end) ///////////////////
522 //Return to previous values
526 /////////////////////////// Search of DF1_v derivative (begin) ///////////////////
527 if(myV-myVinfium < 2*aDeltaV)
530 math_Vector UV2(1,2), UV3(1,2);
534 UV3(2)=myV+2*aDeltaV;
536 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
539 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
540 cout << "There are many points close to singularity points "
541 "and which have zero-derivative." << endl;
542 cout << "Try to decrease aDelta variable's value. ---" << endl;
544 return Standard_False;
549 Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
550 }//if(myV-myVinfium < 2*aDeltaV)
554 math_Vector UV2(1,2), UV1(1,2);
558 UV1(2)=myV-2*aDeltaV;
559 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
562 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
563 cout << "There are many points close to singularity points "
564 "and which have zero-derivative." << endl;
565 cout << "Try to decrease aDelta variable's value. ---" << endl;
567 return Standard_False;
573 Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
574 }//else of if(myV-myVinfium < 2*aDeltaV)
575 /////////////////////////// Search of DF1_v derivative (end) ///////////////////
577 //Return to previous values
580 myP1 = myP1_old, myP2 = myP2_old;
581 myDu = myDu_old, myDv = myDv_old;
582 }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
585 const Standard_Real Ndu = myDu.Magnitude();
586 Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
587 Df(1,2) = myDv.Dot(myDu)/Ndu;
588 }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
590 if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
592 //Derivative is approximated by three points
594 math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
595 Standard_Real F1, F2, F3;
597 /////////////////////////// Search of DF2_v derivative (begin) ///////////////////
598 if(myV-myVinfium < 2*aDeltaV)
601 math_Vector UV2(1,2), UV3(1,2);
605 UV3(2)=myV+2*aDeltaV;
607 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
610 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
611 cout << "There are many points close to singularity points "
612 "and which have zero-derivative." << endl;
613 cout << "Try to decrease aDelta variable's value. ---" << endl;
615 return Standard_False;
621 Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
623 }//if(myV-myVinfium < 2*aDeltaV)
627 math_Vector UV2(1,2), UV1(1,2);
631 UV1(2)=myV-2*aDeltaV;
633 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
636 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
637 cout << "There are many points close to singularity points "
638 "and which have zero-derivative." << endl;
639 cout << "Try to decrease aDelta variable's value. ---" << endl;
641 return Standard_False;
647 Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
648 }//else of if(myV-myVinfium < 2*aDeltaV)
649 /////////////////////////// Search of DF2_v derivative (end) ///////////////////
651 //Return to previous values
655 /////////////////////////// Search of DF2_u derivative (begin) ///////////////////
656 if(myU-myUinfium < 2*aDeltaU)
659 math_Vector UV2(1,2), UV3(1,2);
662 UV3(1)=myU+2*aDeltaU;
664 if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
667 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
668 cout << "There are many points close to singularity points "
669 "and which have zero-derivative." << endl;
670 cout << "Try to decrease aDelta variable's value. ---" << endl;
672 return Standard_False;
678 Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
679 }//if(myU-myUinfium < 2*aDelta)
683 math_Vector UV2(1,2), UV1(1,2);
686 UV1(1)=myU-2*aDeltaU;
689 if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
692 cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
693 cout << "There are many points close to singularity points "
694 "and which have zero-derivative." << endl;
695 cout << "Try to decrease aDelta variable's value. ---" << endl;
697 return Standard_False;
703 Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
704 }//else of if(myU-myUinfium < 2*aDeltaU)
705 /////////////////////////// Search of DF2_u derivative (end) ///////////////////
707 //Return to previous values
714 }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
717 Standard_Real Ndv = myDv.Magnitude();
718 Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
719 Df(2,1) = -myDu.Dot(myDv)/Ndv;
720 }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
722 return Standard_True;
725 //=============================================================================
727 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
729 Vec Du (myDu), Dv (myDv);
730 Vec P1P2 (myP1, myP2);
732 Standard_Real mod = Du.Magnitude();
736 mod = Dv.Magnitude();
741 if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) {
742 mySqDist.Append(myP1.SquareDistance(myP2));
743 myPoints.Append(POnC(myU, myP1));
744 myPoints.Append(POnC(myV, myP2));
748 //=============================================================================
750 void Extrema_FuncExtCC::Points (const Standard_Integer N,
754 P1 = myPoints.Value(2*N-1);
755 P2 = myPoints.Value(2*N);
757 //=============================================================================
759 void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
760 const math_Vector& theSupBound)
762 myUinfium = theInfBound(1);
763 myUsupremum = theSupBound(1);
764 myVinfium = theInfBound(2);
765 myVsupremum = theSupBound(2);
767 //=============================================================================