0024921: ShapeAnalysis_Curve::ValidateRange doesn't adjust the range for periodic...
[occt.git] / src / ShapeAnalysis / ShapeAnalysis_Curve.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 3030cf9..8b83b32
@@ -1,3 +1,16 @@
+// 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.
+
 // pdn 04.12.98 Add method using Adaptor_Curve
 //:j8 abv 10.12.98 TR10 r0501_db.stp #9423
 //pdn 25.12.98 private method ProjectAct
@@ -38,7 +51,7 @@
 #include <Geom2dInt_Geom2dCurveTool.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
 #include <Geom_Circle.hxx>
-
+#include <Extrema_LocateExtPC.hxx>
 
 //=======================================================================
 //function : ProjectOnSegments
@@ -53,14 +66,17 @@ static void ProjectOnSegments (const Adaptor3d_Curve& AC, const gp_Pnt& P3D,
   //  On considere <nbseg> points sur [uMin,uMax]
   //  Quel est le plus proche. Et quel est le nouvel intervalle
   //  (il ne peut pas deborder l ancien)
-  Standard_Real u, dist, delta = (nbseg == 0)? 0 : (uMax-uMin)/nbseg; //szv#4:S4163:12Mar99 anti-exception
+  Standard_Real u, dist2, delta = (nbseg == 0)? 0 : (uMax-uMin)/nbseg; //szv#4:S4163:12Mar99 anti-exception
+  Standard_Real  distmin2 = distmin * distmin;
+  Standard_Boolean aHasChanged = Standard_False;
   for (Standard_Integer i = 0; i <= nbseg; i ++) {
     u = uMin + (delta * i);
     gp_Pnt PU = AC.Value (u);
-    dist = PU.Distance (P3D);
-    if (dist < distmin)  {  distmin = dist;  proj = PU;  param = u;  }
+    dist2 = PU.SquareDistance (P3D);
+    if (dist2 < distmin2)  {  distmin2 = dist2;  proj = PU;  param = u; aHasChanged = Standard_True;  }
   }
-
+  if (aHasChanged)
+    distmin = Sqrt (distmin2);
 #ifdef DEBUG
   cout<<"ShapeAnalysis_Geom:Project, param="<<param<<" -> distmin="<<distmin<<endl;
 #endif
@@ -70,58 +86,6 @@ static void ProjectOnSegments (const Adaptor3d_Curve& AC, const gp_Pnt& P3D,
 }
 
 
-//=======================================================================
-//function : CurveNewton
-//purpose  : Newton algo on curve (study S4030)
-//=======================================================================
-
-static Standard_Boolean CurveNewton(const Standard_Real paramPrev,
-                                   const Adaptor3d_Curve& Ad,
-                                   const gp_Pnt& P3D,
-                                   const Standard_Real /*preci*/,
-                                   Standard_Real& param,
-                                   const Standard_Real cf,
-                                   const Standard_Real cl)
-{
-  //szv#4:S4163:12Mar99 optimized
-  Standard_Real uMin = cf, uMax = cl;
-
-  Standard_Real Tol = Precision::Confusion();
-  Standard_Real Tol2 = Tol*Tol, rs2p = 1.e+10;
-  Standard_Real X = paramPrev;
-
-  for ( Standard_Integer i=0; i<20; i++) {
-    gp_Vec v1, v2;
-    gp_Pnt pnt;
-    Ad.D2 (X, pnt, v1, v2);
-    Standard_Real nv1 = v1.SquareMagnitude();
-    if (nv1 < 1e-10) break;
-    
-    gp_Vec rs (P3D, pnt);
-    Standard_Real D = nv1 + (rs * v2);
-    if ( fabs( D ) <1e-10) break;
-    
-    Standard_Real dX = -(rs *v1)/D;
-    X += dX;
-    
-    Standard_Real rs2 = rs.SquareMagnitude();
-    if (rs2 > 4.*rs2p) break;
-    rs2p = rs2;
-    
-    if ( fabs(dX) > 1e-12) continue;
-     
-    if (X < uMin || X > uMax) break;
-    
-    Standard_Real rsn = rs * v1;
-    if (rsn*rsn/nv1 > Tol2) break;
-    param = X;
-    return Standard_True;
-  }
-
-  // cout <<"NewtonC: failed after " << i+1 << " iterations (fail " << fail << " )" << endl; 
-  return Standard_False;
-}
-
 //=======================================================================
 //function : Project
 //purpose  : 
@@ -208,27 +172,46 @@ Standard_Real ShapeAnalysis_Curve::Project(const Adaptor3d_Curve& C3D,
                                           const Standard_Boolean AdjustToEnds) const
 
 {
+
   Standard_Real uMin = C3D.FirstParameter();
   Standard_Real uMax = C3D.LastParameter();
-  Standard_Real distmin;
-  if (!Precision::IsInfinite(uMin)||!Precision::IsInfinite(uMax)) {
-    Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
-    gp_Pnt LowBound = C3D.Value(uMin);
-    gp_Pnt HigBound = C3D.Value(uMax);
-    distmin = LowBound.Distance(P3D);
-    if (distmin <= prec) {
-      param = uMin;
-      proj  = LowBound;
-      return distmin;
-    }
-    distmin = HigBound.Distance(P3D);
-    if (distmin <= prec) {
-      param = uMax;
-      proj  = HigBound;
-      return distmin;
-    } 
+  
+  if (Precision::IsInfinite(uMin) && Precision::IsInfinite(uMax))
+    return ProjectAct(C3D, P3D, preci, proj, param);
+
+  Standard_Real distmin_L = Precision::Infinite(), distmin_H = Precision::Infinite();
+  Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
+  gp_Pnt LowBound = C3D.Value(uMin);
+  gp_Pnt HigBound = C3D.Value(uMax);
+  distmin_L = LowBound.Distance(P3D);
+  distmin_H = HigBound.Distance(P3D);
+
+  if (distmin_L <= prec) {
+    param = uMin;
+    proj  = LowBound;
+    return distmin_L;
   }
-  return ProjectAct(C3D, P3D, preci, proj, param);
+
+  if (distmin_H <= prec) {
+    param = uMax;
+    proj  = HigBound;
+    return distmin_H;
+  } 
+
+  Standard_Real distProj = ProjectAct(C3D, P3D, preci, proj, param);
+  if(  distProj < distmin_L +  Precision::Confusion() && distProj < distmin_H +  Precision::Confusion())
+    return distProj;
+
+  if( distmin_L < distmin_H)
+  {
+    param = uMin;
+    proj  = LowBound;
+    return distmin_L;
+  }
+  param = uMax;
+  proj  = HigBound;
+  return distmin_H;
+
 }
 
 //=======================================================================
@@ -243,56 +226,35 @@ Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
                                              Standard_Real& param) const
        
 {
-  Standard_Boolean useExtrema = Standard_True; //:i5
-
-#ifdef WNT
-  //:i5 abv 11 Sep 98: UKI60195.igs on NT: hanging on null-length bsplines
-  if (C3D.IsClosed() && C3D.GetType()==GeomAbs_BSplineCurve) {
-    Handle(Geom_BSplineCurve) bspl = C3D.BSpline();
-    Standard_Real prec2 = gp::Resolution();
-    prec2 *= prec2;
-    gp_Pnt p0 = bspl->Pole(1), p1;
-    for ( Standard_Integer i=2; i <= bspl->NbPoles(); i++, p0 = p1 ) {
-      p1 = bspl->Pole(i);
-      if ( p0.SquareDistance ( p1 ) <= prec2 ) { 
-       useExtrema = Standard_False;
-       break;
-      }
-    }
-  }
-#endif
-
   Standard_Boolean OK = Standard_False;
-  if ( useExtrema ) { //:i5 //szv#4:S4163:12Mar99 try moved into if
-    try {
-      OCC_CATCH_SIGNALS
-      Extrema_ExtPC myExtPC(P3D,C3D);
-      //Standard_Boolean done = myExtPC.IsDone() && ( myExtPC.NbExt() > 0); //szv#4:S4163:12Mar99 not needed
-      if ( myExtPC.IsDone() && ( myExtPC.NbExt() > 0) ) {
-       Standard_Real dist2, dist2Min = myExtPC.SquareDistance(1);
-       Standard_Integer index = 1;
-       for ( Standard_Integer i = 2; i <= myExtPC.NbExt(); i++) {
-         dist2 = myExtPC.SquareDistance(i);
-         if ( dist2 < dist2Min) { dist2Min = dist2; index = i; }
-       }
-       param = (myExtPC.Point(index)).Parameter();
-       proj  = (myExtPC.Point(index)).Value();
-       OK = Standard_True;
+  param = 0.;
+  try {
+    OCC_CATCH_SIGNALS
+    Extrema_ExtPC myExtPC(P3D,C3D);
+    if ( myExtPC.IsDone() && ( myExtPC.NbExt() > 0) ) {
+      Standard_Real dist2, dist2Min = myExtPC.SquareDistance(1);
+      Standard_Integer index = 1;
+      for ( Standard_Integer i = 2; i <= myExtPC.NbExt(); i++) {
+        dist2 = myExtPC.SquareDistance(i);
+        if ( dist2 < dist2Min) { dist2Min = dist2; index = i; }
       }
+      param = (myExtPC.Point(index)).Parameter();
+      proj  = (myExtPC.Point(index)).Value();
+      OK = Standard_True;
     }
-    catch(Standard_Failure) {
-      OK = Standard_False;
+  }
+  catch(Standard_Failure) {
+    OK = Standard_False;
 #ifdef DEB //:s5
-      cout << "\nWarning: ShapeAnalysis_Curve::ProjectAct(): Exception in Extrema_ExtPC: "; 
-      Standard_Failure::Caught()->Print(cout); cout << endl;
+    cout << "\nWarning: ShapeAnalysis_Curve::ProjectAct(): Exception in Extrema_ExtPC: "; 
+    Standard_Failure::Caught()->Print(cout); cout << endl;
 #endif
-    }
   }
-
+  
   //szv#4:S4163:12Mar99 moved
   Standard_Real uMin = C3D.FirstParameter(), uMax = C3D.LastParameter();
   Standard_Boolean closed = Standard_False;  // si on franchit les bornes ...
-  Standard_Real distmin = RealLast(), valclosed = 0.;
+  Standard_Real distmin = Precision::Infinite(), valclosed = 0.;
   Standard_Real aModParam = param;
   Standard_Real aModMin = distmin;
   
@@ -338,7 +300,7 @@ Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
           proj  = ElCLib::Value(param, aCirc);
         }
        closed = Standard_True;
-       valclosed = 2.*PI;
+       valclosed = 2.*M_PI;
       }
       break;
     case GeomAbs_Hyperbola:
@@ -364,7 +326,7 @@ Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
        param = ElCLib::Parameter(C3D.Ellipse(), P3D);
        proj  = ElCLib::Value(param, C3D.Ellipse());
        closed = Standard_True;
-       valclosed = 2.*PI;
+       valclosed = 2.*M_PI;
 
       }
       break;
@@ -374,13 +336,17 @@ Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
         //  on tente ceci : 21 points sur la courbe, quel est le plus proche
        distmin = Precision::Infinite();
        ProjectOnSegments (C3D,P3D,25,uMin,uMax,distmin,proj,param);
-       if (distmin <= preci) return distmin;
-       if (CurveNewton(param, C3D, P3D, preci, param, uMin, uMax)) { //:S4030
-         C3D.D0(param, proj);
+       if (distmin <= preci) 
+          return distmin;
+        Extrema_LocateExtPC aProjector (P3D, C3D, param/*U0*/, uMin, uMax, preci/*TolU*/);
+        if (aProjector.IsDone())
+        {
+          param = aProjector.Point().Parameter();
+          proj = aProjector.Point().Value();
           Standard_Real aDistNewton = P3D.Distance(proj);
-          if(aDistNewton < aModMin)
+          if (aDistNewton < aModMin)
             return aDistNewton;
-       }
+        }
         // cout <<"newton failed"<<endl;    
        ProjectOnSegments (C3D,P3D,40,uMin,uMax,distmin,proj,param);
        if (distmin <= preci) return distmin;
@@ -436,7 +402,7 @@ Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
 {
   Standard_Real uMin = (cf < cl ? cf : cl);
   Standard_Real uMax = (cf < cl ? cl : cf);
-  Standard_Real distmin;
+  Standard_Real distmin = Precision::Infinite();
   if (C3D->IsKind(STANDARD_TYPE(Geom_BoundedCurve))) {
     Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
     gp_Pnt LowBound = C3D->Value(uMin);
@@ -488,12 +454,13 @@ Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
   Standard_Real uMin = C3D.FirstParameter();
   Standard_Real uMax = C3D.LastParameter();
   
-  if (CurveNewton(paramPrev, C3D, P3D, preci, param, uMin, uMax)) { 
-    C3D.D0(param, proj);
+  Extrema_LocateExtPC aProjector (P3D, C3D, paramPrev/*U0*/, uMin, uMax, preci/*TolU*/);
+  if (aProjector.IsDone()){
+    param = aProjector.Point().Parameter();
+    proj = aProjector.Point().Value();
     return P3D.Distance(proj);
   }
-
-  return Project(C3D, P3D, preci, proj, param);
+  return Project(C3D, P3D, preci, proj, param, Standard_False);
 }
 
 //=======================================================================
@@ -539,11 +506,13 @@ Standard_Boolean ShapeAnalysis_Curve::ValidateRange (const Handle(Geom_Curve)& t
     }
   }
 
-  if (First < Last) return Standard_True;
-
   // 15.11.2002 PTV OCC966
-  if (ShapeAnalysis_Curve::IsPeriodic(theCurve)) 
+  if (ShapeAnalysis_Curve::IsPeriodic(theCurve)) {
     ElCLib::AdjustPeriodic(cf,cl,Precision::PConfusion(),First,Last); //:a7 abv 11 Feb 98: preci -> PConfusion()
+  }
+  else if (First < Last) {
+    // nothing to fix
+  }
   else if (theCurve->IsClosed()) {
     // l'un des points projecte se trouve sur l'origine du parametrage
     // de la courbe 3D. L algo a donne cl +- preci au lieu de cf ou vice-versa
@@ -802,8 +771,8 @@ static void AppendControlPoles (TColgp_SequenceOfPnt& seq,
     seq.Append(curve->Value(1));
   } else if ( curve->IsKind(STANDARD_TYPE(Geom_Conic))) {
     seq.Append(curve->Value(0));
-    seq.Append(curve->Value(PI/2));
-    seq.Append(curve->Value(PI));
+    seq.Append(curve->Value(M_PI/2));
+    seq.Append(curve->Value(M_PI));
   } else if ( curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
     //DeclareAndCast(Geom_TrimmedCurve, Trimmed, curve);
     Handle(Geom_TrimmedCurve) Trimmed = *((Handle(Geom_TrimmedCurve) *) &curve);
@@ -963,7 +932,7 @@ Standard_Boolean ShapeAnalysis_Curve::IsPlanar (const TColgp_Array1OfPnt& pnts,
       return Standard_True;
     }
     gp_XYZ aVecMul = N1^Normal;
-    return aVecMul.SquareModulus() < Precision::Confusion()*Precision::Confusion();
+    return aVecMul.SquareModulus() < Precision::SquareConfusion();
   }
 
   if (curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
@@ -1087,7 +1056,7 @@ Standard_Boolean ShapeAnalysis_Curve::GetSamplePoints (const Handle(Geom2d_Curve
     return Standard_True;
   }
   else if(curve->IsKind(STANDARD_TYPE(Geom2d_Conic))) {
-    step = Min ( PI, last-first ) / 19; //:abv 05.06.02 TUBE.stp #19209...: PI/16
+    step = Min ( M_PI, last-first ) / 19; //:abv 05.06.02 TUBE.stp #19209...: M_PI/16
 //     if( step>(last-first) ) {
 //       seq.Append(curve->Value(first));
 //       seq.Append(curve->Value((last+first)/2));