[Teem-commit] SF.net SVN: teem:[7197] teem/trunk/src/limn/splineFit .c

[<kin...@us...>]

Revision: 7197
          http://sourceforge.net/p/teem/code/7197
Author:   kindlmann
Date:     2024-07-11 12:11:37 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
refactored Newton-restep code to be less dumb and ugly

Modified Paths:
--------------
    teem/trunk/src/limn/splineFit.c

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-10 23:58:51 UTC (rev 7196)
+++ teem/trunk/src/limn/splineFit.c	2024-07-11 12:11:37 UTC (rev 7197)
@@ -982,18 +982,20 @@
 
 /* Paper pg 621 eqs (7) and (8): the change in the spline parameter to try to improve how
 Q(t) approximates vertex i; naming the parameter "t" since that's what the paper does.
-Note that the Newton step in (8), of the function that seems to be named f(t),
-is not itself the distance (or squared distance) between Q(t) and P, but the
-dot product betweeen (Q(t) - P) and Q'(t).  This function makes sure that the
-Newton step not increase the distance (it tries very hard to find a step that decreases
-the distance); this check is not described in the paper, and is not in the author's
-code, but it does matter in cases where the spline is a poor fit. */
+Note that the Newton step in paper's eq (8), i.e. the Newton step of the function that
+seems to be named f(t), is not itself the distance (or squared distance) between Q(t) and
+P, but instead the dot product betweeen (Q(t) - P) and Q'(t).  This function makes sure
+that the Newton step not increase the distance (it tries very hard to find a step that
+decreases the distance); this caution is not described in the paper, and is not in the
+author's code, but it does matter in cases where the spline is a poor fit, or where
+there are a small number of points */
 static double
 delta_t(double t0, double cap, const double P[2], const double V0[2], const double V1[2],
         const double V2[2], const double V3[2]) {
-  double Q0[2], Q1[2], Q2[2], QmP[2], ww[4], denom, delt = 0, sqd0;
+  double Q0[2], Q1[2], Q2[2], QmP[2], ww[4], denom, delt, sqd0, absdelt, chop = 0.25;
+  uint try;
 
-  /* SQD(T) = squared distance between Q(T) and P */
+  /* SQD(T) = squared distance between Q(T) and P, with side-effects! */
 #define SQD(TT)                                                                         \
   (CBD0(Q0, V0, V1, V2, V3, (TT), ww), /* */                                            \
    ELL_2V_SUB(QmP, Q0, P),             /* */                                            \
@@ -1002,35 +1004,28 @@
   CBD1(Q1, V0, V1, V2, V3, t0, ww);
   CBD2(Q2, V0, V1, V2, V3, t0, ww);
   denom = ELL_2V_DOT(Q1, Q1) + ELL_2V_DOT(QmP, Q2);
-  if (denom) {
-    double absdelt, chop = 4;
-    delt = -ELL_2V_DOT(QmP, Q1) / denom;
-    absdelt = fabs(delt);
-    if (absdelt > cap) {
-      /* cap Newton step */
-      delt *= cap / absdelt;
+  if (!denom) {
+    /* we can't do anything useful if we're about to have divide-by-zero, so bail */
+    return 0;
+  }
+  /* else we can compute Newton step */
+  delt = -ELL_2V_DOT(QmP, Q1) / denom;
+  absdelt = fabs(delt);
+  if (absdelt > cap) {
+    /* cap Newton step */
+    delt *= cap / absdelt;
+  }
+  try = 1;                        /* start with 1st try (1-based numbering) */
+  while (SQD(t0 + delt) > sqd0) { /* while this is an unhelpful step */
+    if (9 == try) {               /* 4 pairs of chop and flip didn't help; bail */
+      delt = 0;
+      break;
     }
-    if (SQD(t0 + delt) > sqd0) { /* true if squared distance has increased */
-      /* the Newton step took Q(t) farther from P, not closer, try to rescue it */
-      delt /= chop;
-      if (SQD(t0 + delt) > sqd0) {
-        /* the chopped down Newton step was unhelpful, try the other direction? */
-        delt *= -1;
-        if (SQD(t0 + delt) > sqd0) {
-          /* the negated chopped down Newton step was unhelpful! not giving up ... */
-          delt /= -chop;
-          if (SQD(t0 + delt) > sqd0) {
-            /* the positive doubly-chopped step still bad, trying one last flip */
-            /* (GLK tests produced evidence that this can in fact help) */
-            delt *= -1;
-            if (SQD(t0 + delt) > sqd0) {
-              /* ok finally giving up */
-              delt = 0;
-            }
-          }
-        }
-      }
-    }
+    delt *= (try % 2  /* alternate between two ideas on how to adjust delt */
+               ? chop /* on 1st (and subsequent odd) tries: scale down by chop */
+               : -1   /* else: flip sign (GLK tests show this does sometimes help) */
+    );
+    try++;
   }
 #undef SQD
   return delt;
@@ -1776,7 +1771,6 @@
       limnCbfPathNix(subpath);
     }
   }
-
   path->isLoop = lpnt->isLoop;
   return 0;
 }

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.