Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[7189] teem/trunk/src/limn
|
From: <kin...@us...> - 2024-07-09 21:22:41
|
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.
|
