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