0030199: Extrema Point<->Curve gives inaccurate result
authorifv <ifv@opencascade.com>
Mon, 15 Oct 2018 12:38:00 +0000 (15:38 +0300)
committerapn <apn@opencascade.com>
Tue, 23 Oct 2018 17:08:38 +0000 (20:08 +0300)
Special treatment of bspline curve of first degree is implemented in Extrema_GExtPC.gxx
Test case is added
Some test cases are modified according to actual state of algorithm

src/Extrema/Extrema_GExtPC.gxx
tests/bugs/modalg_7/bug30199 [new file with mode: 0644]
tests/de/step_1/I5
tests/de/step_1/R1
tests/de/step_2/T9
tests/de/step_2/U8

index 3e60b75..e9f7610 100644 (file)
@@ -154,108 +154,189 @@ void Extrema_GExtPC::Perform(const ThePoint& P)
 
       mysample = (TheCurveTool::BSpline(aCurve))->Degree() + 1;
 
-      // Fill sample points.
-      Standard_Integer aValIdx = 1;
-      NCollection_Array1<Standard_Real> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
-      NCollection_Array1<Standard_Real> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
-      for(anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
+      if (mysample == 2)
       {
-        Standard_Real aF = aKnots(anIdx) + aPeriodJump,
-                      aL = aKnots(anIdx + 1) + aPeriodJump;
+        //BSpline of first degree, direct seaching extrema for each knot interval
+        ThePoint aPmin;
+        Standard_Real tmin = 0., distmin = RealLast();
+        Standard_Real aMin1 = 0., aMin2 = 0.;
+        myExtPC.Initialize(aCurve);
+        for (anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
+        {
+          Standard_Real aF = aKnots(anIdx) + aPeriodJump,
+            aL = aKnots(anIdx + 1) + aPeriodJump;
+
+          if (anIdx == aFirstUsedKnot)
+          {
+            aF = myuinf;
+          }
+          else
+          if (anIdx == aLastUsedKnot - 1)
+          {
+            aL = myusup;
+          }
 
-        if (anIdx == aFirstUsedKnot)
-          aF = myuinf;
-        if (anIdx == aLastUsedKnot - 1)
-          aL = myusup;
 
-        Standard_Real aStep = (aL - aF) / mysample;
-        for(Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
+          ThePoint aP1, aP2;
+          TheCurveTool::D0(aCurve, aF, aP1);
+          TheCurveTool::D0(aCurve, aL, aP2);
+          TheVector aBase1(P, aP1), aBase2(P, aP2);
+          TheVector aV(aP2, aP1);
+          Standard_Real aVal1 = aV.Dot(aBase1); // Derivative of (C(u) - P)^2
+          Standard_Real aVal2 = aV.Dot(aBase2); // Derivative of (C(u) - P)^2
+          if (anIdx == aFirstUsedKnot)
+          {
+            aMin1 = P.SquareDistance(aP1);
+          }
+          else
+          {
+            aMin1 = aMin2;          
+            if (distmin > aMin1)
+            {
+              distmin = aMin1;
+              tmin = aF;
+              aPmin = aP1;
+            }
+          }
+          aMin2 = P.SquareDistance(aP2);
+          Standard_Real aMinSqDist = Min(aMin1, aMin2);
+          Standard_Real aMinDer = Min(Abs(aVal1), Abs(aVal2));
+          if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
+          {
+            // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
+            if (aVal1 * aVal2 <= 0.0 ||
+                aMinSqDist < 100. * Precision::SquareConfusion() ||
+                2.*aMinDer < Precision::Confusion())
+            {
+              myintuinf = aF;
+              myintusup = aL;
+              IntervalPerform(P);
+            }
+          }
+        }
+        if (!Precision::IsInfinite(distmin))
         {
-          Standard_Real aCurrentParam = aF + aStep * aPntIdx;
-          aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
-          aParam(aValIdx) = aCurrentParam;
-          aValIdx++;
+          Standard_Boolean isToAdd = Standard_True;
+          NbExt = mypoint.Length();
+          for (i = 1; i <= NbExt && isToAdd; i++)
+          {
+            Standard_Real t = mypoint.Value(i).Parameter();
+            isToAdd = (distmin < mySqDist(i)) && (Abs(t - tmin) > mytolu);
+          }
+          if (isToAdd)
+          {
+            ThePOnC PC(tmin, aPmin);
+            mySqDist.Append(distmin);
+            myismin.Append(Standard_True);
+            mypoint.Append(PC);
+          }
         }
       }
-      // Fill last point.
-      aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
-      aParam(aValIdx) = myusup;
-
-       myExtPC.Initialize(aCurve);
-
-      // Find extremas.
-      for(anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
+      else
       {
-        if (aVal(anIdx) <= Precision::SquareConfusion())
+
+        // Fill sample points.
+        Standard_Integer aValIdx = 1;
+        NCollection_Array1<Standard_Real> aVal(1, (mysample)* (aLastUsedKnot - aFirstUsedKnot) + 1);
+        NCollection_Array1<Standard_Real> aParam(1, (mysample)* (aLastUsedKnot - aFirstUsedKnot) + 1);
+        for (anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
         {
-          mySqDist.Append(aVal(anIdx));
-          myismin.Append(Standard_True);
-          mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
+          Standard_Real aF = aKnots(anIdx) + aPeriodJump,
+            aL = aKnots(anIdx + 1) + aPeriodJump;
+
+          if (anIdx == aFirstUsedKnot)
+            aF = myuinf;
+          if (anIdx == aLastUsedKnot - 1)
+            aL = myusup;
+
+          Standard_Real aStep = (aL - aF) / mysample;
+          for (Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
+          {
+            Standard_Real aCurrentParam = aF + aStep * aPntIdx;
+            aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
+            aParam(aValIdx) = aCurrentParam;
+            aValIdx++;
+          }
         }
-        if ((aVal(anIdx) >= aVal(anIdx + 1) &&
-             aVal(anIdx) >= aVal(anIdx - 1)) ||
-            (aVal(anIdx) <= aVal(anIdx + 1) &&
-             aVal(anIdx) <= aVal(anIdx - 1)) )
+        // Fill last point.
+        aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
+        aParam(aValIdx) = myusup;
+
+        myExtPC.Initialize(aCurve);
+
+        // Find extremas.
+        for (anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
         {
-          myintuinf = aParam(anIdx - 1);
-          myintusup = aParam(anIdx + 1);
+          if (aVal(anIdx) <= Precision::SquareConfusion())
+          {
+            mySqDist.Append(aVal(anIdx));
+            myismin.Append(Standard_True);
+            mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
+          }
+          if ((aVal(anIdx) >= aVal(anIdx + 1) &&
+            aVal(anIdx) >= aVal(anIdx - 1)) ||
+            (aVal(anIdx) <= aVal(anIdx + 1) &&
+            aVal(anIdx) <= aVal(anIdx - 1)))
+          {
+            myintuinf = aParam(anIdx - 1);
+            myintusup = aParam(anIdx + 1);
 
-          IntervalPerform(P);
+            IntervalPerform(P);
+          }
         }
-      }
 
-      // Solve on the first and last intervals.
-      if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
-      {
-        ThePoint aP1, aP2;
-        TheVector aV1, aV2;
-        TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()),     aP1, aV1);
-        TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
-        TheVector aBase1(P, aP1), aBase2(P, aP2);
-        Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
-        Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
-        if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
+        // Solve on the first and last intervals.
+        if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
         {
-          // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
-          // Necessary condition - when point lies on curve.
-          // Necessary condition - when derivative of point is too small.
-          if(aVal1 * aVal2 <= 0.0 ||
-            aBase1.Dot(aBase2) <= 0.0 ||
-            2.0 * Abs(aVal1) < Precision::Confusion() )
+          ThePoint aP1, aP2;
+          TheVector aV1, aV2;
+          TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()), aP1, aV1);
+          TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
+          TheVector aBase1(P, aP1), aBase2(P, aP2);
+          Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
+          Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
+          if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
           {
-            myintuinf = aParam(aVal.Lower());
-            myintusup = aParam(aVal.Lower() + 1);
-            IntervalPerform(P);
+            // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
+            // Necessary condition - when point lies on curve.
+            // Necessary condition - when derivative of point is too small.
+            if (aVal1 * aVal2 <= 0.0 ||
+              aBase1.Dot(aBase2) <= 0.0 ||
+              2.0 * Abs(aVal1) < Precision::Confusion())
+            {
+              myintuinf = aParam(aVal.Lower());
+              myintusup = aParam(aVal.Lower() + 1);
+              IntervalPerform(P);
+            }
           }
         }
-      }
 
-      if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
-      {
-        ThePoint aP1, aP2;
-        TheVector aV1, aV2;
-        TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
-        TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()),     aP2, aV2);
-        TheVector aBase1(P, aP1), aBase2(P, aP2);
-        Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
-        Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
-
-        if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
+        if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
         {
-          // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
-          // Necessary condition - when point lies on curve.
-          // Necessary condition - when derivative of point is too small.
-          if(aVal1 * aVal2 <= 0.0 ||
-            aBase1.Dot(aBase2) <= 0.0 ||
-            2.0 * Abs(aVal2) < Precision::Confusion() )
+          ThePoint aP1, aP2;
+          TheVector aV1, aV2;
+          TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
+          TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()), aP2, aV2);
+          TheVector aBase1(P, aP1), aBase2(P, aP2);
+          Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
+          Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
+
+          if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
           {
-            myintuinf = aParam(aVal.Upper() - 1);
-            myintusup = aParam(aVal.Upper());
-            IntervalPerform(P);
+            // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
+            // Necessary condition - when point lies on curve.
+            // Necessary condition - when derivative of point is too small.
+            if (aVal1 * aVal2 <= 0.0 ||
+              aBase1.Dot(aBase2) <= 0.0 ||
+              2.0 * Abs(aVal2) < Precision::Confusion())
+            {
+              myintuinf = aParam(aVal.Upper() - 1);
+              myintusup = aParam(aVal.Upper());
+              IntervalPerform(P);
+            }
           }
         }
       }
-
       mydone = Standard_True;
       break;
     }
diff --git a/tests/bugs/modalg_7/bug30199 b/tests/bugs/modalg_7/bug30199
new file mode 100644 (file)
index 0000000..539e750
--- /dev/null
@@ -0,0 +1,16 @@
+puts "========"
+puts "0030199: Extrema Point<->Curve gives inaccurate result"
+puts "========"
+puts ""
+2dbsplinecurve c2d 1 8 0 2 147.948545311472 1 227.205037349287 1 282.782457537237 1 341.526569446034 1 412.981618615062 1 493.370506013535 1 820.162497637298  2 70.4310501764176 -125.064069601289 1 202.351851590745 -192.039245703947  1 281.504332439341 -187.980144121968  1 310.594560443526 -140.62395898876  1 274.062646205713 -94.6208077364441  1 341.714339238701 -71.6192321052281  1 397.188727525752  -129.799688113598  1 70.4310501764176  -125.064069601289 1
+to3d c3d c2d
+mkedge e c3d
+vertex v 360 -127 0
+distmini dd e v
+set dist [dval dd_val]
+checkreal dist $dist 2.26048366458175 1.e-7 0
+
+vertex v1 199.49416090801449 -197.77438365880749 0
+distmini dd1 e v1
+set dist1 [dval dd1_val]
+checkreal dist1 $dist1 6.4076675475126184 1.e-7 0
index 5ef5631..d640c03 100644 (file)
@@ -10,7 +10,7 @@ TPSTAT      : Faulties = 0  ( 0 )  Warnings = 14  ( 3 )  Summary  = 14  ( 3 )
 CHECKSHAPE  : Wires    = 0  ( 0 )  Faces    = 0  ( 0 )  Shells   = 0  ( 0 )   Solids   = 0 ( 0 )
 NBSHAPES    : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 41  ( 41 ) 
 STATSHAPE   : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 41  ( 41 )   FreeWire = 0  ( 4 ) 
-TOLERANCE   : MaxTol   =  0.07980045861  (  0.07980045861 )  AvgTol   =  0.006103073943  (  0.006647359845 )
+TOLERANCE   : MaxTol   =  0.07980045861  (   0.0832075472 )  AvgTol   =  0.006103073943  (  0.007097682834 )
 LABELS      : N0Labels = 3  ( 3 )  N1Labels = 2  ( 2 )  N2Labels = 0  ( 0 )   TotalLabels = 5  ( 5 )   NameLabels = 5  ( 5 )   ColorLabels = 0  ( 0 )   LayerLabels = 0  ( 0 )
 PROPS       : Centroid = 0  ( 0 )  Volume   = 0  ( 0 )  Area     = 0  ( 0 )
 NCOLORS     : NColors  = 0  ( 0 )
index 1914394..97a620e 100644 (file)
@@ -7,7 +7,7 @@ TPSTAT      : Faulties = 0  ( 0 )  Warnings = 0  ( 0 )  Summary  = 0  ( 0 )
 CHECKSHAPE  : Wires    = 0  ( 0 )  Faces    = 0  ( 0 )  Shells   = 0  ( 0 )   Solids   = 0 ( 0 )
 NBSHAPES    : Solid    = 0  ( 0 )  Shell    = 1  ( 1 )  Face     = 47  ( 47 ) 
 STATSHAPE   : Solid    = 0  ( 0 )  Shell    = 1  ( 1 )  Face     = 47  ( 47 )   FreeWire = 0  ( 0 ) 
-TOLERANCE   : MaxTol   =  0.09974534731  (    0.127442531 )  AvgTol   =   0.03628276746  (   0.04204326369 )
+TOLERANCE   : MaxTol   =  0.09974534731  (    0.127442531 )  AvgTol   =   0.03628276746  (   0.04253031643 )
 LABELS      : N0Labels = 1  ( 1 )  N1Labels = 0  ( 0 )  N2Labels = 0  ( 0 )   TotalLabels = 1  ( 1 )   NameLabels = 1  ( 1 )   ColorLabels = 0  ( 0 )   LayerLabels = 0  ( 0 )
 PROPS       : Centroid = 0  ( 0 )  Volume   = 0  ( 0 )  Area     = 0  ( 0 )
 NCOLORS     : NColors  = 0  ( 0 )
index e80f196..381b5bb 100755 (executable)
@@ -12,7 +12,7 @@ TPSTAT      : Faulties = 0  ( 2 )  Warnings = 4  ( 30 )  Summary  = 4  ( 32 )
 CHECKSHAPE  : Wires    = 1  ( 2 )  Faces    = 1  ( 2 )  Shells   = 0  ( 0 )   Solids   = 0 ( 0 )
 NBSHAPES    : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 416  ( 415 ) 
 STATSHAPE   : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 416  ( 415 )   FreeWire = 0  ( 0 ) 
-TOLERANCE   : MaxTol   =    133200.3972  (   0.9492292985 )  AvgTol   =     320.1208367  (   0.03928716897 )
+TOLERANCE   : MaxTol   =    133200.3972  (   0.9492292985 )  AvgTol   =     320.1208627  (   0.04011190067 )
 LABELS      : N0Labels = 1  ( 1 )  N1Labels = 28  ( 28 )  N2Labels = 0  ( 0 )   TotalLabels = 29  ( 29 )   NameLabels = 1  ( 1 )   ColorLabels = 29  ( 29 )   LayerLabels = 0  ( 0 )
 PROPS       : Centroid = 1  ( 1 )  Volume   = 1  ( 1 )  Area     = 1  ( 1 )
 NCOLORS     : NColors  = 2  ( 2 )
index e72e62a..0e64ce4 100644 (file)
@@ -7,7 +7,7 @@ TPSTAT      : Faulties = 0  ( 0 )  Warnings = 1  ( 27 )  Summary  = 1  ( 27 )
 CHECKSHAPE  : Wires    = 0  ( 0 )  Faces    = 0  ( 0 )  Shells   = 0  ( 0 )   Solids   = 0 ( 0 )
 NBSHAPES    : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 415  ( 415 ) 
 STATSHAPE   : Solid    = 1  ( 1 )  Shell    = 1  ( 1 )  Face     = 415  ( 415 )   FreeWire = 0  ( 0 ) 
-TOLERANCE   : MaxTol   =  0.09895613597  (   0.9492292985 )  AvgTol   =   0.01325689195  (   0.04014390707 )
+TOLERANCE   : MaxTol   =  0.09895613597  (   0.9492292985 )  AvgTol   =   0.01328296919  (   0.04097063947 )
 LABELS      : N0Labels = 1  ( 1 )  N1Labels = 28  ( 28 )  N2Labels = 0  ( 0 )   TotalLabels = 29  ( 29 )   NameLabels = 1  ( 1 )   ColorLabels = 29  ( 29 )   LayerLabels = 0  ( 0 )
 PROPS       : Centroid = 1  ( 1 )  Volume   = 1  ( 1 )  Area     = 1  ( 1 )
 NCOLORS     : NColors  = 2  ( 2 )