path: root/editor/pamrubber.c

/*=============================================================================
                              pamrubber
===============================================================================
  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 "rand.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 randomseedSpec;
    unsigned int randomseed;
};


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->randomseed,
            &cmdlineP->randomseedSpec, 0);
    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randomseed,
            &cmdlineP->randomseedSpec, 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 void
findIntersection(const line *  const l1P,
                 const line *  const l2P,
                 bool *        const theyIntersectP,
                 point *       const intersectionP) {

    bool theyIntersect;

    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 */

        theyIntersect = false;

        if ((l1P->p1.x == l1P->p2.x) && (l2P->p1.x == l2P->p2.x)) {
            /* two vertical lines */
            intersectionP->x = (l1P->p1.x + l2P->p1.x) / 2.0;
            intersectionP->y = 1e10;
        } else if ((l1P->p1.y == l1P->p2.y) && (l2P->p1.y == l2P->p2.y)) {
            /* two horizontal lines */
            intersectionP->x = 1e10;
            intersectionP->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 */
                intersectionP->y = 1e10;
                intersectionP->x =
                    (l1P->p2.x - l1P->p1.x) / (l1P->p2.y - l1P->p1.y) * 1e10;
            } else {
                /* even slope */
                intersectionP->x = 1e10;
                intersectionP->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));

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

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



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 triP,
             point      const p1,
             point      const p2,
             point      const p3) {
    point f, c;
    bool cw;

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

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



static double
tiny(struct pm_randSt * const randStP) {

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



static void
angle(point * const p1P,
      point * const p2P,
      struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
   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(randStP);
    }
    if (p1P->y == p2P->y) { /* horizontal line */
        p2P->y += tiny(randStP);
    }
}



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 */

    triangle tri;

    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, &tri);
    } else {
        quadCornerSized(p2, p3, p0, p1, mid, quadP, &tri);
    }
    if (triP)
        *triP = tri;
}



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.y);

    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->height - 1)
        clippedY = pamP->height - 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,
         bool         const randomseedSpec,
         unsigned int const randomseed) {

/* 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;
    point p0;
    struct pm_randSt randSt;

    pm_randinit(&randSt);
    pm_srand2(&randSt, randomseedSpec, randomseed);

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

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

    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, &randSt);
        angle (&c2p1, &c2p2, &randSt);

        /* 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, &randSt);
        angle(&c1p2, &c1p3, &randSt);
        angle(&c1p3, &c1p1, &randSt);

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

        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, &randSt);
        angle (&c1p2, &c1p3, &randSt);
        angle (&c1p3, &c1p4, &randSt);
        angle (&c1p4, &c1p1, &randSt);
        angle (&c1p1, &c1p3, &randSt);
        angle (&c1p2, &c1p4, &randSt);

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

        /*-------------------------------------------------------------------*/
        /*        -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 */
            line const l1 = makeline(c1p1, c1p4);
            line const l2 = makeline(c1p2, c1p3);
            bool theyIntersect;
            findIntersection(&l1, &l2, &theyIntersect, &p0);
            if (theyIntersect) {
                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);
                }
            }
        }
        {
            line const l1 = makeline(c1p1, c1p3);
            line const l2 = makeline(c1p2, c1p4);
            bool theyIntersect;
            findIntersection(&l1, &l2, &theyIntersect, &p0);
            if (theyIntersect) {
                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);
                }
            }
        }
        {
            line const l1 = makeline(c1p1, c1p2);
            line const l2 = makeline(c1p3, c1p4);
            bool theyIntersect;
            findIntersection(&l1, &l2, &theyIntersect, &p0);
            if (theyIntersect) {
                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;
    }

    pm_randterm(&randSt);
}



static void
prepQuad(void) {

/* order quad control points */

    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
            */

            double const d01 = distance(newCP[0], newCP[1]);
            double const d12 = distance(newCP[1], newCP[2]);
            double const 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 */
            line const l1 = makeline(oldCP[0], oldCP[1]);
            line const l2 = makeline(oldCP[2], oldCP[3]);
            bool theyIntersect;
            point mid;
            findIntersection(&l1, &l2, &theyIntersect, &mid);
            if (theyIntersect)
                quadCorner(oldCP[0], oldCP[1], oldCP[2], oldCP[3],
                           mid, &quad1, NULL);
            else {
                line const l1 = makeline(oldCP[0], oldCP[2]);
                line const l2 = makeline(oldCP[1], oldCP[3]);
                bool theyIntersect;
                point mid;
                findIntersection(&l1, &l2, &theyIntersect, &mid);
                if (theyIntersect)
                    quadCorner(oldCP[0], oldCP[2], oldCP[1], oldCP[3],
                               mid, &quad1, NULL);
                else {
                    line const l1 = makeline(oldCP[0], oldCP[3]);
                    line const l2 = makeline(oldCP[1], oldCP[2]);
                    bool theyIntersect;
                    point mid;
                    findIntersection(&l1, &l2, &theyIntersect, &mid);
                    if (theyIntersect)
                        quadCorner(oldCP[0], oldCP[3],
                                   oldCP[1], oldCP[2], mid, &quad1, NULL);
                    else
                        pm_error("The four old control points don't seem "
                                 "to be corners.");
                }
            }
        }
        {
            /* repeat for the "to-be" control points */
            line const l1 = makeline(newCP[0], newCP[1]);
            line const l2 = makeline(newCP[2], newCP[3]);
            bool theyIntersect;
            point mid;
            findIntersection(&l1, &l2, &theyIntersect, &mid);
            if (theyIntersect)
                quadCorner(newCP[0], newCP[1], newCP[2], newCP[3],
                           mid, &quad2, NULL);
            else {
                line const l1 = makeline(newCP[0], newCP[2]);
                line const l2 = makeline(newCP[1], newCP[3]);
                bool theyIntersect;
                point mid;
                findIntersection(&l1, &l2, &theyIntersect, &mid);
                if (theyIntersect)
                    quadCorner(newCP[0], newCP[2], newCP[1], newCP[3],
                               mid, &quad2, NULL);
                else {
                    line const l1 = makeline(newCP[0], newCP[3]);
                    line const l2 = makeline(newCP[1], newCP[2]);
                    bool theyIntersect;
                    point mid;
                    findIntersection(&l1, &l2, &theyIntersect, &mid);
                    if (theyIntersect)
                        quadCorner(newCP[0], newCP[3],
                                   newCP[1], newCP[2], mid, &quad2, NULL);
                    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;
    unsigned 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 */
        double const d1 =
            fabs (p2.x - tri2s[i].p1.x) + fabs (p2.y - tri2s[i].p1.y);
        double const d2 =
            fabs (p2.x - tri2s[i].p2.x) + fabs (p2.y - tri2s[i].p2.y);
        double const 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 */
            line const lp = makeline(tri2s[i].p1, p2);
            line const le = makeline(tri2s[i].p2, tri2s[i].p3);
            findIntersection(&lp, &le, NULL, &e2p1);
        }

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

        {
            /* line through p3 and p intersecting with edge p1-p2 */
            line const lp = makeline(tri2s[i].p3, p2);
            line const le = makeline(tri2s[i].p1, tri2s[i].p2);
            findIntersection(&lp, &le, NULL, &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 */
            line const l1 = makeline(tri1s[i].p1, e1p1);
            line const l2 = makeline(tri1s[i].p2, e1p2);
            line const l3 = makeline(tri1s[i].p3, e1p3);

            if ((d1 < d2) && (d1 < d3))
                findIntersection(&l2, &l3, NULL, p1P);
            else if (d2 < d3)
                findIntersection(&l1, &l3, NULL, p1P);
            else
                findIntersection(&l1, &l2, NULL, 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;

    c1tl = quad1.tl;
    c1tr = quad1.tr;
    c1bl = quad1.bl;
    c1br = quad1.br;

    c2tl = quad2.tl;
    c2tr = quad2.tr;
    c2bl = quad2.bl;
    c2br = quad2.br;

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

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

        {
            /* find intersections of axes through P with control point box */
            line const lv = makeline(p2, v2);
            line const lh = makeline(p2, h2);

            findIntersection(&l2t, &lv, NULL, &e2t);
            findIntersection(&l2b, &lv, NULL, &e2b);
            findIntersection(&l2l, &lh, NULL, &e2l);
            findIntersection(&l2r, &lh, NULL, &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 */
        line const lv = makeline(e1t, e1b);
        line const lh = makeline(e1l, e1r);
        findIntersection(&lh, &lv, NULL, 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);

    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,
                 cmdline.randomseedSpec, cmdline.randomseed);
    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;
}