warp_mbyte_field_bicubic
Tcl_Obj* imageObj
Tcl_Obj* xvecObj
Tcl_Obj* yvecObj

/*
 * Warp image using the given vector field.
 */

crimp_image* image;
crimp_image* xvector;
crimp_image* yvector;
crimp_image* result;
int x, y, c;
double xf, yf;

crimp_input_any (imageObj, image);
CRIMP_ASSERT_NOTIMGTYPE (image, float);
CRIMP_ASSERT_NOTIMGTYPE (image, grey16);
CRIMP_ASSERT_NOTIMGTYPE (image, grey32);

crimp_input (xvecObj, xvector, float);
crimp_input (yvecObj, yvector, float);

if (!crimp_eq_dim (xvector, yvector)) {
    Tcl_SetResult(interp, "Unable to warp, expected equally-sized coordinate fields", TCL_STATIC);
    return TCL_ERROR;
}

/*
 * Create result and scan through it, sampling the input under the guidance of
 * the coordinate fields.
 */

result = crimp_new_at (image->itype, crimp_x (xvector), crimp_y (xvector), crimp_w (xvector), crimp_h (xvector));

for (y = 0; y < crimp_h (result); y++) {
    for (x = 0; x < crimp_w (result); x++) {
	int ixw, iyw;

	xf = FLOATP (xvector, x, y);
	yf = FLOATP (yvector, x, y);

	/*
	 * Perform bicubic interpolation (1,2) using the nearest 4x4 pixels
	 * around the sampling location.
	 *
	 * (Ad 1) http://en.wikipedia.org/wiki/Bicubic_interpolation
	 * (Ad 2) http://www.paulinternet.nl/?page=bicubic
	 */

        ixw = xf;  xf -= ixw;
        iyw = yf;  yf -= iyw;

	ixw --; xf += 1.; xf /= 3.;
	iyw --; yf += 1.; yf /= 3.;

#undef  SAMPLE
#define SAMPLE(dx,dy) ((((ixw+(dx)) < 0) || ((ixw+(dx)) >= crimp_w (image)) || ((iyw+(dy)) < 0) || ((iyw+(dy)) >= crimp_h (image))) ? BLACK : (CH (image, c, (ixw+(dx)), (iyw+(dy)))))

        for (c = 0; c < 4; ++c) {
	    double p00 = SAMPLE(0,0);
	    double p01 = SAMPLE(0,1);
	    double p02 = SAMPLE(0,2);
	    double p03 = SAMPLE(0,3);
	    double p10 = SAMPLE(1,0);
	    double p11 = SAMPLE(1,1);
	    double p12 = SAMPLE(1,2);
	    double p13 = SAMPLE(1,3);
	    double p20 = SAMPLE(2,0);
	    double p21 = SAMPLE(2,1);
	    double p22 = SAMPLE(2,2);
	    double p23 = SAMPLE(2,3);
	    double p30 = SAMPLE(3,0);
	    double p31 = SAMPLE(3,1);
	    double p32 = SAMPLE(3,2);
	    double p33 = SAMPLE(3,3);

	    double a00 =      p11;
	    double a01 = -.50*p10 +  .50*p12;
	    double a02 =      p10 - 2.50*p11 + 2.00*p12 -  .50*p13;
	    double a03 = -.50*p10 + 1.50*p11 - 1.50*p12 +  .50*p13;
	    double a10 = -.50*p01 +  .50*p21;
	    double a11 =  .25*p00 -  .25*p02 -  .25*p20 +  .25*p22;
	    double a12 = -.50*p00 + 1.25*p01 -      p02 +  .25*p03 +  .50*p20 - 1.25*p21 +      p22 -  .25*p23;
	    double a13 =  .25*p00 -  .75*p01 +  .75*p02 -  .25*p03 -  .25*p20 +  .75*p21 -  .75*p22 +  .25*p23;
	    double a20 =      p01 - 2.50*p11 + 2.00*p21 -  .50*p31;
	    double a21 = -.50*p00 +  .50*p02 + 1.25*p10 - 1.25*p12 -      p20 +      p22 +  .25*p30 -  .25*p32;
	    double a22 =      p00 - 2.50*p01 + 2.00*p02 -  .50*p03 - 2.50*p10 + 6.25*p11 - 5.00*p12 + 1.25*p13 +
		2.00*p20 - 5.00*p21 + 4.00*p22 -     p23 - .50*p30 + 1.25*p31 -     p32 + .25*p33;
	    double a23 = -.50*p00 + 1.50*p01 - 1.50*p02 +  .50*p03 + 1.25*p10 - 3.75*p11 + 3.75*p12 - 1.25*p13 -
		p20 + 3.00*p21 - 3.00*p22 +     p23 + .25*p30 -  .75*p31 + .75*p32 - .25*p33;
	    double a30 = -.50*p01 + 1.50*p11 - 1.50*p21 +  .50*p31;
	    double a31 =  .25*p00 -  .25*p02 -  .75*p10 +  .75*p12 +  .75*p20 -  .75*p22 -  .25*p30 +  .25*p32;
	    double a32 = -.50*p00 + 1.25*p01 -      p02 +  .25*p03 + 1.50*p10 - 3.75*p11 + 3.00*p12 -  .75*p13 -
		1.50*p20 + 3.75*p21 - 3.00*p22 + .75*p23 + .50*p30 - 1.25*p31 +     p32 - .25*p33;
	    double a33 =  .25*p00 -  .75*p01 +  .75*p02 -  .25*p03 -  .75*p10 + 2.25*p11 - 2.25*p12 +  .75*p13 +
		.75*p20 - 2.25*p21 + 2.25*p22 - .75*p23 - .25*p30 +  .75*p31 - .75*p32 + .25*p33;

	    double x2 = xf * xf;
	    double x3 = x2 * xf;
	    double y2 = yf * yf;
	    double y3 = y2 * yf;

	    double r =
		(a00 + a01 * yf + a02 * y2 + a03 * y3)      +
		(a10 + a11 * yf + a12 * y2 + a13 * y3) * xf +
		(a20 + a21 * yf + a22 * y2 + a23 * y3) * x2 +
		(a30 + a31 * yf + a32 * y2 + a33 * y3) * x3;	    

	    CH (result, c, x, y) = CLAMPT (0, int, r, MAXVAL_GREY8);
        }
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;


/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */
