/**************************************************************************
 *{@C
 *      Copyright:      1988-2020 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         ImageProcessing
 *      Filename:       IP_WarpFunct.c
 *
 *      Author:         Paul Obermeier
 *
 *	Description:	Functions for free-form image warping and blending.
 *
 *      Additional documentation:
 *                      None.
 *
 *	Exported functions:
 *			IP_WarpFunct
 *			IP_BlendFunct
 *
 **************************************************************************/

#include <stdio.h>
#include <math.h>

#include "UT_Compat.h"
#include "UT_Macros.h"
#include "UT_Error.h"
#include "UT_Memory.h"
#include "UT_Parallel.h"

#define USE_IMG_MACROS
#include "FF_Image.h"

#include "IP_Image.h"
#include "IP_ImagePrivate.h"

/***************************************************************************
 *
 *	Macros and type definitions
 *
 ***************************************************************************/

/* A structure to store a rectangular grid of pixel
displacement and sampling area size data */

typedef struct {
	Int32	numrow,		/* Number of rows in the grid */
		numcol;		/* Number of columns in the grid */
	Float32	*xdata,		/* Displacement, x component */
		*ydata,		/* Displacement, y component */
		*sdata;		/* Sampling area size */
	} xyinterp;



/* A structure to store a rectangular grid of pixel blending factors */

typedef struct {
	Int32	numrow,		/* Number of rows in the grid */
		numcol;		/* Number of columns in the grid */
	Float32	*fdata;		/* Blending factors */
	} finterp;



/* Parameters passed to the "conc_warp" function */

typedef struct {
	IP_ImageId	srcimg,		/* Source image */
			destimg;	/* Destination image */
	IP_FillModeType	oor;		/* Border treatment mode */
	pyramid		*srcpyr;	/* Source image pyramid */
	xyinterp	xyitp;		/* X-Y interpolation data */
	Float32		mx,		/* Coordinate scaling factors */
			my,
			msize;
	} warpparms;



/* Parameters passed to the "conc_blend_1" and "conc_blend_2" functions */

typedef struct {
	IP_ImageId	srcimg,		/* Source image */
			destimg;	/* Destination image */
	IP_FillModeType	oor;		/* Border treatment mode */
	pyramid		*srcpyr;	/* Source image pyramid */
	xyinterp	xyitp;		/* X-Y interpolation data */
	finterp		fitp;		/* F interpolation data */
	Float32		mx,		/* Coordinate scaling factors */
			my,
			msize;
	} blendparms;



/* Amount of memory reserved for the concurrently
invoked functions' auto variable stacks */

#define stacksize	(1 << 12)

/***************************************************************************
 *
 *	Image warping and blending functions
 *
 ***************************************************************************/

/* Find a pixel displacement, (xout, yout), and a sampling patch size, "size",
by sampling the data stored in "xyitp" at position (xin, yin). */

static void sample_xyinterp
	(Float32 xin, Float32 yin,
	 const xyinterp *xyitp,
	 Float32 *xout, Float32 *yout, Float32 *size)
{
    Int32		col0, col1, row0, row1,
		LL_offset, LR_offset, UL_offset, UR_offset;
    Float32	clipx, clipy, fcol, frow, u0, u1, v0, v1;
    Float32	LL, LR, UL, UR;

    clipx = UT_MAX (0.0, UT_MIN (xin, 1.0));
    clipy = UT_MAX (0.0, UT_MIN (yin, 1.0));
    fcol = clipx * (xyitp->numcol - 1);
    col0 = (Int32)fcol;
    col1 = UT_MIN (col0 + 1, xyitp->numcol - 1);
    u0 = fcol - col0;
    u1 = 1.0 - u0;
    frow = clipy * (xyitp->numrow - 1);
    row0 = (Int32)frow;
    row1 = UT_MIN (row0 + 1, xyitp->numrow - 1);
    v0 = frow - row0;
    v1 = 1.0 - v0;
    LL_offset = row0 * xyitp->numcol + col0;
    LR_offset = row0 * xyitp->numcol + col1;
    UL_offset = row1 * xyitp->numcol + col0;
    UR_offset = row1 * xyitp->numcol + col1;
    LL = xyitp->xdata[LL_offset];
    LR = xyitp->xdata[LR_offset];
    UL = xyitp->xdata[UL_offset];
    UR = xyitp->xdata[UR_offset];
    *xout = u1 * (v1 * LL + v0 * UL) + u0 * (v1 * LR + v0 * UR);
    LL = xyitp->ydata[LL_offset];
    LR = xyitp->ydata[LR_offset];
    UL = xyitp->ydata[UL_offset];
    UR = xyitp->ydata[UR_offset];
    *yout = u1 * (v1 * LL + v0 * UL) + u0 * (v1 * LR + v0 * UR);
    LL = xyitp->sdata[LL_offset];
    LR = xyitp->sdata[LR_offset];
    UL = xyitp->sdata[UL_offset];
    UR = xyitp->sdata[UR_offset];
    *size = u1 * (v1 * LL + v0 * UL) + u0 * (v1 * LR + v0 * UR);
    return;
}



/* Find a pixel weighting factor, "fout", by sampling the data stored in
"fitp" at position (xin, yin). */

static void sample_finterp
	(Float32 xin, Float32 yin,
	 const finterp *fitp,
	 Float32 *fout)
{
    Int32		col0, col1, row0, row1;
    Float32	clipx, clipy, fcol, frow, u0, u1, v0, v1;
    Float32	LL, LR, UL, UR;

    clipx = UT_MAX (0.0, UT_MIN (xin, 1.0));
    clipy = UT_MAX (0.0, UT_MIN (yin, 1.0));
    fcol = clipx * (fitp->numcol - 1);
    col0 = (Int32)fcol;
    col1 = UT_MIN (col0 + 1, fitp->numcol - 1);
    u0 = fcol - col0;
    u1 = 1.0 - u0;
    frow = clipy * (fitp->numrow - 1);
    row0 = (Int32)frow;
    row1 = UT_MIN (row0 + 1, fitp->numrow - 1);
    v0 = frow - row0;
    v1 = 1.0 - v0;
    LL = fitp->fdata[row0 * fitp->numcol + col0];
    LR = fitp->fdata[row0 * fitp->numcol + col1];
    UL = fitp->fdata[row1 * fitp->numcol + col0];
    UR = fitp->fdata[row1 * fitp->numcol + col1];
    *fout = u1 * (v1 * LL + v0 * UL) + u0 * (v1 * LR + v0 * UR);
    return;
}



/* Concurrent main loop for function "warp". */

static void conc_warp (void *ptr, Int32 n_conc, Int32 i_conc)
{
    IP_ImageId		srcimg, destimg;
    IP_FillModeType	oor;
    pyramid		*srcpyr;
    const xyinterp	*xyitp;
    Float32		mx, my, msize, xs, ys, size, cdata[FF_NumImgChanTypes];
    Int32		i, xd, yd, ymin, ymax;

    /* Initialization */

    srcimg = ((warpparms *) ptr)->srcimg;
    destimg = ((warpparms *) ptr)->destimg;
    oor = ((warpparms *) ptr)->oor;
    srcpyr = ((warpparms *) ptr)->srcpyr;
    xyitp = &((warpparms *) ptr)->xyitp;
    mx = ((warpparms *) ptr)->mx;
    my = ((warpparms *) ptr)->my;
    msize = ((warpparms *) ptr)->msize;
    n_conc = (n_conc == 0? 1: n_conc);
    ymin = (destimg->pimg.desc.height * i_conc) / n_conc;
    ymax = (destimg->pimg.desc.height * (i_conc + 1)) / n_conc;
    for (i = 0; i < FF_NumImgChanTypes; ++i) {
	cdata[i] = 0.0;
    }

    /* Main loop */

    for (yd = ymin; yd < ymax; ++yd) {
	for (xd = 0; xd < destimg->pimg.desc.width; ++xd) {
	    sample_xyinterp (xd * mx, yd * my, xyitp, &xs, &ys, &size);
	    xs *= srcimg->pimg.desc.width - 1;
	    ys *= srcimg->pimg.desc.height - 1;
	    size *= msize;
	    IP_SamplePyramid (srcpyr, xs, ys, size, cdata, oor);
	    for (i = 0; i < FF_NumImgChanTypes; ++i) {
		switch (destimg->pimg.channel[i]) {
		    case FF_ImgFmtTypeNone: {
			break;
		    }
		    case FF_ImgFmtTypeUByte: {
			*FF_ImgUBytePixel (&destimg->pimg, i, xd, yd) =
			    (UInt8) (255.0 * cdata[i]);
			break;
		    }
		    case FF_ImgFmtTypeFloat: {
			*FF_ImgFloatPixel (&destimg->pimg, i, xd, yd) =
			    cdata[i];
			break;
		    }
		    default: {
			UT_ErrFatal ("conc_warp", "invalid channel type");
		    }
		}
	    }
	}
    }
    return;
}



/* Concurrent main loop for function "blend", first part. */

static void conc_blend_1 (void *ptr, Int32 n_conc, Int32 i_conc)
{
    IP_ImageId		srcimg, destimg;
    IP_FillModeType	oor;
    pyramid		*srcpyr;
    const xyinterp	*xyitp;
    const finterp	*fitp;
    Float32		mx, my, msize, xt, yt, xs, ys,
			size, fs, cdata[FF_NumImgChanTypes];
    Int32		i, xd, yd, ymin, ymax;

    /* Initialization */

    srcimg = ((blendparms *) ptr)->srcimg;
    destimg = ((blendparms *) ptr)->destimg;
    oor = ((blendparms *) ptr)->oor;
    srcpyr = ((blendparms *) ptr)->srcpyr;
    xyitp = &((blendparms *) ptr)->xyitp;
    fitp = &((blendparms *) ptr)->fitp;
    mx = ((blendparms *) ptr)->mx;
    my = ((blendparms *) ptr)->my;
    msize = ((blendparms *) ptr)->msize;
    n_conc = (n_conc == 0? 1: n_conc);
    ymin = (destimg->pimg.desc.height * i_conc) / n_conc;
    ymax = (destimg->pimg.desc.height * (i_conc + 1)) / n_conc;
    for (i = 0; i < FF_NumImgChanTypes; ++i) {
	cdata[i] = 0.0;
    }

    /* Main loop */

    for (yd = ymin; yd < ymax; ++yd) {
	for (xd = 0; xd < destimg->pimg.desc.width; ++xd) {
	    xt = xd * mx;
	    yt = yd * my;
	    sample_finterp (xt, yt, fitp, &fs);
	    sample_xyinterp (xt, yt, xyitp, &xs, &ys, &size);
	    xs *= srcimg->pimg.desc.width - 1;
	    ys *= srcimg->pimg.desc.height - 1;
	    size *= msize;
	    IP_SamplePyramid (srcpyr, xs, ys, size, cdata, oor);
	    for (i = 0; i < FF_NumImgChanTypes; ++i) {
		switch (destimg->pimg.channel[i]) {
		    case FF_ImgFmtTypeNone: {
			break;
		    }
		    case FF_ImgFmtTypeUByte: {
			*FF_ImgUBytePixel (&destimg->pimg, i, xd, yd) =
			    (UInt8) (255.0 * cdata[i] * fs);
			break;
		    }
		    case FF_ImgFmtTypeFloat: {
			*FF_ImgFloatPixel (&destimg->pimg, i, xd, yd) =
			    cdata[i] * fs;
			break;
		    }
		    default: {
			UT_ErrFatal ("conc_blend_1", "invalid channel type");
		    }
		}
	    }
	}
    }
    return;
}



/* Concurrent main loop for function "blend", second part. */

static void conc_blend_2 (void *ptr, Int32 n_conc, Int32 i_conc)
{
    IP_ImageId		srcimg, destimg;
    IP_FillModeType	oor;
    pyramid		*srcpyr;
    const xyinterp	*xyitp;
    const finterp	*fitp;
    Float32		mx, my, msize, xt, yt, xs, ys,
			size, fs, cdata[FF_NumImgChanTypes];
    Int32		i, xd, yd, ymin, ymax;

    /* Initialization */

    srcimg = ((blendparms *) ptr)->srcimg;
    destimg = ((blendparms *) ptr)->destimg;
    oor = ((blendparms *) ptr)->oor;
    srcpyr = ((blendparms *) ptr)->srcpyr;
    xyitp = &((blendparms *) ptr)->xyitp;
    fitp = &((blendparms *) ptr)->fitp;
    mx = ((blendparms *) ptr)->mx;
    my = ((blendparms *) ptr)->my;
    msize = ((blendparms *) ptr)->msize;
    n_conc = (n_conc == 0? 1: n_conc);
    ymin = (destimg->pimg.desc.height * i_conc) / n_conc;
    ymax = (destimg->pimg.desc.height * (i_conc + 1)) / n_conc;
    for (i = 0; i < FF_NumImgChanTypes; ++i) {
	cdata[i] = 0.0;
    }

    /* Main loop */

    for (yd = ymin; yd < ymax; ++yd) {
	for (xd = 0; xd < destimg->pimg.desc.width; ++xd) {
	    xt = xd * mx;
	    yt = yd * my;
	    sample_finterp (xt, yt, fitp, &fs);
	    fs = 1.0 - fs;
	    sample_xyinterp (xt, yt, xyitp, &xs, &ys, &size);
	    xs *= srcimg->pimg.desc.width - 1;
	    ys *= srcimg->pimg.desc.height - 1;
	    size *= msize;
	    IP_SamplePyramid (srcpyr, xs, ys, size, cdata, oor);
	    for (i = 0; i < FF_NumImgChanTypes; ++i) {
		switch (destimg->pimg.channel[i]) {
		    case FF_ImgFmtTypeNone: {
			break;
		    }
		    case FF_ImgFmtTypeUByte: {
			*FF_ImgUBytePixel (&destimg->pimg, i, xd, yd) +=
			    (UInt8) (255.0 * cdata[i] * fs);
			break;
		    }
		    case FF_ImgFmtTypeFloat: {
			*FF_ImgFloatPixel (&destimg->pimg, i, xd, yd) +=
			    cdata[i] * fs;
			break;
		    }
		    default: {
			UT_ErrFatal ("conc_blend_2", "invalid channel type");
		    }
		}
	    }
	}
    }
    return;
}



/* Compute pixel displacements and sampling area sizes for image warping
according to "dfunct" and "dderiv".  Store the result in an "xyinterp"
structure, "xyitp". */

static UT_Bool make_xyinterp
	(Int32 istep,
	 Int32 width, Int32 height,
	 UT_Bool (* dfunct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* dderiv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 xyinterp *xyitp)
{
    Int32		row, col, offset, ndata;
    Float32	mx, my, x, y, u, v, ux, uy, vx, vy, s;

    xyitp->numcol = width / istep;
    if (xyitp->numcol < 2)
	xyitp->numcol = 2;
    xyitp->numrow = height / istep;
    if (xyitp->numrow < 2)
	xyitp->numrow = 2;
    mx = 1.0 / (Float32) (xyitp->numcol - 1);
    my = 1.0 / (Float32) (xyitp->numrow - 1);

    ndata = xyitp->numcol * xyitp->numrow;
    if (!(xyitp->xdata = UT_MemTempArray (ndata, Float32)) ||
	!(xyitp->ydata = UT_MemTempArray (ndata, Float32)) ||
	!(xyitp->sdata = UT_MemTempArray (ndata, Float32)))
	return UT_False;

    for (row = 0; row < xyitp->numrow; ++row)
	for (col = 0; col < xyitp->numcol; ++col)
	{
	    offset = row * xyitp->numcol + col;
	    x = col * mx;
	    y = row * my;
	    if (!(* dfunct) (x, y, &u,  &v) ||
		!(* dderiv) (x, y, &ux, &uy, &vx, &vy))
		return UT_False;
	    if (ux < 0)
		ux = -ux;
	    if (uy < 0)
		uy = -uy;
	    if (vx < 0)
		vx = -vx;
	    if (vy < 0)
		vy = -vy;
	    if (ux > uy)
		s = ux;
	    else
		s = uy;
	    if (vx > s)
		s = vx;
	    if (vy > s)
		s = vy;
	    xyitp->xdata[offset] = u;
	    xyitp->ydata[offset] = v;
	    xyitp->sdata[offset] = s;
	}

    return UT_True;
}



/* Compute pixel weighting factors for image blending according to "bfunct".
Store the result in an "finterp" structure, "fitp". */

static UT_Bool make_finterp
	(Int32 istep, Int32 width, Int32 height,
	 UT_Bool (* bfunct) (Float32, Float32, Float32 *),
	 finterp *fitp)
{
    Int32		row, col, offset, ndata;
    Float32	mx, my, x, y, f;

    fitp->numcol = width / istep;
    if (fitp->numcol < 2)
	fitp->numcol = 2;
    fitp->numrow = height / istep;
    if (fitp->numrow < 2)
	fitp->numrow = 2;
    mx = 1.0 / (Float32) (fitp->numcol - 1);
    my = 1.0 / (Float32) (fitp->numrow - 1);

    ndata = fitp->numcol * fitp->numrow;
    if (!(fitp->fdata = UT_MemTempArray (ndata, Float32)))
	return UT_False;

    for (row = 0; row < fitp->numrow; ++row)
	for (col = 0; col < fitp->numcol; ++col)
	{
	    offset = row * fitp->numcol + col;
	    x = col * mx;
	    y = row * my;
	    if (!(* bfunct) (x, y, &f))
		return UT_False;
	    fitp->fdata[offset] = UT_MAX (0.0, UT_MIN (f, 1.0));
	}

    return UT_True;
}



/* Perform free-form deformations on an image. See the comments
on function IP_WarpFunct for an explanation of the parameters. */

static UT_Bool warp
	(IP_ImageId srcimg, IP_ImageId destimg,
	 UT_Bool (* dfunct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* dderiv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 Int32 istep, IP_FillModeType oor)
{
    warpparms parms;

    if (!(parms.srcpyr = IP_GenPyramid (srcimg)) ||
	!make_xyinterp (istep,
			srcimg->pimg.desc.width,
			srcimg->pimg.desc.height,
			dfunct, dderiv,
			&parms.xyitp))
    {
	return UT_False;
    }

    parms.mx =     
	(Float32) (srcimg->pimg.desc.width - 1) /
	(Float32) (destimg->pimg.desc.width - 1);
    parms.my =     
	(Float32) (srcimg->pimg.desc.height - 1) /
	(Float32) (destimg->pimg.desc.height - 1);
    parms.msize = UT_MAX (parms.mx, parms.my);
    parms.mx = 1.0 / (Float32) (destimg->pimg.desc.width - 1);
    parms.my = 1.0 / (Float32) (destimg->pimg.desc.height - 1);
    parms.srcimg = srcimg;
    parms.destimg = destimg;
    parms.oor = oor;

    return UT_THREAD_EXEC (conc_warp, &parms, IP_StateNumThreads, stacksize);
}



/* Blend two images. See the comments on function
IP_BlendFunct for an explanation of the parameters. */

static UT_Bool blend
	(IP_ImageId srcimg1,
	 IP_ImageId srcimg2,
	 IP_ImageId destimg,
	 UT_Bool (* d1funct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* d1deriv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 UT_Bool (* d2funct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* d2deriv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 UT_Bool (* bfunct) (Float32, Float32, Float32 *),
	 Int32 istep,
	 IP_FillModeType oor)
{
    UT_MemState	memstate;
    blendparms	parms;

    if (!make_finterp (istep,
		       srcimg1->pimg.desc.width,
		       srcimg1->pimg.desc.height,
		       bfunct,
		       &parms.fitp))
    {
	return UT_False;
    }

    memstate = UT_MemRemember ();

    if (!(parms.srcpyr = IP_GenPyramid (srcimg1)) ||
	!make_xyinterp (istep,
			srcimg1->pimg.desc.width,
			srcimg1->pimg.desc.height,
			d1funct, d1deriv,
			&parms.xyitp))
    {
	return UT_False;
    }

    parms.mx =
	(Float32) (srcimg1->pimg.desc.width - 1) /
	(Float32) (destimg->pimg.desc.width - 1);
    parms.my =
	(Float32) (srcimg1->pimg.desc.height - 1) /
	(Float32) (destimg->pimg.desc.height - 1);
    parms.msize = UT_MAX (parms.mx, parms.my);
    parms.mx = 1.0 / (Float32) (destimg->pimg.desc.width - 1);
    parms.my = 1.0 / (Float32) (destimg->pimg.desc.height - 1);
    parms.srcimg = srcimg1;
    parms.destimg = destimg;
    parms.oor = oor;

    if (!UT_THREAD_EXEC (conc_blend_1, &parms,
			      IP_StateNumThreads, stacksize))
    {
	return UT_False;
    }

    UT_MemRestore (memstate);

    if (!(parms.srcpyr = IP_GenPyramid (srcimg2)) ||
	!make_xyinterp (istep,
			srcimg2->pimg.desc.width,
			srcimg2->pimg.desc.height,
			d2funct, d2deriv,
			&parms.xyitp))
    {
	return UT_False;
    }

    parms.mx =
	(Float32) (srcimg2->pimg.desc.width - 1) /
	(Float32) (destimg->pimg.desc.width - 1);
    parms.my =
	(Float32) (srcimg2->pimg.desc.height - 1) /
	(Float32) (destimg->pimg.desc.height - 1);
    parms.msize = UT_MAX (parms.mx, parms.my);
    parms.mx = 1.0 / (Float32) (destimg->pimg.desc.width - 1);
    parms.my = 1.0 / (Float32) (destimg->pimg.desc.height - 1);
    parms.srcimg = srcimg2;

    return UT_THREAD_EXEC (conc_blend_2, &parms,
				IP_StateNumThreads, stacksize);
}



/***************************************************************************
 *[@e
 *	Name:		IP_WarpFunct
 *
 *	Usage:		Perform free-form deformations on an image.
 *
 *	Synopsis:	UT_Bool IP_WarpFunct
 *				(IP_ImageId srcimg,
 *				 IP_ImageId destimg,
 *				 UT_Bool (* dfunct)
 *					(Float32, Float32, Float32 *, Float32 *),
 *				 UT_Bool (* dderiv)
 *					(Float32, Float32,
 *					 Float32 *, Float32 *, Float32 *, Float32 *),
 *				 Int32 istep,
 *				 IP_FillModeType oor)
 *
 *	Description:	Image "srcimg" is deformed.  The result is stored
 *			in image "destimg". The deformation is defined by
 *			an application-supplied displacement function,
 *			"dfunct", and its derivatives, "dderiv", which
 *			should be defined as follows:
 *
 *			  UT_Bool dfunct
 *				(Float32 x, Float32 y,
 *				 Float32 *u,  Float32 *v)
 *
 *			  UT_Bool dderiv
 *				(Float32 x, Float32 y,
 *				 Float32 *ux, Float32 *uy,
 *				 Float32 *vx, Float32 *vy)
 *
 *			Given a position (x, y) in the destination image,
 *			"dfunct" computes the corresonding position (u, v)
 *			in the source image.
 *			"Dderiv" computes the four partial derivatives of
 *			"dfunct"; given a position (x, y), "dderiv" returns
 *
 *			    ux = du (x, y) / dx,
 *			    uy = du (x, y) / dy,
 *			    vx = dv (x, y) / dx,
 *			    vy = dv (x, y) / dy.
 *
 *			Evaluating "dfunct" and "dderiv" for every pixel can
 *			be quite time-consuming.  In order to save some time,
 *			"dfunct" and "dderiv" are sampled only every "istep"
 *			pixels; in-between values are interpolated linearly.
 *
 *			"Oor" selects how colors are assigned to output image
 *			pixels which are not covered by pixels in the input
 *			image. The following values for "oor" are accepted:
 *
 *			    IP_FillModeFill 	Pixels not covered by the input
 *						image are filled with the
 *						current drawing color.
 *
 *			    IP_FillModeWrap	The input image is repeated in
 *						horizontal and vertical
 *						direction so that the output
 *						image is tiled with copies of
 *						the input image.
 *
 *			    IP_FillModeClip	Input pixel coordinates are
 *						clipped so that regions outside
 *						the input image are filled with
 *						the color of the nearest input
 *						image pixel.
 *
 *			Notes:
 *
 *			- Right-handed coordinate systems are used for
 *			  positions in the source and destination images.
 *			  The lower left corner of "srcimg" and "destimg"
 *			  is at (0, 0); the upper right corner is at (1, 1).
 *
 *			- The return value of "dfunct" and "dderiv" should
 *			  normally be UT_True.  Image warping is aborted if
 *			  either "dfunct" or "dderiv" returns UT_False.
 *
 *			- Parts of IP_WarpFunct can be run concurrently on multiple
 *			  processors; see IP_SetNumThreads.
 *
 *	Return value:	UT_True if successful, else UT_False.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool IP_WarpFunct
	(IP_ImageId srcimg, IP_ImageId destimg,
	 UT_Bool (* dfunct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* dderiv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 Int32 istep, IP_FillModeType oor)
{
    UT_MemState	memstate;
    UT_Bool	success;

    if (istep < 1)
    {
	UT_ErrSetNum (UT_ErrParamInvalid, str_inotless, str_istep, 1);
	return UT_False;
    }
    if (oor < 0 || oor >= IP_NumFillModeTypes)
    {
	UT_ErrSetNum (UT_ErrParamInvalid, str_oormode, oor);
	return UT_False;
    }
    memstate = UT_MemRemember ();
    success = warp (srcimg, destimg, dfunct, dderiv, istep, oor);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *	Name:		IP_BlendFunct
 *
 *	Usage:		Blend two images, producing an "in-between" image.
 *
 *	Synopsis:	UT_Bool IP_BlendFunct
 *				(IP_ImageId srcimg1,
 *				 IP_ImageId srcimg2,
 *				 IP_ImageId destimg,
 *				 UT_Bool (* d1funct)
 *					(Float32, Float32, Float32 *, Float32 *),
 *				 UT_Bool (* d1deriv)
 *					(Float32, Float32,
 *					 Float32 *, Float32 *, Float32 *, Float32 *),
 *				 UT_Bool (* d2funct)
 *					(Float32, Float32, Float32 *, Float32 *),
 *				 UT_Bool (* d2deriv)
 *					(Float32, Float32,
 *					 Float32 *, Float32 *, Float32 *, Float32 *),
 *				 UT_Bool (* bfunct)
 *					(Float32, Float32, Float32 *),
 *				 Int32 istep,
 *				 IP_FillModeType oor)
 *
 *	Description:	IP_BlendFunct interpolates between two source images,
 *			"srcimg1" and "srcimg2", by displacing pixels and
 *			blending the colors of the displaced pixels. The
 *			result is stored in image "destimg". The interpola-
 *			tion is defined by a set of functions, which should
 *			be defined as follows:
 *
 *			  UT_Bool d1funct
 *				(Float32 x, Float32 y,
 *				 Float32 *u,  Float32 *v)
 *
 *			  UT_Bool d1deriv
 *				(Float32 x, Float32 y,
 *				 Float32 *ux, Float32 *uy,
 *				 Float32 *vx, Float32 *vy)
 *
 *			  UT_Bool d2funct
 *				(Float32 x, Float32 y,
 *				 Float32 *u,  Float32 *v)
 *
 *			  UT_Bool d2deriv
 *				(Float32 x, Float32 y,
 *				 Float32 *ux, Float32 *uy,
 *				 Float32 *vx, Float32 *vy)
 *
 *			  UT_Bool bfunct
 *				(Float32 x, Float32 y,
 *				 Float32 *f)
 *
 *			"D1funct" and "d1deriv" define how pixels from
 *			"srcimg1" are displaced:
 *			Given a position (x, y) in "destimg", "d1funct"
 *			finds the corresponding position (u, v) in "srcimg1".
 *			"D1deriv" computes the four partial derivatives of
 *			"d1funct"; given a position (x, y), "dderiv1" returns
 *
 *			    ux = du (x, y) / dx,
 *			    uy = du (x, y) / dy,
 *			    vx = dv (x, y) / dx,
 *			    vy = dv (x, y) / dy.
 *			    Dv = dv (x, y) / dy.
 *
 *			"D2funct" and "d2deriv" define how pixels from
 *			"srcimg2" are displaced, analogous to "d1funct"
 *			and "d1deriv".
 *
 *			"Bfunct" defines how the colors of the displaced pixels
 *			from the source images are combined into colors for
 *			the pixels in the destination image.  Given a position
 *			(x, y) in "destimg", "bfunct" computes a blending
 *			factor "f".   The color in "destimg" is "f" times the
 *			color from "srcimg1" plus (1-f) times the color from
 *			"srcimg2".
 *
 *			Evaluating the five interpolation functions for every
 *			pixel can be quite time-consuming.  In order to save
 *			time, the functions are sampled only every "istep"
 *			pixels; in-between values are interpolated linearly.
 *
 *			"Oor" selects how colors are assigned to output image
 *			pixels which are not covered by pixels in the input
 *			images.  The following values for "oor" are accepted:
 *
 *			    IP_FillModeFill 	Pixels not covered by the input
 *						images are filled with the
 *						current drawing color.
 *
 *			    IP_FillModeWrap	The input images are repeated
 *						in horizontal and vertical
 *						direction so that the output
 *						image is tiled with copies of
 *						the input images.
 *
 *			    IP_FillModeClip	Input pixel coordinates are
 *						clipped so that regions outside
 *						the input images are filled
 *						with the color of the nearest
 *						input image pixel.
 *
 *			Notes:
 *
 *			- Right-handed coordinate systems are used for
 *			  positions in the source and destination images.
 *			  The lower left corner of each image is at (0, 0);
 *			  the upper right corner is at (1, 1).
 *
 *			- Blending factors returned by "bfunct" should be
 *			  in the range from 0.0 to 1.0.  Values less than
 *			  0.0 are treated as 0.0; values above 1.0 are
 *			  treated as 1.0.
 *
 *			- The return value of all five interpolation functions
 *			  should normally be UT_True.  Image blending is aborted
 *			  as soon as one of the interpolation functions returns
 *			  UT_False.
 *
 *			- Parts of IP_BlendFunct can be run concurrently on
 *			  multiple processors; see IP_SetNumThreads.
 *
 *	Return value:	UT_True if successful, else UT_False.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool IP_BlendFunct
	(IP_ImageId srcimg1,
	 IP_ImageId srcimg2,
	 IP_ImageId destimg,
	 UT_Bool (* d1funct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* d1deriv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 UT_Bool (* d2funct) (Float32, Float32, Float32 *, Float32 *),
	 UT_Bool (* d2deriv) (Float32, Float32, Float32 *, Float32 *, Float32 *, Float32 *),
	 UT_Bool (* bfunct) (Float32, Float32, Float32 *),
	 Int32 istep,
	 IP_FillModeType oor)
{
    UT_MemState	memstate;
    UT_Bool	success;

    if (istep < 1)
    {
	UT_ErrSetNum (UT_ErrParamInvalid, str_inotless, str_istep, 1);
	return UT_False;
    }
    if (oor < 0 || oor >= IP_NumFillModeTypes)
    {
	UT_ErrSetNum (UT_ErrParamInvalid, str_oormode, oor);
	return UT_False;
    }
    memstate = UT_MemRemember ();
    success = blend (srcimg1, srcimg2, destimg, d1funct, d1deriv,
                     d2funct, d2deriv, bfunct, istep, oor);
    UT_MemRestore (memstate);
    return success;
}
