about summary refs log tree commit diff
path: root/editor/pamrubber.c
diff options
context:
space:
mode:
Diffstat (limited to 'editor/pamrubber.c')
-rw-r--r--editor/pamrubber.c1475
1 files changed, 1475 insertions, 0 deletions
diff --git a/editor/pamrubber.c b/editor/pamrubber.c
new file mode 100644
index 00000000..4378c340
--- /dev/null
+++ b/editor/pamrubber.c
@@ -0,0 +1,1475 @@
+/*----------------------------------------------------------------------------*/
+
+/* 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 : time(NULL));
+
+    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;
+}
+
+
+