#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 */ }