Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[7166] teem/trunk/src/limn
|
From: <kin...@us...> - 2024-06-21 07:43:19
|
Revision: 7166
http://sourceforge.net/p/teem/code/7166
Author: kindlmann
Date: 2024-06-21 07:43:16 +0000 (Fri, 21 Jun 2024)
Log Message:
-----------
finally re-started progress on limnCBF functions
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-20 21:03:47 UTC (rev 7165)
+++ teem/trunk/src/limn/limn.h 2024-06-21 07:43:16 UTC (rev 7166)
@@ -514,11 +514,12 @@
******** limnCBFSeg
**
** 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 )
*/
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 */
+ 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 */
@@ -537,13 +538,9 @@
} limnCBFPath;
/*
-******** limnCBFContext
+******** limnCBFCtx
**
-** The bag of state for limnCBF functions. Callers of limnCBF functions do not
-** need to worry about the dynamically allocated things within (so: no
-** limnCBFContextNew or limnCBFContextNix), but a limnCBFContext variable
-** should be initialized with limnCBFContextInit() in order to set default
-** parameters, before passing to limnCBF functions.
+** The bag of state for limnCBF functions.
**
** note: "nrp" = Newton-based Re-Parameterization of where the given points
** fall along the spline, the iterative process inside limnCBFSingle
@@ -554,38 +551,31 @@
cornNMS; /* 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 */
- 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 that 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, /* absolute value of determinant of 2x2 matrix
- to invert can't 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 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. */
/* ----------- internal --------- */
double *uu, /* buffer used for nrp */
*vw, /* weights for endpoint vertex calculation */
@@ -605,13 +595,13 @@
2: DM < dist <= fD
3: fD < dist
where
+ nD = nrpDistScl*distMin,
DM = distMin,
- nD = nrpDistScl*distMin,
fD = nrpPsi*distMin: */
-} limnCBFContext;
+} limnCBFCtx;
/*
-******** limnPoints
+******** limnCBFPoints
**
** 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
@@ -618,9 +608,11 @@
** 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
+** 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.
*/
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 */
@@ -627,7 +619,7 @@
int isLoop; /* points form a loop: logical indices into coord
array are . . . num-2, num-1, 0, 1, . . .
and index 0 is effectively arbitrary */
-} limnPoints;
+} limnCBFPoints;
/* defaultsLimn.c */
LIMN_EXPORT const int limnPresent;
@@ -905,27 +897,33 @@
double maxT);
/* splineFit.c */
-LIMN_EXPORT limnPoints *limnPointsNew(const double *pp, unsigned int nn, int isLoop);
-LIMN_EXPORT limnPoints *limnPointsNix(limnPoints *lpnt);
+LIMN_EXPORT limnCBFPoints *limnCBFPointsNew(const double *pp, unsigned int nn,
+ int isLoop);
+LIMN_EXPORT limnCBFPoints *limnCBFPointsNix(limnCBFPoints *lpnt);
+LIMN_EXPORT int limnCBFPointsCheck(const limnCBFPoints *lpnt);
+LIMN_EXPORT limnCBFCtx *limnCBFCtxNew(unsigned int pointNum, double scale);
+LIMN_EXPORT limnCBFCtx *limnCBFCtxNix(limnCBFCtx *fctx);
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 pNum,
+LIMN_EXPORT void limnCBFPathSample(double *xy, unsigned int pointNum,
const limnCBFPath *path);
-LIMN_EXPORT void limnCBFContextInit(limnCBFContext *fctx, int outputOnly);
-LIMN_EXPORT int limnCBFCheck(const limnCBFContext *fctx, const limnPoints *lpnt);
-LIMN_EXPORT int limnCBFitSingle(double alpha[2], limnCBFContext *fctx,
- const double vv0[2], const double tt1[2],
- const double tt2[2], const double vv3[2],
- const double *xy, unsigned int pNum, int isLoop);
-LIMN_EXPORT int limnCBFMulti(limnCBFPath *path, limnCBFContext *fctx,
- const double vv0[2], const double tt1[2],
- const double tt2[2], const double vv3[2],
- const limnPoints *lpnt, unsigned int loi, unsigned int hii);
+LIMN_EXPORT 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);
+LIMN_EXPORT int limnCBFCtxCheck(const limnCBFCtx *fctx, const limnCBFPoints *lpnt);
+LIMN_EXPORT 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);
+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,
- limnCBFContext *fctx, const limnPoints *lpnt);
-LIMN_EXPORT int limnCBFit(limnCBFPath *path, limnCBFContext *fctx, const double *xy,
- unsigned int pNum, int isLoop);
+ limnCBFCtx *fctx, const limnCBFPoints *lpnt);
+LIMN_EXPORT int limnCBFit(limnCBFPath *path, limnCBFCtx *fctx, const double *xy,
+ unsigned int pointNum, int isLoop);
/* 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-20 21:03:47 UTC (rev 7165)
+++ teem/trunk/src/limn/lpu_cbfit.c 2024-06-21 07:43:16 UTC (rev 7166)
@@ -39,7 +39,7 @@
unsigned int ii, pNum, iterMax;
int loop, petc, verbose, synth, nofit;
char *synthOut;
- limnCBFContext fctx;
+ limnCBFCtx *fctx;
limnCBFPath *path;
hestOptAdd_1_Other(&hopt, "i", "input", &_nin, NULL, "input xy points", nrrdHestNrrd);
@@ -175,15 +175,14 @@
}
path = limnCBFPathNew();
airMopAdd(mop, path, (airMopper)limnCBFPathNix, airMopAlways);
- limnCBFContextInit(&fctx, AIR_FALSE);
- fctx.nrpIterMax = iterMax;
- fctx.nrpDeltaMin = deltaMin;
- fctx.distMin = distMin;
- fctx.nrpDistScl = distScl;
- fctx.verbose = verbose;
- fctx.nrpPsi = psi;
- fctx.cornAngle = cangle;
- fctx.scale = scale;
+ fctx = limnCBFCtxNew(pNum, scale);
+ fctx->nrpIterMax = iterMax;
+ fctx->nrpDeltaMin = deltaMin;
+ fctx->distMin = distMin;
+ fctx->nrpDistScl = distScl;
+ fctx->verbose = verbose;
+ fctx->nrpPsi = psi;
+ fctx->cornAngle = cangle;
time0 = airTime();
if (petc) {
fprintf(stderr, "%s: Press Enter to Continue ... ", me);
@@ -190,7 +189,7 @@
fflush(stderr);
getchar();
}
- if (limnCBFit(path, &fctx, xy, pNum, loop)) {
+ if (limnCBFit(path, fctx, xy, pNum, loop)) {
airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
fprintf(stderr, "%s: trouble:\n%s", me, err);
airMopError(mop);
@@ -198,7 +197,7 @@
}
dtime = (airTime() - time0) * 1000;
printf("%s: time= %g ms;iterDone= %u ;deltaDone=%g, dist=%g (@%u)\n", me, dtime,
- fctx.nrpIterDone, fctx.nrpDeltaDone, fctx.dist, fctx.distIdx);
+ fctx->nrpIterDone, fctx->nrpDeltaDone, fctx->dist, fctx->distIdx);
{
unsigned int si;
printf("%s: path has %u segments:\n", me, path->segNum);
Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c 2024-06-20 21:03:47 UTC (rev 7165)
+++ teem/trunk/src/limn/splineFit.c 2024-06-21 07:43:16 UTC (rev 7166)
@@ -31,18 +31,18 @@
The author's code is here:
http://www.realtimerendering.com/resources/GraphicsGems/gems/FitCurves.c
- The functions below do not use any existing limnSpline structs or functions;
- those were written a long time ago, and reflect GLK's ignorance about
- splines at the time. Hopefully this will be revisited and re-organized in a
- later version of Teem, at which point the code below can be integrated with
- the rest of limn, but this too will benefit from ongoing scrutiny and
- re-writing; ignorance persists.
+ The functions below do not use any other limnSpline structs or functions, since those
+ were written a long time ago and thus they reflect GLK's ignorance about splines at the
+ time. Hopefully this will be revisited and re-organized in a later version of Teem, at
+ which point the code below can be integrated with the rest of limn, but this too will
+ benefit from ongoing scrutiny and re-writing; ignorance persists.
*/
-limnPoints * /* Biff: nope */
-limnPointsNew(const double *pp, uint nn, int isLoop) {
- limnPoints *lpnt;
- lpnt = AIR_CALLOC(1, limnPoints);
+limnCBFPoints * /* Biff: nope */
+limnCBFPointsNew(const double *pp, uint nn, int isLoop) {
+ /* implicit: DIM=2 */
+ limnCBFPoints *lpnt;
+ lpnt = AIR_CALLOC(1, limnCBFPoints);
assert(lpnt);
if (pp) {
/* we are wrapping around a given pre-allocated buffer */
@@ -52,7 +52,7 @@
/* we are allocating our own buffer */
lpnt->pp = NULL;
lpnt->ppOwn = AIR_CALLOC(nn, double);
- assert(lpnt->pp);
+ assert(lpnt->ppOwn);
}
lpnt->num = nn;
lpnt->isLoop = isLoop;
@@ -59,8 +59,8 @@
return lpnt;
}
-limnPoints * /* Biff: nope */
-limnPointsNix(limnPoints *lpnt) {
+limnCBFPoints * /* Biff: nope */
+limnCBFPointsNix(limnCBFPoints *lpnt) {
if (lpnt) {
/* don't touch lpnt->pp */
if (lpnt->ppOwn) free(lpnt->ppOwn);
@@ -69,9 +69,9 @@
return NULL;
}
-static int /* Biff: 1 */
-pointsCheck(const limnPoints *lpnt) {
- static const char me[] = "pointsCheck";
+int /* Biff: 1 */
+limnCBFPointsCheck(const limnCBFPoints *lpnt) {
+ static const char me[] = "limnCBFPointsCheck";
uint pnmin;
int have;
@@ -81,7 +81,7 @@
}
pnmin = lpnt->isLoop ? 3 : 2;
if (!(lpnt->num >= pnmin)) {
- biffAddf(LIMN, "%s: need %u or more points in limnPoints (not %u)%s", me, pnmin,
+ biffAddf(LIMN, "%s: need %u or more points in limnCBFPoints (not %u)%s", me, pnmin,
lpnt->num, lpnt->isLoop ? " for loop" : "");
return 1;
}
@@ -93,11 +93,12 @@
return 0;
}
+/* cheesy macro to access either pp or ppOwn, useful since they differ in const-ness */
#define PP(lpnt) ((lpnt)->pp ? (lpnt)->pp : (lpnt)->ppOwn)
/* number of points between low,high indices loi,hii */
static uint
-pntNum(const limnPoints *lpnt, uint loi, uint hii) {
+pntNum(const limnCBFPoints *lpnt, uint loi, uint hii) {
if (hii < loi) {
assert(lpnt->isLoop);
hii += lpnt->num;
@@ -107,13 +108,155 @@
/* coordinates of point with index loi+ii */
static const double *
-pntCrd(const limnPoints *lpnt, uint loi, uint ii) {
+pntCrd(const limnCBFPoints *lpnt, uint loi, uint ii) {
uint jj = loi + ii;
while (jj >= lpnt->num)
jj -= lpnt->num;
- return PP(lpnt) + 2 * jj;
+ return PP(lpnt) + 2 * jj; /* DIM=2 */
}
+static void
+ctxInit(limnCBFCtx *fctx) {
+ if (!fctx) return;
+ /* defaults for input parameters to various CBF functions */
+ fctx->verbose = 0;
+ fctx->cornNMS = AIR_TRUE;
+ fctx->nrpIterMax = 10;
+ fctx->distMin = 0;
+ fctx->nrpDeltaMax = 3.0;
+ fctx->nrpDistScl = 0.8;
+ fctx->nrpPsi = 6;
+ fctx->nrpDeltaMin = 0.001;
+ fctx->alphaMin = 0.001;
+ fctx->detMin = 0.01;
+ fctx->cornAngle = 100.0; /* degrees */
+ /* internal */
+ fctx->uu = fctx->vw = fctx->tw = NULL;
+ fctx->mine = NULL;
+ fctx->wLen = 0;
+ fctx->lenF2L = AIR_NAN;
+ /* initialize outputs to bogus valus */
+ fctx->nrpIterDone = (uint)(-1);
+ fctx->distIdx = (uint)(-1);
+ fctx->dist = AIR_POS_INF;
+ fctx->nrpDeltaDone = AIR_POS_INF;
+ fctx->alphaDet = 0;
+ fctx->distBig = 0;
+ return;
+}
+
+/*
+** ctxBuffersSet: allocates in fctx:
+** uu, vw, tw
+*/
+static int /* Biff: 1 */
+ctxBuffersSet(limnCBFCtx *fctx, uint pNum, double scl) {
+ static const char me[] = "ctxBuffersSet";
+ double kparm[2], vsum, tsum;
+ /* one: what value in summing kernel weights should count as 1.0. This
+ should probably be a parm in fctx, but not very interesting to
+ change */
+ double one = 0.999;
+ uint ii, len;
+
+ if (fctx->uu) free(fctx->uu);
+ fctx->uu = AIR_CALLOC(pNum * 2, double); /* DIM=2 */
+ if (!fctx->uu) {
+ biffAddf(LIMN, "%s: failed to allocate parameter buffer", me);
+ return 1;
+ }
+ if (0 == scl) {
+ /* will do simplest possible finite differences; we're done */
+ fctx->vw = fctx->tw = NULL;
+ return 0;
+ }
+ /* else need to allocate and set vw and tw buffers */
+ kparm[0] = scl;
+ kparm[1] = 1000; /* effectively, no limit on scale; sanity check comes later */
+ ii = 0;
+ vsum = 0;
+ do {
+ double kw;
+ kw = nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
+ vsum += (!ii ? 1 : 2) * kw;
+ ii++;
+ } while (vsum < one);
+ /* len = intended length of blurring kernel weight vectors */
+ len = ii + 1;
+ if (len > 128) {
+ biffAddf(LIMN,
+ "%s: weight buffer length %u (from scale %g) seems "
+ "unreasonable",
+ me, len, scl);
+ return 1;
+ }
+ if (fctx->vw) free(fctx->vw);
+ if (fctx->tw) free(fctx->tw);
+ fctx->vw = AIR_CALLOC(len, double);
+ fctx->tw = AIR_CALLOC(len, double);
+ if (!(fctx->vw && fctx->tw)) {
+ biffAddf(LIMN, "%s: couldn't allocate weight buffers (len %u)", me, len);
+ return 1;
+ }
+ fctx->wLen = len;
+ /* normalization intent:
+ 1 = sum_i(vw[|i|]) for i=-(len-1)...len-1
+ 1 = sum_i(tw[i]) for i=0...len-1
+ */
+ vsum = tsum = 0;
+ for (ii = 0; ii < len; ii++) {
+ double kw;
+ kw = nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
+ vsum += (!ii ? 1 : 2) * (fctx->vw[ii] = kw);
+ tsum += (fctx->tw[ii] = ii * kw);
+ }
+ for (ii = 0; ii < len; ii++) {
+ fctx->vw[ii] /= vsum;
+ fctx->tw[ii] /= tsum;
+ /* printf("!%s: %u %g %g\n", me, ii, fctx->vw[ii], fctx->tw[ii]); */
+ }
+ return 0;
+}
+
+limnCBFCtx * /* Biff: NULL */
+limnCBFCtxNew(unsigned int pointNum, double scale) {
+ static const char me[] = "limnCBFCtxNew";
+ limnCBFCtx *ret;
+ if (!(pointNum >= 4)) {
+ biffAddf(LIMN, "%s: got tiny #points %u", me, pointNum);
+ return NULL;
+ }
+ if (!(scale >= 0)) {
+ biffAddf(LIMN, "%s: need scale >= 0 (not %g)", me, scale);
+ return NULL;
+ }
+ ret = AIR_CALLOC(1, limnCBFCtx);
+ if (!ret) {
+ biffAddf(LIMN, "%s: allocation failure?", me);
+ return NULL;
+ }
+ ctxInit(ret);
+ if (ctxBuffersSet(ret, pointNum, scale)) {
+ biffAddf(LIMN, "%s: trouble allocating buffers", me);
+ free(ret);
+ return NULL;
+ }
+ ret->scale = scale;
+
+ return ret;
+}
+
+limnCBFCtx * /* Biff: nope */
+limnCBFCtxNix(limnCBFCtx *fctx) {
+ if (fctx) {
+ if (fctx->uu) free(fctx->uu);
+ if (fctx->vw) free(fctx->vw);
+ if (fctx->tw) free(fctx->tw);
+ free(fctx);
+ }
+ return NULL;
+}
+
/* CB0, CB1, CB2, CB3 = degree 3 Bernstein polynomials, for *C*ubic
*B*ezier curves, and their derivatives D0, D1, D2 (not using any
nice recursion properties for evaluation, oh well) */
@@ -123,12 +266,12 @@
#define CB3D0(T) ((T) * (T) * (T))
#define CB0D1(T) (-3 * (1 - (T)) * (1 - (T)))
-#define CB1D1(T) (3 * ((T)-1) * (3 * (T)-1))
+#define CB1D1(T) (3 * ((T) - 1) * (3 * (T) - 1))
#define CB2D1(T) (3 * (T) * (2 - 3 * (T)))
#define CB3D1(T) (3 * (T) * (T))
#define CB0D2(T) (6 * (1 - (T)))
-#define CB1D2(T) (6 * (3 * (T)-2))
+#define CB1D2(T) (6 * (3 * (T) - 2))
#define CB2D2(T) (6 * (1 - 3 * (T)))
#define CB3D2(T) (6 * (T))
@@ -145,6 +288,7 @@
#define CBDI(P, CB, V0, V1, V2, V3, T, W) \
(CB(W, T), \
ELL_2V_SCALE_ADD4(P, (W)[0], (V0), (W)[1], (V1), (W)[2], (V2), (W)[3], (V3)))
+/* _2V_ above: DIM=2 */
#define CBD0(P, V0, V1, V2, V3, T, W) CBDI(P, VCBD0, V0, V1, V2, V3, T, W)
#define CBD1(P, V0, V1, V2, V3, T, W) CBDI(P, VCBD1, V0, V1, V2, V3, T, W)
#define CBD2(P, V0, V1, V2, V3, T, W) CBDI(P, VCBD2, V0, V1, V2, V3, T, W)
@@ -158,7 +302,7 @@
limnCBFSegEval(double *vv, const limnCBFSeg *seg, double tt) {
double ww[4];
const double *xy = seg->xy;
- CBD0(vv, xy + 0, xy + 2, xy + 4, xy + 6, tt, ww);
+ 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",
@@ -185,7 +329,7 @@
const limnCBFSeg *seg = path->seg + segi;
double tmpf = AIR_AFFINE(0, ii, pNum - 1, 0, sNum);
double tt = tmpf - segi;
- limnCBFSegEval(xy + 2 * ii, seg, tt);
+ limnCBFSegEval(xy + 2 * ii, seg, tt); /* DIM=2 */
/*
fprintf(stderr, "!%s: %u -> %u (%g) %g -> (%g,%g)\n",
"limnCBFPathSample", ii, segi, tmpf, tt,
@@ -196,13 +340,13 @@
}
/*
-** Find endpoint vertex vv and tangent tt (constraints for spline fitting)
-** from the given points lpnt at coord index ii within index range [loi,hoi]
-** (e.g. ii=1 means looking at lpnt coord index loi+1). The tangent direction
+** limnCBFFindVT: Find endpoint vertex vv and tangent tt (constraints for spline fitting)
+** from the given points lpnt at offset index ofi within index range [loi,hii]
+** (e.g. ofi=1 means looking at lpnt coord index loi+1). The tangent direction
** dir controls which points are looked at:
-** >0: considering only ii and higher-index vertices,
-** 0: for tangent centered at ii, using lower- and higher-index vertices
-** <0: considering only ii and lower-index vertices
+** >0: considering only ofi and higher-index vertices,
+** 0: for tangent centered at ofi, using lower- and higher-index vertices
+** <0: considering only ofi and lower-index vertices
** For >0 and 0: the tangent points towards the positions of higher-
** index vertices. For <0, it points the other way.
** The only point indices accessed will be in [loi,hii]; this is what
@@ -209,38 +353,79 @@
** enforces the possible corner-ness of those indices (which prevents
** vertices past corners influencing how vv or tt are found)
*/
-static void
-findVT(double vv[2], double tt[2], const limnCBFContext *fctx, const limnPoints *lpnt,
- uint loi, uint hii, uint ii, int dir) {
- /* static const char me[] = "findVT"; */
+int /* Biff: 1 */
+limnCBFFindVT(double vv[2], double tt[2], const limnCBFCtx *fctx,
+ const limnCBFPoints *lpnt, uint _loi, uint _hii, uint _ofi, int dir) {
+ static const char me[] = "limnCBFFindVT";
double len;
- uint pNum, /* total number of points in lpnts */
- sgsz; /* segment size: number of points in [loi,hii] */
+ /* we use here (signed) int for things that might seem better as uint, but that's
+ because of the frequent need to handle how indices loop around in loops */
+ int loi, hii, ofi, /* int versions of given args */
+ pNum, /* total number of points in lpnts */
+ sgsz, /* "segment" size: number of points in [loi,hii]
+ (but not in sense of the limnCBFSeg specifically) */
+ icent, iplus, imnus; /* icent is the actual data index corresponding to _loi + _ofi;
+ it is used for both the scale==0 and scale>0 cases;
+ iplus and imnus are only needed with scale==0, but too annoying
+ to break those out into that specific branch */
+ if (!(/* vv can be NULL */ tt && fctx && lpnt)) {
+ biffAddf(LIMN, "%s: got NULL pointer", me);
+ return 1;
+ }
+ if (!(_loi < lpnt->num && _hii < lpnt->num)) {
+ biffAddf(LIMN, "%s: _loi %u or _hii %u too high for %u points", me, _loi, _hii,
+ lpnt->num);
+ }
+ /* now both _loi and _hii are valid indices in [0,..,lpnt->num-1] */
+ if (_loi == _hii) {
+ biffAddf(LIMN, "%s: got _loi==_hii %u", me, _loi);
+ return 1;
+ }
+ if (_hii < _loi && !lpnt->isLoop) {
+ biffAddf(LIMN, "%s: _hii %u < _loi %u sensible only in point loop", me, _hii, _loi);
+ return 1;
+ }
+ /* we proceed with either _loi < _hii or (_hii < _loi and) lpnt->isLoop
+ i.e. now _hii < _loi (and soon hii < loi) implies lpnt->isLoop */
dir = airSgn(dir);
- pNum = lpnt->num;
+ pNum = AIR_INT(lpnt->num);
+ loi = AIR_INT(_loi);
+ hii = AIR_INT(_hii);
+ sgsz = (hii < loi ? pNum : 0) + hii - loi + 1;
+ if (sgsz <= 1) {
+ biffAddf(LIMN, "%s: sgsz %d <= 1", me, sgsz);
+ return 1;
+ }
+ /* now sgsz >= 2 */
+ if (!(_ofi < AIR_UINT(sgsz))) {
+ biffAddf(LIMN, "%s: _ofi %u too high for segment size %d", me, _ofi, sgsz);
+ return 1;
+ }
+ /* now _ofi is a valid index in [0,..,sgsz-1] */
+ ofi = AIR_INT(_ofi);
+ icent = loi + ofi;
+ iplus = icent + 1;
+ imnus = icent - 1;
if (lpnt->isLoop) {
- sgsz = (hii < loi ? pNum : 0) + hii - loi + 1;
+ icent = AIR_MOD(icent, pNum);
+ iplus = AIR_MOD(iplus, pNum);
+ imnus = AIR_MOD(imnus, pNum);
} else {
- sgsz = hii - loi + 1;
+ icent = AIR_CLAMP(loi, icent, hii);
+ iplus = AIR_CLAMP(loi, iplus, hii);
+ imnus = AIR_CLAMP(loi, imnus, hii);
}
if (0 == fctx->scale) {
- uint mi, pi, iplus, imnus;
const double *xy, *xyP, *xyM;
- if (lpnt->isLoop) {
- iplus = (loi + ii + 1) % pNum;
- imnus = (uint)AIR_MOD((int)(loi + ii) - 1, (int)pNum);
- } else {
- /* regardless of lpnt->isLoop, we only look in [loi,hii] */
- iplus = loi + AIR_MIN(ii + 1, sgsz - 1);
- imnus = loi + AIR_MAX(1, ii) - 1;
+ int mi, ci, pi;
+ if (vv) {
+ ELL_2V_COPY(vv, PP(lpnt) + 2 * icent); /* DIM=2 */
}
- xy = pntCrd(lpnt, loi, ii);
- if (vv) ELL_2V_COPY(vv, xy);
switch (dir) {
case 1:
pi = iplus;
- mi = ii;
+ mi = icent;
break;
case 0:
pi = iplus;
@@ -248,120 +433,131 @@
break;
case -1:
/* mi and pi switched to point other way */
- mi = ii;
+ mi = icent;
pi = imnus;
break;
}
- /* if (with !isLoop) ii=0 and dir=-1, or, ii=pNum-1 and dir=+1
- ==> mi=pi ==> tt will be (nan,nan), which is appropriate */
- xyP = pntCrd(lpnt, loi, pi);
- xyM = pntCrd(lpnt, loi, mi);
+ if (pi == mi) {
+ biffAddf(LIMN, "%s: imnus,icent,iplus=%d,%d,%d --dir=%d--> mi == pi %d", me, imnus,
+ icent, iplus, dir, mi);
+ return 1;
+ }
+ /* printf("!%s: iplus=%u imnus=%u xy=%g %g\n", me, iplus, imnus, xy[0], xy[1]); */
+ xyP = PP(lpnt) + 2 * pi; /* DIM=2 and following */
+ xyM = PP(lpnt) + 2 * mi;
ELL_2V_SUB(tt, xyP, xyM);
ELL_2V_NORM(tt, tt, len);
+ /* printf("!%s: xyP=%g %g xyM=%g %g len=%g tt=%g %g\n", me, xyP[0], xyP[1],
+ xyM[0], xyM[1], len, tt[0], tt[1]); */
} else {
-#if 0
/* using scale>0 for endpoint and tangent estimation */
+ /* regardless of dir, we compute average positions for points
+ centered around loi + ofi (posC), and for lower,higher indices (posM,posP) */
+ double posP[2] = {0, 0}, posC[2] = {0, 0}, posM[2] = {0, 0};
const double *vw = fctx->vw;
const double *tw = fctx->tw;
+ int smax = (int)fctx->wLen - 1, /* bound of loop index */
+ ci; /* loop through [-smax,smax] */
+ if (!(vw && tw)) {
+ biffAddf(LIMN, "%s: fctx internal buffers vw and tw not both allocated", me);
+ return 1;
+ }
/* various signed indices */
- int sii=(int)ii, /* we compute around vertex ii */
- smax=(int)fctx->wLen - 1, /* bound of loop index */
- sj; /* loop through [-smax,smax] */
- if (vv) ELL_2V_SET(vv, 0, 0);
- ELL_2V_SET(tt, 0, 0);
- /* printf("!%s: ii = %u, dir=%d\n", me, ii, dir); */
- /* j indices are for the local looping */
- for (sj=-smax; sj<=smax; sj++) {
- uint xj, /* eventual index into data */
- asj = (uint)AIR_ABS(sj); /* index into vw, tw */
- int sgn=1,
- sxj = sii + sj; /* signed (tmp) j idx into data */
- double ttw;
- /* printf("!%s[sj=%d,asj=%u]: sxj0 = %d\n", me, sj, asj, sxj); */
- switch (dir) {
- case 1:
- sxj = AIR_MAX(sxj, sii);
- break;
- case -1:
- sgn=-1;
- sxj = AIR_MIN(sxj, sii);
- break;
- }
- /* sxj = sii+sj, but capped at sii according to dir */
- /* printf("!%s[sj=%d]: dir=%d -> sxj1 = %d\n", me, sj, dir, sxj); */
+ /* printf("!%s: ofi = %d, dir=%d\n", me, ofi, dir); */
+ for (ci = -smax; ci <= smax; ci++) {
+ uint wi = abs(ci); /* weight index into vw, tw */
+ int di = loi + ofi + ci; /* signed index into data */
+ const double *xy;
if (lpnt->isLoop) {
- sxj = AIR_MOD(sxj, (int)pNum);
+ di = AIR_MOD(di, pNum);
} else {
- sxj = AIR_CLAMP(0, sxj, (int)pNum-1);
+ di = AIR_CLAMP(loi, di, hii);
}
- xj = (uint)sxj;
- /* printf("!%s[sj=%d]: isLoop=%d -> sxj2 = %d -> xj = %u\n", me, sj, lpnt->isLoop, sxj, xj); */
- if (vv) ELL_2V_SCALE_INCR(vv, vw[asj], xy + 2*xj);
- /* printf("!%s[sj=%d]: vv += %g*(%g,%g) -> (%g,%g)\n", me, sj, vw[asj], (xy + 2*xj)[0], (xy + 2*xj)[1], vv[0], vv[1]); */
- ttw = sgn*airSgn(sj)*tw[asj];
- /* printf("!%s[sj=%d]: %d * %d * %g = %g\n", me, sj, sgn, airSgn(sj), tw[asj], ttw); */
- ELL_2V_SCALE_INCR(tt, ttw, xy + 2*xj);
- /* printf("!%s[sj=%d]: tt += %g*(%g,%g) -> (%g,%g)\n", me, sj, ttw, (xy + 2*xj)[0], (xy + 2*xj)[1], tt[0], tt[1]); */
+ xy = PP(lpnt) + 2 * di;
+ ELL_2V_SCALE_INCR(posC, vw[wi], xy);
+ /* printf("!%s: (ci=%d/wi=%u) posC += %g*(%g %g) --> %g %g\n", me, ci, wi, vw[wi],
+ xy[0], xy[1], posC[0], posC[1]); */
+ if (ci < 0) {
+ ELL_2V_SCALE_INCR(posM, tw[wi], xy);
+ /* printf("!%s: (ci=%d/wi=%u) posM += %g*(%g %g) --> %g %g\n", me, ci, wi,
+ tw[wi], xy[0], xy[1], posM[0], posM[1]); */
+ }
+ if (ci > 0) {
+ ELL_2V_SCALE_INCR(posP, tw[wi], xy);
+ /* printf("!%s: (ci=%d/wi=%u) posP += %g*(%g %g) --> %g %g\n", me, ci, wi,
+ tw[wi], xy[0], xy[1], posP[0], posP[1]); */
+ }
}
+ switch (dir) {
+ case 1:
+ ELL_2V_SUB(tt, posP, posC);
+ break;
+ case 0:
+ ELL_2V_SUB(tt, posP, posM);
+ break;
+ case -1:
+ ELL_2V_SUB(tt, posM, posC);
+ break;
+ }
ELL_2V_NORM(tt, tt, len);
- /* fix the boundary conditions as a post-process */
- if ( 0==ii && -1==dir) ELL_2V_SET(tt, AIR_NAN, AIR_NAN);
- if (pNum-1==ii && +1==dir) ELL_2V_SET(tt, AIR_NAN, AIR_NAN);
if (vv) {
+ ELL_2V_COPY(vv, posC);
/* some post-proceessing of computed spline endpoint */
- double off[2], pp[2], operp;
- ELL_2V_SET(pp, tt[1], -tt[0]); /* pp is perpendicular to tt */
- ELL_2V_SUB(off, vv, xy + 2*ii);
+ double off[2], pp[2], operp, okoff;
+ const double *xy;
+ /* DIM=2 */
+ xy = PP(lpnt) + 2 * icent; /* center vertex in given data */
+ ELL_2V_SET(pp, tt[1], -tt[0]); /* pp is perpendicular to computed tt */
+ ELL_2V_SUB(off, vv, xy); /* off = vv - xy, from given to computed */
operp = ELL_2V_DOT(off, pp);
/* limit distance from chosen (x,y) datapoint to spline endpoint to be
(HEY harcoded) 95% of fctx->distMin. Being allowed to be further away
can cause annoyances */
- operp = AIR_MIN(0.95*fctx->distMin, operp);
+ okoff = 0.95 * fctx->distMin;
+ operp = AIR_CLAMP(-okoff, operp, okoff);
/* constrain difference between chosen (x,y) datapoint and spline
endpoint to be perpendicular to estimated tangent */
- ELL_2V_SCALE_ADD2(vv, 1, xy + 2*ii, operp, pp);
+ ELL_2V_SCALE_ADD2(vv, 1, xy, operp, pp);
}
-#endif
}
- return;
+ return 0;
}
-static int /* Biff: 1 */
-setVTTV(int *given, double vv0[2], double tt1[2], double tt2[2], double vv3[2],
+static int /* Biff: 1 */
+setVTTV(int *given, /* */
+ double vv0[2], double tt1[2], double tt2[2], double vv3[2], /* */
const double _vv0[2], const double _tt1[2], const double _tt2[2],
- const double _vv3[2], const limnCBFContext *fctx, const limnPoints *lpnt,
- uint loi, uint hii) {
+ const double _vv3[2], /* */
+ const limnCBFCtx *fctx, const limnCBFPoints *lpnt, uint loi, uint hii) {
static const char me[] = "setVTTV";
- /* either all the _vv0, _tt1, _tt2, _vv3 can be NULL, or none */
- if (!(_vv0 && _tt1 && _tt2 && _vv3)) {
+ /* either: all the _vv0, _tt1, _tt2, _vv3 can be NULL, or none are NULL */
+ if (_vv0 && _tt1 && _tt2 && _vv3) {
+ /* copy the given endpoint geometry */
+ ELL_2V_COPY(vv0, _vv0);
+ ELL_2V_COPY(tt1, _tt1);
+ ELL_2V_COPY(tt2, _tt2);
+ ELL_2V_COPY(vv3, _vv3);
+ if (given) {
+ *given = AIR_TRUE;
+ }
+ } else {
+ /* not given v, t, t, v values, so we compute them */
if (_vv0 || _tt1 || _tt2 || _vv3) {
- biffAddf(LIMN,
- "%s: either all or none of vv0,tt1,tt2,vv3 "
- "should be NULL",
- me);
+ biffAddf(LIMN, "%s: either all or none of _vv0,_tt1,_tt2,_vv3 should be NULL", me);
return 1;
}
if (lpnt->isLoop) {
- findVT(vv0, tt1, fctx, lpnt, loi, hii, loi, 0);
+ limnCBFFindVT(vv0, tt1, fctx, lpnt, loi, hii, 0, 0);
ELL_2V_COPY(vv3, vv0);
ELL_2V_SCALE(tt2, -1, tt1);
} else {
- findVT(vv0, tt1, fctx, lpnt, loi, hii, loi, +1);
- findVT(vv0, tt1, fctx, lpnt, loi, hii, hii, -1);
+ limnCBFFindVT(vv0, tt1, fctx, lpnt, loi, hii, 0, +1);
+ limnCBFFindVT(vv3, tt2, fctx, lpnt, loi, hii, hii - loi, -1);
}
if (given) {
*given = AIR_FALSE;
}
- } else {
- /* copy the given endpoint geometry */
- ELL_2V_COPY(vv0, _vv0);
- ELL_2V_COPY(tt1, _tt1);
- ELL_2V_COPY(tt2, _tt2);
- ELL_2V_COPY(vv3, _vv3);
- if (given) {
- *given = AIR_TRUE;
- }
}
return 0;
}
@@ -388,9 +584,9 @@
** (possibly unstable?)
*/
static void
-findalpha(double alpha[2], limnCBFContext *fctx, /* must be non-NULL */
+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 limnPoints *lpnt, uint loi, uint hii) {
+ const double vv3[2], const limnCBFPoints *lpnt, uint loi, uint hii) {
static const char me[] = "findalpha";
uint ii, pNum;
double det;
@@ -397,7 +593,7 @@
pNum = pntNum(lpnt, loi, hii);
if (pNum > 2) {
- double xx[2], m11, m12, m22, MM[4], MI[4];
+ double xx[2], m11, m12, m22, MM[4], MI[4]; /* DIM=2 throughout this */
const double *uu = fctx->uu;
xx[0] = xx[1] = m11 = m12 = m22 = 0;
for (ii = 0; ii < pNum; ii++) {
@@ -457,9 +653,9 @@
** to evaluate the spline in order to match the given points xy
*/
static double
-reparm(const limnCBFContext *fctx, /* must be non-NULL */
+reparm(const limnCBFCtx *fctx, /* must be non-NULL */
const double alpha[2], const double vv0[2], const double tt1[2],
- const double tt2[2], const double vv3[2], const limnPoints *lpnt, uint loi,
+ const double tt2[2], const double vv3[2], const limnCBFPoints *lpnt, uint loi,
uint hii) {
static const char me[] = "reparm";
uint ii, pNum;
@@ -506,12 +702,12 @@
return delta;
}
-/* sets fctx->dist to max distance to spline, at point fctx->distIdx,
- and then sets fctx->distBig accordingly */
+/* (assuming current parameterization in fctx->uu) sets fctx->dist to max distance to
+ spline, at point fctx->distIdx, and then sets fctx->distBig accordingly */
static void
-finddist(limnCBFContext *fctx, const double alpha[2], const double vv0[2],
+finddist(limnCBFCtx *fctx, const double alpha[2], const double vv0[2],
const double tt1[2], const double tt2[2], const double vv3[2],
- const limnPoints *lpnt, uint loi, uint hii) {
+ const limnCBFPoints *lpnt, uint loi, uint hii) {
uint ii, distI, pNum;
double vv1[2], vv2[2], dist;
const double *uu = fctx->uu;
@@ -518,14 +714,13 @@
pNum = pntNum(lpnt, loi, hii);
assert(pNum >= 3);
- ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1);
+ ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1); /* DIM=2 everywhere here */
ELL_2V_SCALE_ADD2(vv2, 1, vv3, alpha[1], tt2);
dist = AIR_NAN;
- /* NOTE that the first and last points are actually not part of the max
- distance calculation, which motivates ensuring that the endpoints
- generated by findVT are actually sufficiently close to the first and last
- points (or else the fit spline won't meet the expected accuracy
- threshold) */
+ /* NOTE that the first and last points are actually not part of the max distance
+ calculation, which motivates ensuring that the endpoints generated by limnCBFFindVT
+ are actually sufficiently close to the first and last points (or else the fit spline
+ won't meet the expected accuracy threshold) */
for (ii = 1; ii < pNum - 1; ii++) {
double len, Q[2], df[2], ww[4];
const double *xy;
@@ -547,53 +742,20 @@
return;
}
-void
-limnCBFContextInit(limnCBFContext *fctx, int outputOnly) {
- if (!fctx) return;
- if (!outputOnly) {
- /* defaults for input parameters to various CBF functions */
- fctx->verbose = 0;
- fctx->cornNMS = AIR_TRUE;
- fctx->nrpIterMax = 10;
- fctx->scale = 0;
- fctx->distMin = 0;
- fctx->nrpDeltaMax = 3.0;
- fctx->nrpDistScl = 0.8;
- fctx->nrpPsi = 6;
- fctx->nrpDeltaMin = 0.001;
- fctx->alphaMin = 0.001;
- fctx->detMin = 0.01;
- fctx->cornAngle = 100.0; /* degrees */
- }
- /* internal */
- fctx->uu = fctx->vw = fctx->tw = NULL;
- fctx->mine = NULL;
- fctx->wLen = 0;
- fctx->lenF2L = AIR_NAN;
- /* initialize outputs to bogus valus */
- fctx->nrpIterDone = (uint)(-1);
- fctx->distIdx = (uint)(-1);
- fctx->dist = AIR_POS_INF;
- fctx->nrpDeltaDone = AIR_POS_INF;
- fctx->alphaDet = 0;
- fctx->distBig = 0;
- return;
-}
-
/*
-******** limnCBFCheck
+******** limnCBFCtxCheck
**
** checks the things that are going to be passed around a lot
*/
int /* Biff: 1 */
-limnCBFCheck(const limnCBFContext *fctx, const limnPoints *lpnt) {
- static const char me[] = "limnCBFCheck";
+limnCBFCtxCheck(const limnCBFCtx *fctx, const limnCBFPoints *lpnt) {
+ static const char me[] = "limnCBFCtxCheck";
if (!(fctx && lpnt)) {
biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
- if (pointsCheck(lpnt)) {
+ if (limnCBFPointsCheck(lpnt)) {
biffAddf(LIMN, "%s: problem with points", me);
return 1;
}
@@ -631,7 +793,7 @@
** fitSingle: fits a single cubic Bezier spline, w/out error checking,
** limnCBFSingle is a wrapper around this.
**
-** The given points coordinates are in limnPoints lpnt, between low/high
+** 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 initial endpoint vv0, initial tangent tt1, final endpoint vv3
** and final tangent tt2 (pointing backwards), this function finds alpha such
@@ -652,9 +814,9 @@
** fctx.
*/
static void
-fitSingle(double alpha[2], limnCBFContext *fctx, const double vv0[2],
- const double tt1[2], const double tt2[2], const double vv3[2],
- const limnPoints *lpnt, uint loi, uint hii) {
+fitSingle(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[] = "fitSingle";
uint iter, pNum;
const double *xy;
@@ -754,147 +916,50 @@
}
/*
-** buffersNew: allocates in fctx:
-** uu, vw, tw
-*/
-static int /* Biff: 1 */
-buffersNew(limnCBFContext *fctx, uint pNum) {
- static const char me[] = "buffersNew";
- double kw, kparm[2], vsum, tsum, scl = fctx->scale;
- /* one: what value in summing kernel weights should count as 1.0. This
- should probably be a parm in fctx, but not very interesting to
- change */
- double one = 0.999;
- uint ii, len;
-
- fctx->uu = AIR_CALLOC(pNum * 2, double);
- if (!fctx->uu) {
- biffAddf(LIMN, "%s: failed to allocate parameter buffer", me);
- return 1;
- }
- if (0 == scl) {
- /* will do simplest possible finite differences; we're done */
- fctx->vw = fctx->tw = NULL;
- return 0;
- }
- /* else need to allocate and set vw and tw buffers */
- kparm[0] = scl;
- kparm[1] = 1000000; /* effectively no cut-off */
- ii = 0;
- vsum = 0;
- do {
- kw = nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
- vsum += (!ii ? 1 : 2) * kw;
- ii++;
- } while (vsum < one);
- /* intended length of vectors */
- len = ii + 1;
- if (len > 128) {
- biffAddf(LIMN,
- "%s: weight buffer length %u (from scale %g) seems "
- "unreasonable",
- me, len, scl);
- return 1;
- }
- fctx->vw = AIR_CALLOC(len, double);
- fctx->tw = AIR_CALLOC(len, double);
- if (!(fctx->vw && fctx->tw)) {
- biffAddf(LIMN, "%s: couldn't allocate weight buffers (len %u)", me, len);
- return 1;
- }
- fctx->wLen = len;
- /* normalization intent:
- 1 = sum_i(vw[|i|]) for i=-(len-1)...len-1
- 1 = sum_i(tw[i]) for i=0...len-1
- */
- vsum = tsum = 0;
- for (ii = 0; ii < len; ii++) {
- kw = nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
- vsum += (!ii ? 1 : 2) * (fctx->vw[ii] = kw);
- tsum += (fctx->tw[ii] = ii * kw);
- }
- for (ii = 0; ii < len; ii++) {
- fctx->vw[ii] /= vsum;
- fctx->tw[ii] /= tsum;
- /* printf("!%s: %u %g %g\n", me, ii, fctx->vw[ii], fctx->tw[ii]); */
- }
- return 0;
-}
-
-/* returning a pointer so compatible with an airMopper */
-static void *
-buffersNix(limnCBFContext *fctx) {
- fctx->uu = (double *)airFree(fctx->uu);
- fctx->vw = (double *)airFree(fctx->vw);
- fctx->tw = (double *)airFree(fctx->tw);
- return NULL;
-}
-
-/* macros to manage the heap-allocated things inside limnCBFContext; working
- with the idea that each caller passes an OWN variable on their stack, so
- the NIX macro only frees thing when the address of OWN matches that passed
- to the NEW. Nothing else in Teem uses this strategy; it may be exploring
- the clever/stupid boundary that David and Nigel famously identified. */
-#define BUFFERS_NEW(FCTX, NN, OWN) \
- if (!(FCTX)->uu) { \
- if (buffersNew((FCTX), (NN))) { \
- biffAddf(LIMN, "%s: failed to allocate buffers", me); \
- return 1; \
- } \
- (FCTX)->mine = &(OWN); \
- }
-
-#define BUFFERS_NIX(FCTX, OWN) \
- if ((FCTX)->mine == &(OWN)) { \
- buffersNix(FCTX); \
- (FCTX)->mine = NULL; \
- }
-
-/*
******** limnCBFitSingle
**
-** builds a limnPoints around given xy, determines spline
+** builds a limnCBFPoints around given xy, determines spline
** constraints if necessary, and calls fitSingle
*/
int /* Biff: 1 */
-limnCBFitSingle(double alpha[2], limnCBFContext *_fctx, const double _vv0[2],
+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, uint pNum, int isLoop) {
static const char me[] = "limnCBFitSingle";
- double own, vv0[2], tt1[2], tt2[2], vv3[2];
+ double vv0[2], tt1[2], tt2[2], vv3[2];
uint loi, hii;
- limnCBFContext *fctx, myfctx;
- limnPoints *lpnt;
+ limnCBFCtx *fctx;
+ limnCBFPoints *lpnt;
+ airArray *mop;
if (!(alpha && xy && pNum)) {
biffAddf(LIMN, "%s: got NULL pointer or 0 points", me);
return 1;
}
- lpnt = limnPointsNew(xy, pNum, isLoop);
+ mop = airMopNew();
+ lpnt = limnCBFPointsNew(xy, pNum, isLoop);
+ airMopAdd(mop, lpnt, (airMopper)limnCBFPointsNix, airMopAlways);
loi = 0;
hii = pNum - 1;
if (_fctx) {
fctx = _fctx; /* caller has supplied info */
- if (limnCBFCheck(fctx, lpnt)) {
- biffAddf(LIMN, "%s: problem with fctx", me);
- limnPointsNix(lpnt);
+ if (limnCBFCtxCheck(fctx, lpnt)) {
+ biffAddf(LIMN, "%s: problem with fctx given points lpnt", me);
+ airMopError(mop);
return 1;
}
- limnCBFContextInit(fctx, AIR_TRUE /* outputOnly */);
} else {
- fctx = &myfctx; /* caller supplied nothing: use defaults */
- limnCBFContextInit(fctx, AIR_FALSE /* outputOnly */);
+ fctx = limnCBFCtxNew(pNum, 0.0 /* scale */);
+ airMopAdd(mop, fctx, (airMopper)limnCBFCtxNix, airMopAlways);
}
- BUFFERS_NEW(fctx, pNum, own);
if (setVTTV(NULL, vv0, tt1, tt2, vv3, _vv0, _tt1, _tt2, _vv3, fctx, lpnt, loi, hii)) {
biffAddf(LIMN, "%s: trouble", me);
- limnPointsNix(lpnt);
+ airMopError(mop);
return 1;
}
fitSingle(alpha, fctx, vv0, tt1, tt2, vv3, lpnt, loi, hii);
- BUFFERS_NIX(fctx, own);
- limnPointsNix(lpnt);
+ airMopOkay(mop);
return 0;
}
@@ -901,7 +966,7 @@
static void
segInit(void *_seg) {
limnCBFSeg *seg = (limnCBFSeg *)_seg;
- ELL_2V_NAN_SET(seg->xy + 0);
+ ELL_2V_NAN_SET(seg->xy + 0); /* DIM=2 */
ELL_2V_NAN_SET(seg->xy + 2);
ELL_2V_NAN_SET(seg->xy + 4);
ELL_2V_NAN_SET(seg->xy + 6);
@@ -948,9 +1013,9 @@
** left and right sides around points with the highest error from fitSingle.
*/
int /* Biff: 1 */
-limnCBFMulti(limnCBFPath *path, limnCBFContext *fctx, const double _vv0[2],
+limnCBFMulti(limnCBFPath *path, limnCBFCtx *fctx, const double _vv0[2],
const double _tt1[2], const double _tt2[2], const double _vv3[2],
- const limnPoints *lpnt, uint loi, uint hii) {
+ const limnCBFPoints *lpnt, uint loi, uint hii) {
static const char me[] = "limnCBFMulti";
double vv0[2], tt1[2], tt2[2], vv3[2], alpha[2];
/* &ownbuff determines who frees buffers inside fctx, since each
@@ -960,7 +1025,7 @@
uint pNum;
/* need non-NULL fctx in order to know fctx->distMin */
- if (limnCBFCheck(fctx, lpnt)) {
+ if (limnCBFCtxCheck(fctx, lpnt)) {
biffAddf(LIMN, "%s: got bad args", me);
return 1;
}
@@ -977,7 +1042,6 @@
return 1;
}
pNum = pntNum(lpnt, loi, hii);
- BUFFERS_NEW(fctx, pNum, ownbuff);
if (setVTTV(&geomGiven, vv0, tt1, tt2, vv3, _vv0, _tt1, _tt2, _vv3, fctx, lpnt, loi,
hii)) {
biffAddf(LIMN, "%s: trouble", me);
@@ -1013,14 +1077,14 @@
uint mi = fctx->distIdx;
double ttL[2], mid[2], ttR[2];
limnCBFPath *prth = limnCBFPathNew(); /* right path */
- limnCBFContext fctxL, fctxR;
- memcpy(&fctxL, fctx, sizeof(limnCBFContext));
- memcpy(&fctxR, fctx, sizeof(limnCBFContext));
+ limnCBFCtx fctxL, fctxR;
+ 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, fctx->dist,
fctx->distBig, mi);
}
- findVT(mid, ttR, fctx, lpnt, loi, hii, mi, 0);
+ limnCBFFindVT(mid, ttR, fctx, lpnt, loi, hii, mi, 0);
ELL_2V_SCALE(ttL, -1, ttR);
/* on recursion, can't be a loop, so isLoop is AIR_FALSE */
if (limnCBFMulti(path, &fctxL, vv0, tt1, ttL, mid, lpnt, loi, mi)
@@ -1045,13 +1109,12 @@
fctx->alphaDet = AIR_MIN(fctxL.alphaDet, fctxR.alphaDet);
}
- BUFFERS_NIX(fctx, ownbuff);
return 0;
}
int /* Biff: 1 */
-limnCBFCorners(uint **cornIdx, uint *cornNum, limnCBFContext *fctx,
- const limnPoints *lpnt) {
+limnCBFCorners(uint **cornIdx, uint *cornNum, limnCBFCtx *fctx,
+ const limnCBFPoints *lpnt) {
static const char me[] = "limnCBFCorners";
airArray *mop, *cornArr;
double ownbuff, *angle;
@@ -1062,7 +1125,7 @@
biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
- if (limnCBFCheck(fctx, lpnt)) {
+ if (limnCBFCtxCheck(fctx, lpnt)) {
biffAddf(LIMN, "%s: got bad args", me);
return 1;
}
@@ -1086,11 +1149,10 @@
cornArr = airArrayNew(AIR_CAST(void **, cornIdx), cornNum, sizeof(uint), 32);
/* free with Nix, not Nuke, because we are managing the given pointers */
airMopAdd(mop, cornArr, (airMopper)airArrayNix, airMopAlways);
- BUFFERS_NEW(fctx, pNum, ownbuff);
for (ii = 0; ii < pNum; ii++) {
double LT[2], RT[2];
- findVT(NULL, LT, fctx, lpnt, loi, hii, ii, -1);
- findVT(NULL, RT, fctx, lpnt, loi, hii, ii, +1);
+ limnCBFFindVT(NULL, LT, fctx, lpnt, loi, hii, ii, -1);
+ limnCBFFindVT(NULL, RT, fctx, lpnt, loi, hii, ii, +1);
angle[ii] = 180 * ell_2v_angle_d(LT, RT) / AIR_PI;
corn[ii] = (angle[ii] < fctx->cornAngle);
}
@@ -1109,7 +1171,6 @@
ci = airArrayLenIncr(cornArr, 1);
(*cornIdx)[ci] = ii;
}
- BUFFERS_NIX(fctx, ownbuff);
airMopOkay(mop);
return 0;
}
@@ -1120,12 +1181,11 @@
** top-level function for fitting cubic beziers to given points
*/
int /* Biff: 1 */
-limnCBFit(limnCBFPath *path, limnCBFContext *fctx, const double *xy, uint pNum,
- int isLoop) {
+limnCBFit(limnCBFPath *path, limnCBFCtx *fctx, const double *xy, uint pNum, int isLoop) {
static const char me[] = "limnCBFit";
uint *cornIdx = NULL, cornNum = 0, cii, loi, hii;
limnCBFPath *rpth;
- limnPoints *lpnt;
+ limnCBFPoints *lpnt;
int ret;
airArray *mop;
@@ -1133,25 +1193,22 @@
biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
- lpnt = limnPointsNew(xy, pNum, isLoop);
+ lpnt = limnCBFPointsNew(xy, pNum, isLoop);
mop = airMopNew();
- airMopAdd(mop, lpnt, (airMopper)limnPointsNix, airMopAlways);
- if (limnCBFCheck(fctx, lpnt)) {
+ airMopAdd(mop, lpnt, (airMopper)limnCBFPointsNix, airMopAlways);
+ if (limnCBFCtxCheck(fctx, lpnt)) {
biffAddf(LIMN, "%s: got bad args", me);
airMopError(mop);
return 1;
}
if (fctx->uu) {
- biffAddf(LIMN, "%s: not expecting limnCBFContext buffers allocated", me);
+ biffAddf(LIMN,
+ "%s: not expecting limnCBFCtx internal buffers to already "
+ "be allocated",
+ me);
airMopError(mop);
return 1;
}
- if (buffersNew(fctx, pNum)) {
- biffAddf(LIMN, "%s: failed to allocate buffers", me);
- airMopError(mop);
- return 1;
- }
- airMopAdd(mop, fctx, (airMopper)buffersNix, airMopAlways);
if (limnCBFCorners(&cornIdx, &cornNum, fctx, lpnt)) {
biffAddf(LIMN, "%s: trouble finding corners", me);
@@ -1204,12 +1261,6 @@
/*
TODO:
-rewrite things to use limnPointList, with first and last indices,
-naturally handling the case that last < first, with isLoop
-and new logic: isLoop does NOT depend on duplicate 1st,last coords
-and subtlty that if (with isLoop) hii = (loi-1 % #points) then using
-all points, with no notion of corner possible
-
testing corners: corners at start==stop of isLoop
corners not at start or stop of isLoop: do spline wrap around from last to first index?
@@ -1221,4 +1272,9 @@
(may want to pay in time for more economical representation)
valgrind everything
+
+remove big debugging comment blocks
+
+(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
*/
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
