teem-commit — gets email notification of all Teem svn commits

[Teem-commit] SF.net SVN: teem:[7211] teem/trunk/src/limn

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

Revision: 7211
          http://sourceforge.net/p/teem/code/7211
Author:   kindlmann
Date:     2024-07-16 23:00:09 +0000 (Tue, 16 Jul 2024)
Log Message:
-----------
it may actually be working ...

Modified Paths:
--------------
    teem/trunk/src/limn/lpu_cbfit.c
    teem/trunk/src/limn/test/02-test-fs.sh
    teem/trunk/src/limn/test/03-single.sh
    teem/trunk/src/limn/test/04-multi.sh
    teem/trunk/src/limn/test/05-go.sh

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-16 22:15:26 UTC (rev 7210)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-16 23:00:09 UTC (rev 7211)
@@ -137,7 +137,9 @@
   sprintf(buff, "%.17g", fctx->nrpPsi);
   hestOptAdd_1_Double(&hopt, "psi", "psi", &psi, buff, "psi, of course");
   sprintf(buff, "%.17g", fctx->cornAngle);
-  hestOptAdd_1_Double(&hopt, "ca", "angle", &cangle, buff, "angle indicating a corner");
+  hestOptAdd_1_Double(&hopt, "ca", "angle", &cangle, buff,
+                      "angle indicating a corner, or, "
+                      "0 to say that no corner finding should be done");
   sprintf(buff, "%.17g", fctx->scale);
   hestOptAdd_1_Double(&hopt, "scl", "scale", &scale, buff,
                       "scale for geometry estimation");
@@ -324,7 +326,12 @@
   fctx->nrpDeltaThresh = deltaThresh;
   fctx->nrpIota = nrpIota;
   fctx->nrpPsi = psi;
-  fctx->cornAngle = cangle;
+  if (cangle) {
+    fctx->cornAngle = cangle;
+    fctx->cornerFind = AIR_TRUE;
+  } else {
+    fctx->cornerFind = AIR_FALSE;
+  }
 
   if (tvt[3] >= 0) { /* here just to call limnCbfTVT once */
     double lt[2], vv[2], rt[2];

Modified: teem/trunk/src/limn/test/02-test-fs.sh
===================================================================
--- teem/trunk/src/limn/test/02-test-fs.sh	2024-07-16 22:15:26 UTC (rev 7210)
+++ teem/trunk/src/limn/test/02-test-fs.sh	2024-07-16 23:00:09 UTC (rev 7211)
@@ -32,8 +32,8 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-#IN=circ.txt
-IN=pointy.txt
+IN=circ.txt
+#IN=pointy.txt
 N=$(cat $IN | wc -l | xargs echo)
 BIN=900
 
@@ -40,7 +40,7 @@
 MMB="-min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN"
 
 unu jhisto -i $IN    $MMB -t float |
-    unu resample -s x1 x1 -k gauss:1.3,6 | unu quantize -b 8 | unu 2op x - 0.8 -o in.png
+    unu resample -s x1 x1 -k gauss:1.5,6 | unu quantize -b 8 | unu 2op x - 1 -o in.png
 junk in.png
 
 rm -f test-??.png
@@ -47,15 +47,16 @@
 
 for LO in $(seq 0 $((N-1))); do
     echo $LO
-    HI=$((LO+10))
+    HI=$((LO+7))
     LOO=$(printf %02d $LO)
-    CMD="./lpu cbfit -i $IN -loop -scl 2 -psi 1000 -fs $LO $HI -v 0 -eps 0.01"
+    # CMD="./lpu cbfit -i $IN -loop -scl 0.5 -psi 1000 -fs $LO $HI -v 0 -eps 0.5"
+      CMD="./lpu cbfit -i $IN -loop -scl  0  -psi 1000 -fs $LO $HI -v 0 -eps 0.01"
     echo "==================== $LO $HI --> test-$LOO.png : $CMD"
     eval $CMD 2>&1 > log-$LOO.txt
     # cat log-$LOO.txt
     junk log-$LOO.txt
     tail -n 4 log-$LOO.txt |
-        ./lpu cbfit -i - -synthn $N -syntho out.txt 2>&1 | grep -v _hestDefaults
+        ./lpu cbfit -i - -synthn 200 -syntho out.txt 2>&1 | grep -v _hestDefaults
     # lots of extraneous printfs thwart piping
     unu jhisto -i out.txt $MMB | unu quantize -b 8 -max 1 -o out.png
     #  in: green

Modified: teem/trunk/src/limn/test/03-single.sh
===================================================================
--- teem/trunk/src/limn/test/03-single.sh	2024-07-16 22:15:26 UTC (rev 7210)
+++ teem/trunk/src/limn/test/03-single.sh	2024-07-16 23:00:09 UTC (rev 7211)
@@ -37,9 +37,9 @@
 #    because the nrp parms that make sense for a small number of points don't work great
 #    for a huge number of points
 
-# # Good debugging test case: N=18 is a bad fit, N=19 is a perfect fit
-# BUT THEN the improved delta_t fixed this; so both are equally good!
-# # likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
+## Good debugging test case: N=18 is a bad fit, N=19 is a perfect fit
+## BUT THEN the improved delta_t fixed this; so both are equally good!
+## likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
 #N=18
 #echo "-0.5 0.5
 # 2.0  0.5
@@ -53,7 +53,6 @@
 # 1 1.5
 #-1 1.5
 # 1 -1" | ./lpu cbfit -i - -synthn $N -sup 3 -syntho xy-inn-$N.txt
-#junk xy-inn-$N.txt
 
 # This is demo of why we might want a step that optimizes tangent directions
 # after the parameterization has been found
@@ -62,8 +61,9 @@
  1 3
 -3 -1
  1 -1" | ./lpu cbfit -i - -synthn $N -v 0 -syntho - | unu crop -min 0 0 -max M M-1 -o xy-inn-$N.txt
-# junk xy-inn-$N.txt
 
+junk xy-inn-$N.txt
+
 CMD="./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 1 -psi 1000000000 -nim 4000 -deltathr 0.000000000001"
 echo "====== $CMD"
 eval $CMD > log.txt

Modified: teem/trunk/src/limn/test/04-multi.sh
===================================================================
--- teem/trunk/src/limn/test/04-multi.sh	2024-07-16 22:15:26 UTC (rev 7210)
+++ teem/trunk/src/limn/test/04-multi.sh	2024-07-16 23:00:09 UTC (rev 7211)
@@ -45,11 +45,11 @@
 cat 0.txt 1.txt | uniq > xy-inn.txt
 junk {0,1}.txt xy-inn.txt
 
-# IN=xy-inn.txt
+IN=xy-inn.txt
 #IN=pointy.txt
-IN=circ.txt
+#IN=circ.txt
 
-CMD="./lpu cbfit -i $IN -scl 0 -fm 0 -1 -v 2 -psi 3 -eps 0.01"
+CMD="./lpu cbfit -i $IN -scl 0.5 -fm 0 -1 -v 2 -psi 3 -eps 0.01"
 echo "====== $CMD"
 eval $CMD > log.txt
 cat log.txt # ; junk log.txt

Modified: teem/trunk/src/limn/test/05-go.sh
===================================================================
--- teem/trunk/src/limn/test/05-go.sh	2024-07-16 22:15:26 UTC (rev 7210)
+++ teem/trunk/src/limn/test/05-go.sh	2024-07-16 23:00:09 UTC (rev 7211)
@@ -54,8 +54,21 @@
 
 rm -f tmp.png
 
+SCL=3
+EPS=0.08
 rm -f log.txt
-CMD="./lpu cbfit -i $IN -loop -scl 0 -v 0 -eps 0.03 -roll 10"
+# CMD="./lpu cbfit -i $IN             -scl   0  -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN       -ca 0 -scl   0  -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN             -scl $SCL -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN       -ca 0 -scl $SCL -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN -loop       -scl   0  -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN -loop -ca 0 -scl   0  -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN -loop       -scl   0  -v 0 -eps $EPS -roll 5"
+# CMD="./lpu cbfit -i $IN -loop -ca 0 -scl   0  -v 0 -eps $EPS -roll 5"
+# CMD="./lpu cbfit -i $IN -loop       -scl $SCL -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN -loop -ca 0 -scl $SCL -v 0 -eps $EPS -roll 0"
+# CMD="./lpu cbfit -i $IN -loop       -scl $SCL -v 0 -eps $EPS -roll 5"
+  CMD="./lpu cbfit -i $IN -loop -ca 0 -scl $SCL -v 0 -eps $EPS -roll 5"
 eval $CMD  2>&1 > log.txt ||:
 echo "====== $CMD"
 cat log.txt # ; junk log.txt
@@ -64,11 +77,11 @@
 echo "====== RESULTS: --> $OUT"
 grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9
 grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9 |
-    ./lpu cbfit -i - -loop -synthn 300 -sup 1 -syntho $OUT
+    ./lpu cbfit -i - -loop -synthn 900 -sup 1 -syntho $OUT
 junk $OUT
 
 
-BIN=900
+BIN=660
 MM="-min -1.1 1.1 -max 1.1 -1.1"
 unu jhisto -i  $IN $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
 unu jhisto -i $OUT $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png

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

[Teem-commit] SF.net SVN: teem:[7210] teem/trunk/src/limn

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

Revision: 7210
          http://sourceforge.net/p/teem/code/7210
Author:   kindlmann
Date:     2024-07-16 22:15:26 +0000 (Tue, 16 Jul 2024)
Log Message:
-----------
tweaks to new flexiblity in loi==hii, and in its testing

Modified Paths:
--------------
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/00-data.sh
    teem/trunk/src/limn/test/01-test-tvt.sh

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-16 20:26:18 UTC (rev 7209)
+++ teem/trunk/src/limn/splineFit.c	2024-07-16 22:15:26 UTC (rev 7210)
@@ -650,8 +650,8 @@
 
 /* idxLift: error-checked index lifting (from "actual" to "lifted" indices, explained
 above) for limnCbfTVT and maybe others. If no error, *loiP will be same as given gloi,
-but in loops *vviP and *hiiP may be lifted (relative to gvvi and ghii), and outside loops
-*hiiP may be changed to #points-1.
+but in loops *vviP and/or *hiiP may be lifted (relative to gvvi and ghii), and outside
+loops *hiiP may be changed to #points-1.
 
 That sounds like nothing fancy, but this is messy because of the flexibility in how we
 handle points: might not be a loop or might not, and, consideration of vertices should
@@ -665,7 +665,7 @@
   static const char me[] = "idxLift";
   uint pnum = lpnt->num, loi, hii, vvi;
 
-  *loiP = *hiiP = *vviP = UINT_MAX; /* initialize to bogus indices */
+  *loiP = *hiiP = *vviP = UINT_MAX; /* ini  tialize to bogus indices */
   if (!(pnum < (1U << 29))) {
     /* UB = undefined behavior */
     biffAddf(LIMN, "%s: # points %u seems too big (to stay well clear of UB)", me, pnum);
@@ -682,8 +682,8 @@
   hii = ghii;
   vvi = gvvi;
   if (lpnt->isLoop) {
-    if (!(gloi < ghii)) hii += pnum;
-    /* now loi < hii */
+    if (!(gloi <= ghii)) hii += pnum;
+    /* now loi <= hii ; loi==hii possible for any valid index */
     if (gvvi < gloi) vvi += pnum;
     /* now loi <= vvi */
   } else {
@@ -698,7 +698,7 @@
       /* else gloi==ghii==0 */
       hii = pnum - 1;
       /* now loi < hii */
-      /* because of uint type 0 == loi <= vvi */
+      /* because of uint type, 0 == loi <= vvi */
     } else { /* gloi != ghii */
       if (!(gloi < ghii)) {
         biffAddf(LIMN, "%s: in non-loop, need loi (%u) < hii (%u)", me, gloi, ghii);
@@ -713,18 +713,21 @@
       /* now loi <= vii */
     }
   }
-  /* now, in any case: must have loi < hii and loi <= vvi */
+  /* now, in any case: must have loi <= hii (and loi==hii only possible in loop)
+     and loi <= vvi */
   if (verbose) {
     printf("%s: given loi,hii,vvi %u %u %u --> lifted %u %u %u\n", me, gloi, ghii, gvvi,
            loi, hii, vvi);
   }
-  /* make sure that vvi is below upper bound hii */
-  if (!(vvi <= hii)) {
-    biffAddf(LIMN, "%s: vvi %u->%u not in [%u,%u]->[%u,%u] span", me, gvvi, vvi, /* */
-             gloi, ghii, loi, hii);
-    return 1;
+  if (loi < hii) {
+    /* if have consequential bounds, check that vvi is in closed-interval [loi,hii];
+       (the test w.r.t. loi is actually redundant with above, but here for clarity) */
+    if (!(loi <= vvi && vvi <= hii)) {
+      biffAddf(LIMN, "%s: vvi %u->%u not in [%u,%u]->[%u,%u] span", me, gvvi, vvi, /* */
+               gloi, ghii, loi, hii);
+      return 1;
+    }
   }
-  /* now, in any case: loi <= vvi <= hii */
 
   /* all's well, set output values */
   *loiP = loi;
@@ -866,7 +869,7 @@
                 : sui;
       const double *xy = PPlowerI(lpnt, sbi); /* coords at sbi */
       ELL_2V_SCALE_INCR(posC, vw, xy);
-      if (fctx->verbose > 1) {
+      if (fctx->verbose > 2) {
         printf("%s: ci=%d (in [%d,%d]) idx %d --[%d,%d]--> %d;  v,t w %g,%g on "
                "xy=(%g,%g)\n",
                me, ci, -lim, lim, sui, slo, shi, sbi, vw, tw, xy[0], xy[1]);
@@ -874,13 +877,13 @@
       }
       if (ci < 0) {
         ELL_2V_SCALE_INCR(posM, tw, xy);
-        if (fctx->verbose > 1) {
+        if (fctx->verbose > 2) {
           printf("%s:   ---> posM=(%g,%g)\n", me, posM[0], posM[1]);
         }
       }
       if (ci > 0) {
         ELL_2V_SCALE_INCR(posP, tw, xy);
-        if (fctx->verbose > 1) {
+        if (fctx->verbose > 2) {
           printf("%s:   ---> posP=(%g,%g)\n", me, posP[0], posP[1]);
         }
       }
@@ -1430,7 +1433,7 @@
     ELL_2V_COPY(vttv[3], vv3);
     if (givenP) *givenP = AIR_TRUE;
   } else {
-    double v0c[2], t1c[2], t2c[2], v3c[3]; /* locally computed info */
+    double v0c[2], t1c[2], t2c[2], v3c[3]; /* to store locally computed info */
     if (vv0 || tt1 || tt2 || vv3) {
       biffAddf(LIMN,
                "%s: should either give all vv0,tt1,tt2,vv3 or none (not %p,%p,%p,%p)",
@@ -1438,14 +1441,14 @@
                (const void *)vv3);
       return 1;
     }
-    if (loi == hii && hii != 0) {
+    /* if (loi == hii && hii != 0) {    (no longer true)
       biffAddf(LIMN, "%s: can only handle loi=hii=%u if both 0", me, hii);
       return 1;
-    }
-    if (loi == hii) { /* and hence both 0, due to above */
-      /* e.g., we're doing the first fit on of a circular (corner-less) path */
+    } */
+    if (loi == hii) {
+      /* e.g., we're doing the first fit on a circular path with 0 or 1 corners */
       if (!lpnt->isLoop) {
-        biffAddf(LIMN, "%s: can handle loi=hii=0 only with point loop", me);
+        biffAddf(LIMN, "%s: can handle loi=hii only with point loop", me);
         return 1;
       }
       if (limnCbfTVT(t2c, v0c, t1c,             /* */
@@ -1456,7 +1459,7 @@
       }
       ELL_2V_COPY(v3c, v0c);
     } else {
-      /* do not have geometry info; must find it all */
+      /* loi < hii */
       if (limnCbfTVT(NULL, v0c, t1c,            /* */
                      fctx, lpnt, loi, hii, loi, /* */
                      AIR_TRUE /* yes oneSided */)

Modified: teem/trunk/src/limn/test/00-data.sh
===================================================================
--- teem/trunk/src/limn/test/00-data.sh	2024-07-16 20:26:18 UTC (rev 7209)
+++ teem/trunk/src/limn/test/00-data.sh	2024-07-16 22:15:26 UTC (rev 7210)
@@ -32,7 +32,7 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-N=50
+N=24
 
 echo 0 1 |
     unu reshape -s 2 |

Modified: teem/trunk/src/limn/test/01-test-tvt.sh
===================================================================
--- teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-16 20:26:18 UTC (rev 7209)
+++ teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-16 22:15:26 UTC (rev 7210)
@@ -40,24 +40,60 @@
 RTOUT=out-lt.txt
 BIN=670
 
-rm -f $VVOUT; touch $VVOUT; junk $VVOUT
-rm -f $LTOUT; touch $LTOUT; junk $LTOUT
-rm -f $RTOUT; touch $RTOUT; junk $RTOUT
+EPS=0.2
+SCL=1.1
+VERB=2
+# initialize outputs to be non-empty but also non-consequential
+rm -f $VVOUT; echo "nan nan" > $VVOUT; junk $VVOUT
+rm -f $LTOUT; echo "nan nan" > $LTOUT; junk $LTOUT
+rm -f $RTOUT; echo "nan nan" > $RTOUT; junk $RTOUT
+echo "for I in seq 0 $((N-1)) ..."
 for I in $(seq 0 $((N-1))); do
-    # 16-fold (!) TEST:
-    # 8: without -loop, versus with -loop
-    # 4: LO=HI=0 (or LO=HI=5 in loop), versus some interval around I
-    # 2: oneside (4th arg to -tvt) 0 versus 1
-    # 1: -scl 0 versus >0
-    LO=0 # $((I-4))
-    HI=0 # $((I+4))
-    CMD="./lpu cbfit -i $IN -tvt $LO $HI $I 0 -scl 0 -eps 1 -v 0"
-    echo $CMD
+    ## 16-fold, and then some (!) TEST:
+    ## 8: without -loop, versus with -loop
+    ## 4: LO=HI=0 (or LO=HI=5 in loop), versus some interval around I (centered or not)
+    ## 2: -scl 0 versus >0
+    ## 1: oneside (4th arg to -tvt) 0 versus 1
+      CMD="./lpu cbfit -i $IN       -tvt     0        0    $I 0 -scl   0  -eps $EPS -v $VERB" # 0
+    # CMD="./lpu cbfit -i $IN       -tvt     0        0    $I 1 -scl   0  -eps $EPS -v $VERB" # 1
+    # CMD="./lpu cbfit -i $IN       -tvt     0        0    $I 0 -scl $SCL -eps $EPS -v $VERB" # 2
+    # CMD="./lpu cbfit -i $IN       -tvt     0        0    $I 1 -scl $SCL -eps $EPS -v $VERB" # 3
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-4)) $((I+4)) $I 0 -scl   0  -eps $EPS -v $VERB" # 4
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-1)) $((I+7)) $I 0 -scl   0  -eps $EPS -v $VERB" # 4.5
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-4)) $((I+4)) $I 1 -scl   0  -eps $EPS -v $VERB" # 5
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-1)) $((I+7)) $I 1 -scl   0  -eps $EPS -v $VERB" # 5.5
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-4)) $((I+4)) $I 0 -scl $SCL -eps $EPS -v $VERB" # 6
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-1)) $((I+7)) $I 0 -scl $SCL -eps $EPS -v $VERB" # 6.5
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-4)) $((I+4)) $I 1 -scl $SCL -eps $EPS -v $VERB" # 7
+    # CMD="./lpu cbfit -i $IN       -tvt $((I-1)) $((I+7)) $I 1 -scl $SCL -eps $EPS -v $VERB" # 7.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt     0        0    $I 0 -scl   0  -eps $EPS -v $VERB" # 8
+    # CMD="./lpu cbfit -i $IN -loop -tvt     4        4    $I 0 -scl   0  -eps $EPS -v $VERB" # 8.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt     0        0    $I 1 -scl   0  -eps $EPS -v $VERB" # 9
+    # CMD="./lpu cbfit -i $IN -loop -tvt     4        4    $I 1 -scl   0  -eps $EPS -v $VERB" # 9.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt     0        0    $I 0 -scl $SCL -eps $EPS -v $VERB" # 10
+    # CMD="./lpu cbfit -i $IN -loop -tvt     4        4    $I 0 -scl $SCL -eps $EPS -v $VERB" # 10.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt     0        0    $I 1 -scl $SCL -eps $EPS -v $VERB" # 11
+    # CMD="./lpu cbfit -i $IN -loop -tvt     4        4    $I 1 -scl $SCL -eps $EPS -v $VERB" # 11.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-4)) $((I+4)) $I 0 -scl   0  -eps $EPS -v $VERB" # 12
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-1)) $((I+7)) $I 0 -scl   0  -eps $EPS -v $VERB" # 12.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-4)) $((I+4)) $I 1 -scl   0  -eps $EPS -v $VERB" # 13
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-1)) $((I+7)) $I 1 -scl   0  -eps $EPS -v $VERB" # 13.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-4)) $((I+4)) $I 0 -scl $SCL -eps $EPS -v $VERB" # 14
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-1)) $((I+7)) $I 0 -scl $SCL -eps $EPS -v $VERB" # 14.5
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-4)) $((I+4)) $I 1 -scl $SCL -eps $EPS -v $VERB" # 15
+    # CMD="./lpu cbfit -i $IN -loop -tvt $((I-1)) $((I+7)) $I 1 -scl $SCL -eps $EPS -v $VERB" # 15.5
+    echo "$I/$N: $CMD"
     rm -f log.txt
+    set +o errexit
     (eval $CMD 2>&1) > log.txt
-    grep "lt =" log.txt | cut -d' ' -f 3,4 >> $LTOUT
-    grep "vv =" log.txt | cut -d' ' -f 3,4 >> $VVOUT
-    grep "rt =" log.txt | cut -d' ' -f 3,4 >> $RTOUT
+    RET=$?
+    set -o errexit
+    if [ 0 -eq $RET ]; then
+        grep "^limnCbfTVT" log.txt ||:
+        grep "lt =" log.txt | cut -d' ' -f 3,4 >> $LTOUT
+        grep "vv =" log.txt | cut -d' ' -f 3,4 >> $VVOUT
+        grep "rt =" log.txt | cut -d' ' -f 3,4 >> $RTOUT
+    fi
 done
 rm -f log.txt
 
@@ -68,7 +104,7 @@
 unu jhisto -i $VVOUT $MMB | unu quantize -b 8 -max 1 -o vv.png
 junk vv.png
 
-TSCL=0.02
+TSCL=0.08
 unu 2op x $LTOUT $TSCL | unu 2op + - $VVOUT | unu jhisto $MMB | unu quantize -b 8 -max 1 -o lt.png
 unu 2op x $RTOUT $TSCL | unu 2op + - $VVOUT | unu jhisto $MMB | unu quantize -b 8 -max 1 -o rt.png
 junk {l,r}t.png

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

[Teem-commit] SF.net SVN: teem:[7209] teem/trunk/src/limn

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

Revision: 7209
          http://sourceforge.net/p/teem/code/7209
Author:   kindlmann
Date:     2024-07-16 20:26:18 +0000 (Tue, 16 Jul 2024)
Log Message:
-----------
tweaks to new flexiblity in loi==hii, and in its testing

Modified Paths:
--------------
    teem/trunk/src/limn/lpu_cbfit.c
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/00-data.sh
    teem/trunk/src/limn/test/01-test-tvt.sh
    teem/trunk/src/limn/test/05-go.sh

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-16 18:50:21 UTC (rev 7208)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-16 20:26:18 UTC (rev 7209)
@@ -118,7 +118,8 @@
                       "when synthesizing data on a single segment, warp U parameters "
                       "by raising to this power.");
   hestOptAdd_4_Int(&hopt, "tvt", "loi hii vvi 1s", tvt, "0 0 0 -1",
-                   "if last value is >= 0: make single call to limnCbfTVT and quit");
+                   "if last value is >= 0: quit after making single call to limnCbfTVT, "
+                   "with final 1s value being passed for oneSided");
   sprintf(buff, "%u", fctx->nrpIterMax);
   hestOptAdd_1_UInt(&hopt, "nim", "max", &nrpIterMax, buff,
                     "max # nrp iterations to run");

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-16 18:50:21 UTC (rev 7208)
+++ teem/trunk/src/limn/splineFit.c	2024-07-16 20:26:18 UTC (rev 7209)
@@ -129,7 +129,7 @@
 
 limnCbfGo -- limnCbfCtxPrep
           -- limnCbfCorners
-          -- limnCbfMulti (either with loi==hii==0 or with loi,hii at corners)
+          -- limnCbfMulti (either with loi==hii or with loi,hii at corners)
           -- limnCbfPathNew, limnCbfPathJoin, limnCbfPathNix
 */
 
@@ -655,8 +655,10 @@
 
 That sounds like nothing fancy, but this is messy because of the flexibility in how we
 handle points: might not be a loop or might not, and, consideration of vertices should
-either be bounded in a specific [loi,hii] or be "unbounded" loi==hii==0 (truly unbounded
-in a loop, or bounded only as much as needed in non-loop data). */
+either be bounded in a specific [loi,hii] or be "unbounded" loi==hii.  "unbounded" is
+requested in different ways: in a loop, loi==hii and they can be any valid index
+(e.g. a loop with a single corner at index 5), or, in a non-loop, loi==hii==0 which
+is a way of saying loi==0,hii=pnum-1. */
 static int /* Biff: 1 */
 idxLift(uint *loiP, uint *hiiP, uint *vviP, int verbose, const limnCbfPoints *lpnt,
         uint gloi, uint ghii, uint gvvi) {
@@ -675,56 +677,54 @@
     return 1;
   }
   /* now all of gloi, gvvi, ghii are all valid actual indices */
-  if (gloi == ghii && ghii != 0) {
-    biffAddf(LIMN,
-             "%s: can only have gloi == ghii if both 0 (not %u), "
-             "to signify unbounded vertex consideration",
-             me, gloi);
-    return 1;
-  }
   /* initialize values to return */
   loi = gloi;
   hii = ghii;
   vvi = gvvi;
   if (lpnt->isLoop) {
-    if (gloi != ghii) { /* implies both == 0 because of test above */
-      if (gloi > ghii) hii += pnum;
-      if (gloi > gvvi) vvi += pnum;
-    }
+    if (!(gloi < ghii)) hii += pnum;
+    /* now loi < hii */
+    if (gvvi < gloi) vvi += pnum;
+    /* now loi <= vvi */
   } else {
-    if (gloi == ghii) { /* (implies both == 0, again) */
-      /* we allow loi==hii==0 in non-loop to say: only bounded by data itself
-      loi is already 0, but hii needs fixing */
-      hii = pnum - 1;
-    } else {
-      if (gloi > ghii) {
+    if (gloi == ghii) {
+      if (ghii != 0) {
         biffAddf(LIMN,
-                 "%s: if loi != hii, need loi (%u) < hii (%u) since not in a "
-                 "point loop",
-                 me, gloi, ghii);
+                 "%s: in non-loop, can only have gloi == ghii if both 0 (not %u) "
+                 "to signify vertex consideration only bounded by valid indices",
+                 me, gloi);
         return 1;
       }
-      if (gloi > gvvi) {
-        biffAddf(LIMN, "%s: need given loi (%u) < vvi (%u) since not in point loop", me,
-                 gloi, gvvi);
+      /* else gloi==ghii==0 */
+      hii = pnum - 1;
+      /* now loi < hii */
+      /* because of uint type 0 == loi <= vvi */
+    } else { /* gloi != ghii */
+      if (!(gloi < ghii)) {
+        biffAddf(LIMN, "%s: in non-loop, need loi (%u) < hii (%u)", me, gloi, ghii);
         return 1;
       }
+      /* now loi < hii */
+      if (gvvi < gloi) {
+        biffAddf(LIMN, "%s: in non-loop, need given loi (%u) <= vvi (%u)", me, gloi,
+                 gvvi);
+        return 1;
+      }
+      /* now loi <= vii */
     }
   }
-  /* now: must have loi <= vvi and loi <= hii */
+  /* now, in any case: must have loi < hii and loi <= vvi */
   if (verbose) {
     printf("%s: given loi,hii,vvi %u %u %u --> lifted %u %u %u\n", me, gloi, ghii, gvvi,
            loi, hii, vvi);
   }
-  if (loi < hii) {
-    /* need to check that vvi is inside consequential bounds [loi,hii] */
-    if (vvi > hii) {
-      biffAddf(LIMN, "%s: vvi %u->%u not in [%u,%u]->[%u,%u] span", me, gvvi, vvi, gloi,
-               ghii, loi, hii);
-      return 1;
-    }
-    /* now (if bounded) have vvi <= hii */
+  /* make sure that vvi is below upper bound hii */
+  if (!(vvi <= hii)) {
+    biffAddf(LIMN, "%s: vvi %u->%u not in [%u,%u]->[%u,%u] span", me, gvvi, vvi, /* */
+             gloi, ghii, loi, hii);
+    return 1;
   }
+  /* now, in any case: loi <= vvi <= hii */
 
   /* all's well, set output values */
   *loiP = loi;
@@ -752,10 +752,10 @@
 
 /* utility function for counting how many vertices are in (actual) index span [loi,hii]
 inclusive. It is not our job here to care about lpnt->isLoop; we just assume that if
-we're faced with hii<loi, it must be because of a loop */
+we're faced with hii<=loi, it must be because of a loop */
 static uint
 spanLength(const limnCbfPoints *lpnt, uint loi, uint hii) {
-  uint topi = hii + (hii < loi) * (lpnt->num);
+  uint topi = hii + (hii <= loi) * (lpnt->num);
   return topi - loi + 1;
 }
 
@@ -1132,11 +1132,6 @@
   const double *uu = fctx->uu;
 
   spanlen = spanLength(lpnt, loi, hii);
-  if (0 == loi && 0 == hii && lpnt->isLoop) {
-    /* in the very specific (and rare) case that we've fit a *single* spline
-    to a point loop, starting and ending at 0, then we have to lift spanlen */
-    spanlen += lpnt->num;
-  }
   if (!(spanlen >= 3)) {
     biffAddf(LIMN, "%s: [loi,hii] [%u,%u] -> spanlen %u too small", me, loi, hii,
              spanlen);
@@ -1265,16 +1260,17 @@
     biffAddf(LIMN, "%s: got NULL pointer", me);
     return 1;
   }
-  spanlen = spanLength(lpnt, loi, hii);
   if (!(fctx->uu)) {
     biffAddf(LIMN, "%s: fcgtx->uu NULL; was limnCbfCtxPrep called?", me);
     return 1;
   }
+  spanlen = spanLength(lpnt, loi, hii);
   /* DIM=2 pretty much everywhere here */
   if (fctx->verbose) {
-    printf("%s[%d,%d]: hello, vv0=(%g,%g), tt1=(%g,%g), "
+    printf("%s[%d,%d(spanlen=%u)]: hello, vv0=(%g,%g), tt1=(%g,%g), "
            "tt2=(%g,%g), vv3=(%g,%g)\n",
-           me, loi, hii, vv0[0], vv0[1], tt1[0], tt1[1], tt2[0], tt2[1], vv3[0], vv3[1]);
+           me, loi, hii, spanlen, vv0[0], vv0[1], tt1[0], tt1[1], tt2[0], tt2[1], vv3[0],
+           vv3[1]);
   }
   if (2 == spanlen) {
     /* relying on code in findAlpha() that handles spanlen==2; return should be 1 */

Modified: teem/trunk/src/limn/test/00-data.sh
===================================================================
--- teem/trunk/src/limn/test/00-data.sh	2024-07-16 18:50:21 UTC (rev 7208)
+++ teem/trunk/src/limn/test/00-data.sh	2024-07-16 20:26:18 UTC (rev 7209)
@@ -32,7 +32,7 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-N=160
+N=50
 
 echo 0 1 |
     unu reshape -s 2 |

Modified: teem/trunk/src/limn/test/01-test-tvt.sh
===================================================================
--- teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-16 18:50:21 UTC (rev 7208)
+++ teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-16 20:26:18 UTC (rev 7209)
@@ -44,14 +44,14 @@
 rm -f $LTOUT; touch $LTOUT; junk $LTOUT
 rm -f $RTOUT; touch $RTOUT; junk $RTOUT
 for I in $(seq 0 $((N-1))); do
-    LO=$((I-4))
-    HI=$((I+4))
     # 16-fold (!) TEST:
-    # * without -loop and with -loop
-    # * -scl 0 and >0
-    # * LO=HI=0  versus something around I
-    # * oneside (4th arg to -tvt) 0 and 1
-    CMD="./lpu cbfit -i $IN -loop -scl 2 -tvt $LO $HI $I 1 -eps 1 -v 0"
+    # 8: without -loop, versus with -loop
+    # 4: LO=HI=0 (or LO=HI=5 in loop), versus some interval around I
+    # 2: oneside (4th arg to -tvt) 0 versus 1
+    # 1: -scl 0 versus >0
+    LO=0 # $((I-4))
+    HI=0 # $((I+4))
+    CMD="./lpu cbfit -i $IN -tvt $LO $HI $I 0 -scl 0 -eps 1 -v 0"
     echo $CMD
     rm -f log.txt
     (eval $CMD 2>&1) > log.txt

Modified: teem/trunk/src/limn/test/05-go.sh
===================================================================
--- teem/trunk/src/limn/test/05-go.sh	2024-07-16 18:50:21 UTC (rev 7208)
+++ teem/trunk/src/limn/test/05-go.sh	2024-07-16 20:26:18 UTC (rev 7209)
@@ -52,9 +52,12 @@
 #IN=pointy.txt
 #IN=circ.txt
 
-CMD="./lpu cbfit -i $IN -loop -scl 0.5 -v 3 -psi 3 -eps 0.01 -roll 1"
+rm -f tmp.png
+
+rm -f log.txt
+CMD="./lpu cbfit -i $IN -loop -scl 0 -v 0 -eps 0.03 -roll 10"
+eval $CMD  2>&1 > log.txt ||:
 echo "====== $CMD"
-eval $CMD > log.txt
 cat log.txt # ; junk log.txt
 
 OUT=xy-out.txt
@@ -64,12 +67,12 @@
     ./lpu cbfit -i - -loop -synthn 300 -sup 1 -syntho $OUT
 junk $OUT
 
+
 BIN=900
 MM="-min -1.1 1.1 -max 1.1 -1.1"
 unu jhisto -i  $IN $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
 unu jhisto -i $OUT $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
 
-rm -f tmp.png
 
 unu join -i xy-{out,inn,out}.png -a 0 -incr |
     unu resample -s = x2 x2 -k box -o tmp.png

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

[Teem-commit] SF.net SVN: teem:[7208] teem/trunk/src/limn

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

Revision: 7208
          http://sourceforge.net/p/teem/code/7208
Author:   kindlmann
Date:     2024-07-16 18:50:21 +0000 (Tue, 16 Jul 2024)
Log Message:
-----------
still debugging

Modified Paths:
--------------
    teem/trunk/src/limn/lpu_cbfit.c
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/03-single.sh
    teem/trunk/src/limn/test/05-go.sh

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-12 19:23:37 UTC (rev 7207)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-16 18:50:21 UTC (rev 7208)
@@ -262,13 +262,14 @@
     airMopAdd(mop, xy, airFree, airMopAlways);
     if (2 == size0) {
       unsigned int ci;
-      printf("%s: synthetically sampling single spline with %u points\n", me, synthNum);
+      fprintf(stderr, "%s: synthetically sampling single spline with %u points\n", me,
+              synthNum);
       for (ci = 0; ci < 8; ci++) {
         seg.xy[ci] = cpt[ci];
       }
-      printf("%s: synth seg: (%g,%g) -- (%g,%g) -- (%g,%g) -- (%g,%g)\n", me, seg.xy[0],
-             seg.xy[1], seg.xy[2], seg.xy[3], seg.xy[4], seg.xy[5], seg.xy[6],
-             seg.xy[7]);
+      fprintf(stderr, "%s: synth seg: (%g,%g) -- (%g,%g) -- (%g,%g) -- (%g,%g)\n", me,
+              seg.xy[0], seg.xy[1], seg.xy[2], seg.xy[3], seg.xy[4], seg.xy[5],
+              seg.xy[6], seg.xy[7]);
       for (ii = 0; ii < synthNum; ii++) {
         double uu = AIR_AFFINE(0, ii, synthNum - 1, 0, 1);
         uu = pow(uu, synthPow);
@@ -278,8 +279,8 @@
       unsigned int ci, si;
       limnCbfPath *spath = limnCbfPathNew(size1);
       airMopAdd(mop, spath, (airMopper)limnCbfPathNix, airMopAlways);
-      printf("%s: synthetically sampling %u splines with %u points\n", me, size1,
-             synthNum);
+      fprintf(stderr, "%s: synthetically sampling %u splines with %u points\n", me,
+              size1, synthNum);
       /* copy in control point data */
       for (si = 0; si < size1; si++) {
         for (ci = 0; ci < 8; ci++) {
@@ -435,6 +436,7 @@
     getchar();
   }
 
+  fprintf(stderr, "%s: calling limnCbfGo (verbose=%d)\n", me, fctx->verbose);
   if (limnCbfGo(path, fctx, lpnt)) {
     airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
     fprintf(stderr, "%s: trouble doing fitting:\n%s", me, err);

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-12 19:23:37 UTC (rev 7207)
+++ teem/trunk/src/limn/splineFit.c	2024-07-16 18:50:21 UTC (rev 7208)
@@ -1717,9 +1717,10 @@
         od = fctx->ctvt + 6 * ci;
         ELL_6V_COPY(od, id);
         if (fctx->verbose) {
-          printf("    corner %u is vertex %u\n"
+          printf("    corner %u is vertex %u%s\n"
                  "        T,V,T = (%g,%g)  (%g,%g)  (%g,%g)\n",
-                 ci, vi, od[0], od[1], od[2], od[3], od[4], od[5]);
+                 ci, vi, vi == pnum - 1 ? " (last vert)" : "", od[0], od[1], od[2],
+                 od[3], od[4], od[5]);
         }
         ci++;
       }
@@ -1773,10 +1774,11 @@
 
   /* first try fitting a single spline */
   if (fctx->verbose) {
-    printf("%s[%u,%u]_%u: trying single fit on all points\n", me, loi, hii, recDepth);
+    printf("%s[%u,%u]_%u: trying initial fitSingle on all points\n", me, loi, hii,
+           recDepth);
   }
   if (fitSingle(alpha, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
-    biffAddf(LIMN, "%s: fitSingle failed", me);
+    biffAddf(LIMN, "%s: initial fitSingle failed", me);
     return 1;
   }
   if (fctx->distBig <= 1) {

Modified: teem/trunk/src/limn/test/03-single.sh
===================================================================
--- teem/trunk/src/limn/test/03-single.sh	2024-07-12 19:23:37 UTC (rev 7207)
+++ teem/trunk/src/limn/test/03-single.sh	2024-07-16 18:50:21 UTC (rev 7208)
@@ -44,16 +44,25 @@
 #echo "-0.5 0.5
 # 2.0  0.5
 #-0.5  0.0
-# 0.5 -0.5" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
+# 0.5 -0.5" |./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
 
+## This is demo of why we might want a step that optimizes tangent directions
+## after the parameterization has been found
+#N=42
+#echo "-1 -1
+# 1 1.5
+#-1 1.5
+# 1 -1" | ./lpu cbfit -i - -synthn $N -sup 3 -syntho xy-inn-$N.txt
+#junk xy-inn-$N.txt
+
 # This is demo of why we might want a step that optimizes tangent directions
 # after the parameterization has been found
-N=42
-echo "-1 -1
- 1 1.5
--1 1.5
- 1 -1" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 3 -syntho xy-inn-$N.txt
-junk xy-inn-$N.txt
+N=60
+echo "1 -1
+ 1 3
+-3 -1
+ 1 -1" | ./lpu cbfit -i - -synthn $N -v 0 -syntho - | unu crop -min 0 0 -max M M-1 -o xy-inn-$N.txt
+# junk xy-inn-$N.txt
 
 CMD="./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 1 -psi 1000000000 -nim 4000 -deltathr 0.000000000001"
 echo "====== $CMD"

Modified: teem/trunk/src/limn/test/05-go.sh
===================================================================
--- teem/trunk/src/limn/test/05-go.sh	2024-07-12 19:23:37 UTC (rev 7207)
+++ teem/trunk/src/limn/test/05-go.sh	2024-07-16 18:50:21 UTC (rev 7208)
@@ -33,11 +33,12 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
+if false; then
 N=199
 echo "-1 -1
 -1 1
- 1 1
- 0 0" | ./lpu cbfit -i - -synthn $N  -syntho 0.txt
+1 1
+0 0" | ./lpu cbfit -i - -synthn $N  -syntho 0.txt
 echo "0 0
 -0.5 -0.5
 1 -1
@@ -44,12 +45,14 @@
 1 1" | ./lpu cbfit -i - -synthn $N  -syntho 1.txt
 cat 0.txt 1.txt | uniq > xy-inn.txt
 junk {0,1}.txt xy-inn.txt
+IN=xy-inn.txt
+fi
 
-IN=xy-inn.txt
+IN=xy-inn-60.txt
 #IN=pointy.txt
 #IN=circ.txt
 
-CMD="./lpu cbfit -i $IN -scl 0.5 -v 0 -psi 3 -eps 0.003 -roll 0"
+CMD="./lpu cbfit -i $IN -loop -scl 0.5 -v 3 -psi 3 -eps 0.01 -roll 1"
 echo "====== $CMD"
 eval $CMD > log.txt
 cat log.txt # ; junk log.txt
@@ -58,7 +61,7 @@
 echo "====== RESULTS: --> $OUT"
 grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9
 grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9 |
-    ./lpu cbfit -i - -synthn $((5*N)) -sup 1 -syntho $OUT
+    ./lpu cbfit -i - -loop -synthn 300 -sup 1 -syntho $OUT
 junk $OUT
 
 BIN=900

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

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

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

Revision: 7207
          http://sourceforge.net/p/teem/code/7207
Author:   kindlmann
Date:     2024-07-12 19:23:37 +0000 (Fri, 12 Jul 2024)
Log Message:
-----------
recognizing that having one corner in path is possible

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

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-12 19:15:08 UTC (rev 7206)
+++ teem/trunk/src/limn/splineFit.c	2024-07-12 19:23:37 UTC (rev 7207)
@@ -1892,13 +1892,12 @@
       return 1;
     }
   } else {
-    assert(fctx->cnum >= 2);
-    /* we do have corners: find segments between corners. The corner vertex part of two
-    segments: the last point in segment I and first point in segment I+1. How many
-    segments we analyze depends on whether they're in a loop: if they're in a loop
-    then we do one more */
-    uint cii, ciLast = fctx->cnum - 2 + !!lpnt->isLoop;
-    for (cii = 0; cii <= ciLast; cii++) {
+    /* we do have corners: find segments between corners, but maybe as few as one corner.
+    The corner vertex part of two segments: the last point in segment I and first point
+    in segment I+1. How many segments we analyze depends on whether they're in a loop: if
+    they're in a loop then we do one more */
+    uint cii, ciNum = fctx->cnum - 1 + !!lpnt->isLoop;
+    for (cii = 0; cii < ciNum; cii++) {
       uint cjj = (cii + 1) % fctx->cnum;
       limnCbfPath *subpath = limnCbfPathNew(0);
       /* 0: left tangent   2: vertex   4: right tangent */

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

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

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

Revision: 7206
          http://sourceforge.net/p/teem/code/7206
Author:   kindlmann
Date:     2024-07-12 19:15:08 +0000 (Fri, 12 Jul 2024)
Log Message:
-----------
now actually accepting more that type double

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

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-12 19:10:32 UTC (rev 7205)
+++ teem/trunk/src/limn/splineFit.c	2024-07-12 19:15:08 UTC (rev 7206)
@@ -179,9 +179,10 @@
     biffAddf(LIMN, "%s: point data type %d not valid", me, ptype);
     return NULL;
   }
-  if (ptype != nrrdTypeDouble) {
-    biffAddf(LIMN, "%s: sorry, only %s-type data implemented now (not %s)", me,
-             airEnumStr(nrrdType, nrrdTypeDouble), airEnumStr(nrrdType, ptype));
+  if (!(nrrdTypeDouble == ptype || nrrdTypeFloat == ptype)) {
+    biffAddf(LIMN, "%s: sorry, only %s or %s type data implemented now (not %s)", me,
+             airEnumStr(nrrdType, nrrdTypeDouble), airEnumStr(nrrdType, nrrdTypeFloat),
+             airEnumStr(nrrdType, ptype));
     return NULL;
   }
   if (2 != dim) {

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

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

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

Revision: 7205
          http://sourceforge.net/p/teem/code/7205
Author:   kindlmann
Date:     2024-07-12 19:10:32 +0000 (Fri, 12 Jul 2024)
Log Message:
-----------
usefully acting on float-type data

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

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-11 22:19:19 UTC (rev 7204)
+++ teem/trunk/src/limn/splineFit.c	2024-07-12 19:10:32 UTC (rev 7205)
@@ -197,8 +197,9 @@
     biffAddf(LIMN, "%s: couldn't allocate point container", me);
     return NULL;
   }
-  if (pdata) {
-    /* we are wrapping around a given pre-allocated buffer */
+  if (pdata && nrrdTypeDouble == ptype) {
+    /* we are wrapping around a given pre-allocated buffer
+       and it just happens to already be type double */
     lpnt->pp = pdata;
     lpnt->ppOwn = NULL;
   } else {
@@ -206,8 +207,24 @@
     lpnt->pp = NULL;
     if (!(lpnt->ppOwn = AIR_CALLOC(dim * pnum, double))) {
       biffAddf(LIMN, "%s: couldn't allocate %u %u-D point", me, pnum, dim);
+      free(lpnt);
       return NULL;
     }
+    if (pdata) {
+      uint ii, nn = dim * pnum;
+      if (nrrdTypeFloat == ptype) {
+        const float *fdata = (const float *)pdata;
+        for (ii = 0; ii < nn; ii++) {
+          lpnt->ppOwn[ii] = (double)fdata[ii];
+        }
+      } else {
+        biffAddf(LIMN, "%s: sorry, can't currently handle given %s type data", me,
+                 airEnumStr(nrrdType, ptype));
+        free(lpnt->ppOwn);
+        free(lpnt);
+        return NULL;
+      }
+    }
   }
   lpnt->num = pnum;
   lpnt->dim = dim; /* but really DIM=2 because of above */

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

[Teem-commit] SF.net SVN: teem:[7204] teem/trunk/python/cffi

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

Revision: 7204
          http://sourceforge.net/p/teem/code/7204
Author:   kindlmann
Date:     2024-07-11 22:19:19 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
synching with Teem source

Modified Paths:
--------------
    teem/trunk/python/cffi/biffdata/limn.csv
    teem/trunk/python/cffi/cdef/cdef_air.h
    teem/trunk/python/cffi/cdef/cdef_limn.h
    teem/trunk/python/cffi/teem.py

Modified: teem/trunk/python/cffi/biffdata/limn.csv
===================================================================
--- teem/trunk/python/cffi/biffdata/limn.csv	2024-07-11 22:07:25 UTC (rev 7203)
+++ teem/trunk/python/cffi/biffdata/limn.csv	2024-07-11 22:19:19 UTC (rev 7204)
@@ -60,14 +60,14 @@
 limnSplineUpdate,int,1,0,limn,limn/splineMethods.c:422
 limnSplineTypeSpecParse,limnSplineTypeSpec *,NULL,0,limn,limn/splineMisc.c:222
 limnSplineParse,limnSpline *,NULL,0,limn,limn/splineMisc.c:278
-limnCBFPointsCheck,int,1,0,limn,limn/splineFit.c:73
-limnCBFCtxNew,limnCBFCtx *,NULL,0,limn,limn/splineFit.c:222
-limnCBFFindVT,int,1,0,limn,limn/splineFit.c:357
-limnCBFCtxCheck,int,1,0,limn,limn/splineFit.c:751
-limnCBFitSingle,int,1,0,limn,limn/splineFit.c:925
-limnCBFMulti,int,1,0,limn,limn/splineFit.c:1016
-limnCBFCorners,int,1,0,limn,limn/splineFit.c:1116
-limnCBFit,int,1,0,limn,limn/splineFit.c:1184
+limnCbfPointsNew,limnCbfPoints *,NULL,0,limn,limn/splineFit.c:175
+limnCbfPointsCheck,int,1,0,limn,limn/splineFit.c:229
+limnCbfCtxPrep,int,1,0,limn,limn/splineFit.c:497
+limnCbfTVT,int,1,0,limn,limn/splineFit.c:762
+limnCbfSingle,int,1,0,limn,limn/splineFit.c:1494
+limnCbfCorners,int,1,0,limn,limn/splineFit.c:1545
+limnCbfMulti,int,1,0,limn,limn/splineFit.c:1727
+limnCbfGo,int,1,0,limn,limn/splineFit.c:1844
 limnObjectWorldHomog,int,1,0,limn,limn/transform.c:25
 limnObjectFaceNormals,int,1,0,limn,limn/transform.c:47
 limnObjectSpaceTransform,int,1,0,limn,limn/transform.c:210

Modified: teem/trunk/python/cffi/cdef/cdef_air.h
===================================================================
--- teem/trunk/python/cffi/cdef/cdef_air.h	2024-07-11 22:07:25 UTC (rev 7203)
+++ teem/trunk/python/cffi/cdef/cdef_air.h	2024-07-11 22:19:19 UTC (rev 7204)
@@ -764,13 +764,14 @@
 /*
 ******** AIR_MOD(i, N)
 **
-** returns that integer in [0, N-1] which is i plus a multiple of N. It
-** may be unfortunate that the expression (i)%(N) appears three times;
-** this should be inlined.  Or perhaps the compiler's optimizations
-** (common sub-expression elimination) will save us.
+** returns that integer in [0, N-1] which is i plus a multiple of N. It may be
+** unfortunate that the expression (i)%(N) appears three times; we can hope the compiler
+** optimizes the common sub-expression.
 **
-** Note: integer divisions are not very fast on some modern chips;
-** don't go silly using this one.
+** NOTE: !! Results will not be as expected if i and N differ in sign-ed-ness !!
+**
+** Note: integer divisions are not very fast on some modern chips; don't go silly using
+** this one.
 */
 /*
 ******** AIR_LERP(w, a, b)

Modified: teem/trunk/python/cffi/cdef/cdef_limn.h
===================================================================
--- teem/trunk/python/cffi/cdef/cdef_limn.h	2024-07-11 22:07:25 UTC (rev 7203)
+++ teem/trunk/python/cffi/cdef/cdef_limn.h	2024-07-11 22:19:19 UTC (rev 7204)
@@ -448,112 +448,131 @@
   double B, C; /* B,C values for BC-splines */
 } limnSplineTypeSpec;
 /*
-******** limnCBFSeg
+******** limnCbfSeg
 **
-** how one cubic Bezier spline segment is represented for limnCBF functions
+** how one cubic Bezier spline segment is represented for limnCbf functions
 ** (using DIM=2 to mark places where the 2D-ness of the code surfaces )
+**
+** No constructor (New) or destructor (Nix) functions since it is so simple
 */
 typedef struct {
-  double xy[8];      /* four control points of cubic Bezier:
-                        x0, y0,   x1, y1,   x2, y2,   x3, y3
-                        0   1     2   3     4   5     6  7   DIM=2 */
-  int corner[2];     /* corner[0,1] non-zero if xy[0,3] are corner vertices;
-                        segments otherwise assumed geometrically continuous */
-  unsigned int pNum; /* (if non-zero) this segment approximates pNum points */
-} limnCBFSeg;
+  double xy[8];  /* four control points of cubic Bezier:
+                    x0, y0,   x1, y1,   x2, y2,   x3, y3
+                    0   1     2   3     4   5     6  7   DIM=2 */
+  int corner[2]; /* corner[0,1] non-zero if xy[0,3] are either corner vertices
+                    or path-ending vertices, i.e. reasons to not have geometric
+                    continuity here, either because we intend to have a corner,
+                    or because there's nothing else to be continuous with */
+                 /* how many points does this represent */
+  unsigned int pointNum;
+} limnCbfSeg;
 /*
-******** limnCBFPath
+******** limnCbfPath
 **
-** a multi-segment path in the context of cubic Bezier fitting
+** a multi-spline path in the context of cubic Bezier fitting
 */
 typedef struct {
-  limnCBFSeg *seg;     /* array of limnCBFSeg */
+  limnCbfSeg *seg;     /* array of limnCbfSeg structs (NOT pointers to them) */
   unsigned int segNum; /* length of seg array */
   airArray *segArr;    /* manages seg and segNum */
   int isLoop;          /* path is closed loop */
-} limnCBFPath;
+} limnCbfPath;
 /*
-******** limnCBFCtx
+******** limnCbfCtx
 **
-** The bag of state for limnCBF functions.
+** The complete bag of input parameters and state for limnCbf functions, especially the
+** top-level limnCbfGo function.
 **
 ** note: "nrp" = Newton-based Re-Parameterization of where the given points
-** fall along the spline, the iterative process inside limnCBFSingle
+** fall along the spline, the iterative process inside limnCbfSingle
 */
 typedef struct {
   /* ----------- input ---------- */
-  int verbose,             /* verbosity level */
-    cornNMS;               /* non-minimal-suppression of corners: accept as
-                              corners only those with locally minimal angle */
+  int verbose,  /* verbosity level */
+    cornerFind, /* do first search for corners: places where the path is not
+                   geometrically continuous (between corners, the path is geometrically
+                   continuous between multiple spline segments) */
+    cornerNMS;  /* (if cornerFind) non-minimal-suppression of corners: accept as
+                   corners only those with locally minimal angle */
   unsigned int nrpIterMax; /* max # iters of nrp */
-  double scale,  /* scale (in sense of nrrdKernelDiscreteGaussian) at which to estimate
-                    spline endpoints and tangents; scale=0 means the endpoints are
-                    exactly on vertices, and tangents are from the smallest-support
-                    finite differences. This is the ONLY floating point that should be set
-                    by a method (limnCBFScaleSet); the rest can be set directly. */
-    distMin,     /* min distance to given points: this controls both splitting done by
-                    limnCBFMulti, and nrp within limnCBFSingle */
-    nrpDeltaMax, /* in nrp, capping parameterization change to this scaling of average
-                    u[i+1]-u[i]. This wasn't in author's original code (so their idea of
-                    doing at most ~5 iters of nrp may no longer hold), but it can help
-                    stabilize things */
-    nrpDistScl,  /* scaling on distMin to use when testing distance during nrp; setting
-                    this < 1 means that nrp tries to be more stringent than the overall
-                    fitting, but with the benefit of sometimes being smarter about where
-                    to split, when that is needed */
-    nrpPsi, /* don't even try nrp if max dist is bigger than nrpPsi*distMin, instead just
-               subdivide */
-    nrpDeltaMin, /* min total parameterization change by nrp */
-    alphaMin,    /* alpha can't be negative, and we enforce distinct positivity to ensure
-                    that spline doesn't slow down too much near endpoints */
-    detMin,      /* abs(determinant) of 2x2 matrix to invert can't go below this  */
-    cornAngle; /* angle, in degrees, between (one-sided) incoming and outgoing tangents,
-                  *below* which a vertex should be considered a corner. Vertices in a
-                  straight line have an angle of 180 degrees. Or, if 0, no effort is made
-                  to detect corners. */
+  double
+    epsilon, /* error threshold on min distance from spline (as currently parameterized)
+                to given points: this affects both splitting done by limnCbfMulti, and
+                nrp within limnCbfSingle. Fitting has successfully finished when spline
+                path never strays further than this from input points */
+    scale,   /* scale (in sense of nrrdKernelDiscreteGaussian) at which to estimate
+                spline endpoints and tangents; scale=0 means the endpoints are
+                exactly on vertices, and tangents are from the smallest-support
+                finite differences */
+    nrpCap,  /* in nrp, cap parameterization change to this scaling of average
+                u[i+1]-u[i]. This wasn't in author's original code (so their idea of doing
+                at most ~5 iters of nrp no longer holds), but it can help stabilize things
+                with gnarly inputs */
+    nrpIota, /* (also not in author's original code) scaling on epsilon to use when
+                testing distance during nrp; setting this < 1 means that nrp tries to be
+                more stringent than the overall fitting, but with the benefit of
+                sometimes being smarter about where to split, when that is needed */
+    nrpPsi,  /* don't even try nrp if max dist is bigger than nrpPsi*epsilon, instead
+                just subdivide */
+    nrpDeltaThresh, /* finish npr when mean parameterization change fall below this */
+    alphaMin,   /* alpha can't be negative, and we enforce distinct positivity to ensure
+                   that spline doesn't slow down too much near endpoints */
+    detMin,     /* abs(determinant) of 2x2 matrix to invert can't go below this  */
+    cornAngle,  /* interior angle, in degrees, between (one-sided) incoming and outgoing
+                   tangents, *below* which a vertex should be considered a corner.
+                   Vertices in a straight line have an angle of 180 degrees. */
+    wackyAngle; /* in cases where we are only looking at three points: a spline can
+                   always be fit through the middle point, even with constraints on
+                   position and tangent at first and last points, but the spline looks
+                   wacky if its tangent at the middle point is wildly different than
+                   the (two-sided) tangent that would have been estimated at that point
+                   for the purpose of splitting. If the angle (in degrees) between the
+                   two tangents exceeds this, then fitting will generate the simple
+                   (punted) arc, which will likely trigger splitting. */
   /* ----------- internal --------- */
-  double *uu,        /* buffer used for nrp */
+  double *uu,        /* used for nrp: buffer of parameterizations in [0,1] of point along
+                        currently considered spline segment */
     *vw,             /* weights for endpoint vertex calculation */
     *tw,             /* weights for endpoint tangent calculation */
-    *mine;           /* helps remember who allocated the above */
-  unsigned int wLen; /* how long are vw, tw */
-  double lenF2L;     /* length of segment from first to last */
+    *ctvt;           /* corner info: 6-by-cNum (DIM=2) of tangent,vertex,tangent */
+  unsigned int ulen, /* how long is uu */
+    wlen,            /* how long are vw, tw */
+    *cidx,           /* indices (into lpnt point data) of corners */
+    cnum;            /* number of corners described by ctvt and cidx */
   /* ----------- output --------- */
   unsigned int nrpIterDone, /* number of nrp iters taken */
-    distIdx;                /* which point had distance distDone */
-  double dist,              /* max distance to given points */
-    nrpDeltaDone,           /* latest total parameterization change by nrp */
+    distMaxIdx,             /* which point had distance distMax */
+    nrpPuntFlop;            /* # times that single-spline fit flip-flopped between a
+                               well-computed spline vs a punted one */
+  double distMax,           /* max distance to given points */
+    nrpDeltaDone,           /* latest mean parameterization change by nrp */
     alphaDet;               /* min det of matrix inverted to find alpha */
-  int distBig;              /* how big dist (above) is:
-                               0: dist <= nD
-                               1: nD < dist <= DM
-                               2: DM < dist <= fD
-                               3: fD < dist
-                               where
-                               nD = nrpDistScl*distMin,
-                               DM = distMin,
-                               fD = nrpPsi*distMin: */
-} limnCBFCtx;
+  int distBig;              /* how big distMax (above) is:
+                               0: distMax <= wee        (wee = nrpIota*epsilon)
+                               1: wee < distMax <= eps  (eps = epsilon)
+                               2: eps < distMax <= far  (far = nrpPsi*epsilon)
+                               3: far < distMax */
+} limnCbfCtx;
 /*
-******** limnCBFPoints
+******** limnCbfPoints
 **
-** a container for 1D array of points; currently used for limnCBF functions
+** a container for 1D array of points; currently used for limnCbf functions
 ** Both pp and ppOwn can point to the array of point locations, but exactly
 ** one of pp and ppOwn can be non-NULL.
 **
 ** NOTE: For now, point data is only double (not float), and only in 2D (not
 ** 3D), but if this becomes more general, that generality will be inside here.
-** For time being DIM=2 tags locations where 2D-ness is explicit in code.
+** For time being "DIM=2" tags locations where 2D-ness is explicit in code.
 */
 typedef struct {
-  /* assuming DIM=2: 2 values per logical element pp */
   const double *pp; /* point coords, we do not own buffer */
   double *ppOwn;    /* point coords, we DO own buffer */
-  unsigned int num; /* how many points */
-  int isLoop;       /* points form a loop: logical indices into coord
-                       array are . . . num-2, num-1, 0, 1, . . .
-                       and index 0 is effectively arbitrary */
-} limnCBFPoints;
+  unsigned int num, /* how many points */
+    dim;      /* points live in what dimension (currently only dim=2 implemented) */
+  int isLoop; /* points form a loop: logical indices into coord
+                 array are . . . num-2, num-1, 0, 1, . . .
+                 and index 0 is effectively arbitrary */
+} limnCbfPoints;
 /* defaultsLimn.c */
 extern const int limnPresent;
 extern const char *const limnBiffKey;
@@ -811,33 +830,36 @@
 extern int limnSplineSample(Nrrd *nout, limnSpline *spline, double minT, size_t M,
                                  double maxT);
 /* splineFit.c */
-extern limnCBFPoints *limnCBFPointsNew(const double *pp, unsigned int nn,
+extern limnCbfPoints *limnCbfPointsNew(const void *pdata, int ptype,
+                                            unsigned int dim, unsigned int pnum,
                                             int isLoop);
-extern limnCBFPoints *limnCBFPointsNix(limnCBFPoints *lpnt);
-extern int limnCBFPointsCheck(const limnCBFPoints *lpnt);
-extern limnCBFCtx *limnCBFCtxNew(unsigned int pointNum, double scale);
-extern limnCBFCtx *limnCBFCtxNix(limnCBFCtx *fctx);
-extern void limnCBFSegEval(double *xy, const limnCBFSeg *seg, double tt);
-extern limnCBFPath *limnCBFPathNew(void);
-extern limnCBFPath *limnCBFPathNix(limnCBFPath *path);
-extern void limnCBFPathSample(double *xy, unsigned int pointNum,
-                                   const limnCBFPath *path);
-extern int limnCBFFindVT(double vv[2], double tt[2], const limnCBFCtx *fctx,
-                              const limnCBFPoints *lpnt, unsigned int loi,
-                              unsigned int hii, unsigned int ofi, int dir);
-extern int limnCBFCtxCheck(const limnCBFCtx *fctx, const limnCBFPoints *lpnt);
-extern int limnCBFitSingle(double alpha[2], limnCBFCtx *fctx, const double vv0[2],
-                                const double tt1[2], const double tt2[2],
-                                const double vv3[2], const double *xy,
-                                unsigned int pointNum, int isLoop);
-extern int limnCBFMulti(limnCBFPath *path, limnCBFCtx *fctx, const double vv0[2],
-                             const double tt1[2], const double tt2[2],
-                             const double vv3[2], const limnCBFPoints *lpnt,
-                             unsigned int loi, unsigned int hii);
-extern int limnCBFCorners(unsigned int **cornIdx, unsigned int *cornNum,
-                               limnCBFCtx *fctx, const limnCBFPoints *lpnt);
-extern int limnCBFit(limnCBFPath *path, limnCBFCtx *fctx, const double *xy,
-                          unsigned int pointNum, int isLoop);
+extern limnCbfPoints *limnCbfPointsNix(limnCbfPoints *lpnt);
+extern int limnCbfPointsCheck(const limnCbfPoints *lpnt);
+extern limnCbfPath *limnCbfPathNew(unsigned segNum);
+extern limnCbfPath *limnCbfPathNix(limnCbfPath *path);
+extern void limnCbfPathJoin(limnCbfPath *dst, const limnCbfPath *src);
+extern limnCbfCtx *limnCbfCtxNew();
+extern limnCbfCtx *limnCbfCtxNix(limnCbfCtx *fctx);
+extern int limnCbfCtxPrep(limnCbfCtx *fctx, const limnCbfPoints *lpnt);
+extern void limnCbfSegEval(double *xy, const limnCbfSeg *seg, double tt);
+extern void limnCbfPathSample(double *xy, unsigned int pointNum,
+                                   const limnCbfPath *path);
+extern int limnCbfTVT(double lt[2], double vv[2], double rt[2],
+                           const limnCbfCtx *fctx, const limnCbfPoints *lpnt,
+                           unsigned int loi, unsigned int hii, unsigned int vvi,
+                           int oneSided);
+extern int limnCbfSingle(limnCbfSeg *seg, const double vv0[2], const double tt1[2],
+                              const double tt2[2], const double vv3[2], limnCbfCtx *fctx,
+                              const limnCbfPoints *lpnt, unsigned int loi,
+                              unsigned int hii);
+extern int limnCbfCorners(limnCbfCtx *fctx, const limnCbfPoints *lpnt);
+extern int limnCbfMulti(limnCbfPath *path, const double vv0[2], const double tt1[2],
+                             const double tt2[2], const double vv3[2],
+                             unsigned int recurseDepth, limnCbfCtx *fctx,
+                             const limnCbfPoints *lpnt, unsigned int loi,
+                             unsigned int hii);
+extern int limnCbfGo(limnCbfPath *path, limnCbfCtx *fctx,
+                          const limnCbfPoints *lpnt);
 /* lpu{Flotsam,. . .}.c */
 /* F(clip) \ */
 /* F(vwflip) \ */

Modified: teem/trunk/python/cffi/teem.py
===================================================================
--- teem/trunk/python/cffi/teem.py	2024-07-11 22:07:25 UTC (rev 7203)
+++ teem/trunk/python/cffi/teem.py	2024-07-11 22:19:19 UTC (rev 7204)
@@ -293,6 +293,14 @@
     'mossSamplerSample': (_equals_one, 0, b'moss', 'moss/sampler.c:195'),
     'mossLinearTransform': (_equals_one, 0, b'moss', 'moss/xform.c:140'),
     'mossFourPointTransform': (_equals_one, 0, b'moss', 'moss/xform.c:219'),
+    'alanUpdate': (_equals_one, 0, b'alan', 'alan/coreAlan.c:60'),
+    'alanInit': (_equals_one, 0, b'alan', 'alan/coreAlan.c:99'),
+    'alanRun': (_equals_one, 0, b'alan', 'alan/coreAlan.c:453'),
+    'alanDimensionSet': (_equals_one, 0, b'alan', 'alan/methodsAlan.c:104'),
+    'alan2DSizeSet': (_equals_one, 0, b'alan', 'alan/methodsAlan.c:119'),
+    'alan3DSizeSet': (_equals_one, 0, b'alan', 'alan/methodsAlan.c:139'),
+    'alanTensorSet': (_equals_one, 0, b'alan', 'alan/methodsAlan.c:161'),
+    'alanParmSet': (_equals_one, 0, b'alan', 'alan/methodsAlan.c:208'),
     'gageContextCopy': (_equals_null, 0, b'gage', 'gage/ctx.c:88'),
     'gageKernelSet': (_equals_one, 0, b'gage', 'gage/ctx.c:199'),
     'gagePerVolumeAttach': (_equals_one, 0, b'gage', 'gage/ctx.c:398'),
@@ -339,6 +347,31 @@
     'gageOptimSigErrorPlotSliding': (_equals_one, 0, b'gage', 'gage/optimsig.c:1253'),
     'dyeConvert': (_equals_one, 0, b'dye', 'dye/convertDye.c:351'),
     'dyeColorParse': (_equals_one, 0, b'dye', 'dye/methodsDye.c:185'),
+    'baneClipNew': (_equals_null, 0, b'bane', 'bane/clip.c:102'),
+    'baneClipAnswer': (_equals_one, 0, b'bane', 'bane/clip.c:152'),
+    'baneClipCopy': (_equals_null, 0, b'bane', 'bane/clip.c:167'),
+    'baneFindInclusion': (_equals_one, 0, b'bane', 'bane/hvol.c:87'),
+    'baneMakeHVol': (_equals_one, 0, b'bane', 'bane/hvol.c:248'),
+    'baneGKMSHVol': (_equals_null, 0, b'bane', 'bane/hvol.c:447'),
+    'baneIncNew': (_equals_null, 0, b'bane', 'bane/inc.c:251'),
+    'baneIncAnswer': (_equals_one, 0, b'bane', 'bane/inc.c:360'),
+    'baneIncCopy': (_equals_null, 0, b'bane', 'bane/inc.c:375'),
+    'baneMeasrNew': (_equals_null, 0, b'bane', 'bane/measr.c:33'),
+    'baneMeasrCopy': (_equals_null, 0, b'bane', 'bane/measr.c:149'),
+    'baneRangeNew': (_equals_null, 0, b'bane', 'bane/rangeBane.c:89'),
+    'baneRangeCopy': (_equals_null, 0, b'bane', 'bane/rangeBane.c:130'),
+    'baneRangeAnswer': (_equals_one, 0, b'bane', 'bane/rangeBane.c:144'),
+    'baneRawScatterplots': (_equals_one, 0, b'bane', 'bane/scat.c:26'),
+    'baneOpacInfo': (_equals_one, 0, b'bane', 'bane/trnsf.c:29'),
+    'bane1DOpacInfoFrom2D': (_equals_one, 0, b'bane', 'bane/trnsf.c:144'),
+    'baneSigmaCalc': (_equals_one, 0, b'bane', 'bane/trnsf.c:222'),
+    'banePosCalc': (_equals_one, 0, b'bane', 'bane/trnsf.c:253'),
+    'baneOpacCalc': (_equals_one, 0, b'bane', 'bane/trnsf.c:403'),
+    'baneInputCheck': (_equals_one, 0, b'bane', 'bane/valid.c:26'),
+    'baneHVolCheck': (_equals_one, 0, b'bane', 'bane/valid.c:64'),
+    'baneInfoCheck': (_equals_one, 0, b'bane', 'bane/valid.c:106'),
+    'banePosCheck': (_equals_one, 0, b'bane', 'bane/valid.c:144'),
+    'baneBcptsCheck': (_equals_one, 0, b'bane', 'bane/valid.c:179'),
     'limnCameraUpdate': (_equals_one, 0, b'limn', 'limn/cam.c:33'),
     'limnCameraAspectSet': (_equals_one, 0, b'limn', 'limn/cam.c:130'),
     'limnCameraPathMake': (_equals_one, 0, b'limn', 'limn/cam.c:189'),
@@ -400,14 +433,14 @@
     'limnSplineUpdate': (_equals_one, 0, b'limn', 'limn/splineMethods.c:422'),
     'limnSplineTypeSpecParse': (_equals_null, 0, b'limn', 'limn/splineMisc.c:222'),
     'limnSplineParse': (_equals_null, 0, b'limn', 'limn/splineMisc.c:278'),
-    'limnCBFPointsCheck': (_equals_one, 0, b'limn', 'limn/splineFit.c:73'),
-    'limnCBFCtxNew': (_equals_null, 0, b'limn', 'limn/splineFit.c:222'),
-    'limnCBFFindVT': (_equals_one, 0, b'limn', 'limn/splineFit.c:357'),
-    'limnCBFCtxCheck': (_equals_one, 0, b'limn', 'limn/splineFit.c:751'),
-    'limnCBFitSingle': (_equals_one, 0, b'limn', 'limn/splineFit.c:925'),
-    'limnCBFMulti': (_equals_one, 0, b'limn', 'limn/splineFit.c:1016'),
-    'limnCBFCorners': (_equals_one, 0, b'limn', 'limn/splineFit.c:1116'),
-    'limnCBFit': (_equals_one, 0, b'limn', 'limn/splineFit.c:1184'),
+    'limnCbfPointsNew': (_equals_null, 0, b'limn', 'limn/splineFit.c:175'),
+    'limnCbfPointsCheck': (_equals_one, 0, b'limn', 'limn/splineFit.c:229'),
+    'limnCbfCtxPrep': (_equals_one, 0, b'limn', 'limn/splineFit.c:497'),
+    'limnCbfTVT': (_equals_one, 0, b'limn', 'limn/splineFit.c:762'),
+    'limnCbfSingle': (_equals_one, 0, b'limn', 'limn/splineFit.c:1494'),
+    'limnCbfCorners': (_equals_one, 0, b'limn', 'limn/splineFit.c:1545'),
+    'limnCbfMulti': (_equals_one, 0, b'limn', 'limn/splineFit.c:1727'),
+    'limnCbfGo': (_equals_one, 0, b'limn', 'limn/splineFit.c:1844'),
     'limnObjectWorldHomog': (_equals_one, 0, b'limn', 'limn/transform.c:25'),
     'limnObjectFaceNormals': (_equals_one, 0, b'limn', 'limn/transform.c:47'),
     'limnObjectSpaceTransform': (_equals_one, 0, b'limn', 'limn/transform.c:210'),
@@ -582,6 +615,21 @@
     'pullTraceMultiPlotAdd': (_equals_one, 0, b'pull', 'pull/trace.c:704'),
     'pullTraceMultiWrite': (_equals_one, 0, b'pull', 'pull/trace.c:1014'),
     'pullTraceMultiRead': (_equals_one, 0, b'pull', 'pull/trace.c:1119'),
+    'coilStart': (_equals_one, 0, b'coil', 'coil/coreCoil.c:287'),
+    'coilIterate': (_equals_one, 0, b'coil', 'coil/coreCoil.c:362'),
+    'coilFinish': (_equals_one, 0, b'coil', 'coil/coreCoil.c:407'),
+    'coilVolumeCheck': (_equals_one, 0, b'coil', 'coil/methodsCoil.c:25'),
+    'coilContextAllSet': (_equals_one, 0, b'coil', 'coil/methodsCoil.c:69'),
+    'coilOutputGet': (_equals_one, 0, b'coil', 'coil/methodsCoil.c:200'),
+    'pushOutputGet': (_equals_one, 0, b'push', 'push/action.c:71'),
+    'pushBinProcess': (_equals_one, 0, b'push', 'push/action.c:161'),
+    'pushBinPointAdd': (_equals_one, 0, b'push', 'push/binning.c:180'),
+    'pushRebin': (_equals_one, 0, b'push', 'push/binning.c:197'),
+    'pushStart': (_equals_one, 0, b'push', 'push/corePush.c:183'),
+    'pushIterate': (_equals_one, 0, b'push', 'push/corePush.c:233'),
+    'pushRun': (_equals_one, 0, b'push', 'push/corePush.c:306'),
+    'pushFinish': (_equals_one, 0, b'push', 'push/corePush.c:396'),
+    'pushEnergySpecParse': (_equals_one, 0, b'push', 'push/forces.c:304'),
     'miteSample': (_math.isnan, 0, b'mite', 'mite/ray.c:151'),
     'miteRenderBegin': (_equals_one, 0, b'mite', 'mite/renderMite.c:63'),
     'miteShadeSpecParse': (_equals_one, 0, b'mite', 'mite/shade.c:69'),

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

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

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

Revision: 7203
          http://sourceforge.net/p/teem/code/7203
Author:   kindlmann
Date:     2024-07-11 22:07:25 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
some final comments

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

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-11 21:53:07 UTC (rev 7202)
+++ teem/trunk/src/limn/splineFit.c	2024-07-11 22:07:25 UTC (rev 7203)
@@ -30,9 +30,16 @@
 https://dl.acm.org/doi/10.5555/90767.90941
 The author's code is here:
 http://www.realtimerendering.com/resources/GraphicsGems/gems/FitCurves.c
-but the code here was based more on reading the paper, than their code. Also, that code
-does not handle point loops, and does not handle smoothing built into tangent estimation,
-which were important to GLK, but which added significant implementation complexity.
+but the code here was based more on reading the paper, than the author's code.
+Beyond the author's paper and the author's code, this code here:
+- handles closed point loops
+- implements smoothing as part of vertex and tangent estimation
+- the Newton-based ReParameterization ("nrp") is smarter: it never increases
+  the distance between the data and splines, so it is better at handling cases
+  where the chord-length-based parameterization initialization is terrible
+- is smarter about handling single-spline fits to only 3 points
+- has robust error handling and error reporting
+All of this adds implementation complexity (this is ~4 times longer than author's file)
 
 The functions below do not use any other limnSpline structs or functions, since those
 were written a long time ago when GLK was even more ignorant than now about splines.

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

[Teem-commit] SF.net SVN: teem:[7202] teem/trunk/src/limn

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

Revision: 7202
          http://sourceforge.net/p/teem/code/7202
Author:   kindlmann
Date:     2024-07-11 21:53:07 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
added wacky path handling for 3=spanlen, and maybe done

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

Modified: teem/trunk/src/limn/limn.h
===================================================================
--- teem/trunk/src/limn/limn.h	2024-07-11 21:48:45 UTC (rev 7201)
+++ teem/trunk/src/limn/limn.h	2024-07-11 21:53:07 UTC (rev 7202)
@@ -580,12 +580,20 @@
     nrpPsi,  /* don't even try nrp if max dist is bigger than nrpPsi*epsilon, instead
                 just subdivide */
     nrpDeltaThresh, /* finish npr when mean parameterization change fall below this */
-    alphaMin,  /* alpha can't be negative, and we enforce distinct positivity to ensure
-                  that spline doesn't slow down too much near endpoints */
-    detMin,    /* abs(determinant) of 2x2 matrix to invert can't go below this  */
-    cornAngle; /* interior angle, in degrees, between (one-sided) incoming and outgoing
-                  tangents, *below* which a vertex should be considered a corner.
-                  Vertices in a straight line have an angle of 180 degrees. */
+    alphaMin,   /* alpha can't be negative, and we enforce distinct positivity to ensure
+                   that spline doesn't slow down too much near endpoints */
+    detMin,     /* abs(determinant) of 2x2 matrix to invert can't go below this  */
+    cornAngle,  /* interior angle, in degrees, between (one-sided) incoming and outgoing
+                   tangents, *below* which a vertex should be considered a corner.
+                   Vertices in a straight line have an angle of 180 degrees. */
+    wackyAngle; /* in cases where we are only looking at three points: a spline can
+                   always be fit through the middle point, even with constraints on
+                   position and tangent at first and last points, but the spline looks
+                   wacky if its tangent at the middle point is wildly different than
+                   the (two-sided) tangent that would have been estimated at that point
+                   for the purpose of splitting. If the angle (in degrees) between the
+                   two tangents exceeds this, then fitting will generate the simple
+                   (punted) arc, which will likely trigger splitting. */
   /* ----------- internal --------- */
   double *uu,        /* used for nrp: buffer of parameterizations in [0,1] of point along
                         currently considered spline segment */

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-11 21:48:45 UTC (rev 7201)
+++ teem/trunk/src/limn/splineFit.c	2024-07-11 21:53:07 UTC (rev 7202)
@@ -39,6 +39,24 @@
 Hopefully that older code can be revisited and re-organized for a later version of
 Teem, at which point the code below can be integrated with it.
 
+Terminology as it is used here:
+- points: the input data, a sequence of position coordinate vectors. May be in a (closed)
+  loop (successor of last point is first point), or may in an open path.
+- left,right tangent at a point: the tangents estimated at a point (by looking at
+  neighboring points) pointing away from the point towards points with lower,higher
+  indices, respectively. Tangent estimation can be one-sided (looking at half the
+  neighbors of the point, with only lower or higher indies), or two-sided (looking at
+  neighbors on both sides of the point).
+- corner: a point at which the angle between the one-sided tangents is below some
+  threshold (limnCbfCtx->cornAngle), or, sometimes for completeness, the first and
+  last points in an open path.
+- span: a closed interval of point indices (typically "[loi,hii]"), within which
+  computations are done, often delimited by corners
+- segment: like limnCbfSeg: one cubic Bezier spline. A "single fit" (via fitSingle or
+  limnCbfSingle) finds one segment to approximate some span of points
+- path: a sequence of segments. A "multi fit" (via limnCbfMulti or limnCbfGo) finds
+  a path to approximate some span of points.
+
 NOTE: spline fitting would be useful in 3D (or higher dimensions) too, but currently
 this code only supports 2D. "DIM=2" flags places in code where that is explicit, with
 the hope that this code can later be generalized.
@@ -47,16 +65,17 @@
 initials got annoying with other CamelCase function names.
 
 This code was unusually slow for GLK to write because he struggled to gracefully handle
-the combination of possibilities:
-  - points may trace an open segment, or a closed loop
+the 16-fold combination of possibilities:
+  - points may trace an open path, or a closed loop
   - geometry (spline endpoints and tangents) is computed without or with some smoothing.
+  - tangents can be computed either one-sided or two-sided
   - smoothing may be bounded by the two corners on either side, or smoothing
     may cover all available vertices (which may or may not loop around)
-How to represent an interval of vertex *indices* -- termed a "span", in contrast to the
-single cubic spline "segment" -- is fundamental to all of this. A uint span[2] array was
-tried, but ended up just passing around a lot of loi,hii (low index, high index) pairs to
-represent the span.
 
+How to represent an index interval (a span, in contrast to the single cubic spline
+segment) is fundamental to all of this. A uint span[2] array was tried, but ended up just
+passing around a lot of loi,hii (low index, high index) pairs to represent the span.
+
 In the presence of a loop of N points, we could imagine different kinds of indices:
           0 .. N-1          : "actual" indices
  .. -2 -1 0 .. N-1 N N-1 .. : "lifted" indices, in the sense that the real number line
@@ -109,9 +128,6 @@
 
 /*
 TODO:
-test findDist - is distMaxIdx the correct actual index?
-testing corners: corners near and at start==stop of isLoop
-corners not at start or stop of isLoop: do spline wrap around from last to first index?
 
 valgrind everything
 
@@ -118,6 +134,12 @@
 (DIM=2) explore what would be required to generalize from 2D to 3D,
 perhaps at least at the API level, even if 3D is not yet implemented
 
+use performance tests to explore optimal settings in fctx:
+  nrpIterMax, nrpCap, nrpIota, nrpPsi, nrpDeltaThresh
+evaluated in terms of time and #splines needed for fit
+(may want to pay in time for more economical representation)
+
+Older idea:
 The initialization of uu[] by arc length is especially wrong when the tangents are
 pointing towards each other, and then newton iterations can get stuck in a local minimum.
 As a specific example: with these control points: (-0.5,0.5) (2,0.5) (-0.5,0) (0.5,-0.5)
@@ -127,11 +149,9 @@
 requires exploring what is at least a 4-D space of possible splines (lowered from 8-D
 by symmetries).  The cost of not doing this is less economical representations, because
 we would split these segments needlessly.
-
-use performance tests to explore optimal settings in fctx:
-  nrpIterMax, nrpCap, nrpIota, nrpPsi, nrpDeltaThresh
-evaluated in terms of time and #splines needed for fit
-(may want to pay in time for more economical representation)
+But now that delta_t is smarter:
+Even a bad initial (arc-length) parameterization can be improved, so now with the case
+described here, N=18 and N=19 do equally well (yay!)
 */
 
 #define PNMIN(ISLOOP) ((ISLOOP) ? 4 : 3)
@@ -275,7 +295,7 @@
 static void
 ctxInit(limnCbfCtx *fctx) {
   if (!fctx) return;
-  /* defaults for input parameters to various Cbf functions */
+  /* ----- defaults for input parameters to various Cbf functions */
   fctx->verbose = 0;
   fctx->cornerFind = AIR_TRUE;
   fctx->cornerNMS = AIR_TRUE;
@@ -289,12 +309,13 @@
   fctx->alphaMin = 0.001;
   fctx->detMin = 0.01;
   fctx->cornAngle = 120.0; /* degrees */
-  /* internal state */
+  fctx->wackyAngle = 30.0; /* degrees */
+  /* ----- internal state ----- */
   /* initialize buffer pointers to NULL and buffer lengths to 0 */
   fctx->uu = fctx->vw = fctx->tw = fctx->ctvt = NULL;
   fctx->cidx = NULL;
   fctx->ulen = fctx->wlen = fctx->cnum = 0;
-  /* initialize outputs to bogus valus */
+  /* ----- initialize outputs to bogus valus */
   fctx->nrpIterDone = UINT_MAX;
   fctx->distMaxIdx = UINT_MAX;
   fctx->nrpPuntFlop = UINT_MAX;
@@ -577,15 +598,6 @@
   double ww[4];
   const double *xy = seg->xy;
   CBD0(vv, xy + 0, xy + 2, xy + 4, xy + 6, tt, ww); /* DIM=2 */
-  /*
-  fprintf(stderr, "!%s: tt=%g -> ww={%g,%g,%g,%g} * "
-          "{(%g,%g),(%g,%g),(%g,%g),(%g,%g)} = (%g,%g)\n",
-          "limnCbfSegEval", tt, ww[0], ww[1], ww[2], ww[3],
-          (xy + 0)[0], (xy + 0)[1],
-          (xy + 2)[0], (xy + 2)[1],
-          (xy + 4)[0], (xy + 4)[1],
-          (xy + 6)[0], (xy + 6)[1], vv[0], vv[1]);
-  */
   return;
 }
 
@@ -604,11 +616,6 @@
     double tmpf = AIR_AFFINE(0, ii, pNum - 1, 0, sNum);
     double tt = tmpf - segi;
     limnCbfSegEval(xy + 2 * ii, seg, tt); /* DIM=2 */
-    /*
-    fprintf(stderr, "!%s: %u -> %u (%g) %g -> (%g,%g)\n",
-            "limnCbfPathSample", ii, segi, tmpf, tt,
-            (xy + 2*ii)[0], (xy + 2*ii)[1]);
-    */
   }
   return;
 }
@@ -737,9 +744,9 @@
 
 Given that this is the inner-loop of other things, it would make sense to have a
 non-public version without all the error checking, but given the prolonged birthing pain
-of the code in this file, the error-checking is a useful and welcome safety-net (and
+of the code in this file, the error-checking is a useful and welcome safety net (and
 being a public function permits easier testing), and that is all ok until profiling shows
-that we are a bottleneck.
+that this function is a bottleneck.
 
 NOTE: this assumes that limnCbfCtxPrep(fctx, lpnt) was called without error!
 That (via ctxBuffersSet) allocates things that we depend on here.
@@ -886,6 +893,33 @@
   return 0;
 }
 
+/* weird helper function:
+if given NULL alpha:
+  returns the threshold on a computed alpha below which we should punt
+otherwise:
+  sets given alpha[0,1] to accomplish punting (with some redundant computation, alas)
+*/
+static double /* Biff: nope */
+puntAlpha(double alpha[2], limnCbfCtx *fctx, const double vv0[2], const double tt1[2],
+          const double tt2[2], const double vv3[2]) {
+  double ret, F2L[2], lenF2L;
+
+  /* note that we're using vv0 and vv3, rather that input data at loi and hii, so that
+  we can reflect the consequence of whatever scale > 0 we might be using */
+  ELL_2V_SUB(F2L, vv3, vv0); /* DIM=2 throughout this */
+  lenF2L = ELL_2V_LEN(F2L);
+  if (!alpha) {
+    ret = lenF2L * (fctx->alphaMin);
+  } else {
+    /* set both alphas to 1/3 of distance from first to last point,
+    but also handle non-unit-length tt1 and tt2 */
+    alpha[0] = lenF2L / (3 * ELL_2V_LEN(tt1));
+    alpha[1] = lenF2L / (3 * ELL_2V_LEN(tt2));
+    ret = 0;
+  }
+  return ret;
+}
+
 /*
 (from paper page 620, after "we need only solve"): solves for the alpha[0,1] that
 minimize squared error between xy[i] and Q(uu[i]) where Q(t) is cubic Bezier spline
@@ -914,11 +948,9 @@
   static const char me[] = "findAlpha";
   int ret;
   uint ii, spanlen = spanLength(lpnt, loi, hii);
-  double det, F2L[2], lenF2L;
-  const double *xy = PP(lpnt), *uu = fctx->uu;
+  double det, alphaMin = puntAlpha(NULL, fctx, vv0, tt1, tt2, vv3);
+  const double *uu = fctx->uu;
 
-  ELL_2V_SUB(F2L, xy + 2 * hii, xy + 2 * loi); /* DIM=2 throughout this */
-  lenF2L = ELL_2V_LEN(F2L);
   if (spanlen > 2) {
     /* GLK using "m" and "M" instead of author's "C". Note that Ai1 and Ai2 are scalings
        of (nominally) unit-length tt1 and tt2 by evaluations of the spline basis
@@ -950,24 +982,21 @@
     alpha[0] = alpha[1] = 0;    /* trigger punting */
   }
   /* test if we should return punted arc */
-  if (!(AIR_EXISTS(det) && AIR_ABS(det) > fctx->detMin
-        && alpha[0] > lenF2L * (fctx->alphaMin)
-        && alpha[1] > lenF2L * (fctx->alphaMin))) {
+  if (!(AIR_EXISTS(det) && AIR_ABS(det) > fctx->detMin /* */
+        && alpha[0] > alphaMin                         /* */
+        && alpha[1] > alphaMin)) {
     if (fctx->verbose) {
       if (spanlen > 2) {
-        printf("%s(i%d): bad |det| %g (vs %g) or alpha %g,%g (vs %g*%g) "
+        printf("%s(i%d): bad |det| %g (vs %g) or alpha %g,%g (vs %g) "
                "--> punted arc\n",
-               me, nrpi, AIR_ABS(det), fctx->detMin, alpha[0], alpha[1], lenF2L,
-               fctx->alphaMin);
+               me, nrpi, AIR_ABS(det), fctx->detMin, alpha[0], alpha[1], alphaMin);
       } else {
         printf("%s(i%d): [%u,%u] spanlen %u tiny --> punting\n", me, nrpi, loi, hii,
                spanlen);
       }
     }
-    /* generate punted arc: set both alphas to 1/3 of distance from first to last point,
-    but also handle non-unit-length tt1 and tt2 */
-    alpha[0] = lenF2L / (3 * ELL_2V_LEN(tt1));
-    alpha[1] = lenF2L / (3 * ELL_2V_LEN(tt2));
+    /* set alpha[0,1] for punted arc */
+    puntAlpha(alpha, fctx, vv0, tt1, tt2, vv3);
     ret = 1;
   } else {
     if (fctx->verbose > 1) {
@@ -995,7 +1024,7 @@
 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, sqd0, absdelt, chop = 0.25;
-  uint try;
+  uint try, tryPairsMax = 4; /* HEY hard-coded limit */
 
   /* SQD(T) = squared distance between Q(T) and P, with side-effects! */
 #define SQD(TT)                                                                         \
@@ -1017,9 +1046,9 @@
     /* 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 */
+  try = 1;                            /* start with 1st try (1-based numbering) */
+  while (SQD(t0 + delt) > sqd0) {     /* while this is an unhelpful step */
+    if (try == 1 + 2 * tryPairsMax) { /* tryPairsMax chop,flip pairs didn't help: bail */
       delt = 0;
       break;
     }
@@ -1078,6 +1107,11 @@
   const double *uu = fctx->uu;
 
   spanlen = spanLength(lpnt, loi, hii);
+  if (0 == loi && 0 == hii && lpnt->isLoop) {
+    /* in the very specific (and rare) case that we've fit a *single* spline
+    to a point loop, starting and ending at 0, then we have to lift spanlen */
+    spanlen += lpnt->num;
+  }
   if (!(spanlen >= 3)) {
     biffAddf(LIMN, "%s: [loi,hii] [%u,%u] -> spanlen %u too small", me, loi, hii,
              spanlen);
@@ -1119,6 +1153,58 @@
   return 0;
 }
 
+/* deathToWacky: if the span length is 3, then with the V0, T1, T2, V3 contraints,
+there is exactly one solution for alpha[0,1] that makes the spline exactly hit the
+middle point at index loi+1 == hii-1. But that doesn't mean it's a good solution -
+a wacky spline curve could be going through the middle point in a direction completely
+unrelated to the two-sided tangent through the point that we'd be enforcing at that point
+if we were splitting a multi-fit there. This function (probably) rejects the wacky curve
+and simulates having a low-accuracy fit, to (probably) trigger another split; function
+name is a reference to a good TAL episode. */
+static int /* Biff: 1 */
+deathToWacky(double alpha[2], const double vv0[2], const double tt1[2],
+             const double tt2[2], const double vv3[2], limnCbfCtx *fctx,
+             const limnCbfPoints *lpnt, uint loi, uint hii) {
+  static const char me[] = "deathToWacky";
+  double vv1[2], vv2[2],
+    splTan[2], /* spline's normalized tangent through middle point */
+    midTan[2], /* two-sided tangent estimated at middle point from points */
+    ww[4], um, tmp, angle;
+  uint midi = (loi + 1) % (lpnt->num);
+  int dopunt;
+
+  /* MANY assumptions about the specific context in which we were called */
+  ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1); /* assuming computed alpha[0,1] */
+  ELL_2V_SCALE_ADD2(vv2, 1, vv3, alpha[1], tt2);
+  um = fctx->uu[1]; /* assuming uu[0,1,2] was just used for nrp */
+  CBD1(splTan, vv0, vv1, vv2, vv3, um, ww);
+  ELL_2V_NORM(splTan, splTan, tmp);
+  /* we emulate how tangent estimation would happen if had decided to split at midi */
+  if (limnCbfTVT(NULL, NULL, midTan, fctx, lpnt, loi, hii, midi,
+                 AIR_FALSE /* NOT oneSided */)) {
+    biffAddf(LIMN, "%s: trying to get middle tangent at %u in [%u,%u]", me, midi, loi,
+             hii);
+  }
+  angle = 180 * ell_2v_angle_d(splTan, midTan) / AIR_PI;
+  dopunt = angle > fctx->wackyAngle;
+  if (fctx->verbose > 1) {
+    printf(
+      "%s: alpha[0,1]=%g,%g -> angle %g btw splTan (%g,%g) and midTan (%g,%g) -> %s\n",
+      me, alpha[0], alpha[1], angle, splTan[0], splTan[1], midTan[0], midTan[1],
+      dopunt ? "YES punt" : "no punt, keep result");
+  }
+  if (dopunt) {
+    /* the fitted spline is wacky, so instead we punt;
+       this probably triggers a later split */
+    puntAlpha(alpha, fctx, vv0, tt1, tt2, vv3);
+    if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
+      biffAddf(LIMN, "%s: trouble with punting", me);
+      return 1;
+    }
+  }
+  return 0;
+}
+
 /*
 fitSingle: fits a single cubic Bezier spline, with minimal error checking (limnCbfSingle
 is the error-checking wrapper around this). The given points coordinates are in
@@ -1148,12 +1234,13 @@
           const double vv3[2], limnCbfCtx *fctx, const limnCbfPoints *lpnt, uint loi,
           uint hii) {
   static const char me[] = "fitSingle";
-  uint iter, spanlen = spanLength(lpnt, loi, hii);
+  uint iter, spanlen;
 
   if (!(alpha && vv0 && tt1 && tt2 && vv3 && fctx && lpnt)) {
     biffAddf(LIMN, "%s: got NULL pointer", me);
     return 1;
   }
+  spanlen = spanLength(lpnt, loi, hii);
   if (!(fctx->uu)) {
     biffAddf(LIMN, "%s: fcgtx->uu NULL; was limnCbfCtxPrep called?", me);
     return 1;
@@ -1165,7 +1252,7 @@
            me, loi, hii, vv0[0], vv0[1], tt1[0], tt1[1], tt2[0], tt2[1], vv3[0], vv3[1]);
   }
   if (2 == spanlen) {
-    /* relying on code in findAlpha() that handles slen==2; return should be 1 */
+    /* relying on code in findAlpha() that handles spanlen==2; return should be 1 */
     if (1 != findAlpha(alpha, -2, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
       biffAddf(LIMN, "%s: what? findAlpha should have returned 1 with spanlen=2", me);
       return 1;
@@ -1177,12 +1264,11 @@
     fctx->distMax = fctx->nrpDeltaDone = 0;
     fctx->distMaxIdx = 0;
     fctx->distBig = 0;
-  } else {             /* slen >= 3 */
+  } else {             /* spanlen >= 3 */
     double delta;      /* avg parameterization change of interior points */
     int lastPunt;      /* last return from fitSingle, ==1 if it punted */
     uint puntFlop = 0; /* # times that return from fitSingle changed */
-    {
-      /* initialize uu parameterization to chord length */
+    {                  /* initialize uu parameterization to chord length */
       unsigned int ii;
       double len;
       const double *xyP, *xyM;
@@ -1215,10 +1301,11 @@
       if (fctx->verbose) {
         printf("%s[%d,%d]: initial (chord length) delta = %g\n", me, loi, hii, delta);
       }
-    }
+    } /* (end uu chord-length initialization ) */
     lastPunt = findAlpha(alpha, -1, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
     if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
-      biffAddf(LIMN, "%s: trouble", me);
+      biffAddf(LIMN, "%s: trouble with initial (chord-length-parm) fit (%g,%g)", me,
+               alpha[0], alpha[1]);
       return 1;
     }
     if (fctx->verbose) {
@@ -1251,7 +1338,7 @@
         */
         lastPunt = punt;
         if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
-          biffAddf(LIMN, "%s: trouble", me);
+          biffAddf(LIMN, "%s: trouble after nrp iter %u fit", me, iter);
           return 1;
         }
         if (fctx->verbose > 1) {
@@ -1266,11 +1353,11 @@
           converged = AIR_TRUE;
           break;
         }
-        /* early code contemplated a check here on whether the reparameterization is
-        diverging instead of converging, but now that delta_t makes sure the Newton
-        steps can't make things worse, divergence is impossible. In any case an
-        unimprovably bad fit is not an error, just a reason to split and try fitting
-        sub-segments. */
+        /* early versions of the code contemplated here a check on whether the
+        reparameterization is diverging instead of converging, but now that delta_t makes
+        sure the Newton steps can't make things worse, divergence is impossible. In any
+        case an unimprovably bad fit is not an error, just a reason to split and try
+        fitting sub-paths */
       }
       if (fctx->verbose) {
         printf("%s[%d,%d]: nrp done after %u iters: ", me, loi, hii, iter);
@@ -1286,6 +1373,10 @@
         }
       }
       fctx->nrpIterDone = iter;
+      if (3 == spanlen) {
+        /* avoid generating totally wacky yet interpolating curve */
+        deathToWacky(alpha, vv0, tt1, tt2, vv3, fctx, lpnt, loi, hii);
+      }
     } else {
       /* else fctx->distBig == 3: dist so big that we don't even try nrp */
       fctx->nrpIterDone = 0;
@@ -1326,16 +1417,35 @@
                (const void *)vv3);
       return 1;
     }
-    /* do not have geometry info; must find it all */
-    if (limnCbfTVT(NULL, v0c, t1c,            /* */
-                   fctx, lpnt, loi, hii, loi, /* */
-                   AIR_TRUE /* oneSided */)
-        || limnCbfTVT(t2c, v3c, NULL,            /* */
-                      fctx, lpnt, loi, hii, hii, /* */
-                      AIR_TRUE /* oneSided */)) {
-      biffAddf(LIMN, "%s: trouble finding geometry info", me);
+    if (loi == hii && hii != 0) {
+      biffAddf(LIMN, "%s: can only handle loi=hii=%u if both 0", me, hii);
       return 1;
     }
+    if (loi == hii) { /* and hence both 0, due to above */
+      /* e.g., we're doing the first fit on of a circular (corner-less) path */
+      if (!lpnt->isLoop) {
+        biffAddf(LIMN, "%s: can handle loi=hii=0 only with point loop", me);
+        return 1;
+      }
+      if (limnCbfTVT(t2c, v0c, t1c,             /* */
+                     fctx, lpnt, loi, hii, loi, /* */
+                     AIR_FALSE /* NOT oneSided */)) {
+        biffAddf(LIMN, "%s: trouble finding two-sided geometry info", me);
+        return 1;
+      }
+      ELL_2V_COPY(v3c, v0c);
+    } else {
+      /* do not have geometry info; must find it all */
+      if (limnCbfTVT(NULL, v0c, t1c,            /* */
+                     fctx, lpnt, loi, hii, loi, /* */
+                     AIR_TRUE /* yes oneSided */)
+          || limnCbfTVT(t2c, v3c, NULL,            /* */
+                        fctx, lpnt, loi, hii, hii, /* */
+                        AIR_TRUE /* yes oneSided */)) {
+        biffAddf(LIMN, "%s: trouble finding one-sided geometry info", me);
+        return 1;
+      }
+    }
     if (fctx->verbose) {
       printf("%s[%u,%u]: found geometry (%g,%g) --> (%g,%g) -- (%g,%g) <-- (%g,%g)\n",
              me, loi, hii, v0c[0], v0c[1], t1c[0], t1c[1], t2c[0], t2c[1], v3c[0],
@@ -1431,8 +1541,7 @@
   double *angle, /* angle[i] is angle at vertex i */
     *vtvt;       /* 6-by-pnum array of tangent,vertex,tangent
                     for ALL vertices */
-  int *corny,    /* corny[i] means vertex i seems like a corner */
-    oneSided = AIR_TRUE;
+  int *corny;    /* corny[i] means vertex i seems like a corner */
   uint pnum = lpnt->num, hii, cnum, vi;
 
   if (!(fctx && lpnt)) {
@@ -1460,9 +1569,9 @@
       }
       hii = pnum - 1;
       if (limnCbfTVT(fctx->ctvt + 0, fctx->ctvt + 2, fctx->ctvt + 4, /* */
-                     fctx, lpnt, 0, hii, 0, oneSided)
+                     fctx, lpnt, 0, hii, 0, AIR_TRUE /* yes oneSided */)
           || limnCbfTVT(fctx->ctvt + 6, fctx->ctvt + 8, fctx->ctvt + 10, /* */
-                        fctx, lpnt, 0, hii, hii, oneSided)) {
+                        fctx, lpnt, 0, hii, hii, AIR_TRUE /* yes oneSided */)) {
         biffAddf(LIMN, "%s: trouble with tangents or vertices for endpoints", me);
         return 1;
       }
@@ -1503,7 +1612,8 @@
     revisit the input data (used previously for tangent estimation) to figure out the
     corner vertices. In a non-loop, we know first and last points will be corners, but we
     still need to find the vertex pos and (one-sided) tangent. */
-    if (limnCbfTVT(LT, VV, RT, fctx, lpnt, 0 /* loi */, 0 /* hii */, vi, oneSided)) {
+    if (limnCbfTVT(LT, VV, RT, fctx, lpnt, 0 /* loi */, 0 /* hii */, vi,
+                   AIR_TRUE /* yes oneSided */)) {
       biffAddf(LIMN, "%s: trouble with tangents or vertices for point %u/%u", me, vi,
                pnum);
       airMopError(mop);
@@ -1582,12 +1692,14 @@
         od = fctx->ctvt + 6 * ci;
         ELL_6V_COPY(od, id);
         if (fctx->verbose) {
-          printf("%s: corner %u is vertex %u\n  T,V,T = (%g,%g)  (%g,%g)  (%g,%g)\n", me,
+          printf("    corner %u is vertex %u\n"
+                 "        T,V,T = (%g,%g)  (%g,%g)  (%g,%g)\n",
                  ci, vi, od[0], od[1], od[2], od[3], od[4], od[5]);
         }
         ci++;
       }
     }
+    assert(ci == cnum);
   }
   fctx->cnum = cnum;
 
@@ -1620,10 +1732,17 @@
     biffAddf(LIMN, "%s: problem getting vertex or tangent info", me);
     return 1;
   }
+  if (!(ELL_2V_EXISTS(V0) && ELL_2V_EXISTS(T1) && /* */
+        ELL_2V_EXISTS(T2) && ELL_2V_EXISTS(V3))) {
+    biffAddf(LIMN,
+             "%s: not all %s vecs exist: v0=(%g,%g) t1=(%g,%g) t2=(%g,%g) v3=(%g,%g)",
+             me, geomGiven ? "given" : "computed", V0[0], V0[1], T1[0], T1[1], T2[0],
+             T2[1], V3[0], V3[1]);
+    return 1;
+  }
   if (fctx->verbose) {
-    printf("%s[%u,%u]_%u: hello; %s v0=(%g,%g), t1=(%g,%g), t2=(%g,%g), "
-           "v3=(%g,%g)\n",
-           me, loi, hii, recDepth, geomGiven ? "given" : "computed", V0[0], V0[1], T1[0],
+    printf("%s[%u,%u]_%u: hello; %s v0=(%g,%g) t1=(%g,%g) t2=(%g,%g) v3=(%g,%g)\n", me,
+           loi, hii, recDepth, geomGiven ? "given" : "computed", V0[0], V0[1], T1[0],
            T1[1], T2[0], T2[1], V3[0], V3[1]);
   }
 
@@ -1648,6 +1767,7 @@
     ELL_2V_SCALE_ADD2(path->seg[0].xy + 2, 1, V0, alpha[0], T1);
     ELL_2V_SCALE_ADD2(path->seg[0].xy + 4, 1, V3, alpha[1], T2);
     ELL_2V_COPY(path->seg[0].xy + 6, V3);
+    path->seg[0].pointNum = spanLength(lpnt, loi, hii);
   } else {
     /* need to subdivide at fctx->distMaxIdx and recurse.  But this is NOT an occasion to
     create a new *corner* there (all corners were likely found as pre-process by
@@ -1669,7 +1789,7 @@
     }
     if (limnCbfTVT(TL, VM, TR,                 /* */
                    fctx, lpnt, loi, hii, midi, /* */
-                   AIR_FALSE /* oneSided */)
+                   AIR_FALSE /* NOT oneSided */)
         || limnCbfMulti(path, V0, T1, TL, VM, recDepth + 1, &fctxL, lpnt, loi, midi)
         || limnCbfMulti(pRth, VM, TR, T2, V3, recDepth + 1, &fctxR, lpnt, midi, hii)) {
       biffAddf(LIMN, "%s[%u,%u]_%u: trouble on recursive fit (midvert %u)", me, loi, hii,
@@ -1721,6 +1841,10 @@
     biffAddf(LIMN, "%s: got NULL pointer", me);
     return 1;
   }
+  if (fctx->verbose) {
+    printf("%s: hello, working on all %u points in %sloop\n", me, lpnt->num,
+           lpnt->isLoop ? "" : "NON-");
+  }
   if (limnCbfCtxPrep(fctx, lpnt)) {
     biffAddf(LIMN, "%s: trouble preparing", me);
     return 1;
@@ -1735,7 +1859,7 @@
     assert(lpnt->isLoop);
     /* in a loop and no corners known */
     if (fctx->verbose) {
-      printf("%s: no corners: finding path to fit point loop", me);
+      printf("%s: no corners: finding path to fit point loop\n", me);
     }
     if (limnCbfMulti(path, NULL, NULL, NULL, NULL, 0, fctx, lpnt, 0 /* loi */,
                      0 /* hii */)) {
@@ -1743,19 +1867,22 @@
       return 1;
     }
   } else {
+    assert(fctx->cnum >= 2);
     /* we do have corners: find segments between corners. The corner vertex part of two
-    segments: the last point in segment I and first point in segment I+1 */
-    uint cii;
-    for (cii = 0; cii < fctx->cnum; cii++) {
+    segments: the last point in segment I and first point in segment I+1. How many
+    segments we analyze depends on whether they're in a loop: if they're in a loop
+    then we do one more */
+    uint cii, ciLast = fctx->cnum - 2 + !!lpnt->isLoop;
+    for (cii = 0; cii <= ciLast; cii++) {
       uint cjj = (cii + 1) % fctx->cnum;
       limnCbfPath *subpath = limnCbfPathNew(0);
       /* 0: left tangent   2: vertex   4: right tangent */
-      const double *V0 = fctx->ctvt + 2 * cii, *T1 = fctx->ctvt + 4 * cii,
-                   *T2 = fctx->ctvt + 0 * cjj, *V3 = fctx->ctvt + 2 * cjj;
+      const double *V0 = fctx->ctvt + 2 + 6 * cii, *T1 = fctx->ctvt + 4 + 6 * cii,
+                   *T2 = fctx->ctvt + 0 + 6 * cjj, *V3 = fctx->ctvt + 2 + 6 * cjj;
       uint loi = fctx->cidx[cii];
       uint hii = fctx->cidx[cjj];
       if (fctx->verbose) {
-        printf("%s: finding subpath from between corners [%u,%u] (points [%u,%u])", me,
+        printf("%s: finding subpath from between corners [%u,%u] (points [%u,%u])\n", me,
                cii, cjj, loi, hii);
       }
       if (limnCbfMulti(subpath, V0, T1, T2, V3, 0, fctx, lpnt, loi, hii)) {

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

[Teem-commit] SF.net SVN: teem:[7201] teem/trunk/src/limn/lpu_cbf it.c

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

Revision: 7201
          http://sourceforge.net/p/teem/code/7201
Author:   kindlmann
Date:     2024-07-11 21:48:45 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
-roll options for debugging loops, and now path printing is uniform

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

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-11 21:47:30 UTC (rev 7200)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-11 21:48:45 UTC (rev 7201)
@@ -60,6 +60,21 @@
   return;
 }
 
+static void
+pathPrint(const char *me, const limnCbfPath *path) {
+  unsigned int si;
+  printf("%s: path has %u segments in %sloop:\n", me, path->segNum,
+         path->isLoop ? "" : "NON-");
+  for (si = 0; si < path->segNum; si++) {
+    limnCbfSeg *seg = path->seg + si;
+    printf("seg[%u]      %g %g     %g %g     %g %g     %g %g     %c (%u) %c\n", si,
+           seg->xy[0], seg->xy[1], seg->xy[2], seg->xy[3], seg->xy[4], seg->xy[5],
+           seg->xy[6], seg->xy[7], seg->corner[0] ? 'C' : '-', seg->pointNum,
+           seg->corner[1] ? 'C' : '-');
+  }
+  return;
+}
+
 static int
 limnPu_cbfitMain(int argc, const char **argv, const char *me, hestParm *hparm) {
   hestOpt *hopt = NULL;
@@ -71,7 +86,7 @@
   double *xy, deltaThresh, psi, cangle, epsilon, nrpIota, time0, dtime, scale, nrpCap,
     synthPow, fitTT[4];
   unsigned int size0, size1, ii, synthNum, pNum, nrpIterMax;
-  int loop, petc, verbose, tvt[4], fitSingleLoHi[2], fitMultiLoHi[2], corner2[2];
+  int loop, roll, petc, verbose, tvt[4], fitSingleLoHi[2], fitMultiLoHi[2], corner2[2];
   char *synthOut, buff[AIR_STRLEN_SMALL + 1];
   limnCbfCtx *fctx;
   limnCbfPath *path;
@@ -86,6 +101,10 @@
   hestOptAdd_Flag(&hopt, "loop", &loop,
                   "-i input xy points are actually a loop: the first point logically "
                   "follows the last point");
+  hestOptAdd_1_Int(
+    &hopt, "roll", "n", &roll, "0",
+    "if points are in a loop, then it shouldn't really matter which point has index 0. "
+    "For debugging, roll the input data by this amount prior to doing any work.");
   hestOptAdd_1_Int(&hopt, "v", "verbose", &verbose, "1", "verbosity level");
   hestOptAdd_1_UInt(&hopt, "synthn", "num", &synthNum, "51",
                     "if saving spline sampling to -so, how many sample.");
@@ -187,6 +206,46 @@
     airMopError(mop);
     return 1;
   }
+  if (roll) {
+    Nrrd *ntmp;
+    double *xy, *tmp;
+    int pnum, pi;
+    if (airStrlen(synthOut)) {
+      fprintf(stderr, "%s: can only roll (%d) input XY points, not splines\n", me, roll);
+      airMopError(mop);
+      return 1;
+    }
+    if (!loop) {
+      fprintf(stderr, "%s: can only roll (%d) a point loop (no -loop)\n", me, roll);
+      airMopError(mop);
+      return 1;
+    }
+    if (!(nrrdTypeDouble == nin->type)) {
+      fprintf(stderr, "%s: need type %s (not %s) point data\n", me,
+              airEnumStr(nrrdType, nrrdTypeDouble), airEnumStr(nrrdType, nin->type));
+      airMopError(mop);
+      return 1;
+    }
+    ntmp = nrrdNew();
+    airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
+    if (nrrdCopy(ntmp, nin)) {
+      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
+      fprintf(stderr, "%s: trouble:\n%s", me, err);
+      airMopError(mop);
+      return 1;
+    }
+    pnum = AIR_INT(nin->axis[1].size);
+    xy = nin->data;
+    tmp = ntmp->data;
+    for (pi = 0; pi < pnum; pi++) {
+      int pt = AIR_MOD(pi - roll, pnum);
+      ELL_2V_COPY(xy + 2 * pi, tmp + 2 * pt);
+      if (!pi) {
+        printf("%s: with roll=%d; xy[0] is now original xy[%d]: %g %g\n", me, roll, pt,
+               xy[0], xy[1]);
+      }
+    }
+  }
 
   if (airStrlen(synthOut)) {
     /* synthesize data from control points */
@@ -350,7 +409,7 @@
   }
 
   if (fitMultiLoHi[0] >= 0) { /* here to call limnCbfMulti once */
-    unsigned int loi, hii, segi;
+    unsigned int loi, hii;
     const double *VTTV[4];
     getLoHi(&loi, &hii, lpnt, fitMultiLoHi[0], fitMultiLoHi[1]);
     getVTTV(VTTV, lpnt, fitTT, loi, hii);
@@ -362,19 +421,13 @@
       airMopError(mop);
       return 1;
     }
-    printf("%s: limnCbfMulti results: %u segments in %s\n", me, path->segNum,
-           path->isLoop ? "loop" : "NOT-loop");
-    for (segi = 0; segi < path->segNum; segi++) {
-      const limnCbfSeg *seg = path->seg + segi;
-      const double *xy = seg->xy;
-      printf("        %g %g    %g %g    %g %g    %g %g    %c %c\n", xy[0], xy[1], xy[2],
-             xy[3], xy[4], xy[5], xy[6], xy[7], seg->corner[0] ? 'C' : '-',
-             seg->corner[1] ? 'C' : '-');
-    }
+    printf("%s: limnCbfMulti results:\n", me);
+    pathPrint(me, path);
     airMopOkay(mop);
     return 0;
   }
 
+  /* whoa - we're actually here to call limnCbfGo! */
   time0 = airTime();
   if (petc) {
     fprintf(stderr, "%s: Press Enter to Continue ... ", me);
@@ -390,29 +443,10 @@
   }
 
   dtime = (airTime() - time0) * 1000;
-  printf("%s: time= %g ms;iterDone= %u ;deltaDone=%g, distMax=%g @ %u\n", me, dtime,
+  printf("%s: time=%g ms   iterDone=%u   deltaDone=%g   distMax=%g @ %u\n", me, dtime,
          fctx->nrpIterDone, fctx->nrpDeltaDone, fctx->distMax, fctx->distMaxIdx);
-  {
-    unsigned int si;
-    printf("%s: path has %u segments:\n", me, path->segNum);
-    for (si = 0; si < path->segNum; si++) {
-      limnCbfSeg *seg = path->seg + si;
-      printf("seg %u: (%g,%g) -- (%g,%g) -- (%g,%g) -- (%g,%g)\n", si, seg->xy[0],
-             seg->xy[1], seg->xy[2], seg->xy[3], seg->xy[4], seg->xy[5], seg->xy[6],
-             seg->xy[7]);
-    }
-  }
+  pathPrint(me, path);
 
-  if (1) {
-    unsigned int oNum = pNum * 100;
-    double *pp = AIR_MALLOC(oNum * 2, double);
-    airMopAdd(mop, pp, airFree, airMopAlways);
-    limnCbfPathSample(pp, oNum, path);
-    for (ii = 0; ii < oNum; ii++) {
-      printf("done %u %g %g\n", ii, (pp + 2 * ii)[0], (pp + 2 * ii)[1]);
-    }
-  }
-
   airMopOkay(mop);
   return 0;
 }

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

[Teem-commit] SF.net SVN: teem:[7200] teem/trunk/src/limn/test

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

Revision: 7200
          http://sourceforge.net/p/teem/code/7200
Author:   kindlmann
Date:     2024-07-11 21:47:30 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
debugging scripts, maybe done?

Modified Paths:
--------------
    teem/trunk/src/limn/test/00-data.sh
    teem/trunk/src/limn/test/03-single.sh
    teem/trunk/src/limn/test/04-multi.sh

Added Paths:
-----------
    teem/trunk/src/limn/test/05-go.sh

Modified: teem/trunk/src/limn/test/00-data.sh
===================================================================
--- teem/trunk/src/limn/test/00-data.sh	2024-07-11 18:38:10 UTC (rev 7199)
+++ teem/trunk/src/limn/test/00-data.sh	2024-07-11 21:47:30 UTC (rev 7200)
@@ -45,4 +45,4 @@
 unu join -i - y -a 0 -incr |
     unu crop -min 0 0 -max M M-1 -o circ.txt
 
-unu 2op spow circ.txt 3 -o pointy.txt
+unu 2op spow circ.txt 4 -o pointy.txt

Modified: teem/trunk/src/limn/test/03-single.sh
===================================================================
--- teem/trunk/src/limn/test/03-single.sh	2024-07-11 18:38:10 UTC (rev 7199)
+++ teem/trunk/src/limn/test/03-single.sh	2024-07-11 21:47:30 UTC (rev 7200)
@@ -37,23 +37,25 @@
 #    because the nrp parms that make sense for a small number of points don't work great
 #    for a huge number of points
 
-# # Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
+# # Good debugging test case: N=18 is a bad fit, N=19 is a perfect fit
+# BUT THEN the improved delta_t fixed this; so both are equally good!
 # # likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
-# N=18
-# # with -ftt 1 0 -0.8944 0.4472
-# echo "-0.5 0.5
-#  2.0  0.5
-# -0.5  0.0
-#  0.5 -0.5" | unu 2op x - 1 | ...
+#N=18
+#echo "-0.5 0.5
+# 2.0  0.5
+#-0.5  0.0
+# 0.5 -0.5" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
 
-N=80
-echo "-0.5 -1
- -1 1
-1 1
-0.5 -1" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
+# This is demo of why we might want a step that optimizes tangent directions
+# after the parameterization has been found
+N=42
+echo "-1 -1
+ 1 1.5
+-1 1.5
+ 1 -1" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 3 -syntho xy-inn-$N.txt
 junk xy-inn-$N.txt
 
-CMD="./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 0 -psi 1000000000 -nim 400 -deltathr 0.000000000001"
+CMD="./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 1 -psi 1000000000 -nim 4000 -deltathr 0.000000000001"
 echo "====== $CMD"
 eval $CMD > log.txt
 cat log.txt; junk log.txt

Modified: teem/trunk/src/limn/test/04-multi.sh
===================================================================
--- teem/trunk/src/limn/test/04-multi.sh	2024-07-11 18:38:10 UTC (rev 7199)
+++ teem/trunk/src/limn/test/04-multi.sh	2024-07-11 21:47:30 UTC (rev 7200)
@@ -33,19 +33,6 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-# why changing N can significantly change the fit:
-#    because the nrp parms that make sense for a small number of points don't work great
-#    for a huge number of points
-
-# # Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
-# # likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
-# N=18
-# # with -ftt 1 0 -0.8944 0.4472
-# echo "-0.5 0.5
-#  2.0  0.5
-# -0.5  0.0
-#  0.5 -0.5" | ./lpu cbfit -i - -synth $N ...
-
 N=199
 echo "-1 -1
 -1 1
@@ -62,8 +49,6 @@
 #IN=pointy.txt
 IN=circ.txt
 
-# TODO: figure out why this is generating NON-geometrically continuous segments
-
 CMD="./lpu cbfit -i $IN -scl 0 -fm 0 -1 -v 2 -psi 3 -eps 0.01"
 echo "====== $CMD"
 eval $CMD > log.txt
@@ -71,8 +56,8 @@
 
 OUT=xy-out.txt
 echo "====== RESULTS: --> $OUT"
-grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8
-grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8 |
+grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9
+grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9 |
     ./lpu cbfit -i - -synthn $((5*N)) -sup 1 -syntho $OUT
 junk $OUT
 
@@ -81,6 +66,7 @@
 unu jhisto -i  $IN $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
 unu jhisto -i $OUT $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
 
+rm -f tmp.png
 
 unu join -i xy-{out,inn,out}.png -a 0 -incr |
     unu resample -s = x2 x2 -k box -o tmp.png

Added: teem/trunk/src/limn/test/05-go.sh
===================================================================
--- teem/trunk/src/limn/test/05-go.sh	                        (rev 0)
+++ teem/trunk/src/limn/test/05-go.sh	2024-07-11 21:47:30 UTC (rev 7200)
@@ -0,0 +1,76 @@
+#!/usr/bin/env bash
+#
+# Teem: Tools to process and visualize scientific data and images
+# Copyright (C) 2009--2024  University of Chicago
+# Copyright (C) 2005--2008  Gordon Kindlmann
+# Copyright (C) 1998--2004  University of Utah
+#
+# This library is free software; you can redistribute it and/or modify it under the terms
+# of the GNU Lesser General Public License (LGPL) as published by the Free Software
+# Foundation; either version 2.1 of the License, or (at your option) any later version.
+# The terms of redistributing and/or modifying this software also include exceptions to
+# the LGPL that facilitate static linking.
+#
+# This library is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+# PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License along with
+# this library; if not, write to Free Software Foundation, Inc., 51 Franklin Street,
+# Fifth Floor, Boston, MA 02110-1301 USA
+#
+
+# 03 = FOURTH stage of testing- another test limnCbfSingle, via lpu cbfit -fs,
+# but with completely free-form input
+
+set -o errexit
+set -o nounset
+shopt -s expand_aliases
+JUNK=""
+function junk { JUNK="$JUNK $@"; }
+function cleanup { rm -rf $JUNK; }
+trap cleanup err exit int term
+# https://devmanual.gentoo.org/tools-reference/bash/
+unset UNRRDU_QUIET_QUIT
+
+N=199
+echo "-1 -1
+-1 1
+ 1 1
+ 0 0" | ./lpu cbfit -i - -synthn $N  -syntho 0.txt
+echo "0 0
+-0.5 -0.5
+1 -1
+1 1" | ./lpu cbfit -i - -synthn $N  -syntho 1.txt
+cat 0.txt 1.txt | uniq > xy-inn.txt
+junk {0,1}.txt xy-inn.txt
+
+IN=xy-inn.txt
+#IN=pointy.txt
+#IN=circ.txt
+
+CMD="./lpu cbfit -i $IN -scl 0.5 -v 0 -psi 3 -eps 0.003 -roll 0"
+echo "====== $CMD"
+eval $CMD > log.txt
+cat log.txt # ; junk log.txt
+
+OUT=xy-out.txt
+echo "====== RESULTS: --> $OUT"
+grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9
+grep "^seg" log.txt | xargs -n 12 echo | cut -d' ' -f 2,3,4,5,6,7,8,9 |
+    ./lpu cbfit -i - -synthn $((5*N)) -sup 1 -syntho $OUT
+junk $OUT
+
+BIN=900
+MM="-min -1.1 1.1 -max 1.1 -1.1"
+unu jhisto -i  $IN $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
+unu jhisto -i $OUT $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
+
+rm -f tmp.png
+
+unu join -i xy-{out,inn,out}.png -a 0 -incr |
+    unu resample -s = x2 x2 -k box -o tmp.png
+
+unu cksum tmp.png # with corners, is it stable w.r.t. -roll ?
+
+open tmp.png


Property changes on: teem/trunk/src/limn/test/05-go.sh
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[Teem-commit] SF.net SVN: teem:[7199] teem/trunk/src/limn/test

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

Revision: 7199
          http://sourceforge.net/p/teem/code/7199
Author:   kindlmann
Date:     2024-07-11 18:38:10 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
test script tweaks

Modified Paths:
--------------
    teem/trunk/src/limn/test/00-data.sh
    teem/trunk/src/limn/test/04-multi.sh

Modified: teem/trunk/src/limn/test/00-data.sh
===================================================================
--- teem/trunk/src/limn/test/00-data.sh	2024-07-11 12:46:50 UTC (rev 7198)
+++ teem/trunk/src/limn/test/00-data.sh	2024-07-11 18:38:10 UTC (rev 7199)
@@ -32,7 +32,7 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-N=300
+N=160
 
 echo 0 1 |
     unu reshape -s 2 |

Modified: teem/trunk/src/limn/test/04-multi.sh
===================================================================
--- teem/trunk/src/limn/test/04-multi.sh	2024-07-11 12:46:50 UTC (rev 7198)
+++ teem/trunk/src/limn/test/04-multi.sh	2024-07-11 18:38:10 UTC (rev 7199)
@@ -59,11 +59,12 @@
 junk {0,1}.txt xy-inn.txt
 
 # IN=xy-inn.txt
-IN=pointy.txt
+#IN=pointy.txt
+IN=circ.txt
 
 # TODO: figure out why this is generating NON-geometrically continuous segments
 
-CMD="./lpu cbfit -i $IN -loop -scl 0 -fm 0 -1 -v 1 -scl 0.5 -eps 0.08"
+CMD="./lpu cbfit -i $IN -scl 0 -fm 0 -1 -v 2 -psi 3 -eps 0.01"
 echo "====== $CMD"
 eval $CMD > log.txt
 cat log.txt # ; junk log.txt
@@ -72,12 +73,13 @@
 echo "====== RESULTS: --> $OUT"
 grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8
 grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8 |
-    ./lpu cbfit -i - -synthn $((2*N)) -sup 1 -syntho $OUT
+    ./lpu cbfit -i - -synthn $((5*N)) -sup 1 -syntho $OUT
 junk $OUT
 
-BIN=500
-unu jhisto -i  $IN -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
-unu jhisto -i $OUT -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
+BIN=900
+MM="-min -1.1 1.1 -max 1.1 -1.1"
+unu jhisto -i  $IN $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
+unu jhisto -i $OUT $MM -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
 
 
 unu join -i xy-{out,inn,out}.png -a 0 -incr |

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

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

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

Revision: 7198
          http://sourceforge.net/p/teem/code/7198
Author:   kindlmann
Date:     2024-07-11 12:46:50 +0000 (Thu, 11 Jul 2024)
Log Message:
-----------
clarifyig comments

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

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-11 12:11:37 UTC (rev 7197)
+++ teem/trunk/src/limn/splineFit.c	2024-07-11 12:46:50 UTC (rev 7198)
@@ -980,15 +980,17 @@
   return ret;
 }
 
-/* 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 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 */
+/* delta_t is part of paper pg 621 eqs (7) and (8): we find 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 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 carefulness is NOT described in
+the paper, nor is it 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]) {
@@ -1264,14 +1266,11 @@
           converged = AIR_TRUE;
           break;
         }
-        /* maybe TODO: add logic here to catch if delta is getting bigger and bigger,
-        i.e. the reparm is diverging instead of converging. A younger GLK seemed to think
-        this could happen with the spline making a loop away from a small number of
-        points, e.g.: 4 points on spline defined by vv0 = (1,1), tt1 = (1,2), tt2 =
-        (1,2), vv3 = (0,1). On the other hand, it's not like we have a strategy for
-        doing a different/smarter reparm to handle that, and if it does happen, our
-        failure to fit will likely (in the context of limnCbfMulti) merely trigger
-        subdivision, which isn't terrible */
+        /* early code contemplated a check here on whether the reparameterization is
+        diverging instead of converging, but now that delta_t makes sure the Newton
+        steps can't make things worse, divergence is impossible. In any case an
+        unimprovably bad fit is not an error, just a reason to split and try fitting
+        sub-segments. */
       }
       if (fctx->verbose) {
         printf("%s[%d,%d]: nrp done after %u iters: ", me, loi, hii, iter);
@@ -1354,7 +1353,7 @@
 /*
 limnCbfSingle
 
-Basically a very error-checked version of fitSingle; in the limn API because it's needed
+Basically a very error-checked version of fitSingle; in the limn API because it's useful
 for testing. Unlike fitSingle, the geometry info vv0, tt1, tt2, vv3 can either be punted
 on (by passing NULL for all) or not, by passing specific vectors for all. The results are
 converted into the fields in the given limnCbfSeg *seg.  Despite misgivings, we set

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

[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.

[Teem-commit] SF.net SVN: teem:[7196] teem/trunk/src/limn

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

Revision: 7196
          http://sourceforge.net/p/teem/code/7196
Author:   kindlmann
Date:     2024-07-10 23:58:51 +0000 (Wed, 10 Jul 2024)
Log Message:
-----------
improved the nrp (Newton-based ReParameterization) so that it never takes bad steps; it did that previously

Modified Paths:
--------------
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/00-data.sh

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-10 22:41:23 UTC (rev 7195)
+++ teem/trunk/src/limn/splineFit.c	2024-07-10 23:58:51 UTC (rev 7196)
@@ -980,22 +980,65 @@
   return ret;
 }
 
-#if 0
-/* Paper pg 621 eqs (7) and (8): the change in the spline parameter to improve how Q(u_i)
-approximates vertex i, but naming the paramter "t" since that's what the paper does.
-This is not just the Newton step, but a check that the Newton step is an improvement
-(a check that is not described in the paper, and is not in the author's code, but which
-does matter in cases where the spine is a poor fit) */
+/* 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. */
 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], denom;
+  double Q0[2], Q1[2], Q2[2], QmP[2], ww[4], denom, delt = 0, sqd0;
+
+  /* SQD(T) = squared distance between Q(T) and P */
+#define SQD(TT)                                                                         \
+  (CBD0(Q0, V0, V1, V2, V3, (TT), ww), /* */                                            \
+   ELL_2V_SUB(QmP, Q0, P),             /* */                                            \
+   ELL_2V_DOT(QmP, QmP))
+  sqd0 = SQD(t0);
+  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 (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;
+            }
+          }
+        }
+      }
+    }
+  }
+#undef SQD
+  return delt;
 }
-#endif
 
 /*
 using Newton iterations to try to find a better places at which to evaluate the spline in
-order to match the given points xy
+order to match the given points xy; these locations are updated (in-place) in fctx->uu
 */
 static double
 reparm(const limnCbfCtx *fctx, /* must be non-NULL */
@@ -1002,7 +1045,7 @@
        const double alpha[2], const double vv0[2], const double tt1[2],
        const double tt2[2], const double vv3[2], const limnCbfPoints *lpnt, uint loi,
        uint hii) {
-  static const char me[] = "reparm";
+  /* static const char me[] = "reparm"; */
   uint ii, spanlen;
   double vv1[2], vv2[2], delta, cap;
   double *uu = fctx->uu;
@@ -1017,35 +1060,10 @@
   /* only changing parameterization of interior points,
      not the first (ii=0) or last (ii=pNum-1) */
   for (ii = 1; ii < spanlen - 1; ii++) {
-    double denom, delu = 0, QmP[2], ww[4], tt, Q0[2], Q1[2], Q2[2];
     const double *P = PPlowerI(lpnt, AIR_INT(loi + ii));
-    tt = uu[ii];
-    CBD0(Q0, vv0, vv1, vv2, vv3, tt, ww);
-    CBD1(Q1, vv0, vv1, vv2, vv3, tt, ww);
-    CBD2(Q2, vv0, vv1, vv2, vv3, tt, ww);
-    ELL_2V_SUB(QmP, Q0, P);
-    denom = ELL_2V_DOT(Q1, Q1) + ELL_2V_DOT(QmP, Q2);
-    if (denom) {
-      double absdelu;
-      delu = ELL_2V_DOT(QmP, Q1) / denom;
-      absdelu = fabs(delu);
-      if (absdelu > cap) {
-        /* cap Newton step */
-        delu *= cap / absdelu;
-      }
-      uu[ii] = tt - delu;
-    }
-    /* delu = delta_t(tt, cap, P, vv0, vv1, vv2, vv3); */
+    double delu = delta_t(uu[ii], cap, P, vv0, vv1, vv2, vv3);
+    uu[ii] += delu;
     delta += fabs(delu);
-    if (fctx->verbose > 1) {
-      double R[2], dR[2]; /* R is the new Q */
-      CBD0(R, vv0, vv1, vv2, vv3, uu[ii], ww);
-      ELL_2V_SUB(dR, R, P);
-      printf("%s[%2u]: %g <-- %g - %g\n", me, ii, uu[ii], tt, delu);
-      printf("     %g=|(%g,%g)-(%g,%g)|   <--   %g=|(%g,%g)-(%g,%g)|\n", ELL_2V_LEN(dR),
-             R[0], R[1], P[0], P[1], /* */
-             ELL_2V_LEN(QmP), Q0[0], Q0[1], P[0], P[1]);
-    }
   }
   delta /= spanlen - 2; /* number of interior points */
   return delta;

Modified: teem/trunk/src/limn/test/00-data.sh
===================================================================
--- teem/trunk/src/limn/test/00-data.sh	2024-07-10 22:41:23 UTC (rev 7195)
+++ teem/trunk/src/limn/test/00-data.sh	2024-07-10 23:58:51 UTC (rev 7196)
@@ -32,7 +32,7 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-N=60
+N=300
 
 echo 0 1 |
     unu reshape -s 2 |

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

[Teem-commit] SF.net SVN: teem:[7195] teem/trunk/src/limn

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

Revision: 7195
          http://sourceforge.net/p/teem/code/7195
Author:   kindlmann
Date:     2024-07-10 22:41:23 +0000 (Wed, 10 Jul 2024)
Log Message:
-----------
renaming simple to punt (and removing incorrect statements of straight path), and some more refined debugging statements

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

Modified: teem/trunk/src/limn/limn.h
===================================================================
--- teem/trunk/src/limn/limn.h	2024-07-09 22:40:47 UTC (rev 7194)
+++ teem/trunk/src/limn/limn.h	2024-07-10 22:41:23 UTC (rev 7195)
@@ -599,8 +599,8 @@
   /* ----------- output --------- */
   unsigned int nrpIterDone, /* number of nrp iters taken */
     distMaxIdx,             /* which point had distance distMax */
-    nrpSimpleFlop;          /* # times that single-spline fit flip-flopped between a
-                               curved vs a "simple" straight arc */
+    nrpPuntFlop;            /* # times that single-spline fit flip-flopped between a
+                               well-computed spline vs a punted one */
   double distMax,           /* max distance to given points */
     nrpDeltaDone,           /* latest mean parameterization change by nrp */
     alphaDet;               /* min det of matrix inverted to find alpha */

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 22:40:47 UTC (rev 7194)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-10 22:41:23 UTC (rev 7195)
@@ -313,8 +313,8 @@
       airMopError(mop);
       return 1;
     }
-    printf("%s: nrpIterDone %u    nrpSimpleFlop %u    distMax %g @ %u/%u (big %d)\n", me,
-           fctx->nrpIterDone, fctx->nrpSimpleFlop, fctx->distMax, fctx->distMaxIdx,
+    printf("%s: nrpIterDone %u    nrpPuntFlop %u    distMax %g @ %u/%u (big %d)\n", me,
+           fctx->nrpIterDone, fctx->nrpPuntFlop, fctx->distMax, fctx->distMaxIdx,
            lpnt->num, fctx->distBig);
     printf("%s: limnCbfSingle spline result:\n", me);
     for (ii = 0; ii < 4; ii++) {

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 22:40:47 UTC (rev 7194)
+++ teem/trunk/src/limn/splineFit.c	2024-07-10 22:41:23 UTC (rev 7195)
@@ -297,7 +297,7 @@
   /* initialize outputs to bogus valus */
   fctx->nrpIterDone = UINT_MAX;
   fctx->distMaxIdx = UINT_MAX;
-  fctx->nrpSimpleFlop = UINT_MAX;
+  fctx->nrpPuntFlop = UINT_MAX;
   fctx->distMax = AIR_POS_INF;
   fctx->nrpDeltaDone = AIR_POS_INF;
   fctx->alphaDet = 0;
@@ -891,9 +891,9 @@
 minimize squared error between xy[i] and Q(uu[i]) where Q(t) is cubic Bezier spline
 through vv0, vv0 + alpha[0]*tt1, vv3 + alpha[1]*tt2, and vv3.
 
-There are various conditions where the generated spline ignores the xy array and
-instead is what one could call a "simple arc" (with control points at 1/3 and 2/3 the
-distance between the end points):
+There are various conditions where the generated spline (implied by the alpha[0,1] set
+here) ignores the xy array and is instead what we could call a "punt", with control
+points at 1/3 and 2/3 the distance between the end points:
  - having only two points (xy contains only the end points)
  - the determinant of the 2x2 matrix that is inverted to solve
    for alpha is too close to zero (this test was not part of the
@@ -900,13 +900,15 @@
    author's code)
  - the solved alphas are not convincingly positive
 
-This function is the only place where the straight-line "simple arc" is generated, and
-then we return 1, otherwise 0 (when finding alpha[0,1] for a curved segment). But
-generating the simple arc is not actually an error or problem: if the simple arc is bad
-at fitting the data (as determined by findDist) then it may be subdivided, and that's ok.
+This function is *the* place where the punted spline is generated, in which case we
+return 1, otherwise (when we set alpha[0,1] computed via the 2x2 matrix solve) we return
+0. But generating the punted arc is not actually an error or problem: if the punted arc
+is bad at fitting the data (as determined by findDist) then it may be subdivided, and
+that's ok.
 */
 static int /* Biff: nope */
-findAlpha(double alpha[2], limnCbfCtx *fctx, const double vv0[2], const double tt1[2],
+findAlpha(double alpha[2], int nrpi /* which nrp iter we're on, just for debugging */,
+          limnCbfCtx *fctx, const double vv0[2], const double tt1[2],
           const double tt2[2], const double vv3[2], const limnCbfPoints *lpnt, uint loi,
           uint hii) {
   static const char me[] = "findAlpha";
@@ -945,23 +947,24 @@
     ELL_2MV_MUL(alpha, MI, xx); /* solve for alpha[0,1] */
   } else {                      /* spanlen <= 2 */
     det = 1;                    /* bogus but harmless */
-    alpha[0] = alpha[1] = 0;    /* trigger simple arc code */
+    alpha[0] = alpha[1] = 0;    /* trigger punting */
   }
-  /* test if we should return simple arc */
+  /* test if we should return punted arc */
   if (!(AIR_EXISTS(det) && AIR_ABS(det) > fctx->detMin
         && alpha[0] > lenF2L * (fctx->alphaMin)
         && alpha[1] > lenF2L * (fctx->alphaMin))) {
     if (fctx->verbose) {
       if (spanlen > 2) {
-        printf("%s: bad |det| %g (vs %g) or alpha %g,%g (vs %g*%g) "
-               "--> simple arc\n",
-               me, AIR_ABS(det), fctx->detMin, alpha[0], alpha[1], lenF2L,
+        printf("%s(i%d): bad |det| %g (vs %g) or alpha %g,%g (vs %g*%g) "
+               "--> punted arc\n",
+               me, nrpi, AIR_ABS(det), fctx->detMin, alpha[0], alpha[1], lenF2L,
                fctx->alphaMin);
       } else {
-        printf("%s: [%u,%u] spanlen %u tiny --> simple arc\n", me, loi, hii, spanlen);
+        printf("%s(i%d): [%u,%u] spanlen %u tiny --> punting\n", me, nrpi, loi, hii,
+               spanlen);
       }
     }
-    /* generate simple arc: set both alphas to 1/3 of distance from first to last point,
+    /* generate punted arc: set both alphas to 1/3 of distance from first to last point,
     but also handle non-unit-length tt1 and tt2 */
     alpha[0] = lenF2L / (3 * ELL_2V_LEN(tt1));
     alpha[1] = lenF2L / (3 * ELL_2V_LEN(tt2));
@@ -968,7 +971,8 @@
     ret = 1;
   } else {
     if (fctx->verbose > 1) {
-      printf("%s: all good: det %g, alpha %g,%g\n", me, det, alpha[0], alpha[1]);
+      printf("%s(i%d): all good: det %g, alpha %g,%g\n", me, nrpi, det, alpha[0],
+             alpha[1]);
     }
     ret = 0;
   }
@@ -976,6 +980,19 @@
   return ret;
 }
 
+#if 0
+/* Paper pg 621 eqs (7) and (8): the change in the spline parameter to improve how Q(u_i)
+approximates vertex i, but naming the paramter "t" since that's what the paper does.
+This is not just the Newton step, but a check that the Newton step is an improvement
+(a check that is not described in the paper, and is not in the author's code, but which
+does matter in cases where the spine is a poor fit) */
+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], denom;
+}
+#endif
+
 /*
 using Newton iterations to try to find a better places at which to evaluate the spline in
 order to match the given points xy
@@ -987,13 +1004,13 @@
        uint hii) {
   static const char me[] = "reparm";
   uint ii, spanlen;
-  double vv1[2], vv2[2], delta, maxdelu;
+  double vv1[2], vv2[2], delta, cap;
   double *uu = fctx->uu;
 
   spanlen = spanLength(lpnt, loi, hii);
   assert(spanlen >= 3);
-  /* average interior u[i+1]-u[i] is 1/(spanlen-2) */
-  maxdelu = fctx->nrpCap / (spanlen - 2);
+  /* average u[i+1]-u[i] is 1/(spanlen-1) */
+  cap = fctx->nrpCap / (spanlen - 1);
   ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1);
   ELL_2V_SCALE_ADD2(vv2, 1, vv3, alpha[1], tt2);
   delta = 0;
@@ -1000,26 +1017,34 @@
   /* only changing parameterization of interior points,
      not the first (ii=0) or last (ii=pNum-1) */
   for (ii = 1; ii < spanlen - 1; ii++) {
-    double numer, denom, delu, absdelu, df[2], ww[4], tt, Q[2], QD[2], QDD[2];
-    const double *xy;
+    double denom, delu = 0, QmP[2], ww[4], tt, Q0[2], Q1[2], Q2[2];
+    const double *P = PPlowerI(lpnt, AIR_INT(loi + ii));
     tt = uu[ii];
-    CBD0(Q, vv0, vv1, vv2, vv3, tt, ww);
-    CBD1(QD, vv0, vv1, vv2, vv3, tt, ww);
-    CBD2(QDD, vv0, vv1, vv2, vv3, tt, ww);
-    xy = PPlowerI(lpnt, AIR_INT(loi + ii));
-    ELL_2V_SUB(df, Q, xy);
-    numer = ELL_2V_DOT(df, QD);
-    denom = ELL_2V_DOT(QD, QD) + ELL_2V_DOT(df, QDD);
-    delu = numer / denom;
-    absdelu = fabs(delu);
-    if (absdelu > maxdelu) {
-      /* cap Newton step */
-      delu *= maxdelu / absdelu;
+    CBD0(Q0, vv0, vv1, vv2, vv3, tt, ww);
+    CBD1(Q1, vv0, vv1, vv2, vv3, tt, ww);
+    CBD2(Q2, vv0, vv1, vv2, vv3, tt, ww);
+    ELL_2V_SUB(QmP, Q0, P);
+    denom = ELL_2V_DOT(Q1, Q1) + ELL_2V_DOT(QmP, Q2);
+    if (denom) {
+      double absdelu;
+      delu = ELL_2V_DOT(QmP, Q1) / denom;
+      absdelu = fabs(delu);
+      if (absdelu > cap) {
+        /* cap Newton step */
+        delu *= cap / absdelu;
+      }
+      uu[ii] = tt - delu;
     }
-    uu[ii] = tt - delu;
+    /* delu = delta_t(tt, cap, P, vv0, vv1, vv2, vv3); */
     delta += fabs(delu);
     if (fctx->verbose > 1) {
+      double R[2], dR[2]; /* R is the new Q */
+      CBD0(R, vv0, vv1, vv2, vv3, uu[ii], ww);
+      ELL_2V_SUB(dR, R, P);
       printf("%s[%2u]: %g <-- %g - %g\n", me, ii, uu[ii], tt, delu);
+      printf("     %g=|(%g,%g)-(%g,%g)|   <--   %g=|(%g,%g)-(%g,%g)|\n", ELL_2V_LEN(dR),
+             R[0], R[1], P[0], P[1], /* */
+             ELL_2V_LEN(QmP), Q0[0], Q0[1], P[0], P[1]);
     }
   }
   delta /= spanlen - 2; /* number of interior points */
@@ -1126,21 +1151,21 @@
   }
   if (2 == spanlen) {
     /* relying on code in findAlpha() that handles slen==2; return should be 1 */
-    if (1 != findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
+    if (1 != findAlpha(alpha, -2, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
       biffAddf(LIMN, "%s: what? findAlpha should have returned 1 with spanlen=2", me);
       return 1;
     }
     /* nrp is moot */
     fctx->nrpIterDone = 0;
-    fctx->nrpSimpleFlop = 0;
+    fctx->nrpPuntFlop = 0;
     /* emmulate results of calling findDist() */
     fctx->distMax = fctx->nrpDeltaDone = 0;
     fctx->distMaxIdx = 0;
     fctx->distBig = 0;
-  } else {               /* slen >= 3 */
-    double delta;        /* avg parameterization change of interior points */
-    int lastSimple;      /* last return from fitSingle, ==1 if it found "simple" arc */
-    uint simpleFlop = 0; /* # times that return from fitSingle changed */
+  } else {             /* slen >= 3 */
+    double delta;      /* avg parameterization change of interior points */
+    int lastPunt;      /* last return from fitSingle, ==1 if it punted */
+    uint puntFlop = 0; /* # times that return from fitSingle changed */
     {
       /* initialize uu parameterization to chord length */
       unsigned int ii;
@@ -1176,22 +1201,22 @@
         printf("%s[%d,%d]: initial (chord length) delta = %g\n", me, loi, hii, delta);
       }
     }
-    lastSimple = findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+    lastPunt = findAlpha(alpha, -1, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
     if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
       biffAddf(LIMN, "%s: trouble", me);
       return 1;
     }
     if (fctx->verbose) {
-      printf(
-        "%s[%d,%d]: found alpha %g %g, maxdist %g @ %u (big %d) (%u max nrp iters)\n",
-        me, loi, hii, alpha[0], alpha[1], fctx->distMax, fctx->distMaxIdx, fctx->distBig,
-        fctx->nrpIterMax);
+      printf("%s[%d,%d]: found (%s) alpha %g %g, maxdist %g @ %u (big %d) (%u max nrp "
+             "iters)\n",
+             me, loi, hii, lastPunt ? "punt" : "calc", alpha[0], alpha[1], fctx->distMax,
+             fctx->distMaxIdx, fctx->distBig, fctx->nrpIterMax);
     }
     if (fctx->distBig < 3) {
       /* initial fit isn't awful; try making it better with nrp */
       int converged = AIR_FALSE;
       for (iter = 0; fctx->distBig && iter < fctx->nrpIterMax; iter++) {
-        int simple;
+        int punt;
         if (fctx->verbose > 1) {
           printf("%s[%d,%d]: nrp iter %u starting with alpha %g,%g (det %g) (big %d)\n",
                  me, loi, hii, iter, alpha[0], alpha[1], fctx->alphaDet, fctx->distBig);
@@ -1198,18 +1223,18 @@
         }
         /* *this* is where nrp = Newton-based ReParameterization happens */
         delta = reparm(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii);
-        simple = findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
-        simpleFlop += simple != lastSimple;
-        /* seems like a good idea, but not clearly needed, and has some false positives
-        if (simpleFlop > 3 && simpleFlop - 1 > iter / 2) {
+        punt = findAlpha(alpha, iter, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+        puntFlop += punt != lastPunt;
+        /* seems like a good idea, but not obviously needed, and has some false positives
+        if (puntFlop > 3 && puntFlop - 1 > iter / 2) {
           biffAddf(LIMN,
                    "%s: (nrp iter %u) previous findAlpha return %d but now %d, and has "
-                   "changed %d times => finding simple arc too unstable under nrp",
-                   me, iter, lastSimple, simple, simpleFlop);
+                   "changed %d times => whether to punt too unstable under nrp",
+                   me, iter, lastPunt, punt, puntFlop);
           return 1;
         }
         */
-        lastSimple = simple;
+        lastPunt = punt;
         if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
           biffAddf(LIMN, "%s: trouble", me);
           return 1;
@@ -1258,7 +1283,7 @@
       }
     }
     fctx->nrpDeltaDone = delta;
-    fctx->nrpSimpleFlop = simpleFlop;
+    fctx->nrpPuntFlop = puntFlop;
   }
   if (fctx->verbose) {
     printf("%s[%d,%d]: leaving with alpha %g %g\n", me, loi, hii, alpha[0], alpha[1]);
@@ -1300,8 +1325,9 @@
       return 1;
     }
     if (fctx->verbose) {
-      printf("%s: found geometry (%g,%g) --> (%g,%g) -- (%g,%g) <-- (%g,%g)\n", me,
-             v0c[0], v0c[1], t1c[0], t1c[1], t2c[0], t2c[1], v3c[0], v3c[1]);
+      printf("%s[%u,%u]: found geometry (%g,%g) --> (%g,%g) -- (%g,%g) <-- (%g,%g)\n",
+             me, loi, hii, v0c[0], v0c[1], t1c[0], t1c[1], t2c[0], t2c[1], v3c[0],
+             v3c[1]);
     }
     ELL_2V_COPY(vttv[0], v0c);
     ELL_2V_COPY(vttv[1], t1c);
@@ -1616,6 +1642,7 @@
     limnCbfCorners, called by limnCbfGo). We find a 2-sided tangent and vertex */
     uint midi = fctx->distMaxIdx;
     double TL[2], VM[2], TR[2];
+    const double *pp = PP(lpnt);
     limnCbfCtx fctxL, fctxR;
     limnCbfPath *pRth = limnCbfPathNew(0); /* path on right side of new middle */
     /* holy moly sneakiness: we make two shallow copies of the context, because we
@@ -1624,8 +1651,9 @@
     memcpy(&fctxL, fctx, sizeof(limnCbfCtx));
     memcpy(&fctxR, fctx, sizeof(limnCbfCtx));
     if (fctx->verbose) {
-      printf("%s[%u,%u]_%u: dist %g (big %d) --> split at %u\n", me, loi, hii, recDepth,
-             fctx->distMax, fctx->distBig, midi);
+      printf("%s[%u,%u]_%u: dist %g (big %d) --> split at %u (%g,%g)\n", me, loi, hii,
+             recDepth, fctx->distMax, fctx->distBig, midi, (pp + 2 * midi)[0],
+             (pp + 2 * midi)[1]);
     }
     if (limnCbfTVT(TL, VM, TR,                 /* */
                    fctx, lpnt, loi, hii, midi, /* */

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

[Teem-commit] SF.net SVN: teem:[7194] teem/trunk/src/limn

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

Revision: 7194
          http://sourceforge.net/p/teem/code/7194
Author:   kindlmann
Date:     2024-07-09 22:40:47 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
segfault fixed, debugging continues

Modified Paths:
--------------
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/04-multi.sh

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 22:29:06 UTC (rev 7193)
+++ teem/trunk/src/limn/splineFit.c	2024-07-09 22:40:47 UTC (rev 7194)
@@ -69,10 +69,10 @@
 you can only use actual indices for any memory access. Early iterations of the code used
 lifted indices in lots of places, but GLK got confused, so the following summary of who
 calls whom ("who -- whom") was made to trace where [loi,hii] spans originate, and to
-ensure that all the functions here only take actual indices as parameters, including
-limnCbfTVT (which computes a tangent,vertex,tangent triple at a given point). Nearly all
-of the cleverness with lifted indices happens in limnCbfTVT, and its idxLift helper
-function ensures it gets actual indices.
+ensure that all the functions here only take actual (not-lifted) indices as parameters,
+including limnCbfTVT (which computes a tangent,vertex,tangent triple at a given point).
+Nearly all of the cleverness with lifted indices happens in limnCbfTVT, and its idxLift
+helper function ensures it gets actual indices.
 
 limnCbfTVT -- idxLift to convert given actual indices into lifted,
               and does computations with lifted indices (with signed offsets)
@@ -713,7 +713,7 @@
 that we go from lifted index to actual index */
 static const double *
 PPlowerI(const limnCbfPoints *lpnt, int ssi) {
-  int pnum = lpnt->num;
+  int pnum = AIR_INT(lpnt->num);
   ssi = AIR_MOD(ssi, pnum);
   return PP(lpnt) + 2 * ssi; /* DIM=2 */
 }

Modified: teem/trunk/src/limn/test/04-multi.sh
===================================================================
--- teem/trunk/src/limn/test/04-multi.sh	2024-07-09 22:29:06 UTC (rev 7193)
+++ teem/trunk/src/limn/test/04-multi.sh	2024-07-09 22:40:47 UTC (rev 7194)
@@ -59,9 +59,11 @@
 junk {0,1}.txt xy-inn.txt
 
 # IN=xy-inn.txt
-IN=circ.txt
+IN=pointy.txt
 
-CMD="./lpu cbfit -i $IN -loop -scl 0 -fm 0 -1 -v 1"
+# TODO: figure out why this is generating NON-geometrically continuous segments
+
+CMD="./lpu cbfit -i $IN -loop -scl 0 -fm 0 -1 -v 1 -scl 0.5 -eps 0.08"
 echo "====== $CMD"
 eval $CMD > log.txt
 cat log.txt # ; junk log.txt

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

[Teem-commit] SF.net SVN: teem:[7193] teem/trunk/src/limn

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

Revision: 7193
          http://sourceforge.net/p/teem/code/7193
Author:   kindlmann
Date:     2024-07-09 22:29:06 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
adding memory of recursion depth to limnCbfMulti

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

Modified: teem/trunk/src/limn/limn.h
===================================================================
--- teem/trunk/src/limn/limn.h	2024-07-09 22:28:21 UTC (rev 7192)
+++ teem/trunk/src/limn/limn.h	2024-07-09 22:29:06 UTC (rev 7193)
@@ -932,7 +932,8 @@
                               unsigned int hii);
 LIMN_EXPORT int limnCbfCorners(limnCbfCtx *fctx, const limnCbfPoints *lpnt);
 LIMN_EXPORT int limnCbfMulti(limnCbfPath *path, const double vv0[2], const double tt1[2],
-                             const double tt2[2], const double vv3[2], limnCbfCtx *fctx,
+                             const double tt2[2], const double vv3[2],
+                             unsigned int recurseDepth, limnCbfCtx *fctx,
                              const limnCbfPoints *lpnt, unsigned int loi,
                              unsigned int hii);
 LIMN_EXPORT int limnCbfGo(limnCbfPath *path, limnCbfCtx *fctx,

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 22:28:21 UTC (rev 7192)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 22:29:06 UTC (rev 7193)
@@ -355,7 +355,7 @@
     getLoHi(&loi, &hii, lpnt, fitMultiLoHi[0], fitMultiLoHi[1]);
     getVTTV(VTTV, lpnt, fitTT, loi, hii);
     if (limnCbfCtxPrep(fctx, lpnt)
-        || limnCbfMulti(path, VTTV[0], VTTV[1], VTTV[2], VTTV[3], fctx, lpnt, loi,
+        || limnCbfMulti(path, VTTV[0], VTTV[1], VTTV[2], VTTV[3], 0, fctx, lpnt, loi,
                         hii)) {
       airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
       fprintf(stderr, "%s: trouble doing multi fit:\n%s", me, err);

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 22:28:21 UTC (rev 7192)
+++ teem/trunk/src/limn/splineFit.c	2024-07-09 22:29:06 UTC (rev 7193)
@@ -1568,8 +1568,8 @@
 */
 int /* Biff: 1 */
 limnCbfMulti(limnCbfPath *path, const double vv0[2], const double tt1[2],
-             const double tt2[2], const double vv3[2], limnCbfCtx *fctx,
-             const limnCbfPoints *lpnt, uint loi, uint hii) {
+             const double tt2[2], const double vv3[2], unsigned int recDepth,
+             limnCbfCtx *fctx, const limnCbfPoints *lpnt, uint loi, uint hii) {
   static const char me[] = "limnCbfMulti";
   double alpha[2], V0[2], T1[2], T2[2], V3[2], *vttv[4] = {V0, T1, T2, V3};
   int geomGiven;
@@ -1583,15 +1583,15 @@
     return 1;
   }
   if (fctx->verbose) {
-    printf("%s[%u,%u]: hello; %s v0=(%g,%g), t1=(%g,%g), t2=(%g,%g), "
+    printf("%s[%u,%u]_%u: hello; %s v0=(%g,%g), t1=(%g,%g), t2=(%g,%g), "
            "v3=(%g,%g)\n",
-           me, loi, hii, geomGiven ? "given" : "computed", V0[0], V0[1], T1[0], T1[1],
-           T2[0], T2[1], V3[0], V3[1]);
+           me, loi, hii, recDepth, geomGiven ? "given" : "computed", V0[0], V0[1], T1[0],
+           T1[1], T2[0], T2[1], V3[0], V3[1]);
   }
 
   /* first try fitting a single spline */
   if (fctx->verbose) {
-    printf("%s[%u,%u]: trying single fit on all points\n", me, loi, hii);
+    printf("%s[%u,%u]_%u: trying single fit on all points\n", me, loi, hii, recDepth);
   }
   if (fitSingle(alpha, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
     biffAddf(LIMN, "%s: fitSingle failed", me);
@@ -1600,9 +1600,9 @@
   if (fctx->distBig <= 1) {
     /* max dist was <= fctx->epsilon: single fit was good enough */
     if (fctx->verbose) {
-      printf("%s[%u,%u]: single fit good! nrpi=%u; maxdist=%g @ %u <= %g; "
+      printf("%s[%u,%u]_%u: single fit good! nrpi=%u; maxdist=%g @ %u <= %g; "
              "big=%d det=%g alpha=%g,%g\n",
-             me, loi, hii, fctx->nrpIterDone, fctx->distMax, fctx->distMaxIdx,
+             me, loi, hii, recDepth, fctx->nrpIterDone, fctx->distMax, fctx->distMaxIdx,
              fctx->epsilon, fctx->distBig, fctx->alphaDet, alpha[0], alpha[1]);
     }
     airArrayLenSet(path->segArr, 1);
@@ -1624,16 +1624,16 @@
     memcpy(&fctxL, fctx, sizeof(limnCbfCtx));
     memcpy(&fctxR, fctx, sizeof(limnCbfCtx));
     if (fctx->verbose) {
-      printf("%s[%u,%u]: dist %g (big %d) --> split at %u\n", me, loi, hii,
+      printf("%s[%u,%u]_%u: dist %g (big %d) --> split at %u\n", me, loi, hii, recDepth,
              fctx->distMax, fctx->distBig, midi);
     }
     if (limnCbfTVT(TL, VM, TR,                 /* */
                    fctx, lpnt, loi, hii, midi, /* */
                    AIR_FALSE /* oneSided */)
-        || limnCbfMulti(path, V0, T1, TL, VM, &fctxL, lpnt, loi, midi)
-        || limnCbfMulti(pRth, VM, TR, T2, V3, &fctxR, lpnt, midi, hii)) {
-      biffAddf(LIMN, "%s[%u,%u]: trouble on recursive fit (midvert %u)", me, loi, hii,
-               midi);
+        || limnCbfMulti(path, V0, T1, TL, VM, recDepth + 1, &fctxL, lpnt, loi, midi)
+        || limnCbfMulti(pRth, VM, TR, T2, V3, recDepth + 1, &fctxR, lpnt, midi, hii)) {
+      biffAddf(LIMN, "%s[%u,%u]_%u: trouble on recursive fit (midvert %u)", me, loi, hii,
+               recDepth, midi);
       limnCbfPathNix(pRth);
       return 1;
     }
@@ -1697,7 +1697,7 @@
     if (fctx->verbose) {
       printf("%s: no corners: finding path to fit point loop", me);
     }
-    if (limnCbfMulti(path, NULL, NULL, NULL, NULL, fctx, lpnt, 0 /* loi */,
+    if (limnCbfMulti(path, NULL, NULL, NULL, NULL, 0, fctx, lpnt, 0 /* loi */,
                      0 /* hii */)) {
       biffAddf(LIMN, "%s: trouble fitting point loop", me);
       return 1;
@@ -1718,7 +1718,7 @@
         printf("%s: finding subpath from between corners [%u,%u] (points [%u,%u])", me,
                cii, cjj, loi, hii);
       }
-      if (limnCbfMulti(subpath, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
+      if (limnCbfMulti(subpath, V0, T1, T2, V3, 0, fctx, lpnt, loi, hii)) {
         biffAddf(LIMN, "%s: trouble with corners [%u,%u] (points [%u,%u])", me, cii, cjj,
                  loi, hii);
         limnCbfPathNix(subpath);

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

[Teem-commit] SF.net SVN: teem:[7192] teem/trunk/src/limn/test/04-m ulti.sh

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

Revision: 7192
          http://sourceforge.net/p/teem/code/7192
Author:   kindlmann
Date:     2024-07-09 22:28:21 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
another test script

Added Paths:
-----------
    teem/trunk/src/limn/test/04-multi.sh

Added: teem/trunk/src/limn/test/04-multi.sh
===================================================================
--- teem/trunk/src/limn/test/04-multi.sh	                        (rev 0)
+++ teem/trunk/src/limn/test/04-multi.sh	2024-07-09 22:28:21 UTC (rev 7192)
@@ -0,0 +1,84 @@
+#!/usr/bin/env bash
+#
+# Teem: Tools to process and visualize scientific data and images
+# Copyright (C) 2009--2024  University of Chicago
+# Copyright (C) 2005--2008  Gordon Kindlmann
+# Copyright (C) 1998--2004  University of Utah
+#
+# This library is free software; you can redistribute it and/or modify it under the terms
+# of the GNU Lesser General Public License (LGPL) as published by the Free Software
+# Foundation; either version 2.1 of the License, or (at your option) any later version.
+# The terms of redistributing and/or modifying this software also include exceptions to
+# the LGPL that facilitate static linking.
+#
+# This library is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+# PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License along with
+# this library; if not, write to Free Software Foundation, Inc., 51 Franklin Street,
+# Fifth Floor, Boston, MA 02110-1301 USA
+#
+
+# 03 = FOURTH stage of testing- another test limnCbfSingle, via lpu cbfit -fs,
+# but with completely free-form input
+
+set -o errexit
+set -o nounset
+shopt -s expand_aliases
+JUNK=""
+function junk { JUNK="$JUNK $@"; }
+function cleanup { rm -rf $JUNK; }
+trap cleanup err exit int term
+# https://devmanual.gentoo.org/tools-reference/bash/
+unset UNRRDU_QUIET_QUIT
+
+# why changing N can significantly change the fit:
+#    because the nrp parms that make sense for a small number of points don't work great
+#    for a huge number of points
+
+# # Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
+# # likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
+# N=18
+# # with -ftt 1 0 -0.8944 0.4472
+# echo "-0.5 0.5
+#  2.0  0.5
+# -0.5  0.0
+#  0.5 -0.5" | ./lpu cbfit -i - -synth $N ...
+
+N=199
+echo "-1 -1
+-1 1
+ 1 1
+ 0 0" | ./lpu cbfit -i - -synthn $N  -syntho 0.txt
+echo "0 0
+-0.5 -0.5
+1 -1
+1 1" | ./lpu cbfit -i - -synthn $N  -syntho 1.txt
+cat 0.txt 1.txt | uniq > xy-inn.txt
+junk {0,1}.txt xy-inn.txt
+
+# IN=xy-inn.txt
+IN=circ.txt
+
+CMD="./lpu cbfit -i $IN -loop -scl 0 -fm 0 -1 -v 1"
+echo "====== $CMD"
+eval $CMD > log.txt
+cat log.txt # ; junk log.txt
+
+OUT=xy-out.txt
+echo "====== RESULTS: --> $OUT"
+grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8
+grep "^        " log.txt | xargs -n 10 echo | cut -d' ' -f 1,2,3,4,5,6,7,8 |
+    ./lpu cbfit -i - -synthn $((2*N)) -sup 1 -syntho $OUT
+junk $OUT
+
+BIN=500
+unu jhisto -i  $IN -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-inn.png
+unu jhisto -i $OUT -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN | unu quantize -b 8 -max 1 -o xy-out.png
+
+
+unu join -i xy-{out,inn,out}.png -a 0 -incr |
+    unu resample -s = x2 x2 -k box -o tmp.png
+
+open tmp.png


Property changes on: teem/trunk/src/limn/test/04-multi.sh
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[Teem-commit] SF.net SVN: teem:[7191] teem/trunk/src/limn

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

Revision: 7191
          http://sourceforge.net/p/teem/code/7191
Author:   kindlmann
Date:     2024-07-09 22:18:49 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
multi-fit may be working now

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

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 21:34:52 UTC (rev 7190)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 22:18:49 UTC (rev 7191)
@@ -219,7 +219,7 @@
       unsigned int ci, si;
       limnCbfPath *spath = limnCbfPathNew(size1);
       airMopAdd(mop, spath, (airMopper)limnCbfPathNix, airMopAlways);
-      printf("%s: synthetically sampling %u splines with %u points", me, size1,
+      printf("%s: synthetically sampling %u splines with %u points\n", me, size1,
              synthNum);
       /* copy in control point data */
       for (si = 0; si < size1; si++) {
@@ -354,14 +354,16 @@
     const double *VTTV[4];
     getLoHi(&loi, &hii, lpnt, fitMultiLoHi[0], fitMultiLoHi[1]);
     getVTTV(VTTV, lpnt, fitTT, loi, hii);
-    if (limnCbfMulti(path, VTTV[0], VTTV[1], VTTV[2], VTTV[3], fctx, lpnt, loi, hii)) {
+    if (limnCbfCtxPrep(fctx, lpnt)
+        || limnCbfMulti(path, VTTV[0], VTTV[1], VTTV[2], VTTV[3], fctx, lpnt, loi,
+                        hii)) {
       airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
       fprintf(stderr, "%s: trouble doing multi fit:\n%s", me, err);
       airMopError(mop);
       return 1;
     }
-    printf("%s: limnCbfMulti results: %u segments %s\n", me, path->segNum,
-           path->isLoop ? "in loop" : "NOT in loop");
+    printf("%s: limnCbfMulti results: %u segments in %s\n", me, path->segNum,
+           path->isLoop ? "loop" : "NOT-loop");
     for (segi = 0; segi < path->segNum; segi++) {
       const limnCbfSeg *seg = path->seg + segi;
       const double *xy = seg->xy;

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 21:34:52 UTC (rev 7190)
+++ teem/trunk/src/limn/splineFit.c	2024-07-09 22:18:49 UTC (rev 7191)
@@ -1098,6 +1098,10 @@
    points falls below fctx->nrpIota * fctx->epsilon, or,
  - parameterization change falls below fctx->nrpDeltaThresh
 Information about the results of this process are set in the given fctx.
+
+This assumes that limnCbfCtxPrep(fctx, lpnt) was called without error!
+That (via ctxBuffersSet) allocates fctx->uu that we depend on here
+(and we fail via biff if it seems like that buffer was not set)
 */
 static int /* Biff: 1 */
 fitSingle(double alpha[2], const double vv0[2], const double tt1[2], const double tt2[2],
@@ -1106,6 +1110,14 @@
   static const char me[] = "fitSingle";
   uint iter, spanlen = spanLength(lpnt, loi, hii);
 
+  if (!(alpha && vv0 && tt1 && tt2 && vv3 && fctx && lpnt)) {
+    biffAddf(LIMN, "%s: got NULL pointer", me);
+    return 1;
+  }
+  if (!(fctx->uu)) {
+    biffAddf(LIMN, "%s: fcgtx->uu NULL; was limnCbfCtxPrep called?", me);
+    return 1;
+  }
   /* DIM=2 pretty much everywhere here */
   if (fctx->verbose) {
     printf("%s[%d,%d]: hello, vv0=(%g,%g), tt1=(%g,%g), "
@@ -1180,7 +1192,7 @@
       int converged = AIR_FALSE;
       for (iter = 0; fctx->distBig && iter < fctx->nrpIterMax; iter++) {
         int simple;
-        if (fctx->verbose) {
+        if (fctx->verbose > 1) {
           printf("%s[%d,%d]: nrp iter %u starting with alpha %g,%g (det %g) (big %d)\n",
                  me, loi, hii, iter, alpha[0], alpha[1], fctx->alphaDet, fctx->distBig);
         }
@@ -1202,7 +1214,7 @@
           biffAddf(LIMN, "%s: trouble", me);
           return 1;
         }
-        if (fctx->verbose) {
+        if (fctx->verbose > 1) {
           printf("%s[%d,%d]: nrp iter %u (reparm) delta = %g (big %d)\n", me, loi, hii,
                  iter, delta, fctx->distBig);
         }

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

[Teem-commit] SF.net SVN: teem:[7190] teem/trunk/src/limn

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

Revision: 7190
          http://sourceforge.net/p/teem/code/7190
Author:   kindlmann
Date:     2024-07-09 21:34:52 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
fixing compilation/valgrind issues

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

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 21:22:33 UTC (rev 7189)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 21:34:52 UTC (rev 7190)
@@ -154,6 +154,7 @@
 
   USAGE(myinfo);
   PARSE(myinfo);
+  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
   airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
 
   if (!(2 == _nin->dim)) {

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 21:22:33 UTC (rev 7189)
+++ teem/trunk/src/limn/splineFit.c	2024-07-09 21:34:52 UTC (rev 7190)
@@ -1033,7 +1033,7 @@
          const double tt1[2], const double tt2[2], const double vv3[2],
          const limnCbfPoints *lpnt, uint loi, uint hii) {
   static const char me[] = "findDist";
-  uint ii, distMaxIdx, spanlen;
+  uint ii, distMaxIdx = UINT_MAX, spanlen;
   double vv1[2], vv2[2], distMax;
   const double *uu = fctx->uu;
 

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

[Teem-commit] SF.net SVN: teem:[7189] teem/trunk/src/limn

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

Revision: 7189
          http://sourceforge.net/p/teem/code/7189
Author:   kindlmann
Date:     2024-07-09 21:22:33 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
making more progress, but now debugging segfault

Modified Paths:
--------------
    teem/trunk/src/limn/limn.h
    teem/trunk/src/limn/lpu_cbfit.c
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/01-test-tvt.sh
    teem/trunk/src/limn/test/02-test-fs.sh
    teem/trunk/src/limn/test/03-single.sh

Modified: teem/trunk/src/limn/limn.h
===================================================================
--- teem/trunk/src/limn/limn.h	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/limn.h	2024-07-09 21:22:33 UTC (rev 7189)
@@ -598,7 +598,9 @@
     cnum;            /* number of corners described by ctvt and cidx */
   /* ----------- output --------- */
   unsigned int nrpIterDone, /* number of nrp iters taken */
-    distMaxIdx;             /* which point had distance distMax */
+    distMaxIdx,             /* which point had distance distMax */
+    nrpSimpleFlop;          /* # times that single-spline fit flip-flopped between a
+                               curved vs a "simple" straight arc */
   double distMax,           /* max distance to given points */
     nrpDeltaDone,           /* latest mean parameterization change by nrp */
     alphaDet;               /* min det of matrix inverted to find alpha */

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-09 21:22:33 UTC (rev 7189)
@@ -26,6 +26,40 @@
 static const char *myinfo
   = (INFO ". \"nrp\" == Newton-based ReParameterization of spline domain");
 
+static void
+getLoHi(unsigned int *loiP, unsigned int *hiiP, const limnCbfPoints *lpnt, int slo,
+        int shi) {
+  unsigned int loi, hii;
+  int pnum = AIR_INT(lpnt->num);
+  if (lpnt->isLoop || (0 == slo && -1 == shi)) {
+    loi = AIR_UINT(AIR_MOD(slo, pnum));
+    hii = AIR_UINT(AIR_MOD(shi, pnum));
+  } else {
+    loi = AIR_UINT(AIR_CLAMP(0, slo, pnum - 1));
+    hii = AIR_UINT(AIR_CLAMP(0, shi, pnum - 1));
+  }
+  *loiP = loi;
+  *hiiP = hii;
+  return;
+}
+
+static void
+getVTTV(const double *VTTV[4], const limnCbfPoints *lpnt, const double fitTT[4],
+        unsigned int loi, unsigned int hii) {
+
+  if (ELL_4V_LEN(fitTT)) {
+    /* help out limnCbfSingle with specific V,T,T,V */
+    VTTV[0] = lpnt->pp + 0 + 2 * loi;
+    VTTV[1] = fitTT + 0;
+    VTTV[2] = fitTT + 2;
+    VTTV[3] = lpnt->pp + 0 + 2 * hii;
+  } else {
+    /* no help will be given */
+    VTTV[0] = VTTV[1] = VTTV[2] = VTTV[3] = NULL;
+  }
+  return;
+}
+
 static int
 limnPu_cbfitMain(int argc, const char **argv, const char *me, hestParm *hparm) {
   hestOpt *hopt = NULL;
@@ -35,9 +69,9 @@
 
   Nrrd *_nin, *nin;
   double *xy, deltaThresh, psi, cangle, epsilon, nrpIota, time0, dtime, scale, nrpCap,
-    synthPow, fitSingleTT[4];
+    synthPow, fitTT[4];
   unsigned int size0, size1, ii, synthNum, pNum, nrpIterMax;
-  int loop, petc, verbose, tvt[4], fitSingleLoHi[2];
+  int loop, petc, verbose, tvt[4], fitSingleLoHi[2], fitMultiLoHi[2], corner2[2];
   char *synthOut, buff[AIR_STRLEN_SMALL + 1];
   limnCbfCtx *fctx;
   limnCbfPath *path;
@@ -96,13 +130,23 @@
                    "the loi and hii indices given here. A negative hii will be "
                    "incremented by the number of points, so -1 works to indicate "
                    "the last point.");
-  hestOptAdd_4_Double(&hopt, "fstt", "T1x T1y T2x T2y", fitSingleTT, "0 0 0 0",
-                      "(if non-zero): help out call to limnCbfSingle by giving these "
-                      "vectors for T1 (outgoing from V0) and T2 (incoming to V3) "
-                      "tangents, so they are not estimated from the data. If this is "
-                      "used; V0 and V3 are set as the first and last points (there "
-                      "is currently no ability to set only some of the 4 vector "
-                      "args to limnCbfSingle)");
+  hestOptAdd_4_Double(&hopt, "ftt", "T1x T1y T2x T2y", fitTT, "0 0 0 0",
+                      "(if non-zero): help out call to either -fs limnCbfSingle or -fm "
+                      "limnCbfMulti by giving these vectors for T1 (outgoing from V0) "
+                      "and T2 (incoming to V3) tangents, so they are not estimated from "
+                      "the data. If this is used; V0 and V3 are set as the first and "
+                      "last points (there is currently no ability to set only some of "
+                      "the 4 vector args to limnCbfSingle or limnCbfMulti)");
+  hestOptAdd_2_Int(&hopt, "fm", "loi hii", fitMultiLoHi, "-1 -1",
+                   "(if loi is >= 0) just do a single call to limnCbfMulti and "
+                   "quit, using the -i input points, fitting a multi-spline path "
+                   "between the loi and hii indices given here. A negative hii will be "
+                   "incremented by the number of points, so -1 works to indicate "
+                   "the last point.");
+  hestOptAdd_2_Int(&hopt, "corn", "val nms", corner2, "0 0",
+                   "if 1st val non-zero: call limnCbfCorners and quit. fctx->cornerFind "
+                   "is set to true if 1st value given here is positive, fctx->cornerNMS "
+                   "is set to !! of second value");
   hestOptAdd_Flag(&hopt, "petc", &petc, "(Press Enter To Continue) ");
   /*
   hestOptAdd_1_String(&hopt, NULL, "output", &out, NULL, "output nrrd filename");
@@ -227,15 +271,12 @@
        sign-ed-ness */
     unsigned int loi, hii, vvi;
     int E, oneSided = !!tvt[3];
-    if (loop) {
-      loi = AIR_UINT(AIR_MOD(tvt[0], pnum));
-      hii = AIR_UINT(AIR_MOD(tvt[1], pnum));
+    if (lpnt->isLoop) {
       vvi = AIR_UINT(AIR_MOD(tvt[2], pnum));
     } else {
-      loi = AIR_UINT(AIR_CLAMP(0, tvt[0], pnum - 1));
-      hii = AIR_UINT(AIR_CLAMP(0, tvt[1], pnum - 1));
       vvi = AIR_UINT(AIR_CLAMP(0, tvt[2], pnum - 1));
     }
+    getLoHi(&loi, &hii, lpnt, tvt[0], tvt[1]);
     E = 0;
     if (!E && fctx->verbose)
       printf("%s: int %d in [%d,%d] --> uint %u in [%u,%u]\n", me, /* */
@@ -259,33 +300,22 @@
     return 0;
   }
 
-  if (fitSingleLoHi[0] >= 0) {
+  if (fitSingleLoHi[0] >= 0) { /* here to call limnCbfSingle once */
+    unsigned int loi, hii;
+    const double *VTTV[4];
     limnCbfSeg seg;
-    int pnum = AIR_INT(lpnt->num);
-    /* re-using the logic from the TVT case above */
-    unsigned int loi = AIR_UINT(AIR_MOD(fitSingleLoHi[0], pnum));
-    unsigned int hii = AIR_UINT(AIR_MOD(fitSingleLoHi[1], pnum));
-    double _V0[2], _V3[2];
-    const double *V0, *T1, *T2, *V3;
-    if (ELL_4V_LEN(fitSingleTT)) {
-      /* help out limnCbfSingle with specific V,T,T,V */
-      ELL_2V_COPY(_V0, lpnt->pp + 0 + 2 * 0);
-      V0 = _V0;
-      ELL_2V_COPY(_V3, lpnt->pp + 0 + 2 * (pnum - 1));
-      V3 = _V3;
-      T1 = fitSingleTT + 0;
-      T2 = fitSingleTT + 2;
-    } else {
-      /* no help will be given */
-      V0 = T1 = T2 = V3 = NULL;
-    }
-    if (limnCbfSingle(&seg, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
+    getLoHi(&loi, &hii, lpnt, fitSingleLoHi[0], fitSingleLoHi[1]);
+    getVTTV(VTTV, lpnt, fitTT, loi, hii);
+    if (limnCbfSingle(&seg, VTTV[0], VTTV[1], VTTV[2], VTTV[3], fctx, lpnt, loi, hii)) {
       airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
       fprintf(stderr, "%s: trouble doing single segment fit:\n%s", me, err);
       airMopError(mop);
       return 1;
     }
-    printf("%s: limnCbfSingle results:\n", me);
+    printf("%s: nrpIterDone %u    nrpSimpleFlop %u    distMax %g @ %u/%u (big %d)\n", me,
+           fctx->nrpIterDone, fctx->nrpSimpleFlop, fctx->distMax, fctx->distMaxIdx,
+           lpnt->num, fctx->distBig);
+    printf("%s: limnCbfSingle spline result:\n", me);
     for (ii = 0; ii < 4; ii++) {
       printf("%g %g\n", seg.xy[0 + 2 * ii], seg.xy[1 + 2 * ii]);
     }
@@ -293,6 +323,55 @@
     return 0;
   }
 
+  if (corner2[0]) { /* here to call limnCbfCorners once */
+    fctx->cornerFind = corner2[0] > 0;
+    fctx->cornerNMS = !!corner2[1];
+    if (limnCbfCorners(fctx, lpnt)) {
+      airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
+      fprintf(stderr, "%s: trouble finding corners:\n%s", me, err);
+      airMopError(mop);
+      return 1;
+    }
+    if (!fctx->cnum) {
+      printf("%s: Zero corners found!\n", me);
+    } else {
+      unsigned int ci;
+      const double *tvt;
+      printf("%s: %u corners found:\n", me, fctx->cnum);
+      for (ci = 0; ci < fctx->cnum; ci++) {
+        tvt = fctx->ctvt + 6 * ci;
+        printf("%3u: vi=%3u  lt=(%g,%g)  vv=(%g,%g)  rt=(%g,%g)\n", ci, fctx->cidx[ci],
+               tvt[0], tvt[1], tvt[2], tvt[3], tvt[4], tvt[5]);
+      }
+    }
+    airMopOkay(mop);
+    return 0;
+  }
+
+  if (fitMultiLoHi[0] >= 0) { /* here to call limnCbfMulti once */
+    unsigned int loi, hii, segi;
+    const double *VTTV[4];
+    getLoHi(&loi, &hii, lpnt, fitMultiLoHi[0], fitMultiLoHi[1]);
+    getVTTV(VTTV, lpnt, fitTT, loi, hii);
+    if (limnCbfMulti(path, VTTV[0], VTTV[1], VTTV[2], VTTV[3], fctx, lpnt, loi, hii)) {
+      airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
+      fprintf(stderr, "%s: trouble doing multi fit:\n%s", me, err);
+      airMopError(mop);
+      return 1;
+    }
+    printf("%s: limnCbfMulti results: %u segments %s\n", me, path->segNum,
+           path->isLoop ? "in loop" : "NOT in loop");
+    for (segi = 0; segi < path->segNum; segi++) {
+      const limnCbfSeg *seg = path->seg + segi;
+      const double *xy = seg->xy;
+      printf("        %g %g    %g %g    %g %g    %g %g    %c %c\n", xy[0], xy[1], xy[2],
+             xy[3], xy[4], xy[5], xy[6], xy[7], seg->corner[0] ? 'C' : '-',
+             seg->corner[1] ? 'C' : '-');
+    }
+    airMopOkay(mop);
+    return 0;
+  }
+
   time0 = airTime();
   if (petc) {
     fprintf(stderr, "%s: Press Enter to Continue ... ", me);

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/splineFit.c	2024-07-09 21:22:33 UTC (rev 7189)
@@ -109,26 +109,29 @@
 
 /*
 TODO:
-test fitSingle uu initialization with points that are already exactly uniform
 test findDist - is distMaxIdx the correct actual index?
 testing corners: corners near and at start==stop of isLoop
 corners not at start or stop of isLoop: do spline wrap around from last to first index?
 
+valgrind everything
+
+(DIM=2) explore what would be required to generalize from 2D to 3D,
+perhaps at least at the API level, even if 3D is not yet implemented
+
+The initialization of uu[] by arc length is especially wrong when the tangents are
+pointing towards each other, and then newton iterations can get stuck in a local minimum.
+As a specific example: with these control points: (-0.5,0.5) (2,0.5) (-0.5,0) (0.5,-0.5)
+sampled 19 times (uniformly in u), fitSingle does great, but sampled 18 times it gets
+stuck on a bad answer.  Would be nice to come up with a heuristic for how to warp
+the initial arc-length parameterization to get closer to correct answer, but this
+requires exploring what is at least a 4-D space of possible splines (lowered from 8-D
+by symmetries).  The cost of not doing this is less economical representations, because
+we would split these segments needlessly.
+
 use performance tests to explore optimal settings in fctx:
   nrpIterMax, nrpCap, nrpIota, nrpPsi, nrpDeltaThresh
 evaluated in terms of time and #splines needed for fit
 (may want to pay in time for more economical representation)
-
-"What GLK hasn't thought through is: what is the interaction of nrp iterations and
-findAlpha generating the simple arc on some but not all iterations (possibly
-unstable?)"
-
-valgrind everything
-
-search for HEY (handle two adjacent equally strong corners)
-
-(DIM=2) explore what would be required to generalized from 2D to 3D,
-perhaps at least at the API level, even if 3D is not yet implemented
 */
 
 #define PNMIN(ISLOOP) ((ISLOOP) ? 4 : 3)
@@ -294,6 +297,7 @@
   /* initialize outputs to bogus valus */
   fctx->nrpIterDone = UINT_MAX;
   fctx->distMaxIdx = UINT_MAX;
+  fctx->nrpSimpleFlop = UINT_MAX;
   fctx->distMax = AIR_POS_INF;
   fctx->nrpDeltaDone = AIR_POS_INF;
   fctx->alphaDet = 0;
@@ -771,7 +775,7 @@
   right tangent if hii == vvi, but will avoid NaNs if scale > 0. Because it will be too
   annoying to require being called in different ways depending on scale, we do *not* do
   error-checking to prevent NaN generation. */
-  if (fctx->verbose) {
+  if (fctx->verbose > 1) {
     printf("%s: (post-idxLift) %u in [%u,%u] -> %u in [%u,%u]\n", me, gvvi, gloi, ghii,
            vvi, loi, hii);
   }
@@ -883,9 +887,9 @@
 }
 
 /*
-(from paper page 620) solves for the alpha that minimize squared error between xy[i]
-and Q(uu[i]) where Q(t) is cubic Bezier spline through vv0, vv0 + alpha[0]*tt1,
-vv3 + alpha[1]*tt2, and vv3.
+(from paper page 620, after "we need only solve"): solves for the alpha[0,1] that
+minimize squared error between xy[i] and Q(uu[i]) where Q(t) is cubic Bezier spline
+through vv0, vv0 + alpha[0]*tt1, vv3 + alpha[1]*tt2, and vv3.
 
 There are various conditions where the generated spline ignores the xy array and
 instead is what one could call a "simple arc" (with control points at 1/3 and 2/3 the
@@ -896,52 +900,52 @@
    author's code)
  - the solved alphas are not convincingly positive
 
-This function is the only place where the "simple arc" is generated, and generating the
-simple arc is not actually an error or problem: if it is bad at fitting the data (as
-determined by findDist) then it may be subdivided, and that's ok. What GLK hasn't thought
-through is: what is the interaction of nrp iterations and findAlpha generating the simple
-arc on some but not all iterations (possibly unstable?)
+This function is the only place where the straight-line "simple arc" is generated, and
+then we return 1, otherwise 0 (when finding alpha[0,1] for a curved segment). But
+generating the simple arc is not actually an error or problem: if the simple arc is bad
+at fitting the data (as determined by findDist) then it may be subdivided, and that's ok.
 */
-static void
-findAlpha(double alpha[2], limnCbfCtx *fctx, /* must be non-NULL */
-          const double vv0[2], const double tt1[2], const double tt2[2],
-          const double vv3[2], const limnCbfPoints *lpnt, uint loi, uint hii) {
+static int /* Biff: nope */
+findAlpha(double alpha[2], limnCbfCtx *fctx, const double vv0[2], const double tt1[2],
+          const double tt2[2], const double vv3[2], const limnCbfPoints *lpnt, uint loi,
+          uint hii) {
   static const char me[] = "findAlpha";
+  int ret;
   uint ii, spanlen = spanLength(lpnt, loi, hii);
   double det, F2L[2], lenF2L;
-  const double *xy = PP(lpnt);
+  const double *xy = PP(lpnt), *uu = fctx->uu;
 
   ELL_2V_SUB(F2L, xy + 2 * hii, xy + 2 * loi); /* DIM=2 throughout this */
   lenF2L = ELL_2V_LEN(F2L);
   if (spanlen > 2) {
+    /* GLK using "m" and "M" instead of author's "C". Note that Ai1 and Ai2 are scalings
+       of (nominally) unit-length tt1 and tt2 by evaluations of the spline basis
+       functions, so they (and the M computed from them, and det(M)), are invariant w.r.t
+       over-all rescalings of the data points */
     double xx[2], m11, m12, m22, MM[4], MI[4];
-    const double *uu = fctx->uu;
     xx[0] = xx[1] = m11 = m12 = m22 = 0;
     for (ii = 0; ii < spanlen; ii++) {
-      const double *xy = PPlowerI(lpnt, AIR_INT(loi + ii));
+      const double *xy = PPlowerI(lpnt, AIR_INT(loi + ii)); /* == paper's "d_i" */
+      double ui = uu[ii];
       double bb[4], Ai1[2], Ai2[2], Pi[2], dmP[2];
-      double ui = uu[ii];
       VCBD0(bb, ui);
       ELL_2V_SCALE(Ai1, bb[1], tt1);
       ELL_2V_SCALE(Ai2, bb[2], tt2);
-      /* GLK using "m" and "M" instead of author's "C". Note that Ai1 and Ai2 are
-         scalings of (nominally) unit-length tt1 and tt2 by evaluations of the spline
-         basis functions, so they (and the M computed from them, and det(M)), are
-         invariant w.r.t over-all rescalings of the data points */
       m11 += ELL_2V_DOT(Ai1, Ai1);
       m12 += ELL_2V_DOT(Ai1, Ai2);
       m22 += ELL_2V_DOT(Ai2, Ai2);
+      /* paper doesn't name what we call P */
       ELL_2V_SCALE_ADD2(Pi, bb[0] + bb[1], vv0, bb[2] + bb[3], vv3);
-      ELL_2V_SUB(dmP, xy, Pi);
-      xx[0] += ELL_2V_DOT(dmP, Ai1);
+      ELL_2V_SUB(dmP, xy, Pi);       /* d minus P */
+      xx[0] += ELL_2V_DOT(dmP, Ai1); /* column vector on right-hand side */
       xx[1] += ELL_2V_DOT(dmP, Ai2);
     }
-    ELL_4V_SET(MM, m11, m12, m12, m22);
+    ELL_4V_SET(MM, m11, m12, m12, m22); /* paper's 2x2 [c11 c12; c21 c22]*/
     ELL_2M_INV(MI, MM, det);
-    ELL_2MV_MUL(alpha, MI, xx);
-  } else {                   /* spanlen <= 2 */
-    det = 1;                 /* bogus but harmless */
-    alpha[0] = alpha[1] = 0; /* trigger simple arc code */
+    ELL_2MV_MUL(alpha, MI, xx); /* solve for alpha[0,1] */
+  } else {                      /* spanlen <= 2 */
+    det = 1;                    /* bogus but harmless */
+    alpha[0] = alpha[1] = 0;    /* trigger simple arc code */
   }
   /* test if we should return simple arc */
   if (!(AIR_EXISTS(det) && AIR_ABS(det) > fctx->detMin
@@ -961,13 +965,15 @@
     but also handle non-unit-length tt1 and tt2 */
     alpha[0] = lenF2L / (3 * ELL_2V_LEN(tt1));
     alpha[1] = lenF2L / (3 * ELL_2V_LEN(tt2));
+    ret = 1;
   } else {
     if (fctx->verbose > 1) {
       printf("%s: all good: det %g, alpha %g,%g\n", me, det, alpha[0], alpha[1]);
     }
+    ret = 0;
   }
   fctx->alphaDet = det;
-  return;
+  return ret;
 }
 
 /*
@@ -1022,7 +1028,7 @@
 
 /* (assuming current parameterization in fctx->uu) sets fctx->distMax to max distance
    to spline, at point fctx->distMaxIdx, and then sets fctx->distBig accordingly */
-static void
+static int /* Biff: 1 */
 findDist(limnCbfCtx *fctx, const double alpha[2], const double vv0[2],
          const double tt1[2], const double tt2[2], const double vv3[2],
          const limnCbfPoints *lpnt, uint loi, uint hii) {
@@ -1032,10 +1038,14 @@
   const double *uu = fctx->uu;
 
   spanlen = spanLength(lpnt, loi, hii);
-  assert(spanlen >= 3);
+  if (!(spanlen >= 3)) {
+    biffAddf(LIMN, "%s: [loi,hii] [%u,%u] -> spanlen %u too small", me, loi, hii,
+             spanlen);
+    return 1;
+  }
   ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1); /* DIM=2 everywhere here */
   ELL_2V_SCALE_ADD2(vv2, 1, vv3, alpha[1], tt2);
-  distMax = -1; /* any computed distance will be >= 0 */
+  distMax = -1.0; /* any computed distance will be >= 0 */
   /* NOTE that the first and last points are actually not part of the max distance
      calculation, which motivates ensuring that the endpoints generated by limnCbfTVT
      are actually sufficiently close to the first and last points (or else the fit spline
@@ -1052,8 +1062,8 @@
     }
   }
   fctx->distMax = distMax;
-  /* we could use a lifted index for internal distMaxIdx but upon saving to fctx it needs
-     to be an actual index */
+  /* we could use a lifted index for internal distMaxIdx,
+     but upon saving to fctx it needs to be an actual index */
   fctx->distMaxIdx = distMaxIdx % lpnt->num;
   fctx->distBig = (distMax <= fctx->nrpIota * fctx->epsilon
                      ? 0
@@ -1066,17 +1076,17 @@
     printf("%s[%u,%u]: distMax %g @ %u (big %d)\n", me, loi, hii, fctx->distMax,
            fctx->distMaxIdx, fctx->distBig);
   }
-  return;
+  return 0;
 }
 
 /*
-fitSingle: fits a single cubic Bezier spline, WITHOUT error checking (limnCbfSingle is an
-error-checking wrapper around this). The given points coordinates are in limnCbfPoints
-lpnt, between low/high indices loi/hii (inclusively); hii can be < loi in the case of a
-point loop. From GIVEN initial endpoint vv0, initial tangent tt1, final tangent tt2
-(pointing backwards), and final endpoint vv3, the job of this function is actually just
-to set alpha[0],alpha[1] such that the cubic Bezier spline with control points vv0,
-vv0+alpha[0]*tt1, vv3+alpha[1]*tt2, vv3 approximates all the given points.
+fitSingle: fits a single cubic Bezier spline, with minimal error checking (limnCbfSingle
+is the error-checking wrapper around this). The given points coordinates are in
+limnCbfPoints lpnt, between low/high indices loi/hii (inclusively); hii can be < loi in
+the case of a point loop. From GIVEN initial endpoint vv0, initial tangent tt1, final
+tangent tt2 (pointing backwards), and final endpoint vv3, the job of this function is
+actually just to set alpha[0],alpha[1] such that the cubic Bezier spline with control
+points vv0, vv0+alpha[0]*tt1, vv3+alpha[1]*tt2, vv3 approximates all the given points.
 
 This is an iterative process, in which alpha is solved for multiples times, after taking
 a Newton step to try to optimize the parameterization of the points (in the
@@ -1089,7 +1099,7 @@
  - parameterization change falls below fctx->nrpDeltaThresh
 Information about the results of this process are set in the given fctx.
 */
-static void
+static int /* Biff: 1 */
 fitSingle(double alpha[2], const double vv0[2], const double tt1[2], const double tt2[2],
           const double vv3[2], limnCbfCtx *fctx, const limnCbfPoints *lpnt, uint loi,
           uint hii) {
@@ -1103,16 +1113,22 @@
            me, loi, hii, vv0[0], vv0[1], tt1[0], tt1[1], tt2[0], tt2[1], vv3[0], vv3[1]);
   }
   if (2 == spanlen) {
-    /* relying on code in findAlpha() that handles slen==2 */
-    findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+    /* relying on code in findAlpha() that handles slen==2; return should be 1 */
+    if (1 != findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
+      biffAddf(LIMN, "%s: what? findAlpha should have returned 1 with spanlen=2", me);
+      return 1;
+    }
     /* nrp is moot */
     fctx->nrpIterDone = 0;
+    fctx->nrpSimpleFlop = 0;
     /* emmulate results of calling findDist() */
     fctx->distMax = fctx->nrpDeltaDone = 0;
     fctx->distMaxIdx = 0;
     fctx->distBig = 0;
-  } else {        /* slen >= 3 */
-    double delta; /* avg parameterization change of interior points */
+  } else {               /* slen >= 3 */
+    double delta;        /* avg parameterization change of interior points */
+    int lastSimple;      /* last return from fitSingle, ==1 if it found "simple" arc */
+    uint simpleFlop = 0; /* # times that return from fitSingle changed */
     {
       /* initialize uu parameterization to chord length */
       unsigned int ii;
@@ -1127,7 +1143,7 @@
         len += ELL_2V_LEN(dd);
         fctx->uu[ii] = len;
         xyM = xyP;
-        /* yes on last iter this is set to wrong coord but it's never used */
+        /* yes on last iter this is set to a coord that's not used */
         xyP = PPlowerI(lpnt, AIR_INT(loi + ii + 1));
       }
       delta = 0;
@@ -1136,7 +1152,7 @@
           fctx->uu[ii] /= len;
           delta += fctx->uu[ii];
         } else {
-          /* ii == pNum-1 last vertex */
+          /* ii == spanlen-1 last vertex */
           fctx->uu[ii] = 1;
         }
         if (fctx->verbose > 1) {
@@ -1143,13 +1159,16 @@
           printf("%s[%d,%d]: initial uu[%u] = %g\n", me, loi, hii, ii, fctx->uu[ii]);
         }
       }
-      delta /= spanlen - 2; /* within the pNum verts are pNum-2 interior verts */
+      delta /= spanlen - 2; /* within the spanlen verts are spanlen-2 interior verts */
       if (fctx->verbose) {
         printf("%s[%d,%d]: initial (chord length) delta = %g\n", me, loi, hii, delta);
       }
     }
-    findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
-    findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+    lastSimple = findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+    if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
+      biffAddf(LIMN, "%s: trouble", me);
+      return 1;
+    }
     if (fctx->verbose) {
       printf(
         "%s[%d,%d]: found alpha %g %g, maxdist %g @ %u (big %d) (%u max nrp iters)\n",
@@ -1160,13 +1179,29 @@
       /* initial fit isn't awful; try making it better with nrp */
       int converged = AIR_FALSE;
       for (iter = 0; fctx->distBig && iter < fctx->nrpIterMax; iter++) {
+        int simple;
         if (fctx->verbose) {
           printf("%s[%d,%d]: nrp iter %u starting with alpha %g,%g (det %g) (big %d)\n",
                  me, loi, hii, iter, alpha[0], alpha[1], fctx->alphaDet, fctx->distBig);
         }
+        /* *this* is where nrp = Newton-based ReParameterization happens */
         delta = reparm(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii);
-        findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
-        findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+        simple = findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
+        simpleFlop += simple != lastSimple;
+        /* seems like a good idea, but not clearly needed, and has some false positives
+        if (simpleFlop > 3 && simpleFlop - 1 > iter / 2) {
+          biffAddf(LIMN,
+                   "%s: (nrp iter %u) previous findAlpha return %d but now %d, and has "
+                   "changed %d times => finding simple arc too unstable under nrp",
+                   me, iter, lastSimple, simple, simpleFlop);
+          return 1;
+        }
+        */
+        lastSimple = simple;
+        if (findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii)) {
+          biffAddf(LIMN, "%s: trouble", me);
+          return 1;
+        }
         if (fctx->verbose) {
           printf("%s[%d,%d]: nrp iter %u (reparm) delta = %g (big %d)\n", me, loi, hii,
                  iter, delta, fctx->distBig);
@@ -1211,11 +1246,12 @@
       }
     }
     fctx->nrpDeltaDone = delta;
+    fctx->nrpSimpleFlop = simpleFlop;
   }
   if (fctx->verbose) {
     printf("%s[%d,%d]: leaving with alpha %g %g\n", me, loi, hii, alpha[0], alpha[1]);
   }
-  return;
+  return 0;
 }
 
 static int
@@ -1267,7 +1303,7 @@
 /*
 limnCbfSingle
 
-Basically and error-checking version of fitSingle; in the limn API because it's needed
+Basically a very error-checked version of fitSingle; in the limn API because it's needed
 for testing. Unlike fitSingle, the geometry info vv0, tt1, tt2, vv3 can either be punted
 on (by passing NULL for all) or not, by passing specific vectors for all. The results are
 converted into the fields in the given limnCbfSeg *seg.  Despite misgivings, we set
@@ -1311,7 +1347,11 @@
     return 1;
   }
   /* now actually do the work */
-  fitSingle(alpha, vttv[0], vttv[1], vttv[2], vttv[3], fctx, lpnt, loi, hii);
+  if (fitSingle(alpha, vttv[0], vttv[1], vttv[2], vttv[3], fctx, lpnt, loi, hii)) {
+    biffAddf(LIMN, "%s: fitSingle failed", me);
+    airMopError(mop);
+    return 1;
+  }
 
   /* process the results to generate info in output limnCbfSeg */
   ELL_2V_COPY(seg->xy + 0, V0);
@@ -1339,7 +1379,8 @@
   static const char me[] = "limnCbfCorners";
   airArray *mop;
   double *angle, /* angle[i] is angle at vertex i */
-    *tvt;        /* 6-by-pnum array of tangent,vertex,tangent */
+    *vtvt;       /* 6-by-pnum array of tangent,vertex,tangent
+                    for ALL vertices */
   int *corny,    /* corny[i] means vertex i seems like a corner */
     oneSided = AIR_TRUE;
   uint pnum = lpnt->num, hii, cnum, vi;
@@ -1353,6 +1394,10 @@
   fctx->cidx = airFree(fctx->cidx);
   fctx->cnum = 0;
 
+  if (fctx->verbose) {
+    printf("%s: hello; cornerFind = %d; cornerNMS = %d\n", me, fctx->cornerFind,
+           fctx->cornerNMS);
+  }
   if (!fctx->cornerFind) {
     /* caller not interested in doing computations to discover corners */
     if (!(lpnt->isLoop)) {
@@ -1373,9 +1418,16 @@
       }
       fctx->cidx[0] = 0;
       fctx->cidx[1] = hii;
+      if (fctx->verbose) {
+        printf("%s: leaving with %u \"corners\" at 1st and last point\n", me,
+               fctx->cnum);
+      }
     }
     return 0;
   }
+  if (fctx->verbose) {
+    printf("%s: looking for corners among %u points\n", me, pnum);
+  }
 
   /* else we search for corners */
   mop = airMopNew();
@@ -1384,15 +1436,15 @@
         && (airMopAdd(mop, angle, airFree, airMopAlways), 1)
         && (corny = AIR_CALLOC(pnum, int))
         && (airMopAdd(mop, corny, airFree, airMopAlways), 1)
-        && (tvt = AIR_CALLOC(6 * pnum, double))
-        && (airMopAdd(mop, tvt, airFree, airMopAlways), 1))) {
+        && (vtvt = AIR_CALLOC(6 * pnum, double))
+        && (airMopAdd(mop, vtvt, airFree, airMopAlways), 1))) {
     biffAddf(LIMN, "%s: trouble allocating working buffers for %u points", me, pnum);
     return 1;
   }
   for (vi = 0; vi < pnum; vi++) {
-    double *LT = tvt + 0 + 6 * vi;
-    double *VV = tvt + 2 + 6 * vi;
-    double *RT = tvt + 4 + 6 * vi;
+    double *LT = vtvt + 0 + 6 * vi;
+    double *VV = vtvt + 2 + 6 * vi;
+    double *RT = vtvt + 4 + 6 * vi;
     /* we find TVT for *every* vertex, despite this seeming like computational overkill.
     Why: we don't know which vertex might be corner until we look at the
     tangent-to-tangent angles for EVERY vertex, for which we don't need to know the
@@ -1417,6 +1469,13 @@
       angle[vi] = 180 * ell_2v_angle_d(LT, RT) / AIR_PI;
       corny[vi] = (angle[vi] < fctx->cornAngle);
     }
+    if (fctx->verbose > 1) {
+      printf("%s: vi=%3u   corny %d   angle %g\n", me, vi, corny[vi], angle[vi]);
+      if (corny[vi]) {
+        printf("    (%g,%g)<--(%g,%g)-->(%g,%g)\n", LT[0], LT[1], VV[0], VV[1], RT[0],
+               RT[1]);
+      }
+    }
   }
   if (fctx->cornerNMS) {
     for (vi = 0; vi < pnum; vi++) {
@@ -1424,12 +1483,27 @@
         /* not a loop, and, either at first or last point ==> must stay a "corner" */
         continue;
       }
-      uint iplus, imnus;
-      iplus = (vi < pnum - 1 ? vi + 1 : (lpnt->isLoop ? 0 : pnum - 1));
-      imnus = (vi ? vi - 1 : (lpnt->isLoop ? pnum - 1 : 0));
-      /* stays a corner only if angle smaller than neighbors */
-      corny[vi] &= (angle[vi] < angle[iplus] && angle[vi] < angle[imnus]);
-      /* HEY handle here if two adjacent corners */
+      /* a little weird to be re-implementing here index lowering, but oh well */
+#define PLUS(I) ((I) < pnum - 1 ? (I) + 1 : (lpnt->isLoop ? 0 : pnum - 1))
+#define MNUS(I) ((I) ? (I) - 1 : (lpnt->isLoop ? pnum - 1 : 0))
+      uint ip1 = PLUS(vi);
+      uint ip2 = PLUS(ip1);
+      uint im1 = MNUS(vi);
+      uint im2 = MNUS(im1);
+#undef PLUS
+#undef MNUS
+      corny[vi]
+        &= (/* stays a corner if angle at vi is smaller (pointier) than neighbors */
+            (angle[im1] > angle[vi] && /* */
+             angle[vi] < angle[ip1])
+            || /* OR, we are part of a pair that is surrounded by less pointy angles */
+            (angle[im1] > angle[vi] &&  /* */
+             angle[vi] == angle[ip1] && /* either the lower index of the pair */
+             angle[ip1] < angle[ip2])   /* */
+            ||                          /* */
+            (angle[im2] > angle[im1] && /* */
+             angle[im1] == angle[vi] && /* or the higher index of the pair */
+             angle[vi] < angle[ip1]));
     }
   }
   /* now, corny[vi] iff vert vi really is a corner; so now count # corners */
@@ -1437,6 +1511,9 @@
   for (vi = 0; vi < pnum; vi++) {
     cnum += !!corny[vi];
   }
+  if (fctx->verbose > 1) {
+    printf("%s: final corner count %u\n", me, cnum);
+  }
   if (cnum) {
     unsigned int ci;
     /* can now allocate new corner-related buffers */
@@ -1447,20 +1524,22 @@
     }
     /* now fill in the corner buffers */
     ci = 0;
-    for (vi = 0; vi < cnum; vi++) {
-      double *id, *od;
-      if (!corny[vi]) continue;
-      fctx->cidx[ci] = vi;
-      id = tvt + 6 * vi;
-      od = fctx->ctvt + 6 * ci;
-      ELL_6V_COPY(od, id);
-      if (fctx->verbose) {
-        printf("%s: corner %u is vertex %u\n  T,V,T = (%g,%g)  (%g,%g)  (%g,%g)\n", me,
-               ci, vi, od[0], od[1], od[2], od[3], od[4], od[5]);
+    for (vi = 0; vi < pnum; vi++) {
+      if (corny[vi]) {
+        double *id, *od;
+        fctx->cidx[ci] = vi;
+        id = vtvt + 6 * vi;
+        od = fctx->ctvt + 6 * ci;
+        ELL_6V_COPY(od, id);
+        if (fctx->verbose) {
+          printf("%s: corner %u is vertex %u\n  T,V,T = (%g,%g)  (%g,%g)  (%g,%g)\n", me,
+                 ci, vi, od[0], od[1], od[2], od[3], od[4], od[5]);
+        }
+        ci++;
       }
-      ci++;
     }
   }
+  fctx->cnum = cnum;
 
   airMopOkay(mop);
   return 0;
@@ -1502,7 +1581,10 @@
   if (fctx->verbose) {
     printf("%s[%u,%u]: trying single fit on all points\n", me, loi, hii);
   }
-  fitSingle(alpha, V0, T1, T2, V3, fctx, lpnt, loi, hii);
+  if (fitSingle(alpha, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
+    biffAddf(LIMN, "%s: fitSingle failed", me);
+    return 1;
+  }
   if (fctx->distBig <= 1) {
     /* max dist was <= fctx->epsilon: single fit was good enough */
     if (fctx->verbose) {

Modified: teem/trunk/src/limn/test/01-test-tvt.sh
===================================================================
--- teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/test/01-test-tvt.sh	2024-07-09 21:22:33 UTC (rev 7189)
@@ -46,11 +46,12 @@
 for I in $(seq 0 $((N-1))); do
     LO=$((I-4))
     HI=$((I+4))
-    # 8-fold TEST:
+    # 16-fold (!) TEST:
     # * without -loop and with -loop
     # * -scl 0 and >0
+    # * LO=HI=0  versus something around I
     # * oneside (4th arg to -tvt) 0 and 1
-    CMD="./lpu cbfit -i $IN -scl 0 -tvt $LO $HI $I 1 -eps 1 -v 0"
+    CMD="./lpu cbfit -i $IN -loop -scl 2 -tvt $LO $HI $I 1 -eps 1 -v 0"
     echo $CMD
     rm -f log.txt
     (eval $CMD 2>&1) > log.txt

Modified: teem/trunk/src/limn/test/02-test-fs.sh
===================================================================
--- teem/trunk/src/limn/test/02-test-fs.sh	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/test/02-test-fs.sh	2024-07-09 21:22:33 UTC (rev 7189)
@@ -32,8 +32,8 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-IN=circ.txt
-#IN=pointy.txt
+#IN=circ.txt
+IN=pointy.txt
 N=$(cat $IN | wc -l | xargs echo)
 BIN=900
 
@@ -47,9 +47,9 @@
 
 for LO in $(seq 0 $((N-1))); do
     echo $LO
-    HI=$((LO+20))
+    HI=$((LO+10))
     LOO=$(printf %02d $LO)
-    CMD="./lpu cbfit -i $IN -loop -scl 0 -psi 1000 -fs $LO $HI -v 0 -eps 0.3"
+    CMD="./lpu cbfit -i $IN -loop -scl 2 -psi 1000 -fs $LO $HI -v 0 -eps 0.01"
     echo "==================== $LO $HI --> test-$LOO.png : $CMD"
     eval $CMD 2>&1 > log-$LOO.txt
     # cat log-$LOO.txt
@@ -56,7 +56,7 @@
     junk log-$LOO.txt
     tail -n 4 log-$LOO.txt |
         ./lpu cbfit -i - -synthn $N -syntho out.txt 2>&1 | grep -v _hestDefaults
-    # lots of extraneous printf thwart piping
+    # lots of extraneous printfs thwart piping
     unu jhisto -i out.txt $MMB | unu quantize -b 8 -max 1 -o out.png
     #  in: green
     # out: magenta

Modified: teem/trunk/src/limn/test/03-single.sh
===================================================================
--- teem/trunk/src/limn/test/03-single.sh	2024-07-09 21:09:16 UTC (rev 7188)
+++ teem/trunk/src/limn/test/03-single.sh	2024-07-09 21:22:33 UTC (rev 7189)
@@ -37,17 +37,25 @@
 #    because the nrp parms that make sense for a small number of points don't work great
 #    for a huge number of points
 
-# Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
-N=18
+# # Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
+# # likely due to initial arc-length parameterization being bad, and nrp stuck in local minima
+# N=18
+# # with -ftt 1 0 -0.8944 0.4472
+# echo "-0.5 0.5
+#  2.0  0.5
+# -0.5  0.0
+#  0.5 -0.5" | unu 2op x - 1 | ...
 
-echo "-0.5 0.5
- 2.0  0.5
--0.5  0.0
- 0.5 -0.5" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
-# junk xy-inn-$N.txt
+N=80
+echo "-0.5 -1
+ -1 1
+1 1
+0.5 -1" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
+junk xy-inn-$N.txt
 
-#./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 0 -psi 1000000000 -eps 0.001 -nim 8000 -deltathr 0.0001 -iota 0.01 -cap 10 2>&1 > log.txt
-./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 3 -psi 1000000000 -fstt 1 0 -0.8944 0.4472 -nim 400 -deltathr 0.000000000001 2>&1 > log.txt
+CMD="./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 0 -psi 1000000000 -nim 400 -deltathr 0.000000000001"
+echo "====== $CMD"
+eval $CMD > log.txt
 cat log.txt; junk log.txt
 tail -n 4 log.txt | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-out-$N.txt
 junk xy-out-$N.txt

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

[Teem-commit] SF.net SVN: teem:[7188] teem/trunk/src/limn/test/lpu .c

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

Revision: 7188
          http://sourceforge.net/p/teem/code/7188
Author:   kindlmann
Date:     2024-07-09 21:09:16 +0000 (Tue, 09 Jul 2024)
Log Message:
-----------
number of colummns set by terminal

Modified Paths:
--------------
    teem/trunk/src/limn/test/lpu.c

Modified: teem/trunk/src/limn/test/lpu.c
===================================================================
--- teem/trunk/src/limn/test/lpu.c	2024-07-08 23:37:12 UTC (rev 7187)
+++ teem/trunk/src/limn/test/lpu.c	2024-07-09 21:09:16 UTC (rev 7188)
@@ -73,7 +73,7 @@
   hparm->elideMultipleEmptyStringDefault = AIR_TRUE;
   hparm->cleverPluralizeOtherY = AIR_TRUE;
   hparm->respectDashDashHelp = AIR_TRUE;
-  hparm->columns = 78;
+  hestParmColumnsIoctl(hparm, hestDefaultColumns);
 
   /* if there are no arguments, then we give general usage information */
   if (1 >= argc) {

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

[Teem-commit] SF.net SVN: teem:[7187] teem/trunk/src/limn

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

Revision: 7187
          http://sourceforge.net/p/teem/code/7187
Author:   kindlmann
Date:     2024-07-08 23:37:12 +0000 (Mon, 08 Jul 2024)
Log Message:
-----------
isolated a case where single fit is bad

Modified Paths:
--------------
    teem/trunk/src/limn/lpu_cbfit.c
    teem/trunk/src/limn/splineFit.c
    teem/trunk/src/limn/test/02-test-fs.sh
    teem/trunk/src/limn/test/03-single.sh

Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c	2024-07-08 22:33:25 UTC (rev 7186)
+++ teem/trunk/src/limn/lpu_cbfit.c	2024-07-08 23:37:12 UTC (rev 7187)
@@ -35,7 +35,7 @@
 
   Nrrd *_nin, *nin;
   double *xy, deltaThresh, psi, cangle, epsilon, nrpIota, time0, dtime, scale, nrpCap,
-    synthPow;
+    synthPow, fitSingleTT[4];
   unsigned int size0, size1, ii, synthNum, pNum, nrpIterMax;
   int loop, petc, verbose, tvt[4], fitSingleLoHi[2];
   char *synthOut, buff[AIR_STRLEN_SMALL + 1];
@@ -91,11 +91,18 @@
   hestOptAdd_1_Double(&hopt, "cap", "cap", &nrpCap, buff,
                       "nrp cap parameterization change");
   hestOptAdd_2_Int(&hopt, "fs", "loi hii", fitSingleLoHi, "-1 -1",
-                   "(if loi is >= 0): just do a single call to limnCbfSingle and "
+                   "(if loi is >= 0) just do a single call to limnCbfSingle and "
                    "quit, using the -i input points, and fitting a spline between "
                    "the loi and hii indices given here. A negative hii will be "
                    "incremented by the number of points, so -1 works to indicate "
                    "the last point.");
+  hestOptAdd_4_Double(&hopt, "fstt", "T1x T1y T2x T2y", fitSingleTT, "0 0 0 0",
+                      "(if non-zero): help out call to limnCbfSingle by giving these "
+                      "vectors for T1 (outgoing from V0) and T2 (incoming to V3) "
+                      "tangents, so they are not estimated from the data. If this is "
+                      "used; V0 and V3 are set as the first and last points (there "
+                      "is currently no ability to set only some of the 4 vector "
+                      "args to limnCbfSingle)");
   hestOptAdd_Flag(&hopt, "petc", &petc, "(Press Enter To Continue) ");
   /*
   hestOptAdd_1_String(&hopt, NULL, "output", &out, NULL, "output nrrd filename");
@@ -258,7 +265,21 @@
     /* re-using the logic from the TVT case above */
     unsigned int loi = AIR_UINT(AIR_MOD(fitSingleLoHi[0], pnum));
     unsigned int hii = AIR_UINT(AIR_MOD(fitSingleLoHi[1], pnum));
-    if (limnCbfSingle(&seg, NULL, NULL, NULL, NULL, fctx, lpnt, loi, hii)) {
+    double _V0[2], _V3[2];
+    const double *V0, *T1, *T2, *V3;
+    if (ELL_4V_LEN(fitSingleTT)) {
+      /* help out limnCbfSingle with specific V,T,T,V */
+      ELL_2V_COPY(_V0, lpnt->pp + 0 + 2 * 0);
+      V0 = _V0;
+      ELL_2V_COPY(_V3, lpnt->pp + 0 + 2 * (pnum - 1));
+      V3 = _V3;
+      T1 = fitSingleTT + 0;
+      T2 = fitSingleTT + 2;
+    } else {
+      /* no help will be given */
+      V0 = T1 = T2 = V3 = NULL;
+    }
+    if (limnCbfSingle(&seg, V0, T1, T2, V3, fctx, lpnt, loi, hii)) {
       airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
       fprintf(stderr, "%s: trouble doing single segment fit:\n%s", me, err);
       airMopError(mop);

Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c	2024-07-08 22:33:25 UTC (rev 7186)
+++ teem/trunk/src/limn/splineFit.c	2024-07-08 23:37:12 UTC (rev 7187)
@@ -1194,7 +1194,7 @@
           printf("converged! with maxdist %g @ %u (big %d)\n", fctx->distMax,
                  fctx->distMaxIdx, fctx->distBig);
         } else if (!fctx->distBig) {
-          printf("NICE small dist %g (<%g) @ %u\n", fctx->distMax, fctx->epsilon,
+          printf("NICE small dist %g (<%g) @ %u (big 0)\n", fctx->distMax, fctx->epsilon,
                  fctx->distMaxIdx);
         } else {
           printf("hit nrp itermax %u; maxdist %g @ %u (big %d)\n", fctx->nrpIterMax,

Modified: teem/trunk/src/limn/test/02-test-fs.sh
===================================================================
--- teem/trunk/src/limn/test/02-test-fs.sh	2024-07-08 22:33:25 UTC (rev 7186)
+++ teem/trunk/src/limn/test/02-test-fs.sh	2024-07-08 23:37:12 UTC (rev 7187)
@@ -32,8 +32,8 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-IN=pointy.txt
-#IN=circ.txt
+IN=circ.txt
+#IN=pointy.txt
 N=$(cat $IN | wc -l | xargs echo)
 BIN=900
 
@@ -47,9 +47,9 @@
 
 for LO in $(seq 0 $((N-1))); do
     echo $LO
-    HI=$((LO+7))
+    HI=$((LO+20))
     LOO=$(printf %02d $LO)
-    CMD="./lpu cbfit -i $IN -loop -scl 3 -psi 1000 -fs $LO $HI -v 0 -eps 0.3"
+    CMD="./lpu cbfit -i $IN -loop -scl 0 -psi 1000 -fs $LO $HI -v 0 -eps 0.3"
     echo "==================== $LO $HI --> test-$LOO.png : $CMD"
     eval $CMD 2>&1 > log-$LOO.txt
     # cat log-$LOO.txt

Modified: teem/trunk/src/limn/test/03-single.sh
===================================================================
--- teem/trunk/src/limn/test/03-single.sh	2024-07-08 22:33:25 UTC (rev 7186)
+++ teem/trunk/src/limn/test/03-single.sh	2024-07-08 23:37:12 UTC (rev 7187)
@@ -33,29 +33,32 @@
 # https://devmanual.gentoo.org/tools-reference/bash/
 unset UNRRDU_QUIET_QUIT
 
-# why changing this from 200 to 400 to 800 can significantly change the fit
-# because the nrp parms that make sense for a small number of points don't work great
-# for a huge number of points
-N=100
+# why changing N can significantly change the fit:
+#    because the nrp parms that make sense for a small number of points don't work great
+#    for a huge number of points
 
-echo "-0.7 0.7
-2 0.7
-0 0.7
-0.7 -0.7" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-0.txt
-junk xy-0.txt
+# Good debugging test case, N=18 is a bad fit, N=19 is a perfect fit
+N=18
 
-./lpu cbfit -i xy-0.txt -scl 0 -fs 0 -1 -v 0 -psi 1000000000 -eps 0.001 -nim 8000 -deltathr 0.0001 -iota 0.01 -cap 10 2>&1 > log.txt
+echo "-0.5 0.5
+ 2.0  0.5
+-0.5  0.0
+ 0.5 -0.5" | unu 2op x - 1 | unu 2op + - 0.0 | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-inn-$N.txt
+# junk xy-inn-$N.txt
+
+#./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 0 -psi 1000000000 -eps 0.001 -nim 8000 -deltathr 0.0001 -iota 0.01 -cap 10 2>&1 > log.txt
+./lpu cbfit -i xy-inn-$N.txt -scl 0 -fs 0 -1 -v 3 -psi 1000000000 -fstt 1 0 -0.8944 0.4472 -nim 400 -deltathr 0.000000000001 2>&1 > log.txt
 cat log.txt; junk log.txt
-tail -n 4 log.txt | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-1.txt
-junk xy-1.txt
+tail -n 4 log.txt | ./lpu cbfit -i - -synthn $N -sup 1 -syntho xy-out-$N.txt
+junk xy-out-$N.txt
 
 BIN=500
-for I in 0 1; do
-    unu jhisto -i xy-$I.txt -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN |
-        unu quantize -b 8 -max 1 -o xy-$I.png
+for X in inn out; do
+    unu jhisto -i xy-$X-$N.txt -min -1.1 1.1 -max 1.1 -1.1 -b $BIN $BIN |
+        unu quantize -b 8 -max 1 -o xy-$X.png
 done
 
-unu join -i xy-{1,0,1}.png -a 0 -incr |
+unu join -i xy-{out,inn,out}.png -a 0 -incr |
     unu resample -s = x2 x2 -k box -o tmp.png
 
 open tmp.png

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