Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[7168] teem/trunk/src/limn
|
From: <kin...@us...> - 2024-06-27 02:28:36
|
Revision: 7168
http://sourceforge.net/p/teem/code/7168
Author: kindlmann
Date: 2024-06-27 02:28:35 +0000 (Thu, 27 Jun 2024)
Log Message:
-----------
finally making some progress on cubic spline fitting; still very much in progress; committing work so far to facilitate testing on different platforms
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-21 07:44:59 UTC (rev 7167)
+++ teem/trunk/src/limn/limn.h 2024-06-27 02:28:35 UTC (rev 7168)
@@ -517,21 +517,22 @@
** (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 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 */
+ double xy[8]; /* four control points of cubic Bezier:
+ x0, y0, x1, y1, x2, y2, x3, y3
+ 0 1 2 3 4 5 6 7 DIM=2 */
+ int corner[2]; /* corner[0,1] non-zero if xy[0,3] are either corner vertices
+ or path-ending vertices, i.e. reasons to not have geometric
+ continuity here, either because we intend to have a corner,
+ or because there's nothing else to be continuous with */
} limnCBFSeg;
/*
******** limnCBFPath
**
-** a multi-segment path in the context of cubic Bezier fitting
+** a multi-spline path in the context of cubic Bezier fitting
*/
typedef struct {
- limnCBFSeg *seg; /* array of limnCBFSeg */
+ limnCBFSeg *seg; /* array of limnCBFSeg structs (NOT pointers to them) */
unsigned int segNum; /* length of seg array */
airArray *segArr; /* manages seg and segNum */
int isLoop; /* path is closed loop */
@@ -540,7 +541,7 @@
/*
******** limnCBFCtx
**
-** The bag of state for limnCBF functions.
+** The complete bag of input parameters and state for limnCBF functions.
**
** note: "nrp" = Newton-based Re-Parameterization of where the given points
** fall along the spline, the iterative process inside limnCBFSingle
@@ -547,57 +548,61 @@
*/
typedef struct {
/* ----------- input ---------- */
- int verbose, /* verbosity level */
- cornNMS; /* non-minimal-suppression of corners: accept as
- corners only those with locally minimal angle */
+ int verbose, /* verbosity level */
+ cornerFind, /* do first search for corners: places where the path is not
+ geometrically continuous (between corners, the path is geometrically
+ continuous between multiple spline segments) */
+ cornerNMS; /* (if cornerFind) non-minimal-suppression of corners: accept as
+ corners only those with locally minimal angle */
unsigned int nrpIterMax; /* max # iters of nrp */
- double scale, /* scale (in sense of nrrdKernelDiscreteGaussian) at which to estimate
- spline endpoints and tangents; scale=0 means the endpoints are
- exactly on vertices, and tangents are from the smallest-support
- finite differences. This is the ONLY floating point that should be set
- by a method (limnCBFScaleSet); the rest can be set directly. */
- distMin, /* min distance to given points: this controls both splitting done by
- limnCBFMulti, and nrp within limnCBFSingle */
- nrpDeltaMax, /* in nrp, capping parameterization change to this scaling of average
- u[i+1]-u[i]. This wasn't in author's original code (so their idea of
- doing at most ~5 iters of nrp may no longer hold), but it can help
- stabilize things */
- nrpDistScl, /* scaling on distMin to use when testing distance during nrp; setting
- this < 1 means that nrp tries to be more stringent than the overall
- fitting, but with the benefit of sometimes being smarter about where
- to split, when that is needed */
- nrpPsi, /* don't even try nrp if max dist is bigger than nrpPsi*distMin, instead just
- subdivide */
- nrpDeltaMin, /* min total parameterization change by nrp */
- alphaMin, /* alpha can't be negative, and we enforce distinct positivity to ensure
- that spline doesn't slow down too much near endpoints */
- detMin, /* abs(determinant) of 2x2 matrix to invert can't go below this */
+ double
+ epsilon, /* error threshold on min distance from spline (as currently parameterized)
+ to given points: this affects both splitting done by limnCBFMulti, and
+ nrp within limnCBFSingle. Fitting has successfully finished when spline
+ path never strays further than this from input points */
+ scale, /* scale (in sense of nrrdKernelDiscreteGaussian) at which to estimate
+ spline endpoints and tangents; scale=0 means the endpoints are
+ exactly on vertices, and tangents are from the smallest-support
+ finite differences */
+ nrpCap, /* in nrp, cap parameterization change to this scaling of average
+ u[i+1]-u[i]. This wasn't in author's original code (so their idea of doing
+ at most ~5 iters of nrp no longer holds), but it can help stabilize things
+ with gnarly inputs */
+ nrpIota, /* (also not in author's original code) scaling on epsilon to use when
+ testing distance during nrp; setting this < 1 means that nrp tries to be
+ more stringent than the overall fitting, but with the benefit of
+ sometimes being smarter about where to split, when that is needed */
+ nrpPsi, /* don't even try nrp if max dist is bigger than nrpPsi*epsilon, instead
+ just subdivide */
+ nrpDeltaThresh, /* finish npr when mean parameterization change fall below this */
+ alphaMin, /* alpha can't be negative, and we enforce distinct positivity to ensure
+ that spline doesn't slow down too much near endpoints */
+ detMin, /* abs(determinant) of 2x2 matrix to invert can't go below this */
cornAngle; /* 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. */
+ straight line have an angle of 180 degrees. */
/* ----------- internal --------- */
- double *uu, /* buffer used for nrp */
+ double *uu, /* used for nrp: buffer of parameterizations in [0,1] of point along
+ currently considered spline segment */
*vw, /* weights for endpoint vertex calculation */
*tw, /* weights for endpoint tangent calculation */
- *mine; /* helps remember who allocated the above */
- unsigned int wLen; /* how long are vw, tw */
- double lenF2L; /* length of segment from first to last */
+ *cpp, /* x,y positions of corners */
+ *clt, /* x,y left (incoming) tangents at corners */
+ *crt; /* x,y right (outgoing) tangents at corners */
+ unsigned int uLen, /* how long is uu */
+ wLen, /* how long are vw, tw */
+ cNum; /* number of corners = # positions in cpp = # tangent vecs in clt, crt */
/* ----------- output --------- */
unsigned int nrpIterDone, /* number of nrp iters taken */
- distIdx; /* which point had distance distDone */
- double dist, /* max distance to given points */
- nrpDeltaDone, /* latest total parameterization change by nrp */
+ distMaxIdx; /* which point had distance distMax */
+ double distMax, /* max distance to given points */
+ nrpDeltaDone, /* latest mean parameterization change by nrp */
alphaDet; /* min det of matrix inverted to find alpha */
- int distBig; /* how big dist (above) is:
- 0: dist <= nD
- 1: nD < dist <= DM
- 2: DM < dist <= fD
- 3: fD < dist
- where
- nD = nrpDistScl*distMin,
- DM = distMin,
- fD = nrpPsi*distMin: */
+ int distBig; /* how big distMax (above) is:
+ 0: distMax <= wee (wee = nrpIota*epsilon)
+ 1: wee < distMax <= eps (eps = epsilon)
+ 2: eps < distMax <= far (far = nrpPsi*epsilon)
+ 3: far < distMax */
} limnCBFCtx;
/*
@@ -609,16 +614,16 @@
**
** NOTE: For now, point data is only double (not float), and only in 2D (not
** 3D), but if this becomes more general, that generality will be inside here.
-** For time being DIM=2 tags locations where 2D-ness is explicit in code.
+** For time being "DIM=2" tags locations where 2D-ness is explicit in code.
*/
typedef struct {
- /* assuming DIM=2: 2 values per logical element pp */
const double *pp; /* point coords, we do not own buffer */
double *ppOwn; /* point coords, we DO own buffer */
- unsigned int num; /* how many points */
- int isLoop; /* points form a loop: logical indices into coord
- array are . . . num-2, num-1, 0, 1, . . .
- and index 0 is effectively arbitrary */
+ unsigned int num, /* how many points */
+ dim; /* points live in what dimension (currently only dim=2 implemented) */
+ int isLoop; /* points form a loop: logical indices into coord
+ array are . . . num-2, num-1, 0, 1, . . .
+ and index 0 is effectively arbitrary */
} limnCBFPoints;
/* defaultsLimn.c */
@@ -897,25 +902,28 @@
double maxT);
/* splineFit.c */
-LIMN_EXPORT limnCBFPoints *limnCBFPointsNew(const double *pp, unsigned int nn,
+LIMN_EXPORT limnCBFPoints *limnCBFPointsNew(const void *pdata, int ptype,
+ unsigned int dim, unsigned int pnum,
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 *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 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 limnCBFFindTVT(double lt[2], double vv[2], double rt[2],
+ 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(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);
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,
@@ -922,8 +930,8 @@
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 double *xy,
- unsigned int pointNum, int isLoop);
+LIMN_EXPORT int limnCBFit(limnCBFPath *path, limnCBFCtx *fctx,
+ const limnCBFPoints *lpnt);
/* 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-21 07:44:59 UTC (rev 7167)
+++ teem/trunk/src/limn/lpu_cbfit.c 2024-06-27 02:28:35 UTC (rev 7168)
@@ -34,43 +34,42 @@
int pret;
Nrrd *_nin, *nin;
- double *xy, alpha[2], vv0[2], tt1[2], tt2[2], vv3[2], deltaMin, psi, cangle, distMin,
- distScl, utt1[2], utt2[2], time0, dtime, scale;
+ double *xy, deltaThresh, psi, cangle, epsilon, nrpIota, time0, dtime, scale, supow;
unsigned int ii, pNum, iterMax;
- int loop, petc, verbose, synth, nofit;
+ int loop, petc, verbose, tvt[4];
char *synthOut;
limnCBFCtx *fctx;
limnCBFPath *path;
+ limnCBFPoints *lpnt;
hestOptAdd_1_Other(&hopt, "i", "input", &_nin, NULL, "input xy points", nrrdHestNrrd);
hestOptAdd_1_Int(&hopt, "v", "verbose", &verbose, "1", "verbosity level");
- hestOptAdd_Flag(&hopt, "s", &synth, "synthesize xy points from control points");
hestOptAdd_1_String(&hopt, "so", "synth out", &synthOut, "",
- "if non-empty, filename in which to save synthesized xy pts");
- hestOptAdd_Flag(&hopt, "snf", &nofit,
- "actually do not fit, just save -so synthetic "
- "output and quit");
- hestOptAdd_2_Double(&hopt, "t1", "tan", utt1, "nan nan",
- "if non-nan, the outgoing tangent from the first point");
- hestOptAdd_2_Double(&hopt, "t2", "tan", utt2, "nan nan",
- "if non-nan, the incoming tangent to the last point");
+ "if non-empty, filename in which to save synthesized xy pts, "
+ "and then quit before any fitting.");
+ hestOptAdd_1_Double(&hopt, "sup", "expo", &supow, "1",
+ "when synthesizing data on a single segment, warp U parameters "
+ "by raising to this power.");
+ hestOptAdd_4_Int(&hopt, "tvt", "loi hii ofi 1s", tvt, "-1 -1 -1 -1",
+ "if all values are >= 0: make single call to "
+ "limnCBFFindTVT and quit");
hestOptAdd_1_UInt(&hopt, "im", "max", &iterMax, "0",
"(if non-zero) max # nrp iterations to run");
- hestOptAdd_1_Double(&hopt, "deltam", "delta", &deltaMin, "0.0005",
+ hestOptAdd_1_Double(&hopt, "deltathr", "delta", &deltaThresh, "0.0005",
"(if non-zero) stop nrp when change in spline "
"domain sampling goes below this");
- hestOptAdd_1_Double(&hopt, "distm", "dist", &distMin, "0.01",
+ hestOptAdd_1_Double(&hopt, "eps", "dist", &epsilon, "0.01",
"(if non-zero) stop nrp when distance between spline "
"and points goes below this");
- hestOptAdd_1_Double(&hopt, "dists", "scl", &distScl, "0.25",
- "scaling on nrp distMin check");
+ hestOptAdd_1_Double(&hopt, "iota", "scl", &nrpIota, "0.25",
+ "scaling on nrp epsilon check");
hestOptAdd_1_Double(&hopt, "psi", "psi", &psi, "10", "psi, of course");
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, "loop", &loop,
- "given xy points are actually a loop; BUT "
- "the first and last points need to be the same!");
+ "given xy points are actually a loop: the first point logically "
+ "follows the last point");
hestOptAdd_Flag(&hopt, "petc", &petc, "(Press Enter To Continue) ");
/*
hestOptAdd_1_String(&hopt, NULL, "output", &out, NULL, "output nrrd filename");
@@ -91,7 +90,7 @@
airMopError(mop);
return 1;
}
- if (synth && 6 != _nin->axis[1].size) {
+ if (airStrlen(synthOut) && 6 != _nin->axis[1].size) {
fprintf(stderr, "%s: need 2-by-6 array (not 2-by-%u) for synthetic xy\n", me,
(unsigned int)_nin->axis[1].size);
airMopError(mop);
@@ -107,10 +106,11 @@
return 1;
}
- if (!synth) {
+ if (!airStrlen(synthOut)) {
xy = (double *)nin->data;
pNum = (unsigned int)nin->axis[1].size;
} else {
+ double alpha[2], vv0[2], tt1[2], tt2[2], vv3[2];
/* synthesize data from control points */
double *cpt = (double *)nin->data;
limnCBFSeg seg;
@@ -135,7 +135,9 @@
xy = AIR_MALLOC(2 * pNum, double);
airMopAdd(mop, xy, airFree, airMopAlways);
for (ii = 0; ii < pNum; ii++) {
- limnCBFSegEval(xy + 2 * ii, &seg, AIR_AFFINE(0, ii, pNum - 1, 0, 1));
+ double uu = AIR_AFFINE(0, ii, pNum - 1, 0, 1);
+ uu = pow(uu, supow);
+ limnCBFSegEval(xy + 2 * ii, &seg, uu);
}
if (airStrlen(synthOut)) {
Nrrd *nsyn = nrrdNew();
@@ -147,42 +149,48 @@
airMopError(mop);
return 1;
}
- if (nofit) {
- fprintf(stderr, "%s: got -nf nofit; bye\n", me);
- airMopOkay(mop);
- return 0;
- }
+ printf("%s: saved synthetic output to %s; bye\n", me, synthOut);
+ airMopOkay(mop);
+ return 0;
}
}
- {
- /* set up 2-vector-valued arguments to fitting */
- double len;
- ELL_2V_COPY(vv0, xy);
- ELL_2V_COPY(vv3, xy + 2 * (pNum - 1));
- if (ELL_2V_EXISTS(utt1)) {
- ELL_2V_COPY(tt1, utt1);
- } else {
- /* TODO: better tangent estimation */
- ELL_2V_SUB(tt1, xy + 2, xy);
- }
- if (ELL_2V_EXISTS(utt2)) {
- ELL_2V_COPY(tt2, utt2);
- } else {
- ELL_2V_SUB(tt2, xy + 2 * (pNum - 2), vv3);
- }
- ELL_2V_NORM(tt1, tt1, len);
- ELL_2V_NORM(tt2, tt2, len);
+ if (!(lpnt = limnCBFPointsNew(xy, nrrdTypeDouble, 2, pNum, loop))) {
+ airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
+ fprintf(stderr, "%s: trouble setting up points:\n%s", me, err);
+ airMopError(mop);
+ return 1;
}
path = limnCBFPathNew();
airMopAdd(mop, path, (airMopper)limnCBFPathNix, airMopAlways);
- fctx = limnCBFCtxNew(pNum, scale);
+ fctx = limnCBFCtxNew();
+ fctx->verbose = verbose;
fctx->nrpIterMax = iterMax;
- fctx->nrpDeltaMin = deltaMin;
- fctx->distMin = distMin;
- fctx->nrpDistScl = distScl;
- fctx->verbose = verbose;
+ fctx->scale = scale;
+ fctx->epsilon = epsilon;
+ fctx->nrpDeltaThresh = deltaThresh;
+ fctx->nrpIota = nrpIota;
fctx->nrpPsi = psi;
fctx->cornAngle = cangle;
+ if (tvt[0] >= 0 && tvt[1] >= 0 && tvt[2] >= 0 && tvt[3] >= 0) {
+ double lt[2], vv[2], rt[2];
+ unsigned int loi = AIR_UINT(tvt[0]), hii = AIR_UINT(tvt[1]), ofi = AIR_UINT(tvt[2]);
+ int oneSided = tvt[3];
+ if (limnCBFCtxPrep(fctx, lpnt)
+ || limnCBFFindTVT(lt, vv, rt, fctx, lpnt, loi, hii, ofi, oneSided)) {
+ airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
+ fprintf(stderr, "%s: trouble doing single tangent-vertex-tangent:\n%s", me, err);
+ airMopError(mop);
+ return 1;
+ }
+ printf("%s: loi,hii=[%d,%d] ofi=%d oneSided=%d limnCBFFindTVT:\n", me, loi, hii, ofi,
+ oneSided);
+ printf(" lt = %g %g\n", lt[0], lt[1]);
+ printf(" vv = %g %g\n", vv[0], vv[1]);
+ printf(" rt = %g %g\n", rt[0], rt[1]);
+ printf("(quitting)\n");
+ airMopOkay(mop);
+ return 0;
+ }
time0 = airTime();
if (petc) {
fprintf(stderr, "%s: Press Enter to Continue ... ", me);
@@ -189,23 +197,23 @@
fflush(stderr);
getchar();
}
- if (limnCBFit(path, fctx, xy, pNum, loop)) {
+ if (limnCBFit(path, fctx, lpnt)) {
airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
- fprintf(stderr, "%s: trouble:\n%s", me, err);
+ fprintf(stderr, "%s: trouble doing fitting:\n%s", me, err);
airMopError(mop);
return 1;
}
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);
+ printf("%s: time= %g ms;iterDone= %u ;deltaDone=%g, distMax=%g (@%u)\n", me, dtime,
+ fctx->nrpIterDone, fctx->nrpDeltaDone, fctx->distMax, fctx->distMaxIdx);
{
unsigned int si;
printf("%s: path has %u segments:\n", me, path->segNum);
for (si = 0; si < path->segNum; si++) {
limnCBFSeg *seg = path->seg + si;
- printf("seg %u (%3u): (%g,%g) -- (%g,%g) -- (%g,%g) -- (%g,%g)\n", si, seg->pNum,
- seg->xy[0], seg->xy[1], seg->xy[2], seg->xy[3], seg->xy[4], seg->xy[5],
- seg->xy[6], seg->xy[7]);
+ printf("seg %u: (%g,%g) -- (%g,%g) -- (%g,%g) -- (%g,%g)\n", si, seg->xy[0],
+ seg->xy[1], seg->xy[2], seg->xy[3], seg->xy[4], seg->xy[5], seg->xy[6],
+ seg->xy[7]);
}
}
Modified: teem/trunk/src/limn/splineFit.c
===================================================================
--- teem/trunk/src/limn/splineFit.c 2024-06-21 07:44:59 UTC (rev 7167)
+++ teem/trunk/src/limn/splineFit.c 2024-06-27 02:28:35 UTC (rev 7168)
@@ -1,6 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
- Copyright (C) 2009--2022 University of Chicago
+ Copyright (C) 2009--2024 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -24,37 +24,68 @@
#include <assert.h>
/*
- This file contains GLK's implementation of the curve fitting described in:
- Philip J. Schneider. “An Algorithm for Automatically Fitting Digitized
- Curves”. In Graphics Gems, Academic Press, 1990, pp. 612–626.
- https://dl.acm.org/doi/10.5555/90767.90941
- The author's code is here:
- http://www.realtimerendering.com/resources/GraphicsGems/gems/FitCurves.c
+This file contains GLK's implementation of the curve fitting described in:
+Philip J. Schneider. "An Algorithm for Automatically Fitting Digitized Curves".
+In Graphics Gems, Academic Press, 1990, pp. 612–626.
+https://dl.acm.org/doi/10.5555/90767.90941
+The author's code is here:
+http://www.realtimerendering.com/resources/GraphicsGems/gems/FitCurves.c
- 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.
+The functions below do not use any other limnSpline structs or functions, since those
+were written a long time ago when GLK was even more ignorant than now about splines.
+Hopefully that other code can be revisited and re-organized for a later version of
+Teem, at which point the code below can be integrated with it.
+
+NOTE: spline fitting would be useful in 3D (or higher dimensions) too, but currently
+this code only supports 2D. "DIM=2" flags places in code where that is explicit.
*/
-limnCBFPoints * /* Biff: nope */
-limnCBFPointsNew(const double *pp, uint nn, int isLoop) {
- /* implicit: DIM=2 */
+#define PNMIN(ISLOOP) ((ISLOOP) ? 4 : 3)
+
+/*
+******** limnCBFPointsNew
+**
+** create a point data container, possibly around given pdata pointer. In an aspirational
+** hope of API stability, this is one of the few functions for which the interface itself
+** does not expose the specificity to DIM=2 and type double (though the code inside
+** does (apologetically) enforce that).
+*/
+limnCBFPoints * /* Biff: NULL */
+limnCBFPointsNew(const void *pdata, int ptype, uint dim, uint pnum, int isLoop) {
+ static const char me[] = "limnCBFPointsNew";
limnCBFPoints *lpnt;
+ if (airEnumValCheck(nrrdType, ptype)) {
+ biffAddf(LIMN, "%s: point data type %d not valid", me, ptype);
+ return NULL;
+ }
+ if (ptype != nrrdTypeDouble) {
+ biffAddf(LIMN, "%s: sorry, only %s-type data implemented now (not %s)", me,
+ airEnumStr(nrrdType, nrrdTypeDouble), airEnumStr(nrrdType, ptype));
+ return NULL;
+ }
+ if (2 != dim) {
+ biffAddf(LIMN, "%s: sorry, only 2-D data implemented now (not %u)", me, dim);
+ return NULL;
+ }
+ if (pnum < PNMIN(isLoop)) {
+ biffAddf(LIMN, "%s: need at least %u points in %s (not %u)", me, PNMIN(isLoop),
+ isLoop ? "loop" : "non-loop", pnum);
+ return NULL;
+ }
lpnt = AIR_CALLOC(1, limnCBFPoints);
assert(lpnt);
- if (pp) {
+ if (pdata) {
/* we are wrapping around a given pre-allocated buffer */
- lpnt->pp = pp;
+ lpnt->pp = pdata;
lpnt->ppOwn = NULL;
} else {
/* we are allocating our own buffer */
lpnt->pp = NULL;
- lpnt->ppOwn = AIR_CALLOC(nn, double);
+ lpnt->ppOwn = AIR_CALLOC(pnum, double);
assert(lpnt->ppOwn);
}
- lpnt->num = nn;
+ lpnt->num = pnum;
+ lpnt->dim = dim; /* but really DIM=2 because of above */
lpnt->isLoop = isLoop;
return lpnt;
}
@@ -79,7 +110,7 @@
biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
- pnmin = lpnt->isLoop ? 3 : 2;
+ pnmin = PNMIN(lpnt->isLoop);
if (!(lpnt->num >= pnmin)) {
biffAddf(LIMN, "%s: need %u or more points in limnCBFPoints (not %u)%s", me, pnmin,
lpnt->num, lpnt->isLoop ? " for loop" : "");
@@ -93,52 +124,73 @@
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)
+static void
+segInit(void *_seg) {
+ limnCBFSeg *seg = (limnCBFSeg *)_seg;
+ 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);
+ seg->corner[0] = seg->corner[1] = AIR_FALSE;
+ return;
+}
-/* number of points between low,high indices loi,hii */
-static uint
-pntNum(const limnCBFPoints *lpnt, uint loi, uint hii) {
- if (hii < loi) {
- assert(lpnt->isLoop);
- hii += lpnt->num;
+limnCBFPath * /* Biff: nope */
+limnCBFPathNew() {
+ limnCBFPath *path;
+ path = AIR_MALLOC(1, limnCBFPath);
+ if (path) {
+ path->segArr = airArrayNew((void **)(&path->seg), &path->segNum, sizeof(limnCBFSeg),
+ 128 /* incr */);
+ airArrayStructCB(path->segArr, segInit, NULL);
+ path->isLoop = AIR_FALSE;
}
- return hii - loi + 1;
+ return path;
}
-/* coordinates of point with index loi+ii */
-static const double *
-pntCrd(const limnCBFPoints *lpnt, uint loi, uint ii) {
- uint jj = loi + ii;
- while (jj >= lpnt->num)
- jj -= lpnt->num;
- return PP(lpnt) + 2 * jj; /* DIM=2 */
+limnCBFPath * /* Biff: nope */
+limnCBFPathNix(limnCBFPath *path) {
+ if (path) {
+ airArrayNuke(path->segArr);
+ free(path);
+ }
+ return NULL;
}
static void
+limnCBFPathJoin(limnCBFPath *dst, const limnCBFPath *src) {
+ uint bb = airArrayLenIncr(dst->segArr, src->segNum);
+ memcpy(dst->seg + bb, src->seg, (src->segNum) * sizeof(limnCBFSeg));
+ return;
+}
+
+/* initialize a freshly allocated limnCBFCtx struct;
+ the pointers therein do not point to anything valid */
+static void
ctxInit(limnCBFCtx *fctx) {
if (!fctx) return;
/* defaults for input parameters to various CBF functions */
fctx->verbose = 0;
- fctx->cornNMS = AIR_TRUE;
+ fctx->cornerFind = AIR_TRUE;
+ fctx->cornerNMS = AIR_TRUE;
fctx->nrpIterMax = 10;
- fctx->distMin = 0;
- fctx->nrpDeltaMax = 3.0;
- fctx->nrpDistScl = 0.8;
- fctx->nrpPsi = 6;
- fctx->nrpDeltaMin = 0.001;
+ fctx->epsilon = 0; /* will need to be set to something valid elsewhere */
+ fctx->scale = 0; /* will need to be set to something valid elsewhere */
+ fctx->nrpCap = 3.0;
+ fctx->nrpIota = 0.8;
+ fctx->nrpPsi = 10;
+ fctx->nrpDeltaThresh = 0.01;
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;
+ fctx->cornAngle = 120.0; /* degrees */
+ /* internal state */
+ /* initialize buffer pointers to NULL and buffer lengths to 0 */
+ fctx->uu = fctx->vw = fctx->tw = fctx->cpp = fctx->clt = fctx->crt = NULL;
+ fctx->uLen = fctx->wLen = fctx->cNum = 0;
/* initialize outputs to bogus valus */
fctx->nrpIterDone = (uint)(-1);
- fctx->distIdx = (uint)(-1);
- fctx->dist = AIR_POS_INF;
+ fctx->distMaxIdx = (uint)(-1);
+ fctx->distMax = AIR_POS_INF;
fctx->nrpDeltaDone = AIR_POS_INF;
fctx->alphaDet = 0;
fctx->distBig = 0;
@@ -145,118 +197,215 @@
return;
}
+limnCBFCtx * /* Biff: nope */
+limnCBFCtxNew() {
+ limnCBFCtx *ret;
+
+ ret = AIR_CALLOC(1, limnCBFCtx);
+ assert(ret);
+ ctxInit(ret);
+ 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);
+ if (fctx->cpp) free(fctx->cpp);
+ if (fctx->clt) free(fctx->clt);
+ if (fctx->crt) free(fctx->crt);
+ free(fctx);
+ }
+ return NULL;
+}
+
/*
-** ctxBuffersSet: allocates in fctx:
-** uu, vw, tw
+** ctxBuffersSet: ensures that some buffers in fctx: uu, vw, tw
+** are set up for current #points pNum and measurement scale fctx->scale.
+** The buffers are re-allocated only when necessary.
+** Does NOT touch the corner-related buffers: cpp, clt, crt
*/
static int /* Biff: 1 */
-ctxBuffersSet(limnCBFCtx *fctx, uint pNum, double scl) {
+ctxBuffersSet(limnCBFCtx *fctx, uint pNum) {
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;
+ double scale;
+ uint ulen, ii;
- 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);
+ if (!fctx) {
+ biffAddf(LIMN, "%s: got NULL pointer", me);
return 1;
}
- if (0 == scl) {
- /* will do simplest possible finite differences; we're done */
- fctx->vw = fctx->tw = NULL;
- return 0;
+ scale = fctx->scale;
+ if (!pNum || scale < 0 || !AIR_EXISTS(scale)) {
+ biffAddf(LIMN, "%s: pNum %u or scale %g not valid", me, pNum, scale);
+ return 1;
}
- /* 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);
+
+ /* assuming pNum = lpnt->num for points lpnt in consideration, this allocation size
+ (big enough to parameterize all the points at once) is safe, though it will likely
+ be excessive given how the path may be split at corners into separate segments. */
+ ulen = pNum * 2; /* DIM=2 */
+ if (ulen != fctx->uLen) {
+ airFree(fctx->uu);
+ if (!(fctx->uu = AIR_CALLOC(ulen, double))) {
+ biffAddf(LIMN, "%s: failed to allocate uu buffer (%u doubles)", me, ulen);
+ return 1;
+ }
+ }
+ fctx->uLen = ulen;
+
+ if (0 == scale) {
+ /* will do simplest possible finite differences; no need for weights */
+ fctx->vw = airFree(fctx->vw);
+ fctx->tw = airFree(fctx->tw);
+ fctx->wLen = 0;
+ } else {
+ /* 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; it reflects something
+ about the confidence that the nrrdKernelDiscreteGaussian is working as expected,
+ rather than something about tuning cubic spline fitting */
+ const double one = 0.999;
+ /* if vw and tw are allocated for length wlbig (or bigger) something isn't right */
+ const uint wlbig = 128;
+ /* if the initial weights for the tangent computation sum to smaller than this
+ (they will be later normalized to sum to 1) then something isn't right */
+ const double tinytsum = 1.0 / 64;
+ double kparm[2], vsum, tsum;
+ uint wlen;
+ /* else need to (possibly allocate and) set vw and tw buffers */
+ kparm[0] = scale;
+ kparm[1] = 100000; /* effectively no cut-off; sanity check comes later */
+ ii = 0;
+ vsum = 0;
+ do {
+ vsum += (!ii ? 1 : 2) * nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
+ ii++;
+ } while (vsum < one && ii < wlbig);
+ /* wlen = intended length of blurring kernel weight vectors */
+ wlen = ii;
+ if (wlen > wlbig) {
+ biffAddf(LIMN,
+ "%s: weight buffer length %u (from scale %g) seems "
+ "too large",
+ me, wlen, scale);
+ return 1;
+ }
+ if (wlen > pNum / 2) {
+ biffAddf(LIMN,
+ "%s: weight buffer length %u (from scale %g) seems "
+ "too large compared to #points %u",
+ me, wlen, scale, pNum);
+ return 1;
+ }
+ if (wlen != fctx->wLen) {
+ airFree(fctx->vw);
+ airFree(fctx->tw);
+ if (!((fctx->vw = AIR_CALLOC(wlen, double))
+ && (fctx->tw = AIR_CALLOC(wlen, double)))) {
+ biffAddf(LIMN, "%s: couldn't allocate weight buffers (%u doubles)", me, wlen);
+ return 1;
+ }
+ fctx->wLen = wlen;
+ }
+ /* normalization intent:
+ 1 = sum_i(vw[|i|]) for i=-(len-1)...len-1
+ 1 = sum_i(tw[i]) for i=1...len-1 */
+ vsum = tsum = 0;
+ for (ii = 0; ii < wlen; ii++) {
+ double kw = nrrdKernelDiscreteGaussian->eval1_d(ii, kparm);
+ vsum += (!ii ? 1 : 2) * (fctx->vw[ii] = kw);
+ tsum += (fctx->tw[ii] = ii * kw);
+ }
+ if (tsum < tinytsum) {
+ biffAddf(LIMN,
+ "%s: scale %g led to tiny unnormalized tangent weight sum %g; "
+ "purpose of scale is to do blurring but scale %g won't do that",
+ me, scale, tsum, scale);
+ return 1;
+ }
+ for (ii = 0; ii < wlen; ii++) {
+ fctx->vw[ii] /= vsum;
+ fctx->tw[ii] /= tsum;
+ if (fctx->verbose) {
+ printf("%s: ii=%3u v=%0.17g t=%0.17g\n", me, ii, fctx->vw[ii],
+ fctx->tw[ii]);
+ }
+ }
+ } /* else scale > 0 */
+
+ return 0;
+}
+
+/*
+******** limnCBFCtxPrep
+**
+** checks the things that are going to be passed around a lot,
+** and makes call to initialize buffers inside fctx
+*/
+int /* Biff: 1 */
+limnCBFCtxPrep(limnCBFCtx *fctx, const limnCBFPoints *lpnt) {
+ static const char me[] = "limnCBFCtxPrep";
+
+ if (!(fctx && lpnt)) {
+ biffAddf(LIMN, "%s: got NULL pointer", me);
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);
+ if (limnCBFPointsCheck(lpnt)) {
+ biffAddf(LIMN, "%s: problem with points", me);
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);
+ if (!(fctx->epsilon > 0)) {
+ biffAddf(LIMN, "%s: need positive epsilon (not %g)", me, fctx->epsilon);
+ return 1;
}
- 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]); */
+ if (!(fctx->scale >= 0)) {
+ biffAddf(LIMN, "%s: need non-negative scale (not %g)", me, fctx->scale);
+ return 1;
}
- 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 (!(fctx->nrpCap > 0)) {
+ biffAddf(LIMN, "%s: need positive nrpCap (not %g)", me, fctx->nrpCap);
+ return 1;
}
- if (!(scale >= 0)) {
- biffAddf(LIMN, "%s: need scale >= 0 (not %g)", me, scale);
- return NULL;
+ if (!(0 < fctx->nrpIota && fctx->nrpIota <= 1)) {
+ biffAddf(LIMN, "%s: nrpIota (%g) must be in (0,1]", me, fctx->nrpIota);
+ return 1;
}
- ret = AIR_CALLOC(1, limnCBFCtx);
- if (!ret) {
- biffAddf(LIMN, "%s: allocation failure?", me);
- return NULL;
+ if (!(1 <= fctx->nrpPsi)) {
+ biffAddf(LIMN, "%s: nrpPsi (%g) must be >= 1", me, fctx->nrpPsi);
+ return 1;
}
- ctxInit(ret);
- if (ctxBuffersSet(ret, pointNum, scale)) {
- biffAddf(LIMN, "%s: trouble allocating buffers", me);
- free(ret);
- return NULL;
+ if (!(fctx->nrpDeltaThresh > 0)) {
+ biffAddf(LIMN, "%s: need positive nrpDeltaThresh (not %g) ", me,
+ fctx->nrpDeltaThresh);
+ return 1;
}
- ret->scale = scale;
+ if (!(fctx->alphaMin > 0)) {
+ biffAddf(LIMN, "%s: need positive alphaMin (not %g) ", me, fctx->alphaMin);
+ return 1;
+ }
+ if (!(fctx->detMin > 0)) {
+ biffAddf(LIMN, "%s: need positive detMin (not %g) ", me, fctx->detMin);
+ return 1;
+ }
+ {
+ const double amin = 60;
+ const double amax = 180;
+ if (!(amin <= fctx->cornAngle && fctx->cornAngle <= amax)) {
+ biffAddf(LIMN, "%s: cornAngle (%g) outside sane range [%g,%g]", me,
+ fctx->cornAngle, amin, amax);
+ return 1;
+ }
+ }
+ if (ctxBuffersSet(fctx, lpnt->num)) {
+ biffAddf(LIMN, "%s: trouble setting up buffers", me);
+ return 1;
+ }
- return ret;
+ return 0;
}
-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) */
@@ -339,144 +488,179 @@
return;
}
-/*
-** 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 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
-** enforces the possible corner-ness of those indices (which prevents
-** vertices past corners influencing how vv or tt are found)
-*/
-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;
- /* 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 */
+/* cheesy macro as short-hand to access either pp or ppOwn */
+#define PP(lpnt) ((lpnt)->pp ? (lpnt)->pp : (lpnt)->ppOwn)
- if (!(/* vv can be NULL */ tt && fctx && lpnt)) {
- biffAddf(LIMN, "%s: got NULL pointer", me);
+/* error-checked index processing for limnCBFFindVT and maybe others */
+static int /* Biff: 1 */
+idxPrep(int *sloP, int *shiP, int *loopyP, const limnCBFPoints *lpnt, uint loi,
+ uint hii) {
+ static const char me[] = "idxPrep";
+ int slo, shi, loopy, spanlen;
+
+ *sloP = *shiP = 10 * lpnt->num; /* initialize to bogus indices */
+ if (!(loi < lpnt->num && hii < lpnt->num)) {
+ biffAddf(LIMN, "%s: loi %u, hii %u not both < #points %u", me, loi, hii, lpnt->num);
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);
+ if (loi == hii && hii != 0) {
+ biffAddf(LIMN,
+ "%s: can only have loi == hii if both 0 (not %u), to signify no bounds in "
+ "point loop",
+ me, loi);
+ return 1;
}
- /* 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);
+ slo = AIR_INT(loi);
+ shi = AIR_INT(hii);
+ if (lpnt->isLoop) {
+ loopy = (0 == loi && 0 == hii);
+ if (hii < loi) {
+ shi += lpnt->num;
+ }
+ } else {
+ if (0 == loi && 0 == hii) {
+ biffAddf(LIMN, "%s: can only have loi == hii == 0 with point loop", me);
+ return 1;
+ }
+ if (hii < loi) {
+ biffAddf(LIMN, "%s: can only have hii (%u) < loi (%u) in a point loop", me, hii,
+ loi);
+ return 1;
+ }
+ loopy = AIR_FALSE;
+ }
+ spanlen = shi - slo + 1;
+ if (spanlen <= 1 && !loopy) {
+ biffAddf(LIMN, "%s: [%u,%u]->[%d,%d] span length %d <= 1 but not in loop", me, loi,
+ hii, slo, shi, spanlen);
return 1;
}
- if (_hii < _loi && !lpnt->isLoop) {
- biffAddf(LIMN, "%s: _hii %u < _loi %u sensible only in point loop", me, _hii, _loi);
+ /* all's well, set output values */
+ *sloP = slo;
+ *shiP = shi;
+ *loopyP = loopy;
+ return 0;
+}
+
+static void
+subnorm2(double dir[2], const double aa[2], const double bb[2]) {
+ double len;
+ ELL_2V_SUB(dir, aa, bb); /* dir = aa - bb */
+ ELL_2V_NORM(dir, dir, len); /* normalize(dir) */
+}
+
+/* limnCBFFindTVT: Find constraints for spline fitting: incoming/left tangent lt, center
+or endpoint vertex vv, outgoing/right tangent rt; any but not all can be NULL. These are
+computed from the given points lpnt, at offset index ofi within index range loi, hii
+and only that range: that range is probably delimited by corners, and we have to be blind
+to anything past the corners on either side of us (except if loi==hii==0, in which case
+we can look at all the points in a loop).
+
+NOTE: this assumes that limnCBFCtxPrep(fctx, lpnt) was called without error!
+It (via ctxBuffersSet) allocates things that we depend on here
+*/
+int /* Biff: 1 */
+limnCBFFindTVT(double lt[2], double vv[2], double rt[2], const limnCBFCtx *fctx,
+ const limnCBFPoints *lpnt, uint loi, uint hii, uint ofi, int oneSided) {
+ static const char me[] = "limnCBFFindTVT";
+ /* we use here (signed) int for things that might seem better as uint, but it
+ simplifies implementing arithmetic and comparisons given how indices wrap around in
+ point loops */
+ int slo, /* signed version of loi */
+ shi, /* signed version of hii, but shi = hii + lpnt->num if hii < loi in loop */
+ sof, /* signed versions of ofi */
+ loopy, /* lpnt->isLoop && loi == 0 && hii = 0, i.e. there are no bounds on indices */
+ pnum, /* total number of points in lpnts */
+ spanlen, /* span length: number of points in [loi,hii] */
+ 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 (!((lt || vv || rt) && fctx && lpnt)) {
+ biffAddf(LIMN, "%s: got NULL pointer (or too many NULL pointers)", me);
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 = 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);
+ /* so: each of lt, vv, rt can be NULL (they just can't be all NULL) */
+ if (idxPrep(&slo, &shi, &loopy, lpnt, loi, hii)) {
+ biffAddf(LIMN, "%s: trouble with loi %u or hii %u", me, loi, hii);
return 1;
}
- /* now sgsz >= 2 */
- if (!(_ofi < AIR_UINT(sgsz))) {
- biffAddf(LIMN, "%s: _ofi %u too high for segment size %d", me, _ofi, sgsz);
+ spanlen = shi - slo + 1;
+ pnum = AIR_INT(lpnt->num);
+
+ /* now:
+ 0 == slo == shi implies lpnt->isLoop (and this is called "loopy")
+ slo == shi != 0 is always impossible
+ always: slo < shi (even if given hii < loi), and thus any indices in range [slo,shi]
+ need to be mod'd with pnum before indexing into PP(lpnt)
+ spanlen >= 2 */
+ if (!loopy && !(ofi < AIR_UINT(spanlen))) {
+ biffAddf(LIMN,
+ "%s: ofi %u too high for [%u,%u]->[%d,%d] span length %d (not in loop)", me,
+ ofi, loi, hii, slo, shi, spanlen);
return 1;
}
- /* now _ofi is a valid index in [0,..,sgsz-1] */
- ofi = AIR_INT(_ofi);
- icent = loi + ofi;
+ /* now ofi is a valid index in [0,..,spanlen-1] */
+ sof = AIR_INT(ofi);
+ icent = slo + sof;
iplus = icent + 1;
imnus = icent - 1;
- if (lpnt->isLoop) {
- icent = AIR_MOD(icent, pNum);
- iplus = AIR_MOD(iplus, pNum);
- imnus = AIR_MOD(imnus, pNum);
- } else {
- icent = AIR_CLAMP(loi, icent, hii);
- iplus = AIR_CLAMP(loi, iplus, hii);
- imnus = AIR_CLAMP(loi, imnus, hii);
+ if (!loopy) {
+ /* this is the code that motivated if (hii < loi) shi += pnum;
+ otherwise clamping is too annoying */
+ icent = AIR_CLAMP(slo, icent, shi);
+ iplus = AIR_CLAMP(slo, iplus, shi);
+ imnus = AIR_CLAMP(slo, imnus, shi);
}
+ /* now map to actual indices */
+ icent = AIR_MOD(icent, pnum);
+ iplus = AIR_MOD(iplus, pnum);
+ imnus = AIR_MOD(imnus, pnum);
if (0 == fctx->scale) {
- const double *xy, *xyP, *xyM;
- int mi, ci, pi;
+ /* DIM=2 through-out */
+ const double *xyC = PP(lpnt) + 2 * icent, *xyP, *xyM;
if (vv) {
- ELL_2V_COPY(vv, PP(lpnt) + 2 * icent); /* DIM=2 */
+ ELL_2V_COPY(vv, xyC);
}
- switch (dir) {
- case 1:
- pi = iplus;
- mi = icent;
- break;
- case 0:
- pi = iplus;
- mi = imnus;
- break;
- case -1:
- /* mi and pi switched to point other way */
- mi = icent;
- pi = imnus;
- break;
+ /* printf("!%s: iplus=%u imnus=%u xy=%g %g\n", me, iplus, imnus, xy[0], xy[1]); */
+ xyP = PP(lpnt) + 2 * iplus;
+ xyM = PP(lpnt) + 2 * imnus;
+ if (rt) {
+ subnorm2(rt, xyP, oneSided ? xyC : xyM);
}
- 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;
+ if (lt) {
+ subnorm2(lt, xyM, oneSided ? xyC : xyP);
}
- /* 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 {
/* 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};
+ /* for simplicity: regardless of dir, we compute average positions for points
+ centered around loi + ofi (posC), and for lower/higher indices (posM/posP) */
+ double posM[2] = {0, 0}, posC[2] = {0, 0}, posP[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] */
+ int lim = (int)fctx->wLen - 1, /* limit on loop index */
+ ci; /* loops through [-lim,lim] */
if (!(vw && tw)) {
biffAddf(LIMN, "%s: fctx internal buffers vw and tw not both allocated", me);
return 1;
}
- /* various signed indices */
- /* printf("!%s: ofi = %d, dir=%d\n", me, ofi, dir); */
- for (ci = -smax; ci <= smax; ci++) {
+ if (tw[0] != 0) {
+ biffAddf(LIMN, "%s: first tangent weight fctx->tw[0] %g not zero", me, tw[0]);
+ return 1;
+ }
+ for (ci = -lim; ci <= lim; ci++) {
uint wi = abs(ci); /* weight index into vw, tw */
- int di = loi + ofi + ci; /* signed index into data */
+ int di = slo + sof + ci; /* signed index into data */
const double *xy;
- if (lpnt->isLoop) {
- di = AIR_MOD(di, pNum);
- } else {
- di = AIR_CLAMP(loi, di, hii);
- }
+ if (!loopy) di = AIR_CLAMP(slo, di, shi);
+ di = AIR_MOD(di, pnum);
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]); */
+ /* 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,
@@ -488,789 +672,149 @@
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;
+ {
+ /* limit distance from chosen (x,y) datapoint to posC to be (HEY harcoded) 95% of
+ fctx->epsilon. Being allowed to be further away can cause annoyances (for GLK in
+ some early stage of debugging) */
+ double off[2], offlen, okoff = 0.95 * fctx->epsilon; /* DIM=2 throughout */
+ const double *xy = PP(lpnt) + 2 * icent; /* center vertex in given data */
+ ELL_2V_SUB(off, posC, xy); /* off = posC - xy, from given to computed */
+ ELL_2V_NORM(off, off, offlen);
+ offlen = AIR_MIN(okoff, offlen);
+ /* difference between chosen (x,y) datapoint and spline endpoint
+ can be in any direction, but we limit the length */
+ ELL_2V_SCALE_ADD2(posC, 1, xy, offlen, off);
}
- ELL_2V_NORM(tt, tt, len);
+ if (lt) {
+ subnorm2(lt, posM, oneSided ? posC : posP);
+ }
+ if (rt) {
+ subnorm2(rt, posP, oneSided ? posC : posM);
+ }
if (vv) {
ELL_2V_COPY(vv, posC);
- /* some post-proceessing of computed spline endpoint */
- 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 */
- 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, operp, pp);
}
}
return 0;
}
-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 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 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);
- return 1;
- }
- if (lpnt->isLoop) {
- limnCBFFindVT(vv0, tt1, fctx, lpnt, loi, hii, 0, 0);
- ELL_2V_COPY(vv3, vv0);
- ELL_2V_SCALE(tt2, -1, tt1);
- } else {
- limnCBFFindVT(vv0, tt1, fctx, lpnt, loi, hii, 0, +1);
- limnCBFFindVT(vv3, tt2, fctx, lpnt, loi, hii, hii - loi, -1);
- }
- if (given) {
- *given = AIR_FALSE;
- }
- }
- return 0;
-}
-
/*
-** (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.
+******** limnCBFit
**
-** There are various conditions where the generated spline ignores the
-** xy array and instead is what one could call a "simple arc" (with
-** control points at 1/3 and 2/3 the distance between the end points):
-** - having only two points (xy contains only the end points)
-** - the determinant of the 2x2 matrix that is inverted to solve
-** for alpha is too close to zero (this test was not part of the
-** 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?)
+** top-level function for fitting cubic beziers to given points
*/
-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 const char me[] = "findalpha";
- uint ii, pNum;
- double det;
-
- pNum = pntNum(lpnt, loi, hii);
- if (pNum > 2) {
- 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++) {
- const double *xy;
- 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 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);
- ELL_2V_SCALE_ADD2(Pi, bb[0] + bb[1], vv0, bb[2] + bb[3], vv3);
- xy = pntCrd(lpnt, loi, ii);
- ELL_2V_SUB(dmP, xy, Pi);
- xx[0] += ELL_2V_DOT(dmP, Ai1);
- xx[1] += ELL_2V_DOT(dmP, Ai2);
- }
- ELL_4V_SET(MM, m11, m12, m12, m22);
- ELL_2M_INV(MI, MM, det);
- ELL_2MV_MUL(alpha, MI, xx);
- } else { /* pNum <= 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
- && alpha[0] > (fctx->lenF2L) * (fctx->alphaMin)
- && alpha[1] > (fctx->lenF2L) * (fctx->alphaMin))) {
- if (fctx->verbose) {
- printf("%s: bad |det| %g (vs %g) or alpha %g,%g (vs %g*%g) "
- "--> simple arc\n",
- me, AIR_ABS(det), fctx->detMin, alpha[0], alpha[1], fctx->lenF2L,
- fctx->alphaMin);
- }
- /* generate simple arc: set both alphas to 1/3 of distance from
- first to last point, but also handle non-unit-length tt1 and
- tt2 */
- alpha[0] = fctx->lenF2L / (3 * ELL_2V_LEN(tt1));
- alpha[1] = fctx->lenF2L / (3 * ELL_2V_LEN(tt2));
- } else {
- if (fctx->verbose > 1) {
- printf("%s: all good: det %g, alpha %g,%g\n", me, det, alpha[0], alpha[1]);
- }
- }
- fctx->alphaDet = det;
- return;
-}
-
-/*
-** using Newton iterations to try to find a better places at which
-** to evaluate the spline in order to match the given points xy
-*/
-static double
-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 limnCBFPoints *lpnt, uint loi,
- uint hii) {
- static const char me[] = "reparm";
- uint ii, pNum;
- double vv1[2], vv2[2], delta, maxdelu;
- double *uu = fctx->uu;
-
- pNum = pntNum(lpnt, loi, hii);
- assert(pNum >= 3);
- /* average u[i+1]-u[i] is 1/(pNum-1) */
- maxdelu = fctx->nrpDeltaMax / (pNum - 1);
- ELL_2V_SCALE_ADD2(vv1, 1, vv0, alpha[0], tt1);
- ELL_2V_SCALE_ADD2(vv2, 1, vv3, alpha[1], tt2);
- delta = 0;
- /* only changing parameterization of interior points,
- not the first (ii=0) or last (ii=pNum-1) */
- for (ii = 1; ii < pNum - 1; ii++) {
- double numer, denom, delu, df[2], ww[4], tt, Q[2], QD[2], QDD[2];
- const double *xy;
- tt = uu[ii];
- CBD0(Q, vv0, vv1, vv2, vv3, tt, ww);
- CBD1(QD, vv0, vv1, vv2, vv3, tt, ww);
- CBD2(QDD, vv0, vv1, vv2, vv3, tt, ww);
- xy = pntC...
[truncated message content] |
