/**************************************************************************
 *{@C
 *      Copyright:      1988-2023 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;
}
