/***************************************************************************
 * Func: cubic_interpolation                                               *
 *                                                                         *
 * Desc: takes 16 unsigned characters in an array and outputs the results  *
 *       of a cubic 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 *
 *         a - cubic spline interpolation function constant                *
 *                                                                         *
 * Returns: interpolated value                                             *
 ***************************************************************************/

unsigned char cubic_interpolation(unsigned char image[4][4],
				  float x, float y, double a)
    {
    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 */
    double two_minus_x, two_minus_y; /* 2 - x, 2 - y */
    int i;                           /* loop index */
    double pixel;                    /* newly interpolated value */

    /* Do we really 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, we do need to interpolate rows */
	{
	x_plus_1 = x + 1.0;
	one_minus_x = 1.0 - x;
	two_minus_x = 2.0 - x;

	a0 = ((a * x_plus_1 - 5.0 * a) * x_plus_1 + 8.0 * a) * x_plus_1 - 4.0 * a;
	a1 = ((a + 2.0) * x - (a + 3.0)) * x * x + 1;
	a2 = ((a + 2.0) * one_minus_x - (a + 3.0)) * one_minus_x * one_minus_x + 1;
	a3 = ((a * two_minus_x - 5.0 * a) * two_minus_x + 8.0 * a) * two_minus_x - 4.0 * a;

	if(y == 0.0)   /* can we just process one row */
	    {
	    pixel = a0 * image[1][0] + a1 * image[1][1] +
		    a2 * image[1][2] + a3 * image[1][3];
	    CLIP(pixel, 0.0, 255.0);
	    return (unsigned char) pixel;
	    }
	else
	    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 = ((a * y_plus_1 - 5.0 * a) * y_plus_1 + 8.0 * a) * y_plus_1 - 4.0 * a;
    a1 = ((a + 2.0) * y - (a + 3.0)) * y * y + 1;
    a2 = ((a + 2.0) * one_minus_y - (a + 3.0)) * one_minus_y * one_minus_y + 1;
    a3 = ((a * two_minus_y - 5.0 * a) * two_minus_y + 8.0 * a) * two_minus_y - 4.0 * a;

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