From b0922ff9a12f5aa408cc54c12a5cf1cadcfa37eb Mon Sep 17 00:00:00 2001
From: jfa <jfa@opencascade.com>
Date: Tue, 17 Sep 2024 17:51:36 +0100
Subject: [PATCH] 0033828: Modeling Algorithms - GCPnts_QuasiUniformDeflection
 returns very different results under small change in deflection

---
 src/GCPnts/GCPnts_QuasiUniformDeflection.cxx | 36 +++++++++++++++++---
 tests/bugs/moddata_3/bug33828                | 35 +++++++++++++++++++
 2 files changed, 66 insertions(+), 5 deletions(-)
 create mode 100644 tests/bugs/moddata_3/bug33828

diff --git a/src/GCPnts/GCPnts_QuasiUniformDeflection.cxx b/src/GCPnts/GCPnts_QuasiUniformDeflection.cxx
index fd357bfbc6..78a66eac80 100644
--- a/src/GCPnts/GCPnts_QuasiUniformDeflection.cxx
+++ b/src/GCPnts/GCPnts_QuasiUniformDeflection.cxx
@@ -103,29 +103,55 @@ static void QuasiFleche (const TheCurve& theC,
     aVdelta = theVfin;
   }
 
+  // square length of chord
   const Standard_Real aNorme = gp_Vec (thePdeb, aPdelta).SquareMagnitude();
   Standard_Real aFleche = 0.0;
   Standard_Boolean isFlecheOk = Standard_False;
-  if (aNorme > theEps)
+  if (aNorme > theEps && aNorme > 16. * theDeflection2)
   {
     // Evaluation de la fleche par interpolation . Voir IntWalk_IWalking_5.gxx
     Standard_Real N1 = theVdeb.SquareMagnitude();
     Standard_Real N2 = aVdelta.SquareMagnitude();
     if (N1 > theEps && N2 > theEps)
     {
+      // square distance between ends of two normalized vectors [0; 4]
       Standard_Real aNormediff = (theVdeb.Normalized().XYZ() - aVdelta.Normalized().XYZ()).SquareModulus();
       if (aNormediff > theEps)
       {
         aFleche = aNormediff * aNorme / 64.0;
+        // So, fleche <= (aNorme / 16), independently of Vdeb and Vdelta.
+        // And if (aNorme / 16) < theDeflection2, this approach gives
+        // fleche < theDeflection2 independently of real curve.
+        // That is why we exclude case aNorme < (16. * theDeflection2)
         isFlecheOk = Standard_True;
       }
     }
   }
-  if (!isFlecheOk)
+
+  gp_Pnt aPmid ((thePdeb.XYZ() + aPdelta.XYZ()) * 0.5);
+  gp_Pnt aPverif (Value (theC, theUdeb + aUdelta * 0.5));
+  Standard_Real aFlecheMidMid = aPmid.SquareDistance (aPverif);
+
+  if (isFlecheOk)
   {
-    gp_Pnt aPmid ((thePdeb.XYZ() + aPdelta.XYZ()) * 0.5);
-    gp_Pnt aPverif (Value (theC, theUdeb + aUdelta * 0.5));
-    aFleche = aPmid.SquareDistance (aPverif);
+    // Algorithm, evaluating "fleche" by interpolation,
+    // can give false-positive result.
+    // So we check also distance between Pmid and Pverif (aFlecheMidMid).
+    // But aFlecheMidMid gives worse result in case of non-uniform parameterisation.
+    // Maximum aFlecheMidMid, that seems reasonable, is (chord/2)^2 + theDeflection2
+    //   .---------------.Pverif  .
+    //   |               |        | Deflection
+    //   ._______. ______.        .
+    // Pdeb    Pmid    Pdelta
+    if (aFlecheMidMid > aNorme/4. + theDeflection2)
+    //if (aFlecheMidMid > aNorme/4.)
+    {
+      aFleche = aFlecheMidMid;
+    }
+  }
+  else
+  {
+    aFleche = aFlecheMidMid;
   }
 
   if (aFleche < theDeflection2)
diff --git a/tests/bugs/moddata_3/bug33828 b/tests/bugs/moddata_3/bug33828
new file mode 100644
index 0000000000..9dbdf6a337
--- /dev/null
+++ b/tests/bugs/moddata_3/bug33828
@@ -0,0 +1,35 @@
+puts "==========================================================="
+puts "0033828: Modeling Data - GCPnts_QuasiUniformDeflection"
+puts "returns very different results under small change in deflection"
+puts "==========================================================="
+
+proc check_crvpoints {cc deflection nb_expected} {
+    upvar ${cc} ${cc}
+
+    set str1 "Nb points +: +(\[-0-9.+eE\]+)\n"
+    set str2 "Max defl: +(\[-0-9.+eE\]+) +(\[-0-9.+eE\]+) +(\[-0-9.+eE\]+) +(\[-0-9.+eE\]+)"
+
+    set info [crvpoints r ${cc} ${deflection}]
+    regexp "${str1}${str2}" ${info} full Nb dmax ufmax ulmax i
+
+    if { ${Nb} != ${nb_expected} } {
+        puts "Error : bad value of Nb points = ${Nb}, expected ${nb_expected}"
+    }
+
+    if { ${dmax} > ${deflection} } {
+        puts "Error : bad value of maximum deflection = ${dmax}, expected < ${deflection}"
+    }
+}
+
+bsplinecurve cu 3 7  0 4  0.17 2  0.33 2  0.5 2  0.67 2  0.83 2  1 4   0.163 0.233 0 1  0.158 0.204 0 1  0.139 0.180 0 1  0.086 0.159 0 1  0.055 0.163 0 1  0.009 0.196 0 1  -0.004 0.225 0 1  0.002 0.281 0 1  0.019 0.307 0 1  0.070 0.332 0 1  0.101 0.331 0 1  0.149 0.301 0 1  0.164 0.274 0 1  0.163 0.246 0 1
+
+check_crvpoints cu .5     2
+check_crvpoints cu .1     3
+check_crvpoints cu .05    5
+check_crvpoints cu .025   5
+check_crvpoints cu .007   9
+check_crvpoints cu .005  17
+check_crvpoints cu .0005 33
+check_crvpoints cu .0003 65
+check_crvpoints cu .0002 65
+check_crvpoints cu .0001 73