/***************************************************************************
 * Func: spline_interpolation                                              *
 *                                                                         *
 * Desc: takes 16 unsigned characters in an array and outputs the results  *
 *       of a B-spline interpolation.  The sample points of the            *
 *       interpolation function f(-1-d) and f(-d) are evaluated as f(1+d)  *
 *       and f(d).  The computed constants are a0 = f(1+d)                 *
 *                                             a1 = f(d)                   *
 *                                             a2 = f(1-d)                 *
 *                                             a3 = f(2-d)                 *
 *                                                                         *
 * Params: image - array of 16 image sample points                         *
 *                 (0,0)           (0,3)                                   *
 *                   *---*---*---*---*                                     *
 *                   |   |   |   |   |                                     *
 *                   *---*---*---*---*                                     *
 *                   |   | o |   |   |    o is point to be interpolated    *
 *                   *---*---*---*---*                                     *
 *                   |   |   |   |   |                                     *
 *                   *---*---*---*---*                                     *
 *                 (3,0)           (3,3)                                   *
 *                                                                         *
 *         x_diff - x difference from interpolated address to image sample *
 *         y_diff - y difference from interpolated address to image sample *
 *                                                                         *
 * Returns: interpolated value                                             *
 ***************************************************************************/

unsigned char spline_interpolation(unsigned char image[4][4], float x, float y)
    {
    double column[4];                /* storage for row interpolation results */
    double a0, a1, a2, a3;           /* interpolation coefficients */
    double x_plus_1, y_plus_1;       /* x + 1, y + 1 */
    double one_minus_x, one_minus_y; /* 1 - x, 1 - y */
    int i;                           /* loop index */
    double pixel;                    /* newly interpolated pixel value */

    /* Do we even need to interpolate ? */
    if((x == 0.0) && (y == 0.0))
	return image[1][1];

    if(x == 0.0)                    /* Do we need to interpolate the rows? */
	for(i=0; i<4; i++)          /* NO */
	    column[i] = (double) image[i][1];
    else                            /* YES, interpolate the rows */
	{
	x_plus_1 = x + 1.0;
	one_minus_x = 1.0 - x;
	two_minus_x = 2.0 - x;

	a0 = ((-0.16666666 * x_plus_1 + 1.0) * x_plus_1 - 2.0) * x_plus_1 + 1.33333;
	a1 = (0.5 * x - 1.0) * x * x + 0.666666667;
	a2 = (0.5 * one_minus_x - 1.0) * one_minus_x * one_minus_x + 0.66666667;
	a3 = 1.0 - a0 - a1 - a2;

	if(y == 0.0)
	    return (unsigned char) (a0 * image[1][0] + a1 * image[1][1] +
				    a2 * image[1][2] + a3 * image[1][3]);
	for(i=0; i<4; i++)
	    column[i] = a0 * image[i][0] + a1 * image[i][1] +
			a2 * image[i][2] + a3 * image[i][3];
	}

    /* interpolate the resulting column */

    y_plus_1 = y + 1.0;
    one_minus_y = 1.0 - y;
    two_minus_y = 2.0 - y;

    a0 = ((-0.16666666 * y_plus_1 + 1.0) * y_plus_1 - 2.0) * y_plus_1 + 1.33333;
    a1 = (0.5 * y - 1.0) * y * y + 0.666666667;
    a2 = (0.5 * one_minus_y - 1.0) * one_minus_y * one_minus_y + 0.666666667;
    a3 = 1.0 - a0 - a1 - a2;

    pixel = column[0] * a0 + column[1] * a1 + column[2] * a2 + column[3] * a3;
    return (unsigned char) pixel;
    }