Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[7172] teem/trunk/src/limn
|
From: <kin...@us...> - 2024-06-30 00:59:41
|
Revision: 7172
http://sourceforge.net/p/teem/code/7172
Author: kindlmann
Date: 2024-06-30 00:59:39 +0000 (Sun, 30 Jun 2024)
Log Message:
-----------
limnCBFitSingle seems to work on tests so far; more to come ...
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-06-28 04:44:18 UTC (rev 7171)
+++ teem/trunk/src/limn/limn.h 2024-06-30 00:59:39 UTC (rev 7172)
@@ -911,12 +911,13 @@
int isLoop);
LIMN_EXPORT limnCBFPoints *limnCBFPointsNix(limnCBFPoints *lpnt);
LIMN_EXPORT int limnCBFPointsCheck(const limnCBFPoints *lpnt);
+LIMN_EXPORT limnCBFPath *limnCBFPathNew(void);
+LIMN_EXPORT limnCBFPath *limnCBFPathNix(limnCBFPath *path);
+LIMN_EXPORT void limnCBFPathJoin(limnCBFPath *dst, const limnCBFPath *src);
LIMN_EXPORT limnCBFCtx *limnCBFCtxNew();
LIMN_EXPORT limnCBFCtx *limnCBFCtxNix(limnCBFCtx *fctx);
LIMN_EXPORT int limnCBFCtxPrep(limnCBFCtx *fctx, const limnCBFPoints *lpnt);
LIMN_EXPORT void limnCBFSegEval(double *xy, const limnCBFSeg *seg, double tt);
-LIMN_EXPORT limnCBFPath *limnCBFPathNew(void);
-LIMN_EXPORT limnCBFPath *limnCBFPathNix(limnCBFPath *path);
LIMN_EXPORT void limnCBFPathSample(double *xy, unsigned int pointNum,
const limnCBFPath *path);
LIMN_EXPORT int limnCBFFindTVT(double lt[2], double vv[2], double rt[2],
@@ -923,19 +924,20 @@
const limnCBFCtx *fctx, const limnCBFPoints *lpnt,
unsigned int loi, unsigned int hii, unsigned int ofi,
int oneSided);
-LIMN_EXPORT int limnCBFCtxInit(const limnCBFCtx *fctx, const limnCBFPoints *lpnt);
LIMN_EXPORT int limnCBFitSingle(limnCBFSeg *seg, const double vv0[2],
const double tt1[2], const double tt2[2],
- const double vv3[2], limnCBFCtx *fctx, const double *xy,
- unsigned int pointNum);
+ const double vv3[2], limnCBFCtx *fctx,
+ limnCBFPoints *lpnt, unsigned int loi, unsigned int hii);
+#if 0
+LIMN_EXPORT int limnCBFCorners(unsigned int **cornIdx, unsigned int *cornNum,
+ limnCBFCtx *fctx, const limnCBFPoints *lpnt);
LIMN_EXPORT 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);
-LIMN_EXPORT int limnCBFCorners(unsigned int **cornIdx, unsigned int *cornNum,
- limnCBFCtx *fctx, const limnCBFPoints *lpnt);
LIMN_EXPORT int limnCBFit(limnCBFPath *path, limnCBFCtx *fctx,
const limnCBFPoints *lpnt);
+#endif
/* lpu{Flotsam,. . .}.c */
#define LIMN_DECLARE(C) LIMN_EXPORT const unrrduCmd limnPu_##C##Cmd;
Modified: teem/trunk/src/limn/lpu_cbfit.c
===================================================================
--- teem/trunk/src/limn/lpu_cbfit.c 2024-06-28 04:44:18 UTC (rev 7171)
+++ teem/trunk/src/limn/lpu_cbfit.c 2024-06-30 00:59:39 UTC (rev 7172)
@@ -1,6 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
- Copyright (C) 2009--2023 University of Chicago
+ Copyright (C) 2009--2024 University of Chicago
Copyright (C) 2005--2008 Gordon Kindlmann
Copyright (C) 1998--2004 University of Utah
@@ -36,7 +36,7 @@
Nrrd *_nin, *nin;
double *xy, deltaThresh, psi, cangle, epsilon, nrpIota, time0, dtime, scale, synthPow;
unsigned int ii, synthNum, pNum, nrpIterMax;
- int loop, petc, verbose, tvt[4], fitSingle;
+ int loop, petc, verbose, tvt[4], fitSingleLoHi[2];
char *synthOut;
limnCBFCtx *fctx;
limnCBFPath *path;
@@ -47,10 +47,10 @@
"-i input xy points are actually a loop: the first point logically "
"follows the last point");
hestOptAdd_1_Int(&hopt, "v", "verbose", &verbose, "1", "verbosity level");
- hestOptAdd_1_UInt(&hopt, "sn", "num", &synthNum, "51",
+ hestOptAdd_1_UInt(&hopt, "synthn", "num", &synthNum, "51",
"if saving spline sampling to -so, how many sample.");
hestOptAdd_1_String(
- &hopt, "so", "synth out", &synthOut, "",
+ &hopt, "syntho", "synth out", &synthOut, "",
"if non-empty, input xy points are actually control points for a single spline "
"segment, to be sampled -sn times, and this is the filename into which to save "
"the synthesized xy pts, and then quit without any fitting.");
@@ -75,8 +75,12 @@
hestOptAdd_1_Double(&hopt, "ca", "angle", &cangle, "100", "angle indicating a corner");
hestOptAdd_1_Double(&hopt, "scl", "scale", &scale, "0",
"scale for geometry estimation");
- hestOptAdd_Flag(&hopt, "fs", &fitSingle,
- "just do a single call to limnCBFitSingle and quit");
+ hestOptAdd_2_Int(&hopt, "fs", "loi hii", fitSingleLoHi, "-1 -1",
+ "(if loi is >= 0): just do a single call to limnCBFitSingle 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_Flag(&hopt, "petc", &petc, "(Press Enter To Continue) ");
/*
hestOptAdd_1_String(&hopt, NULL, "output", &out, NULL, "output nrrd filename");
@@ -204,10 +208,13 @@
return 0;
}
- if (fitSingle) {
- int ii;
+ if (fitSingleLoHi[0] >= 0) {
limnCBFSeg seg;
- if (limnCBFitSingle(&seg, NULL, NULL, NULL, NULL, fctx, lpnt->pp, lpnt->num)) {
+ 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));
+ if (limnCBFitSingle(&seg, NULL, NULL, NULL, NULL, 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);
@@ -227,6 +234,7 @@
fflush(stderr);
getchar();
}
+#if 0 /* HEY */
if (limnCBFit(path, fctx, lpnt)) {
airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
fprintf(stderr, "%s: trouble doing fitting:\n%s", me, err);
@@ -233,6 +241,7 @@
airMopError(mop);
return 1;
}
+#endif
dtime = (airTime() - time0) * 1000;
printf("%s: time= %g ms;iterDone= %u ;deltaDone=%g, distMax=%g (@%u)\n", me, dtime,
fctx->nrpIterDone, fctx->nrpDeltaDone, fctx->distMax, fctx->distMaxIdx);
Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c 2024-06-28 04:44:18 UTC (rev 7171)
+++ teem/trunk/src/limn/splineFit.c 2024-06-30 00:59:39 UTC (rev 7172)
@@ -158,7 +158,7 @@
return NULL;
}
-static void
+void
limnCBFPathJoin(limnCBFPath *dst, const limnCBFPath *src) {
uint bb = airArrayLenIncr(dst->segArr, src->segNum);
memcpy(dst->seg + bb, src->seg, (src->segNum) * sizeof(limnCBFSeg));
@@ -508,6 +508,15 @@
return;
}
+/* utility function for counting how many vertices are in 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 */
+static uint
+spanLength(const limnCBFPoints *lpnt, uint loi, uint hii) {
+ uint topi = hii + (hii < loi) * (lpnt->num);
+ return topi - loi + 1;
+}
+
/* cheesy macro as short-hand to access either pp or ppOwn */
#define PP(lpnt) ((lpnt)->pp ? (lpnt)->pp : (lpnt)->ppOwn)
@@ -731,15 +740,6 @@
return 0;
}
-/* utility function for counting how many vertices are in 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 */
-static uint
-spanLength(const limnCBFPoints *lpnt, uint loi, uint hii) {
- uint topi = hii + (hii < loi) * (lpnt->num);
- return topi - loi + 1;
-}
-
/* utility function for getting pointer to coords in lpnt, for point with index loi+ofi,
while accounting for possibility of wrapping */
static const double *
@@ -937,12 +937,12 @@
}
/*
-fitSingle: fits a single cubic Bezier spline, w/out error checking (limnCBFSingle is a
-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,
+fitSingle: fits a single cubic Bezier spline, w/out 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.
This is an iterative process, in which alpha is solved for multiples times, after taking
@@ -1018,9 +1018,10 @@
findAlpha(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
findDist(fctx, alpha, vv0, tt1, tt2, vv3, lpnt, loi, hii);
if (fctx->verbose) {
- printf("%s[%d,%d]: found alpha %g %g, dist %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 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);
}
if (fctx->distBig < 3) {
/* initial fit isn't awful; try making it better with nrp */
@@ -1039,8 +1040,8 @@
}
if (delta <= fctx->nrpDeltaThresh) {
if (fctx->verbose) {
- printf("%s[%d,%d]: nrp iter %u delta %g <= min %g --> break\n", me, loi, hii,
- iter, delta, fctx->nrpDeltaThresh);
+ printf("%s[%d,%d]: nrp iter %u delta %g <= thresh %g --> break\n", me, loi,
+ hii, iter, delta, fctx->nrpDeltaThresh);
}
converged = AIR_TRUE;
break;
@@ -1049,11 +1050,14 @@
if (fctx->verbose) {
printf("%s[%d,%d]: nrp done after %u iters: ", me, loi, hii, iter);
if (converged) {
- printf("converged!\n");
+ printf("converged! with maxdist %g @ %u (big %d)\n", fctx->distMax,
+ fctx->distMaxIdx, fctx->distBig);
} else if (!fctx->distBig) {
- printf("nice small dist %g @ %u\n", fctx->distMax, fctx->distMaxIdx);
+ printf("NICE small dist %g (<%g) @ %u\n", fctx->distMax, fctx->epsilon,
+ fctx->distMaxIdx);
} else {
- printf("hit itermax %u\n", fctx->nrpIterMax);
+ printf("hit nrp itermax %u; maxdist %g @ %u (big %d)\n", fctx->nrpIterMax,
+ fctx->distMax, fctx->distMaxIdx, fctx->distBig);
}
}
fctx->nrpIterDone = iter;
@@ -1076,45 +1080,47 @@
/*
limnCBFitSingle
-builds a limnCBFPoints around given xy, determines spline constraints if necessary, and
-calls fitSingle
+Basically and error-checking 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
+both seg->corner[0,1] to AIR_TRUE.
+
+Perservating on seg->corner[0,1]: we really don't have the information to consistently
+set them with certainty. If not given the geometry vectors, we do assert oneSided when
+estimating the vertices and tangents, so maybe then we can set set->corner[0,1] to true,
+but on the other hand we don't know what to do when the geometry vectors are given. But
+it wouldn't be cool to sometimes set fields in the output struct and sometimes not.
+
+This function used to also be lower-overhead, by not requiring a fctx, and taking a coord
+data pointer instead of a lpnts. But for the sake of testing, there needed to be way of
+passing specific loi and hii in a point loop, so that's what this accepts now. As a
+result, that means this function isn't as convenient a function for one-off single-spline
+fits, and at this point limn doesn't have one. The real entry point for spline fitting
+is limnCBFit().
*/
int /* Biff: 1 */
limnCBFitSingle(limnCBFSeg *seg, const double vv0[2], const double tt1[2],
- const double tt2[2], const double vv3[2], limnCBFCtx *_fctx,
- const double *xy, uint pNum) {
+ const double tt2[2], const double vv3[2], limnCBFCtx *fctx,
+ limnCBFPoints *lpnt, uint loi, uint hii) {
static const char me[] = "limnCBFitSingle";
double v0c[2], t1c[2], t2c[2], v3c[2]; /* locally computed geometry info */
+ const double *v0p, *t1p, *t2p, *v3p; /* pointers for the geometry info */
double alpha[2]; /* recovered by fitting */
- const double *v0p, *t1p, *t2p, *v3p; /* pointers for the geometry info */
- uint loi, hii;
- limnCBFCtx *fctx;
- limnCBFPoints *lpnt;
+ uint spanlen;
airArray *mop;
- if (!(seg && xy && pNum)) {
- biffAddf(LIMN, "%s: got NULL pointer or 0 points", me);
+ if (!(seg && fctx && lpnt)) {
+ biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
mop = airMopNew();
- lpnt = limnCBFPointsNew(xy, nrrdTypeDouble, 2, pNum, AIR_FALSE /* isLoop */);
- airMopAdd(mop, lpnt, (airMopper)limnCBFPointsNix, airMopAlways);
- loi = 0;
- hii = pNum - 1;
- if (_fctx) {
- fctx = _fctx; /* caller has supplied info */
- } else {
- fctx = limnCBFCtxNew();
- airMopAdd(mop, fctx, (airMopper)limnCBFCtxNix, airMopAlways);
- fctx->scale = 0;
- fctx->epsilon = 0.01; /* HEY maybe instead some multiple of stdv? */
- /* (if you don't like these hard-coded defaults then pass your own fctx) */
- }
if (limnCBFCtxPrep(fctx, lpnt)) {
- biffAddf(LIMN, "%s: problem with %s fctx or lpnt", me, _fctx ? "given" : "own");
+ biffAddf(LIMN, "%s: problem with fctx or lpnt", me);
airMopError(mop);
return 1;
}
+ spanlen = spanLength(lpnt, loi, hii);
if (vv0 && tt1 && tt2 && vv3) {
/* point to the given containers of geometry info */
v0p = vv0; /* DIM=2 ? */
@@ -1131,11 +1137,11 @@
airMopError(mop);
return 1;
}
- /* have NO given geometry info; must find it all */
+ /* do not have geometry info; must find it all */
if (limnCBFFindTVT(NULL, v0c, t1c, /* */
fctx, lpnt, loi, hii, 0 /* ofi */, oneSided)
|| limnCBFFindTVT(t2c, v3c, NULL, /* */
- fctx, lpnt, loi, hii, hii /* ofi */, oneSided)) {
+ fctx, lpnt, loi, hii, spanlen - 1 /* ofi */, oneSided)) {
biffAddf(LIMN, "%s: trouble finding geometry info", me);
airMopError(mop);
return 1;
@@ -1149,14 +1155,16 @@
t2p = t2c;
v3p = v3c;
}
+ /* HERE is the call that actually does the work */
fitSingle(alpha, fctx, v0p, t1p, t2p, v3p, lpnt, loi, hii);
+ /* process the results to generate info in output limnCBFSeg */
ELL_2V_COPY(seg->xy + 0, v0p);
ELL_2V_SCALE_ADD2(seg->xy + 2, 1, v0p, alpha[0], t1p);
ELL_2V_SCALE_ADD2(seg->xy + 4, 1, v3p, alpha[1], t2p);
ELL_2V_COPY(seg->xy + 6, v3p);
- /* HEY is it our job to set seg->corner[] ?!? */
- seg->pointNum = pNum;
+ seg->corner[0] = seg->corner[1] = AIR_TRUE; /* misgivings . . . */
+ seg->pointNum = spanlen;
airMopOkay(mop);
return 0;
@@ -1163,91 +1171,6 @@
}
/*
-******** limnCBFit
-**
-** top-level function for fitting cubic beziers to given points
-*/
-int /* Biff: 1 */
-limnCBFit(limnCBFPath *pathWhole, limnCBFCtx *fctx, const limnCBFPoints *lpnt) {
- static const char me[] = "limnCBFit";
- uint *cornIdx, /* (if non-NULL) array of logical indices into PP(lpnt) of corners */
- cornNum, /* length of cornIdx array */
- cii;
- int ret;
- airArray *mop;
-
- if (!(pathWhole && fctx && lpnt)) {
- biffAddf(LIMN, "%s: got NULL pointer", me);
- return 1;
- }
- if (limnCBFCtxPrep(fctx, lpnt)) {
- biffAddf(LIMN, "%s: trouble preparing", me);
- return 1;
- }
- if (fctx->cornerFind) {
- /* HEY RESURRECT ME
- if (limnCBFCorners(&cornIdx, &cornNum, fctx, lpnt)) {
- biffAddf(LIMN, "%s: trouble finding corners", me);
- return 1;
- }
- */
- } else {
- if (lpnt->isLoop) {
- /* there really are no corners */
- cornIdx = NULL;
- cornNum = 0;
- } else {
- /* even without "corners": if not a loop, first and last verts act as corners */
- cornIdx = AIR_CALLOC(2, uint);
- assert(cornIdx);
- cornNum = 2;
- cornIdx[0] = 0;
- cornIdx[1] = lpnt->num - 1;
- }
- }
- mop = airMopNew();
- if (cornIdx) {
- airMopAdd(mop, cornIdx, airFree, airMopAlways);
- }
- airArrayLenSet(pathWhole->segArr, 0);
- if (!cornNum) {
-/* no corners; do everything with one multi call */
-#if 0 /* HEY RESURRECT ME */
- if (limnCBFMulti(pathWhole, fctx, NULL, NULL, NULL, NULL, lpnt, 0 /* loi */,
- 0 /* hii */)) {
- biffAddf(LIMN, "%s: trouble", me);
- airMopError(mop);
- return 1;
- }
-#endif
- } else {
- /* do have corners: split points into segments between corners. The corner vertex is
- both last point in segment I and first point in segment I+1 */
- for (cii = 0; cii < cornNum; cii++) {
- uint cjj = (cii + 1) % cornNum;
- uint loi = cornIdx[cii];
- uint hii = cornIdx[cjj];
- limnCBFPath *subpath = limnCBFPathNew();
-#if 0 /* HEY RESURRECT ME */
- ret = limnCBFMulti(subpath, fctx, NULL, NULL, NULL, NULL, lpnt, loi, hii);
- if (ret) {
- biffAddf(LIMN, "%s: trouble from corners [%u,%u] (points [%u,%u])", me, cii,
- cjj, loi, hii); limnCBFPathNix(subpath); airMopError(mop); return 1;
- }
-#endif
- subpath->seg[0].corner[0] = 1;
- subpath->seg[subpath->segNum - 1].corner[1] = 1;
- limnCBFPathJoin(pathWhole, subpath);
- limnCBFPathNix(subpath);
- }
- }
-
- pathWhole->isLoop = lpnt->isLoop;
- airMopOkay(mop);
- return 0;
-}
-
-/*
TODO:
-- avoid repeated calls to limnCBFCtxPrep
@@ -1267,8 +1190,6 @@
findAlpha generating the simple arc on some but not all iterations (possibly
unstable?)"
-resurrect in segment struct: how many points were represented by a single spline
-
reparm: "HEY TODO: need to make sure that half-way between points,
spline isn't wildly diverging; this can happen with the
spline making a loop away from a small number of points, e.g.:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
