path: root/editor/pamrubber.c

/*----------------------------------------------------------------------------*/

/* pamrubber.c - transform images using Rubber Sheeting algorithm
**               see: http://www.schaik.com/netpbm/rubber/
**
** Copyright (C) 2011 by Willem van Schaik (willem@schaik.com)
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

/*----------------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <time.h>

#include "pm_c_util.h"
#include "mallocvar.h"
#include "shhopt.h"
#include "pam.h"
#include "pamdraw.h"



typedef struct {
  double x;
  double y;
} point;

typedef struct {
  point p1;
  point p2;
} line;

typedef struct {
  point p1;
  point p2;
  point p3;
} triangle;

typedef struct {
    point tl;  /* top left     */
    point tr;  /* top right    */
    point bl;  /* bottom left  */
    point br;  /* bottom right */
} quadrilateral;

struct cmdlineInfo {
    unsigned int nCP;
    point        oldCP[4];
    point        newCP[4];
    const char * fileName;
    unsigned int quad;
    unsigned int tri;
    unsigned int frame;
    unsigned int linear;
    unsigned int verbose;
    unsigned int randseedSpec;
    unsigned int randseed;
};


static void
parseCmdline(int argc, const char ** argv,
             struct cmdlineInfo * const cmdlineP) {

/* parse all parameters from the command line */

    unsigned int option_def_index;
    char * endptr;
    unsigned int i;
    unsigned int nCP;

    /* instructions to optParseOptions3 on how to parse our options. */
    optEntry * option_def;
    optStruct3 opt;
    
    MALLOCARRAY_NOFAIL(option_def, 100);

    option_def_index = 0;   /* incremented by OPTENT3 */
    OPTENT3(0, "quad",     OPT_FLAG, NULL, &cmdlineP->quad,     0);
    OPTENT3(0, "tri",      OPT_FLAG, NULL, &cmdlineP->tri,      0);
    OPTENT3(0, "frame",    OPT_FLAG, NULL, &cmdlineP->frame,    0);
    OPTENT3(0, "linear",   OPT_FLAG, NULL, &cmdlineP->linear,   0);
    OPTENT3(0, "verbose",  OPT_FLAG, NULL, &cmdlineP->verbose,  0);
    OPTENT3(0, "randseed", OPT_UINT, &cmdlineP->randseed,
            &cmdlineP->randseedSpec, 0);
    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randseed,
            &cmdlineP->randseedSpec, 0);

    opt.opt_table = option_def;
    opt.short_allowed = FALSE;  /* we have no short (old-fashioned) options */
    opt.allowNegNum = FALSE;   /* we have no parms that are negative numbers */

    /* uses and sets argc, argv, and some of *cmdlineP and others. */
    pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);

    if (!cmdlineP->tri && !cmdlineP->quad)
        pm_error("You must specify either -tri or -quad");

    if (cmdlineP->tri && cmdlineP->quad)
        pm_error("You may not specify both -tri and -quad");

    /* Parameters are the control points (in qty of 4) and possibly a file name
     */
    nCP = (argc-1) / 4;

    if (nCP > 4)
        pm_error("Too many arguments: %u.  Arguments are "
                 "control point coordinates and an optional file name, "
                 "with a maximum of 4 control points", argc-1);

    cmdlineP->nCP = nCP;

    assert(nCP <= ARRAY_SIZE(cmdlineP->oldCP));
    assert(nCP <= ARRAY_SIZE(cmdlineP->newCP));

    for (i = 0; i < nCP; ++i) {
        cmdlineP->oldCP[i].x = strtol(argv[i * 2 + 1], &endptr, 10);
        cmdlineP->oldCP[i].y = strtol(argv[i * 2 + 2], &endptr, 10);
        cmdlineP->newCP[i].x = strtol(argv[4 * nCP / 2 + i * 2 + 1],
                                      &endptr, 10);
        cmdlineP->newCP[i].y = strtol(argv[4 * nCP / 2 + i * 2 + 2],
                                      &endptr, 10);
    }

    if (argc - 1 == 4 * nCP)
        cmdlineP->fileName = "-";
    else if (argc - 2 == 4 * nCP)
        cmdlineP->fileName = argv[nCP * 4 + 1];
    else
        pm_error("Invalid number of arguments.  Arguments are "
                 "control point coordinates and an optional file name, "
                 "so there must be a multiple of 4 or a multiple of 4 "
                 "plus 1.");
}



/* global variables */

static int nCP;
static point oldCP[4];
static point newCP[4];
static int nTri;
static triangle tri1s[10];
static triangle tri2s[10];
static quadrilateral quad1;
static quadrilateral quad2;
static tuple black;

static point
makepoint(double const x,
          double const y) {

    point retval;

    retval.x = x;
    retval.y = y;

    return retval;
}



static double
distance(point const p1,
         point const p2) {

    return sqrt(SQR(p1.x - p2.x) + SQR(p1.y - p2.y));
}



static line
makeline(point const p1,
         point const p2) {

    line retval;

    retval.p1 = p1;
    retval.p2 = p2;

    return retval;
}



static bool
intersect(line *  const l1P,
          line *  const l2P,
          point * const pP) {

    bool cross;

    if (((l2P->p2.y - l2P->p1.y) * (l1P->p2.x - l1P->p1.x) -
         (l2P->p2.x - l2P->p1.x) * (l1P->p2.y - l1P->p1.y)) == 0) {
        /* parallel lines */

        cross = false;

        if ((l1P->p1.x == l1P->p2.x) && (l2P->p1.x == l2P->p2.x)) {
            /* two vertical lines */
            pP->x = (l1P->p1.x + l2P->p1.x) / 2.0;
            pP->y = 1e10;
        } else if ((l1P->p1.y == l1P->p2.y) && (l2P->p1.y == l2P->p2.y)) {
            /* two horizontal lines */
            pP->x = 1e10;
            pP->y = (l1P->p1.y + l2P->p1.y) / 2.0;
        } else {
            if (fabs(l1P->p2.y - l1P->p1.y) > fabs(l1P->p2.x - l1P->p1.x)) {
                /* steep slope */
                pP->y = 1e10;
                pP->x = (l1P->p2.x - l1P->p1.x) / (l1P->p2.y - l1P->p1.y)
                    * 1e10;
            } else {
                /* even slope */
                pP->x = 1e10;
                pP->y = (l1P->p2.y - l1P->p1.y) / (l1P->p2.x - l1P->p1.x)
                    * 1e10;
            }
        }
    } else {
        /* intersecting lines */
        double const ua =
            ((l2P->p2.x - l2P->p1.x) * (l1P->p1.y - l2P->p1.y)
              - (l2P->p2.y - l2P->p1.y) * (l1P->p1.x - l2P->p1.x))
            / ((l2P->p2.y - l2P->p1.y)
               * (l1P->p2.x - l1P->p1.x) - (l2P->p2.x - l2P->p1.x)
               * (l1P->p2.y - l1P->p1.y));
        double const ub =
            ((l1P->p2.x - l1P->p1.x) * (l1P->p1.y - l2P->p1.y)
              - (l1P->p2.y - l1P->p1.y) * (l1P->p1.x - l2P->p1.x))
            / ((l2P->p2.y - l2P->p1.y)
               * (l1P->p2.x - l1P->p1.x) - (l2P->p2.x - l2P->p1.x)
               * (l1P->p2.y - l1P->p1.y));

        pP->x = l1P->p1.x + ua * (l1P->p2.x - l1P->p1.x);
        pP->y = l1P->p1.y + ua * (l1P->p2.y - l1P->p1.y);

        if ((ua >= 0.0) && (ua <= 1.0) && (ub >= 0.0) && (ub <= 1.0))
            cross = true;
        else
            cross = false;
    }

    return cross;
}



static triangle
maketriangle(point const p1,
             point const p2,
             point const p3) {

    triangle retval;

    retval.p1 = p1;
    retval.p2 = p2;
    retval.p3 = p3;
    
    return retval;
}



static int
insidetri(triangle * const triP,
          point      const p) {

    int cnt;

    cnt = 0;  /* initial value */

    if ((((triP->p1.y <= p.y) && (p.y < triP->p3.y))
         || ((triP->p3.y <= p.y) && (p.y < triP->p1.y)))
        &&
        (p.x < (triP->p3.x - triP->p1.x) * (p.y - triP->p1.y)
         / (triP->p3.y - triP->p1.y) + triP->p1.x))
        cnt = !cnt;

    if ((((triP->p2.y <= p.y) && (p.y < triP->p1.y))
         || ((triP->p1.y <= p.y) && (p.y < triP->p2.y)))
        &&
        (p.x < (triP->p1.x - triP->p2.x) * (p.y - triP->p2.y)
         / (triP->p1.y - triP->p2.y) + triP->p2.x))
        cnt = !cnt;

    if ((((triP->p3.y <= p.y) && (p.y < triP->p2.y))
         || ((triP->p2.y <= p.y) && (p.y < triP->p3.y)))
        &&
        (p.x < (triP->p2.x - triP->p3.x) * (p.y - triP->p3.y)
         / (triP->p2.y - triP->p3.y) + triP->p3.x))
        cnt = !cnt;

    return cnt;
}



static bool
windtriangle(triangle * const tP,
             point      const p1,
             point      const p2,
             point      const p3) {
    point f, c;
    line le, lv;
    bool cw;

    /* find cross of vertical through p3 and the edge p1-p2 */
    f.x = p3.x;
    f.y = -1.0;
    lv = makeline(p3, f);
    le = makeline(p1, p2);
    intersect(&le, &lv, &c);

    /* check for clockwise winding triangle */
    if (((p1.x > p2.x) && (p3.y < c.y)) ||
        ((p1.x < p2.x) && (p3.y > c.y))) {
        *tP = maketriangle(p1, p2, p3);
        cw = true;
    } else { /* p1/2/3 were counter clockwise */
        *tP = maketriangle(p1, p3, p2);
        cw = false;
    }
    return cw;
}



static double
tiny(void) {

    if (rand() % 2)
        return +1E-6 * (double) ((rand() % 90) + 9);
    else
        return -1E-6 * (double) ((rand() % 90) + 9);
}



static void
angle(point * const p1P,
      point * const p2P) {
/*----------------------------------------------------------------------------
   Move *p2P slightly if necessary to make sure the line (*p1P, *p2P)
   is not horizontal or vertical.
-----------------------------------------------------------------------------*/
    if (p1P->x == p2P->x) { /* vertical line */
        p2P->x += tiny();
    }
    if (p1P->y == p2P->y) { /* horizontal line */
        p2P->y += tiny();
    }
}



static void
sideTriangleVerticalEdge(unsigned int const n,
                         triangle *   const trig1P,
                         point        const p11,
                         point        const p12,
                         point        const p13,
                         point        const p14,
                         point        const r11,
                         point        const r12,
                         triangle *   const trig2P,
                         point        const p21,
                         point        const p22,
                         point        const p23,
                         point        const p24,
                         point        const r21,
                         point        const r22) {
                                   
    if (((n >= 4) && (r11.x < p11.x)
         && (p14.x < p13.x) && (p14.x < p12.x)
         && (p14.x < p11.x)) /* left edge */
        || 
        ((n >= 4) && (r11.x > p11.x)
         && (p14.x > p13.x) && (p14.x > p12.x)
         && (p14.x > p11.x))) /* right edge */ {
        *trig1P = maketriangle(r11, r12, p14);
        *trig2P = maketriangle(r21, r22, p24);
    } else if (((n >= 3) && (r11.x < p11.x) && (p13.x < p12.x)
                && (p13.x < p11.x)) /* left edge */
               || 
               ((n >= 3) && (r11.x > p11.x) && (p13.x > p12.x)
                && (p13.x > p11.x))) /* right edge */ {
        *trig1P = maketriangle(r11, r12, p13);
        *trig2P = maketriangle(r21, r22, p23);
    } else if (((n >= 2) && (r11.x < p11.x)
                && (p12.x < p11.x)) /* left edge */
               ||
               ((n >= 2) && (r11.x > p11.x)
                && (p12.x > p11.x))) /* right edge */ { 
        *trig1P = maketriangle(r11, r12, p12);
        *trig2P = maketriangle(r21, r22, p22);
    } else if (n >= 1) {
        *trig1P = maketriangle(r11, r12, p11);
        *trig2P = maketriangle(r21, r22, p21);
    }
}



static void
sideTriangleHorizontalEdge(unsigned int const n,
                           triangle *   const trig1P,
                           point        const p11,
                           point        const p12,
                           point        const p13,
                           point        const p14,
                           point        const r11,
                           point        const r12,
                           triangle *   const trig2P,
                           point        const p21,
                           point        const p22,
                           point        const p23,
                           point        const p24,
                           point        const r21,
                           point        const r22) {
                                   
    if (((n >= 4) && (r11.y < p11.y) && (p14.y < p13.y)
         && (p14.y < p12.y) && (p14.y < p11.y)) /* top edge */
        || 
        ((n >= 4) && (r11.y > p11.y) && (p14.y > p13.y)
         && (p14.y > p12.y) && (p14.y > p11.y))) /* bottom edge */ {
        *trig1P = maketriangle(r11, r12, p14);
        *trig2P = maketriangle(r21, r22, p24);
    } else if (((n >= 3) && (r11.y < p11.y) && (p13.y < p12.y)
                && (p13.y < p11.y)) /* top edge */
               || 
               ((n >= 3) && (r11.y > p11.y) && (p13.y > p12.y)
                && (p13.y > p11.y))) /* bottom edge */ {
        *trig1P = maketriangle(r11, r12, p13);
        *trig2P = maketriangle(r21, r22, p23);
    } else if (((n >= 2) && (r11.y < p11.y)
                && (p12.y < p11.y)) /* top edge */
               || 
               ((n >= 2) && (r11.y > p11.y)
                && (p12.y > p11.y))) /* bottom edge */ {
        *trig1P = maketriangle(r11, r12, p12);
        *trig2P = maketriangle(r21, r22, p22);
    } else if (n >= 1) {
        *trig1P = maketriangle(r11, r12, p11);
        *trig2P = maketriangle(r21, r22, p21);
    }
}



static void
sideTriangle(unsigned int const n,
             triangle *   const trig1P,
             point        const p11,
             point        const p12,
             point        const p13,
             point        const p14,
             point        const r11,
             point        const r12,
             triangle *   const trig2P,
             point        const p21,
             point        const p22,
             point        const p23,
             point        const p24,
             point        const r21,
             point        const r22) {

    if (fabs (r11.x - r12.x) < 1.0)
        sideTriangleVerticalEdge(n,
                                 trig1P, p11, p12, p13, p14, r11, r12,
                                 trig2P, p21, p22, p23, p24, r21, r22);

    else if (fabs (r11.y - r12.y) < 1.0)
        sideTriangleHorizontalEdge(n,
                                 trig1P, p11, p12, p13, p14, r11, r12,
                                 trig2P, p21, p22, p23, p24, r21, r22);
}



static void
edgeTriangle(triangle * const trig1P,
             point      const p11,
             point      const p12,
             point      const tl1,
             point      const tr1,
             point      const bl1,
             point      const br1,
             triangle * const trig2P,
             point      const p21,
             point      const p22,
             point      const tl2,
             point      const tr2,
             point      const bl2,
             point      const br2) {
             
    if ((p11.x < p12.x) && (p11.y < p12.y)) {
        /* up/left to down/right */
        *trig1P = maketriangle(tr1, p12, p11);
        *trig2P = maketriangle(tr2, p22, p21);
    } else if ((p11.x > p12.x) && (p11.y > p12.y)) {
        /* down/right to up/left */
        *trig1P = maketriangle(bl1, p12, p11);
        *trig2P = maketriangle(bl2, p22, p21);
    } else if ((p11.x < p12.x) && (p11.y > p12.y)) {
        /* down/left to up/right */
        *trig1P = maketriangle(tl1, p12, p11);
        *trig2P = maketriangle(tl2, p22, p21);
    } else if ((p11.x > p12.x) && (p11.y < p12.y)) {
        /* up/right to down/left */
        *trig1P = maketriangle(br1, p12, p11);
        *trig2P = maketriangle(br2, p22, p21);
    }
}



static quadrilateral
quadRect(double  const lft,
         double  const rgt,
         double  const top,
         double  const bot) {

    quadrilateral retval;

    retval.tl = makepoint(lft, top);
    retval.tr = makepoint(rgt, top);
    retval.bl = makepoint(lft, bot);
    retval.br = makepoint(rgt, bot);

    return retval;
}



static void
quadCornerSized(point           const p0,
                point           const p1,
                point           const p2,
                point           const p3,
                point           const mid,
                quadrilateral * const quadP,
                triangle *      const triP) {

/* P0-P1 and P2-P3 are the diagonals */
/* P0-P1 are further apart than P2-P3 */

    if ((p0.x < p1.x) && (p0.y < p1.y)) {
        /* p0 is top-left */
        quadP->tl = p0; quadP->br = p1;
        if (windtriangle(triP, p0, p2, p1)) {
            /* p2 is top-right */
            quadP->tr = p2; quadP->bl = p3;
        } else {
            /* p3 is top-right */
            quadP->tr = p3; quadP->bl = p2;
        }
    } else if ((p0.x > p1.x) && (p0.y < p1.y)) {
        /* p0 is top-right */
        quadP->tr = p0; quadP->bl = p1;
        if (windtriangle(triP, p0, p2, p1)) {
            /* p2 is bottom-right */
            quadP->br = p2; quadP->tl = p3;
        } else {
            /* p3 is bottom-right */
            quadP->br = p3; quadP->tl = p2;
        }
    } else if ((p0.x < p1.x) && (p0.y > p1.y)) {
        /* p0 is bottom-left */
        quadP->bl = p0; quadP->tr = p1;
        if (windtriangle(triP, p0, p2, p1)) {
            /* p2 is top-left */
            quadP->tl = p2; quadP->br = p3;
        } else {
            /* p3 is top-left */
            quadP->tl = p3; quadP->br = p2;
        }
    } else if ((p0.x > p1.x) && (p0.y > p1.y)) {
        /* p0 is bottom-right */
        quadP->br = p0; quadP->tl = p1;
        if (windtriangle(triP, p0, p2, p1)) {
            /* p2 is bottom-left */
            quadP->bl = p2; quadP->tr = p3;
        } else {
            /* p3 is bottom-left */
            quadP->bl = p3; quadP->tr = p2;
        }
    }
}



static void
quadCorner(point           const p0,
           point           const p1,
           point           const p2,
           point           const p3,
           point           const mid,
           quadrilateral * const quadP,
           triangle *      const triP) {

    /* p0-p1 and p2-p3 are the diagonals */

    if (fabs(p0.x - p1.x) + fabs(p0.y - p1.y) >=
        fabs(p2.x - p3.x) + fabs(p2.y - p3.y)) {
        quadCornerSized(p0, p1, p2, p3, mid, quadP, triP);
    } else {
        quadCornerSized(p2, p3, p0, p1, mid, quadP, triP);
    }
}



static pamd_drawproc frameDrawproc;

static void
frameDrawproc (tuple **     const tuples,
               unsigned int const cols,
               unsigned int const rows,
               unsigned int const depth,
               sample       const maxval,
               pamd_point   const p,
               const void * const clientdata) {
    
    int yy;

    for (yy = p.y - 1; yy <= p.y + 1; ++yy) {
        int xx;
        for (xx = p.x - 1; xx <= p.x + 1; ++xx) {
            if (xx >= 0 && xx < cols && yy >= 0 && yy < rows) {
                unsigned int i;
                for (i = 0; i < depth; ++i)
                    tuples[yy][xx][i] = (sample) *((tuple *) clientdata + i);
            }
        }
    }
}



static void
drawExtendedLine(const struct pam * const pamP,
                 tuple **           const outTuples,
                 point              const p1,
                 point              const p2) {

    pamd_point const p1ext =
        pamd_makePoint(p1.x - 10 * (p2.x - p1.x),
                       p1.y - 10 * (p2.y - p1.y));

    pamd_point const p2ext =
        pamd_makePoint(p2.x + 10 * (p2.x - p1.x),
                       p2.y + 10 * (p2.y - p1.y));

    pamd_line(outTuples, pamP->width, pamP->height, pamP->depth, pamP->maxval,
              p1ext, p2ext, frameDrawproc, black);
}



static pamd_point
clippedPoint(const struct pam * const pamP,
             point              const p) {

    int const roundedX = ROUND(p.x);
    int const roundedY = ROUND(p.x);

    int clippedX, clippedY;

    assert(pamP->width  >= 2);
    assert(pamP->height >= 2);

    if (roundedX <= 0)
        clippedX = 1;
    else if (roundedX > pamP->width - 1)
        clippedX = pamP->width - 2;
    else
        clippedX = roundedX;
        
    if (roundedY <= 0)
        clippedY = 1;
    else if (roundedY > pamP->width - 1)
        clippedY = pamP->width - 2;
    else
        clippedY = roundedY;
        
    return pamd_makePoint(clippedX, clippedY);
}



static void drawClippedTriangle(const struct pam * const pamP,
                                tuple **           const tuples,
                                triangle           const tri) {

    pamd_point const p1 = clippedPoint(pamP, tri.p1);
    pamd_point const p2 = clippedPoint(pamP, tri.p2);
    pamd_point const p3 = clippedPoint(pamP, tri.p3);

    pamd_line(tuples, pamP->width, pamP->height, pamP->depth, pamP->maxval,
              p1, p2, frameDrawproc, black);
    pamd_line(tuples, pamP->width, pamP->height, pamP->depth, pamP->maxval,
              p2, p3, frameDrawproc, black);
    pamd_line(tuples, pamP->width, pamP->height, pamP->depth, pamP->maxval,
              p3, p1, frameDrawproc, black);
}



static void
prepTrig(int const wd,
         int const ht) {

/* create triangles using control points */

    point rtl1, rtr1, rbl1, rbr1;
    point rtl2, rtr2, rbl2, rbr2;
    point c1p1, c1p2, c1p3, c1p4;
    point c2p1, c2p2, c2p3, c2p4;
    line l1, l2;
    point p0;

    rtl1 = makepoint(0.0 + tiny(),               0.0 + tiny());
    rtr1 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
    rbl1 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
    rbr1 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());

    rtl2 = makepoint(0.0 + tiny(),               0.0 + tiny());
    rtr2 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
    rbl2 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
    rbr2 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());

    if (nCP == 1) {
        c1p1 = oldCP[0];
        c2p1 = newCP[0];

        /* connect control point to all corners to get 4 triangles */
        /* left side triangle */
        sideTriangle(nCP,
                     &tri1s[0], c1p1, p0, p0, p0, rbl1, rtl1, 
                     &tri2s[0], c2p1, p0, p0, p0, rbl2, rtl2);
        /* top side triangle */
        sideTriangle(nCP,
                     &tri1s[1], c1p1, p0, p0, p0, rtl1, rtr1, 
                     &tri2s[1], c2p1, p0, p0, p0, rtl2, rtr2);
        /* right side triangle */
        sideTriangle(nCP,
                     &tri1s[2], c1p1, p0, p0, p0, rtr1, rbr1, 
                     &tri2s[2], c2p1, p0, p0, p0, rtr2, rbr2);
        /* bottom side triangle */
        sideTriangle(nCP,
                     &tri1s[3], c1p1, p0, p0, p0, rbr1, rbl1, 
                     &tri2s[3], c2p1, p0, p0, p0, rbr2, rbl2);

        nTri = 4;
    } else if (nCP == 2) {
        c1p1 = oldCP[0];
        c1p2 = oldCP[1];
        c2p1 = newCP[0];
        c2p2 = newCP[1];

        /* check for hor/ver edges */
        angle (&c1p1, &c1p2);
        angle (&c2p1, &c2p2);

        /* connect two control points to corners to get 6 triangles */
        /* left side */
        sideTriangle(nCP,
                     &tri1s[0], c1p1, c1p2, p0, p0, rbl1, rtl1, 
                     &tri2s[0], c2p1, c2p2, p0, p0, rbl2, rtl2);
        /* top side */
        sideTriangle(nCP, 
                     &tri1s[1], c1p1, c1p2, p0, p0, rtl1, rtr1, 
                     &tri2s[1], c2p1, c2p2, p0, p0, rtl2, rtr2);
        /* right side */
        sideTriangle(nCP, 
                     &tri1s[2], c1p1, c1p2, p0, p0, rtr1, rbr1, 
                     &tri2s[2], c2p1, c2p2, p0, p0, rtr2, rbr2);
        /* bottom side */
        sideTriangle(nCP, 
                     &tri1s[3], c1p1, c1p2, p0, p0, rbr1, rbl1, 
                     &tri2s[3], c2p1, c2p2, p0, p0, rbr2, rbl2);

        /* edge to corner triangles */
        edgeTriangle(&tri1s[4], c1p1, c1p2, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[4], c2p1, c2p2, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[5], c1p2, c1p1, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[5], c2p2, c2p1, rtl2, rtr2, rbl2, rbr2);
        nTri = 6;
    } else if (nCP == 3) {
        c1p1 = oldCP[0];
        c1p2 = oldCP[1];
        c1p3 = oldCP[2];
         
        c2p1 = newCP[0];
        c2p2 = newCP[1];
        c2p3 = newCP[2];

        /* Move vertices slightly if necessary to make sure no edge is
           horizontal or vertical.
        */
        angle(&c1p1, &c1p2);
        angle(&c1p2, &c1p3);
        angle(&c1p3, &c1p1);

        angle(&c2p1, &c2p2);
        angle(&c2p2, &c2p3);
        angle(&c2p3, &c2p1);

        if (windtriangle(&tri1s[0], c1p1, c1p2, c1p3)) {
            tri2s[0] = maketriangle(c2p1, c2p2, c2p3);
        } else {
            tri2s[0] = maketriangle(c2p1, c2p3, c2p2);
        }

        c1p1 = tri1s[0].p1;
        c1p2 = tri1s[0].p2;
        c1p3 = tri1s[0].p3;
         
        c2p1 = tri2s[0].p1;
        c2p2 = tri2s[0].p2;
        c2p3 = tri2s[0].p3;

        /* point to side triangles */
        /* left side */
        sideTriangle(nCP,
                     &tri1s[1], c1p1, c1p2, c1p3, p0, rbl1, rtl1, 
                     &tri2s[1], c2p1, c2p2, c2p3, p0, rbl2, rtl2);
        /* top side */
        sideTriangle(nCP, 
                     &tri1s[2], c1p1, c1p2, c1p3, p0, rtl1, rtr1, 
                     &tri2s[2], c2p1, c2p2, c2p3, p0, rtl2, rtr2);
        /* right side */
        sideTriangle(nCP, 
                     &tri1s[3], c1p1, c1p2, c1p3, p0, rtr1, rbr1, 
                     &tri2s[3], c2p1, c2p2, c2p3, p0, rtr2, rbr2);
        /* bottom side */
        sideTriangle(nCP, 
                     &tri1s[4], c1p1, c1p2, c1p3, p0, rbr1, rbl1, 
                     &tri2s[4], c2p1, c2p2, c2p3, p0, rbr2, rbl2);

        /* edge to corner triangles */
        edgeTriangle(&tri1s[5], c1p1, c1p2, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[5], c2p1, c2p2, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[6], c1p2, c1p3, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[6], c2p2, c2p3, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[7], c1p3, c1p1, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[7], c2p3, c2p1, rtl2, rtr2, rbl2, rbr2);
        nTri = 8;
    } else if (nCP == 4) {
        c1p1 = oldCP[0];
        c1p2 = oldCP[1];
        c1p3 = oldCP[2];
        c1p4 = oldCP[3];
         
        c2p1 = newCP[0];
        c2p2 = newCP[1];
        c2p3 = newCP[2];
        c2p4 = newCP[3];

        /* check for hor/ver edges */
        angle (&c1p1, &c1p2);
        angle (&c1p2, &c1p3);
        angle (&c1p3, &c1p4);
        angle (&c1p4, &c1p1);
        angle (&c1p1, &c1p3);
        angle (&c1p2, &c1p4);

        angle (&c2p1, &c2p2);
        angle (&c2p2, &c2p3);
        angle (&c2p3, &c2p4);
        angle (&c2p4, &c2p1);
        angle (&c2p1, &c2p3);
        angle (&c2p2, &c2p4);

        /*-------------------------------------------------------------------*/
        /*        -1-      -2-        -3-      -4-        -5-      -6-       */
        /*       1   2    1   3      1   2    1   4      1   3    1   4      */
        /*         X        X          X        X          X        X        */
        /*       3   4    2   4      4   3    2   3      4   2    3   2      */
        /*-------------------------------------------------------------------*/

        /* center two triangles */
        l1 = makeline(c1p1, c1p4);
        l2 = makeline(c1p2, c1p3);
        if (intersect(&l1, &l2, &p0)) {
            if (windtriangle(&tri1s[0], c1p1, c1p2, c1p3)) {
                tri1s[1] = maketriangle(c1p4, c1p3, c1p2);
                tri2s[0] = maketriangle(c2p1, c2p2, c2p3);
                tri2s[1] = maketriangle(c2p4, c2p3, c2p2);
            } else {
                tri1s[1] = maketriangle(c1p4, c1p2, c1p3);
                tri2s[0] = maketriangle(c2p1, c2p3, c2p2);
                tri2s[1] = maketriangle(c2p4, c2p2, c2p3);
            }
        }
        l1 = makeline(c1p1, c1p3);
        l2 = makeline(c1p2, c1p4);
        if (intersect(&l1, &l2, &p0)) {
            if (windtriangle(&tri1s[0], c1p1, c1p2, c1p4)) {
                tri1s[1] = maketriangle(c1p3, c1p4, c1p2);
                tri2s[0] = maketriangle(c2p1, c2p2, c2p4);
                tri2s[1] = maketriangle(c2p3, c2p4, c2p2);
            } else {
                tri1s[1] = maketriangle(c1p3, c1p2, c1p4);
                tri2s[0] = maketriangle(c2p1, c2p4, c2p2);
                tri2s[1] = maketriangle(c2p3, c2p2, c2p4);
            }
        }
        l1 = makeline(c1p1, c1p2);
        l2 = makeline(c1p3, c1p4);
        if (intersect(&l1, &l2, &p0)) {
            if (windtriangle(&tri1s[0], c1p1, c1p3, c1p4)) {
                tri1s[1] = maketriangle(c1p2, c1p4, c1p3);
                tri2s[0] = maketriangle(c2p1, c2p3, c2p4);
                tri2s[1] = maketriangle(c2p2, c2p4, c2p3);
            } else {
                tri1s[1] = maketriangle(c1p2, c1p3, c1p4);
                tri2s[0] = maketriangle(c2p1, c2p4, c2p3);
                tri2s[1] = maketriangle(c2p2, c2p3, c2p4);
            }
        }

        /* control points in correct orientation */
        c1p1 = tri1s[0].p1;
        c1p2 = tri1s[0].p2;
        c1p3 = tri1s[0].p3;
        c1p4 = tri1s[1].p1;
        c2p1 = tri2s[0].p1;
        c2p2 = tri2s[0].p2;
        c2p3 = tri2s[0].p3;
        c2p4 = tri2s[1].p1;

        /* triangle from triangle point to side of image */
        /* left side triangle */
        sideTriangle(nCP, 
                     &tri1s[2], c1p1, c1p2, c1p3, c1p4, rbl1, rtl1, 
                     &tri2s[2], c2p1, c2p2, c2p3, c2p4, rbl2, rtl2);
        /* top side triangle */
        sideTriangle(nCP, 
                     &tri1s[3], c1p1, c1p2, c1p3, c1p4, rtl1, rtr1, 
                     &tri2s[3], c2p1, c2p2, c2p3, c2p4, rtl2, rtr2);
        /* right side triangle */
        sideTriangle(nCP, 
                     &tri1s[4], c1p1, c1p2, c1p3, c1p4, rtr1, rbr1, 
                     &tri2s[4], c2p1, c2p2, c2p3, c2p4, rtr2, rbr2);
        /* bottom side triangle */
        sideTriangle(nCP, 
                     &tri1s[5], c1p1, c1p2, c1p3, c1p4, rbr1, rbl1, 
                     &tri2s[5], c2p1, c2p2, c2p3, c2p4, rbr2, rbl2);

        /*-------------------------------------------------------------------*/
        /*        -1-      -2-        -3-      -4-        -5-      -6-       */
        /*       1   2    1   3      1   2    1   4      1   3    1   4      */
        /*         X        X          X        X          X        X        */
        /*       3   4    2   4      4   3    2   3      4   2    3   2      */
        /*-------------------------------------------------------------------*/

        /* edge-corner triangles */
        edgeTriangle(&tri1s[6], c1p1, c1p2, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[6], c2p1, c2p2, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[7], c1p2, c1p4, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[7], c2p2, c2p4, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[8], c1p4, c1p3, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[8], c2p4, c2p3, rtl2, rtr2, rbl2, rbr2);
        edgeTriangle(&tri1s[9], c1p3, c1p1, rtl1, rtr1, rbl1, rbr1,
                     &tri2s[9], c2p3, c2p1, rtl2, rtr2, rbl2, rbr2);
        nTri = 10;
    }
}



static void
prepQuad(void) {

/* order quad control points */

    double d01, d12, d20;
    line l1, l2;
    point mid;
    triangle tri;

    if (nCP == 1) {
        /* create a rectangle from top-left corner of image and control
           point
        */
        quad1 = quadRect(0.0, oldCP[0].x, 0.0, oldCP[0].y);
        quad2 = quadRect(0.0, newCP[0].x, 0.0, newCP[0].y);
    } else if (nCP == 2) {
        /* create a rectangle with the two points as opposite corners */
        if ((oldCP[0].x < oldCP[1].x) && (oldCP[0].y < oldCP[1].y)) {
            /* top-left and bottom-right */
            quad1 = quadRect(oldCP[0].x,oldCP[1].x, oldCP[0].y, oldCP[1].y);
        } else if ((oldCP[0].x > oldCP[1].x) && (oldCP[0].y < oldCP[1].y)) {
            /* top-right and bottom-left */
            quad1 = quadRect(oldCP[1].x, oldCP[0].x, oldCP[0].y, oldCP[1].y);
        } else if ((oldCP[0].x < oldCP[1].x) && (oldCP[0].y < oldCP[1].y)) {
            /* bottom-left and top-right */
            quad1 = quadRect(oldCP[0].x, oldCP[1].x, oldCP[1].y, oldCP[0].y);
        } else if ((oldCP[0].x > oldCP[1].x) && (oldCP[0].y < oldCP[1].y)) {
            /* bottom-right and top-left */
            quad1 = quadRect(oldCP[1].x, oldCP[0].x, oldCP[1].y, oldCP[0].y);
        }
        
        if ((newCP[0].x < newCP[1].x) && (newCP[0].y < newCP[1].y)) {
            /* top-left and bottom-right */
            quad2 = quadRect(newCP[0].x, newCP[1].x, newCP[0].y, newCP[1].y);
        } else if ((newCP[0].x > newCP[1].x) && (newCP[0].y < newCP[1].y)) {
            /* top-right and bottom-left */
            quad2 = quadRect(newCP[1].x, newCP[0].x, newCP[0].y, newCP[1].y);
        } else if ((newCP[0].x < newCP[1].x) && (newCP[0].y < newCP[1].y)) {
            /* bottom-left and top-right */
            quad2 = quadRect(newCP[0].x, newCP[1].x, newCP[1].y, newCP[0].y);
        } else if ((newCP[0].x > newCP[1].x) && (newCP[0].y < newCP[1].y)) {
            /* bottom-right and top-left */
            quad2 = quadRect(newCP[1].x, newCP[0].x, newCP[1].y, newCP[0].y);
        }
    } else {
        if (nCP == 3) {
            /* add the fourth corner of a parallelogram */
            /* diagonal of the parallelogram is the two control points
               furthest apart
            */
            
            d01 = distance(newCP[0], newCP[1]);
            d12 = distance(newCP[1], newCP[2]);
            d20 = distance(newCP[2], newCP[0]);

            if ((d01 > d12) && (d01 > d20)) {
                oldCP[3].x = oldCP[0].x + oldCP[1].x - oldCP[2].x;
                oldCP[3].y = oldCP[0].y + oldCP[1].y - oldCP[2].y;
            } else
                if (d12 > d20) {
                    oldCP[3].x = oldCP[1].x + oldCP[2].x - oldCP[0].x;
                    oldCP[3].y = oldCP[1].y + oldCP[2].y - oldCP[0].y;
                } else {
                    oldCP[3].x = oldCP[2].x + oldCP[0].x - oldCP[1].x;
                    oldCP[3].y = oldCP[2].y + oldCP[0].y - oldCP[1].y;
                }

            if ((d01 > d12) && (d01 > d20)) {
                newCP[3].x = newCP[0].x + newCP[1].x - newCP[2].x;
                newCP[3].y = newCP[0].y + newCP[1].y - newCP[2].y;
            } else
                if (d12 > d20) {
                    newCP[3].x = newCP[1].x + newCP[2].x - newCP[0].x;
                    newCP[3].y = newCP[1].y + newCP[2].y - newCP[0].y;
                } else {
                    newCP[3].x = newCP[2].x + newCP[0].x - newCP[1].x;
                    newCP[3].y = newCP[2].y + newCP[0].y - newCP[1].y;
                }

        } /* end nCP = 3 */

        /* nCP = 3 or 4 */

        /* find the intersection of the diagonals */
        l1 = makeline(oldCP[0], oldCP[1]);
        l2 = makeline(oldCP[2], oldCP[3]);
        if (intersect(&l1, &l2, &mid)) {
            quadCorner(oldCP[0], oldCP[1], oldCP[2], oldCP[3],
                       mid, &quad1, &tri);
        } else {
            l1 = makeline(oldCP[0], oldCP[2]);
            l2 = makeline(oldCP[1], oldCP[3]);
            if (intersect(&l1, &l2, &mid))
                quadCorner(oldCP[0], oldCP[2], oldCP[1], oldCP[3],
                           mid, &quad1, &tri);
            else {
                l1 = makeline(oldCP[0], oldCP[3]);
                l2 = makeline(oldCP[1], oldCP[2]);
                if (intersect(&l1, &l2, &mid))
                    quadCorner(oldCP[0], oldCP[3],
                               oldCP[1], oldCP[2], mid, &quad1, &tri);
                else
                    pm_error("The four old control points don't seem "
                             "to be corners.");
            }
        }

        /* repeat for the "to-be" control points */
        l1 = makeline(newCP[0], newCP[1]);
        l2 = makeline(newCP[2], newCP[3]);
        if (intersect(&l1, &l2, &mid))
            quadCorner(newCP[0], newCP[1], newCP[2], newCP[3],
                       mid, &quad2, &tri);
        else {
            l1 = makeline(newCP[0], newCP[2]);
            l2 = makeline(newCP[1], newCP[3]);
            if (intersect(&l1, &l2, &mid))
                quadCorner(newCP[0], newCP[2], newCP[1], newCP[3],
                           mid, &quad2, &tri);
            else {
                l1 = makeline(newCP[0], newCP[3]);
                l2 = makeline(newCP[1], newCP[2]);
                if (intersect(&l1, &l2, &mid))
                    quadCorner(newCP[0], newCP[3],
                               newCP[1], newCP[2], mid, &quad2, &tri);
                else
                    pm_error("The four new control points don't seem "
                             "to be corners.");
            }
        }
    }
}



static void
warpTrig(point   const p2,
         point * const p1P) {

/* map target to source by triangulation */

    point e1p1, e1p2, e1p3;
    point e2p1, e2p2, e2p3;
    line lp, le;
    line l1, l2, l3;
    double d1, d2, d3;
    int i;

    /* find in which triangle p2 lies */
    for (i = 0; i < nTri; i++) {
        if (insidetri (&tri2s[i], p2))
            break;
    }

    if (i == nTri)
        *p1P = makepoint(0.0, 0.0);
    else {
        /* where in triangle is point */
        d1 = fabs (p2.x - tri2s[i].p1.x) + fabs (p2.y - tri2s[i].p1.y);
        d2 = fabs (p2.x - tri2s[i].p2.x) + fabs (p2.y - tri2s[i].p2.y);
        d3 = fabs (p2.x - tri2s[i].p3.x) + fabs (p2.y - tri2s[i].p3.y);

        /* line through p1 and p intersecting with edge p2-p3 */
        lp = makeline(tri2s[i].p1, p2);
        le = makeline(tri2s[i].p2, tri2s[i].p3);
        intersect (&lp, &le, &e2p1);

        /* line through p2 and p intersecting with edge p3-p1 */
        lp = makeline(tri2s[i].p2, p2);
        le = makeline(tri2s[i].p3, tri2s[i].p1);
        intersect (&lp, &le, &e2p2);

        /* line through p3 and p intersecting with edge p1-p2 */
        lp = makeline(tri2s[i].p3, p2);
        le = makeline(tri2s[i].p1, tri2s[i].p2);
        intersect (&lp, &le, &e2p3);

        /* map target control points to source control points */
        e1p1.x = tri1s[i].p2.x
            + (e2p1.x - tri2s[i].p2.x)/(tri2s[i].p3.x - tri2s[i].p2.x) 
            * (tri1s[i].p3.x - tri1s[i].p2.x);
        e1p1.y = tri1s[i].p2.y
            + (e2p1.y - tri2s[i].p2.y)/(tri2s[i].p3.y - tri2s[i].p2.y)
            * (tri1s[i].p3.y - tri1s[i].p2.y);
        e1p2.x = tri1s[i].p3.x
            + (e2p2.x - tri2s[i].p3.x)/(tri2s[i].p1.x - tri2s[i].p3.x)
            * (tri1s[i].p1.x - tri1s[i].p3.x);
        e1p2.y = tri1s[i].p3.y
            + (e2p2.y - tri2s[i].p3.y)/(tri2s[i].p1.y - tri2s[i].p3.y)
            * (tri1s[i].p1.y - tri1s[i].p3.y);
        e1p3.x = tri1s[i].p1.x
            + (e2p3.x - tri2s[i].p1.x)/(tri2s[i].p2.x - tri2s[i].p1.x)
            * (tri1s[i].p2.x - tri1s[i].p1.x);
        e1p3.y = tri1s[i].p1.y
            + (e2p3.y - tri2s[i].p1.y)/(tri2s[i].p2.y - tri2s[i].p1.y)
            * (tri1s[i].p2.y - tri1s[i].p1.y);

        /* intersect grid lines in source */
        l1 = makeline(tri1s[i].p1, e1p1);
        l2 = makeline(tri1s[i].p2, e1p2);
        l3 = makeline(tri1s[i].p3, e1p3);

        if ((d1 < d2) && (d1 < d3))
            intersect (&l2, &l3, p1P);
        else if (d2 < d3)
            intersect (&l1, &l3, p1P);
        else
            intersect (&l1, &l2, p1P);
    }
}



static void
warpQuad(point   const p2,
         point * const p1P) {

/* map target to source for quad control points */

    point h2, v2;
    point c1tl, c1tr, c1bl, c1br;
    point c2tl, c2tr, c2bl, c2br;
    point e1t, e1b, e1l, e1r;
    point e2t, e2b, e2l, e2r;
    line l2t, l2b, l2l, l2r;
    line lh, lv;

    c1tl = quad1.tl;
    c1tr = quad1.tr;
    c1bl = quad1.bl;
    c1br = quad1.br;
       
    c2tl = quad2.tl;
    c2tr = quad2.tr;
    c2bl = quad2.bl;
    c2br = quad2.br;

    l2t = makeline(c2tl, c2tr);
    l2b = makeline(c2bl, c2br);
    l2l = makeline(c2tl, c2bl);
    l2r = makeline(c2tr, c2br);

    /* find intersections of lines through control points */
    intersect (&l2t, &l2b, &h2);
    intersect (&l2l, &l2r, &v2);

    /* find intersections of axes through P with control point box */
    lv = makeline(p2, v2);
    intersect (&l2t, &lv, &e2t);
    intersect (&l2b, &lv, &e2b);

    lh = makeline(p2, h2);
    intersect (&l2l, &lh, &e2l);
    intersect (&l2r, &lh, &e2r);

    /* map target control points to source control points */
    e1t.x = c1tl.x + (e2t.x - c2tl.x)/(c2tr.x - c2tl.x) * (c1tr.x - c1tl.x);
    if (c1tl.y == c1tr.y)
        e1t.y = c1tl.y;
    else
        e1t.y =
            c1tl.y + (e2t.x - c2tl.x)/(c2tr.x - c2tl.x) * (c1tr.y - c1tl.y);

    e1b.x = c1bl.x + (e2b.x - c2bl.x)/(c2br.x - c2bl.x) * (c1br.x - c1bl.x);
    if (c1bl.y == c1br.y)
        e1b.y = c1bl.y;
    else
        e1b.y =
            c1bl.y + (e2b.x - c2bl.x)/(c2br.x - c2bl.x) * (c1br.y - c1bl.y);

    if (c1tl.x == c1bl.x)
        e1l.x = c1tl.x;
    else
        e1l.x =
            c1tl.x + (e2l.y - c2tl.y)/(c2bl.y - c2tl.y) * (c1bl.x - c1tl.x);
    e1l.y = c1tl.y + (e2l.y - c2tl.y)/(c2bl.y - c2tl.y) * (c1bl.y - c1tl.y);

    if (c1tr.x == c1br.x)
        e1r.x = c1tr.x;
    else
        e1r.x
            = c1tr.x + (e2r.y - c2tr.y)/(c2br.y - c2tr.y) * (c1br.x - c1tr.x);
    e1r.y = c1tr.y + (e2r.y - c2tr.y)/(c2br.y - c2tr.y) * (c1br.y - c1tr.y);

    /* intersect grid lines in source */
    lv = makeline(e1t, e1b);
    lh = makeline(e1l, e1r);
    intersect (&lh, &lv, p1P);
}



static void
setGlobalCP(struct cmdlineInfo const cmdline) {

    unsigned int i;

    for (i = 0; i < cmdline.nCP; ++i) {
        oldCP[i] = cmdline.oldCP[i];
        newCP[i] = cmdline.newCP[i];
    }
    nCP = cmdline.nCP;
    for (i = nCP; i < ARRAY_SIZE(oldCP); ++i)
        oldCP[i] = makepoint(-1.0, -1.0);
    for (i = nCP; i < ARRAY_SIZE(newCP); ++i)
        newCP[i] = makepoint(-1.0, -1.0);
}



static void
createWhiteTuple(const struct pam * const pamP,
                 tuple *            const tupleP) {

    tuple white;
    unsigned int plane;

    white = pnm_allocpamtuple(pamP);

    for (plane = 0; plane < pamP->depth; ++plane)
        white[plane] = pamP->maxval;

    *tupleP = white;
}



static void
makeAllWhite(struct pam * const pamP,
             tuple **     const tuples) {

    tuple white;
    unsigned int row;

    createWhiteTuple(pamP, &white);

    for (row = 0; row < pamP->height; ++row)  {
        unsigned int col;
        for (col = 0; col < pamP->width; ++col)
            pnm_assigntuple(pamP, tuples[row][col], white);
    }

    pnm_freepamtuple(white);
}



static sample
pix(tuple **     const tuples,
    unsigned int const width,
    unsigned int const height,
    point        const p1,
    unsigned int const plane,
    bool         const linear) {

    point p;
    double pix;

    p.x = p1.x + 1E-3;
    p.y = p1.y + 1E-3;

    if (p.x < 0.0)
        p.x = 0.0;
    if (p.x > (double) width - 1.0)
        p.x = (double) width - 1.0;
    if (p.y < 0.0)
        p.y = 0.0;
    if (p.y > (double) height - 1.0)
        p.y = (double) height - 1.0;

    if (!linear) {
        pix = tuples
            [(int) floor(p.y + 0.5)]
            [(int) floor(p.x + 0.5)][plane];
    } else {
        double const rx = p.x - floor(p.x);
        double const ry = p.y - floor(p.y);
        pix = 0.0;
        if (((int) floor(p.x) >= 0) && ((int) floor(p.x) < width) &&
            ((int) floor(p.y) >= 0) && ((int) floor(p.y) < height)) {
            pix += (1.0 - rx) * (1.0 - ry) 
                * tuples[(int) floor(p.y)][(int) floor(p.x)][plane];
        }
        if (((int) floor(p.x) + 1 >= 0) && ((int) floor(p.x) + 1 < width) &&
            ((int) floor(p.y) >= 0) && ((int) floor(p.y) < height)) {
            pix += rx * (1.0 - ry) 
                * tuples[(int) floor(p.y)][(int) floor(p.x) + 1][plane];
        }
        if (((int) floor(p.x) >= 0) && ((int) floor(p.x) < width) &&
            ((int) floor(p.y) + 1 >= 0) && ((int) floor(p.y) + 1 < height)) {
            pix += (1.0 - rx) * ry 
                * tuples[(int) floor(p.y) + 1][(int) floor(p.x)][plane];
        }
        if (((int) floor(p.x) + 1 >= 0) && ((int) floor(p.x) + 1 < width) &&
            ((int) floor(p.y) + 1 >= 0) && ((int) floor(p.y) + 1 < height)) {
            pix += rx * ry 
                * tuples[(int) floor(p.y) + 1][(int) floor(p.x) + 1][plane];
        }
    }

    return (int) floor(pix);
}



int
main(int argc, const char ** const argv) {

    struct cmdlineInfo cmdline;
    FILE * ifP;
    struct pam inpam, outpam;
    tuple ** inTuples;
    tuple ** outTuples;
    unsigned int p2y;
  
    pm_proginit(&argc, argv);

    parseCmdline(argc, argv, &cmdline);

    setGlobalCP(cmdline);

    srand(cmdline.randseedSpec ? cmdline.randseed : pm_randseed());

    ifP = pm_openr(cmdline.fileName);

    inTuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));

    outpam = inpam;  /* initial value */
    outpam.file = stdout;

    outTuples = pnm_allocpamarray(&outpam);

    pnm_createBlackTuple(&outpam, &black);

    makeAllWhite(&outpam, outTuples);

    if (cmdline.tri)
        prepTrig(inpam.width, inpam.height);
    if (cmdline.quad)
        prepQuad();

    for (p2y = 0; p2y < inpam.height; ++p2y) {
        unsigned int p2x;
        for (p2x = 0; p2x < inpam.width; ++p2x) {
            point p1, p2;
            unsigned int plane;

            p2 = makepoint(p2x, p2y);
            if (cmdline.quad)
                warpQuad(p2, &p1);
            if (cmdline.tri)
                warpTrig(p2, &p1);

            for (plane = 0; plane < inpam.depth; ++plane) {
                outTuples[p2y][p2x][plane] =
                    pix(inTuples, inpam.width, inpam.height, p1, plane,
                        cmdline.linear);
            }
        }
    }

    if (cmdline.frame) {
        if (cmdline.quad) {
            drawExtendedLine(&outpam, outTuples, quad2.tl, quad2.tr);
            drawExtendedLine(&outpam, outTuples, quad2.bl, quad2.br);
            drawExtendedLine(&outpam, outTuples, quad2.tl, quad2.bl);
            drawExtendedLine(&outpam, outTuples, quad2.tr, quad2.br);
        }
        if (cmdline.tri) {
            unsigned int i;
            for (i = 0; i < nTri; ++i)
                drawClippedTriangle(&outpam, outTuples, tri2s[i]);
        }
    }

    pnm_writepam(&outpam, outTuples);

    pnm_freepamtuple(black);
    pnm_freepamarray(outTuples, &outpam);
    pnm_freepamarray(inTuples, &inpam);

    pm_close(ifP);
    pm_close(stdout);

    return 0;
}