warp_fpcomplex_projective_bilinear
Tcl_Obj* imageObj
Tcl_Obj* forwardObj

/*
 * Warp image using the given specified transform. The result is made large
 * enough to contain all of the warped image, and will contain meta data about
 * the location of the actual (0,0) origin point relative to the physical top
 * left corner of the result. This last is required because translations in
 * the transform may move pixels to negative positions which we cannot express
 * with the regular memory grid.
 */

crimp_image* image;
crimp_image* forward;
crimp_image* backward;
crimp_image* result;
int x, y, xt, yt, origx, origy;
double xf, yf;

crimp_input (imageObj,   image,   fpcomplex);
crimp_input (forwardObj, forward, fpcomplex);

if (!crimp_require_dim (forward, 3, 3)) {
    Tcl_SetResult(interp, "bad matrix dimensions, expected 3x3", TCL_STATIC);
    return TCL_ERROR;
}

backward = crimp_la_invert_matrix_3x3 (forward);
if (!backward) {
    Tcl_SetResult(interp, "Unable to invert singular matrix", TCL_STATIC);
    return TCL_ERROR;
}

/*
 * Determine size of the result, and the location of the origin point inside
 * based on the four corners of the input image and the forward transformation.
 */

result = crimp_geo_warp_init (image, forward, &origx, &origy);

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

	xf = xt;
	yf = yt;
	crimp_geo_warp_point (backward, &xf, &yf);

	/*
	 * Perform bilinear interpolation (1) using the nearest 2x2 pixels
	 * around the sampling location.
	 *
	 * (Ad 1) http://en.wikipedia.org/wiki/Bilinear_interpolation
	 */

        ixw = xf;
        iyw = yf;

        xf -= ixw;
        yf -= iyw;

	{
            float re = 0;
            float im = 0;
	    int ix, iy;
            for (iy = MAX(iyw, 0); iy < MIN(iyw + 2, crimp_h (image)); iy++) {
                yf = 1 - yf;
                for (ix = MAX(ixw, 0); ix < MIN(ixw + 2, crimp_w (image)); ix++) {
                    xf = 1 - xf;

		    re += RE (image, ix, iy) * yf * xf;
		    im += IM (image, ix, iy) * yf * xf;
                }
            }

	    RE (result, x, y) = re;
	    IM (result, x, y) = im;
        }
    }
}

crimp_del (backward);
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:
 */
