#include "defs.h"
#include "globals.h"
#include "tables.h"
/* &&Module spline.p */
/*
* Procedures below may make free use of the global variables arrayXY [list of control points] pointmatrix [list of spline
* segments] knot [list of spline knots] catrommtx [matrix for Catmull-Rom splines] bsplmtx [matrix for B-splines]
* lastPoint, intervals
*/
/*-----------------------------------------------------*/
int
max(a, b)
int a, b;
{
if (a > b)
return a;
else
return b;
}
/*-----------------------------------------------------*/
int
min(a, b)
int a, b;
{
if (a < b)
return a;
else
return b;
}
/*-----------------------------------------------------*/
static void
matXvector(m, v, result)
double (*m)[4];
double *v, *result;
{
/* IN */
/* IN */
/* OUT */
double t[4];
t[0] = v[0] * m[0][0] + v[1] * m[0][1] + v[2] * m[0][2] + v[3] * m[0][3];
t[1] = v[0] * m[1][0] + v[1] * m[1][1] + v[2] * m[1][2] + v[3] * m[1][3];
t[2] = v[0] * m[2][0] + v[1] * m[2][1] + v[2] * m[2][2] + v[3] * m[2][3];
t[3] = v[0] * m[3][0] + v[1] * m[3][1] + v[2] * m[3][2] + v[3] * m[3][3];
result[0] = t[0];
result[1] = t[1];
result[2] = t[2];
result[3] = t[3];
}
/*-----------------------------------------------------*/
/* actually the dot-product */
static double
vecXvec(v1, v2)
double *v1, *v2;
{
return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] + v1[3] * v2[3]);
}
/*------------------------------------------------------*/
/* basXctl is the pre-computed BasisMatrix times the control-point vector */
static double
splinePosition(basXctl, t)
double *basXctl;
double t;
{
/* IN */
double tvect[4]; /* vector of t values for spline matrix */
tvect[3] = 1.0;
tvect[2] = t;
tvect[1] = t * t;
if (tvect[1] <= MINREAL)/* avoid underflow problems */
tvect[1] = 0.0;
tvect[0] = t * tvect[1];/* t^3 */
return (vecXvec(tvect, basXctl));
}
/*-------------------------------------------------*/
static int
TwoToThe(n)
int n;
{
int i, tmp;
for (i = 0, tmp = 1; i < n; i++)
tmp *= 2;
return (tmp);
}
/*------------------------------------------------------*/
double
distance(x0, y0, x1, y1)
double x0, y0, x1, y1;
{
return (sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)));
}
/*------------------------------------------------------*/
/*
* compute the number of subdivisions for this span. We do this by a quadrature method and a simple linear-distance metric.
* This is not optimal in the number of subdivisions actually required, but is computationally efficient and accurate to the
* nearest power of 2 .
*/
static int32
numsubdivisions(XtimesBas, YtimesBas, resolution)
double *XtimesBas, *YtimesBas;
int32 resolution;
{
int32 n;
double t, x0, y0, xt, yt;
x0 = splinePosition(XtimesBas, 0.0);
y0 = splinePosition(YtimesBas, 0.0);
t = 1.0;
n = 0;
xt = splinePosition(XtimesBas, t);
yt = splinePosition(YtimesBas, t);
while ((int32) floor(distance(x0, y0, xt, yt) + 0.5) > resolution || n < 1)
{
t /= 2.0; /* perform the quadrature */
n++;
xt = splinePosition(XtimesBas, t);
yt = splinePosition(YtimesBas, t);
} /* while */
return (TwoToThe(n));
}
/*------------------------------------------------------------------------*/
/*
* compute new control vertices such that the resulting spline will interpolate through the old control points. This will work
* as int32 as the actual arc length between consecutive nodes does not vary from span to span. The method is noted in Wu and
* Abel's CACM 20(10) Oct 77 paper but the actual working method is from Barsky and Greenberg's paper in CG&IP 14(3) Nov 1980
* pp.203-226
*/
static void
invertsplvertices(numpts, isclosed, xys)
int32 numpts;
int isclosed;
int32(*xys)[2];
{
/* INOUT */
int32 i;
double beta[MAXCTLPTS + 1], Xrprime[MAXCTLPTS + 1], Yrprime[MAXCTLPTS + 1];
ControlPoints tempxys;
/* compute the values of beta */
beta[1] = 0.25;
for (i = 2; i <= numpts + 1; i++)
beta[i] = 1.0 / (4.0 - beta[i - 1]);
/* and the r primes from the original vertices */
Xrprime[1] = beta[1] * xys[1][0] * 5.0;
Yrprime[1] = beta[1] * xys[1][1] * 5.0;
for (i = 2; i < numpts; i++)
{
Xrprime[i] = beta[i] * (6.0 * xys[i][0] - Xrprime[i - 1]);
Yrprime[i] = beta[i] * (6.0 * xys[i][1] - Yrprime[i - 1]);
} /* for */
Xrprime[numpts] = beta[numpts] * (5.0 * xys[numpts][0] - Xrprime[numpts - 1]);
Yrprime[numpts] = beta[numpts] * (5.0 * xys[numpts][1] - Yrprime[numpts - 1]);
/* Now perform the back-substitution from the bottom up */
tempxys[numpts][0] = (int32) floor(Xrprime[numpts] + 0.5);
tempxys[numpts][1] = (int32) floor(Yrprime[numpts] + 0.5);
for (i = numpts - 1; i >= 1; i--)
{
tempxys[i]
[0] = (int32) floor(Xrprime[i] - beta[i] * tempxys[i + 1][0] + 0.5);
tempxys[i]
[1] = (int32) floor(Yrprime[i] - beta[i] * tempxys[i + 1][1] + 0.5);
}
if (isclosed)
{
/*
* at this point, we've probably been through one control-point adjustment, so let's not muck it up
*/
tempxys[numpts + 1][0] = tempxys[1][0];
tempxys[numpts + 1][1] = tempxys[1][1];
tempxys[numpts + 2][0] = tempxys[2][0];
tempxys[numpts + 2][1] = tempxys[2][1];
tempxys[0][0] = tempxys[numpts][0];
tempxys[0][1] = tempxys[numpts][1];
/* copy them back */
for (i = 0; i <= numpts + 2; i++)
{
xys[i][0] = tempxys[i][0];
xys[i][1] = tempxys[i][1];
}
return;
} /* closed */
/* copy back */
for (i = 2; i < numpts; i++)
{
xys[i][0] = tempxys[i][0];
xys[i][1] = tempxys[i][1];
}
/* open */
}
/*-----------------------------------------------------*/
/*
* adjust the list of control points so that we can use it for B-spline interpolation. Add any "phantom" vertices necessary so
* that the end conditions will be correct for interpolation
*/
static void
Bctlptadjust(isclosed, isarc, n, xys, thx)
int isclosed, isarc;
int32 *n;
int32(*xys)[2];
VThickness *thx;
{ /* ctlpt adjust */
/* INOUT */
/* INOUT */
/* INOUT */
int32 j;
ControlPoints tmp;
ThickAryType tmpthx;
int32 FORLIM;
if (isclosed)
{
if (*n == 2)
{
complain(ERRBAD);
fprintf(logfile,
"A closed spline requires more than 2 control points \n");
fprintf(logfile, "making a temporary fix in order to continue...\n");
xys[3][0] = xys[1][0];
xys[3][1] = xys[1][1];
}
FORLIM = *n;
for (j = 1; j <= FORLIM; j++)
{
tmp[j][0] = xys[j][0];
tmp[j][1] = xys[j][1];
tmpthx[j] = thx[j];
}
/* Now take care of the 'phantom' vertices */
tmp[*n + 1][0] = xys[1][0];
tmp[*n + 1][1] = xys[1][1];
tmpthx[*n + 1] = thx[1];
tmp[*n + 2][0] = xys[2][0];
tmp[*n + 2][1] = xys[2][1];
tmpthx[*n + 2] = thx[2];
tmp[*n + 3][0] = xys[3][0];
tmp[*n + 3][1] = xys[3][1];
tmpthx[*n + 3] = thx[3];
if (!isarc)
{
tmp[0][0] = xys[*n][0]; /* wrap around to the real last point */
tmp[0][1] = xys[*n][1];
tmpthx[0] = thx[*n];
}
else
{
tmp[0][0] = xys[0][0];
tmp[0][1] = xys[0][1];
tmpthx[0] = thx[0];
}
(*n)++; /* we supplied the 'last' point for the user */
FORLIM = *n + 2;
for (j = 0; j <= FORLIM; j++)
{
xys[j][0] = tmp[j][0];
xys[j][1] = tmp[j][1];
thx[j] = tmpthx[j];
} /* for */
return;
} /* if closed */
/*
* here, we have to supply the last 'real' point for the user, and add three phantoms-- one before, and two after
*/
if (!isarc)
{
tmp[0][0] = xys[1][0] * 2 - xys[2][0];
tmp[0][1] = xys[1][1] * 2 - xys[2][1];
}
else
{
tmp[0][0] = xys[0][0];
tmp[0][1] = xys[0][1];
}
tmpthx[0] = thx[1];
FORLIM = *n;
for (j = 1; j <= FORLIM; j++)
{
tmp[j][0] = xys[j][0];
tmp[j][1] = xys[j][1];
tmpthx[j] = thx[j];
}
tmp[*n + 1][0] = xys[*n][0] * 2 - xys[*n - 1][0];
tmp[*n + 1][1] = xys[*n][1] * 2 - xys[*n - 1][1];
tmpthx[*n + 1] = thx[*n];
tmp[*n + 2][0] = tmp[*n + 1][0];
tmp[*n + 2][1] = tmp[*n + 1][1];
tmpthx[*n + 2] = thx[*n];
FORLIM = *n + 2;
for (j = 0; j <= FORLIM; j++)
{
xys[j][0] = tmp[j][0];
xys[j][1] = tmp[j][1];
thx[j] = tmpthx[j];
} /* for */
/* OPEN SPLINE */
/* if open */
}
/*-----------------------------------------------------*/
/*
* adjust the list of control points so that we can use it for simple Catmull-Rom spline interpolation. Add any "phantom"
* vertices necessary so that the end conditions will be correct for interpolation
*/
static void
CRctlptadjust(isclosed, isarc, n, xys, thx)
int isclosed, isarc;
int32 *n;
int32(*xys)[2];
VThickness *thx;
{ /* ctlpt adjust */
/* INOUT */
/* INOUT */
/* INOUT */
int32 j;
ControlPoints tmp;
ThickAryType tmpthx;
int32 FORLIM;
if (isclosed)
{
if (*n == 2)
{
complain(ERRBAD);
fprintf(logfile,
"A closed spline requires more than 2 control points \n");
fprintf(logfile, "making a temporary fix in order to continue...\n");
xys[3][0] = xys[1][0];
xys[3][1] = xys[1][1];
}
FORLIM = *n;
for (j = 1; j <= FORLIM; j++)
{
tmp[j][0] = xys[j][0];
tmp[j][1] = xys[j][1];
tmpthx[j] = thx[j];
}
/* the phantom vertices */
tmp[*n + 1][0] = xys[1][0];
tmp[*n + 1][1] = xys[1][1];
tmpthx[*n + 1] = thx[1];
tmp[*n + 2][0] = xys[2][0];
tmp[*n + 2][1] = xys[2][1];
tmpthx[*n + 2] = thx[2];
tmp[*n + 3][0] = xys[3][0];
tmp[*n + 3][1] = xys[3][1];
tmpthx[*n + 3] = thx[3];
if (!isarc)
{
tmp[0][0] = xys[*n][0]; /* wrap around to the real last point */
tmp[0][1] = xys[*n][1];
tmpthx[0] = thx[*n];
}
else
{
tmp[0][0] = xys[0][0];
tmp[0][1] = xys[0][1];
tmpthx[0] = thx[0];
}
(*n)++; /* we supplied the 'last' point for the user */
FORLIM = *n + 2;
for (j = 0; j <= FORLIM; j++)
{
xys[j][0] = tmp[j][0];
xys[j][1] = tmp[j][1];
thx[j] = tmpthx[j];
} /* for */
return;
} /* if closed */
/*
* here, we have to supply the last 'real' point for the user, and add three phantoms-- one before, and two after
*/
if (!isarc)
{
tmp[0][0] = xys[1][0]; /* double the first point */
tmp[0][1] = xys[1][1];
}
else
{
tmp[0][0] = xys[0][0];
tmp[0][1] = xys[0][1];
}
tmpthx[0] = thx[1];
FORLIM = *n;
for (j = 1; j <= FORLIM; j++)
{
tmp[j][0] = xys[j][0];
tmp[j][1] = xys[j][1];
tmpthx[j] = thx[j];
}
tmp[*n + 1][0] = xys[*n][0]; /* and triple the last */
tmp[*n + 1][1] = xys[*n][1];
tmpthx[*n + 1] = thx[*n];
tmp[*n + 2][0] = xys[*n][0];
tmp[*n + 2][1] = xys[*n][1];
tmpthx[*n + 2] = thx[*n];
FORLIM = *n + 2;
for (j = 0; j <= FORLIM; j++)
{
xys[j][0] = tmp[j][0];
xys[j][1] = tmp[j][1];
thx[j] = tmpthx[j];
} /* for */
/* OPEN SPLINE */
/* if open */
} /* ctlpt adjust */
/*----------------------------------------------------------*/
static void
interpsplines(splinetype, isclosed, isanArc, linepatt,
basismatrix, numctls, arrayXY, pointmatrix,
varythicks, thickmatrix, TTmatrix)
SplineKind splinetype;
int isclosed, isanArc;
LineStyle linepatt;
double (*basismatrix)[4];
int32 numctls;
int32(*arrayXY)[2];
int32(*pointmatrix)[2];
int varythicks;
VThickness *thickmatrix, *TTmatrix;
{ /* interp splines */
/* IN */
/* IN */
/* OUT */
/* IN */
/* OUT */
double xctl[4], yctl[4]; /* vectors of x, y posits of control points */
double wctl[4];/* vector of thicknesses at each ctl pt */
double t, incr;
int32 Pi; /* P sub i */
int32 currpt;
ScaledPts theresolution;
if (!isclosed && isanArc) /* lie a little */
numctls++;
switch (splinetype)
{
case BSPL:
Bctlptadjust(isclosed, isanArc, &numctls, arrayXY, thickmatrix);
break;
case CARD:
case CATROM:
CRctlptadjust(isclosed, isanArc, &numctls, arrayXY, thickmatrix);
break;
case INTBSPL:
if (isclosed)
{
Bctlptadjust(true, isanArc, &numctls, arrayXY, thickmatrix);
invertsplvertices(numctls, true, arrayXY);
}
else
{
invertsplvertices(numctls, false, arrayXY);
Bctlptadjust(false, isanArc, &numctls, arrayXY, thickmatrix);
} /* else */
break;
/* Interpolating Bsplines */
}
if (!isclosed && isanArc) /* UN-lie a little */
numctls--;
/*
* this is the scheme: val := t-vector * Basis matrix * point matrix [t^3 t^2 t 1] * [[Ms]] * [Pi-1
* Pi Pi+1 Pi+2] where "Pi-1" is "P sub (i-1)", etc.
*
* But we do this in a round about way: Point matrix * basis then * t-vector will yield the single value
*
* there are certainly faster ways to do this, but this is the easiest to understand
*/
currpt = 1;
switch (linepatt)
{
case solid:
theresolution = MAXVECLENsp;
break;
case dotted:
case dashed:
case dotdash: /* ### */
theresolution = MAXVECLENsp * 3;
break;
}
for (Pi = 1; Pi < numctls; Pi++)
{
xctl[0] = float_(arrayXY[Pi - 1][0]);
xctl[1] = float_(arrayXY[Pi][0]);
xctl[2] = float_(arrayXY[Pi + 1][0]);
xctl[3] = float_(arrayXY[Pi + 2][0]);
yctl[0] = float_(arrayXY[Pi - 1][1]);
yctl[1] = float_(arrayXY[Pi][1]);
yctl[2] = float_(arrayXY[Pi + 1][1]);
yctl[3] = float_(arrayXY[Pi + 2][1]);
matXvector(basismatrix, xctl, xctl);
matXvector(basismatrix, yctl, yctl);
/*
* compute the delta-t increment for this segment based on a metric for subdivision
*/
intervals = numsubdivisions(xctl, yctl, theresolution);
if (linepatt == solid && intervals <= 2)
intervals *= 2;
incr = 1.0 / intervals;
/* avoid over-flowing the "pointmatrix" */
if (currpt + intervals - 1 >= MAXSPLINESEGS)
{
complain(ERRREALBAD);
fprintf(logfile, "error: Too many spline segments required.\n");
fprintf(logfile,
" Reducing the number of control points to get output.\n");
goto _L32;
}
t = 0.0;
while (t < 0.999999999)
{
pointmatrix[currpt - 1][0] = (int32) floor(splinePosition(xctl, t) + 0.5);
pointmatrix[currpt - 1][1] = (int32) floor(splinePosition(yctl, t) + 0.5);
if (varythicks)
{
wctl[0] = float_((int32) thickmatrix[Pi - 1]);
wctl[1] = float_((int32) thickmatrix[Pi]);
wctl[2] = float_((int32) thickmatrix[Pi + 1]);
wctl[3] = float_((int32) thickmatrix[Pi + 2]);
matXvector(catrommtx, wctl, wctl); /* requires using Catmull-Rom */
TTmatrix[currpt] = (int32) floor(splinePosition(wctl, t) + 0.5);
}
t += incr;
currpt++;
} /* while loop */
} /* for loop */
_L32:
/* the END-condtion */
pointmatrix[currpt - 1][0] = (int32) floor(splinePosition(xctl, 1.0) + 0.5);
pointmatrix[currpt - 1][1] = (int32) floor(splinePosition(yctl, 1.0) + 0.5);
if (varythicks)
{
wctl[0] = thickmatrix[numctls - 2];
wctl[1] = thickmatrix[numctls - 1];
wctl[2] = thickmatrix[numctls];
wctl[3] = thickmatrix[numctls + 1];
matXvector(catrommtx, wctl, wctl); /* requires using Catmull-Rom */
TTmatrix[currpt] = (int32) floor(splinePosition(wctl, 1.0) + 0.5);
}
lastPoint = currpt;
} /* interpsplines */
/*----------------------------------------------------------------*/
void
drawSpline(splinetype, isclosed, isanArc, patt, numctls, arrayXY,
pointmatrix, varythicks, thickmatrix, TTmatrix)
SplineKind splinetype;
int isclosed, isanArc;
LineStyle patt;
int32 numctls;
int32(*arrayXY)[2];
int32(*pointmatrix)[2];
int varythicks;
VThickness *thickmatrix, *TTmatrix;
{
/* IN */
/* OUT */
/* IN */
/* OUT */
lastPoint = 0;
switch (splinetype)
{
case CATROM:
interpsplines(splinetype, isclosed, isanArc, patt, catrommtx, numctls,
arrayXY, pointmatrix, varythicks, thickmatrix, TTmatrix);
break;
case CARD:
interpsplines(splinetype, isclosed, isanArc, patt, cardmtx, numctls,
arrayXY, pointmatrix, varythicks, thickmatrix, TTmatrix);
break;
case BSPL:
interpsplines(splinetype, isclosed, isanArc, patt, bsplmtx, numctls,
arrayXY, pointmatrix, varythicks, thickmatrix, TTmatrix);
break;
case INTBSPL:
interpsplines(splinetype, isclosed, isanArc, patt, bsplmtx, numctls,
arrayXY, pointmatrix, varythicks, thickmatrix, TTmatrix);
break;
} /* Case */
}