/*============================================================================= 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 #include #include #include #include #include #include #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; }