From 36b161da4598883d6afb40317de32d70f14e3ff3 Mon Sep 17 00:00:00 2001 From: giraffedata Date: Mon, 4 Jan 2021 00:39:49 +0000 Subject: Add pamhomography git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@4022 9d0c8265-081b-0410-96cb-a4ca84ce46f8 --- editor/pamhomography.c | 677 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 677 insertions(+) create mode 100644 editor/pamhomography.c (limited to 'editor/pamhomography.c') diff --git a/editor/pamhomography.c b/editor/pamhomography.c new file mode 100644 index 00000000..59b59ed7 --- /dev/null +++ b/editor/pamhomography.c @@ -0,0 +1,677 @@ +/* ---------------------------------------------------------------------- + * + * Map one quadrilateral to another + * by Scott Pakin + * + * ---------------------------------------------------------------------- + * + * Copyright (C) 2020 Scott Pakin + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * ---------------------------------------------------------------------- + */ + +#include +#include +#include +#include +#include +#include +#include "pm_c_util.h" +#include "shhopt.h" +#include "mallocvar.h" +#include "pam.h" + +#define MIN4(A, B, C, D) MIN(MIN(A, B), MIN(C, D)) +#define MAX4(A, B, C, D) MAX(MAX(A, B), MAX(C, D)) + +/* A point on the image plane. It may or may not lie within the + bounds of the image itself. */ +typedef struct Point { + int x; + int y; +} Point; + +/* A quadrilateral on the image plane */ +typedef struct Quad { + Point ul; + Point ur; + Point lr; + Point ll; +} Quad; + +/* A user-specified mapping from one quadrilateral to another */ +typedef struct QuadMap { + Quad from; + Quad to; +} QuadMap; + +struct CmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char * inputFilespec; /* "-" if stdin */ + QuadMap qmap; /* Source and target quadrilaterals */ + Quad bbox; /* Bounding box for the target image */ + const char * fillColor; /* Fill color for unused coordinates */ +}; + + + +static unsigned int +parseCoords(const char * const str, + int * const coords) { +/*---------------------------------------------------------------------------- + Parse a list of up to 16 integers. The function returns the number + of integers encountered. +-----------------------------------------------------------------------------*/ + + const char * p; + char * pnext; + unsigned int i; + + for (i = 0, p = str; i < 16; ++i, p = pnext) { + long int val; + + /* Skip punctuation, except "+" and "-", and white space. */ + while (*p != '\0' && *p != '+' && *p != '-' && + (isspace(*p) || ispunct(*p))) + ++p; + + /* Parse the next integer. */ + errno = 0; /* strtol() sets errno on error. */ + val = strtol(p, &pnext, 10); + if (errno == ERANGE) + return i; /* Integer lies out of long int range */ + if (errno != 0 || pnext == p) + return i; /* Too few integers */ + coords[i] = (int)val; + if ((long int)coords[i] != val) + return i; /* Integer lies out of int range */ + } + return i; +} + + + +static void +parseViewString(const char * const str, + Quad * const quad) { +/*---------------------------------------------------------------------------- + Parse a list of four integers in the order {ulx, uly, lrx, lry} into a + quadrilateral. The function aborts on error. +-----------------------------------------------------------------------------*/ + + int coords[16]; + + if (parseCoords(str, coords) != 4) + pm_error("failed to parse \"%s\" as a list of four integers", str); + quad->ul.x = quad->ll.x = coords[0]; + quad->ul.y = quad->ur.y = coords[1]; + quad->lr.x = quad->ur.x = coords[2]; + quad->lr.y = quad->ll.y = coords[3]; +} + + + +static void +parseQuadString(const char * const str, + Quad * const quad) { +/*---------------------------------------------------------------------------- + Parse a list of eight integers in the order {ulx, uly, urx, ury, + lrx, lry, llx, lly} into a quadrilateral. The function aborts on + error. +-----------------------------------------------------------------------------*/ + + int coords[16]; + + if (parseCoords(str, coords) != 8) + pm_error("failed to parse \"%s\" as a list of eight integers", str); + quad->ul.x = coords[0]; + quad->ul.y = coords[1]; + quad->ur.x = coords[2]; + quad->ur.y = coords[3]; + quad->lr.x = coords[4]; + quad->lr.y = coords[5]; + quad->ll.x = coords[6]; + quad->ll.y = coords[7]; +} + + + +static void +readMapFile(const char * const fname, + QuadMap * const qmap) { +/*---------------------------------------------------------------------------- + Read from a file either 16 numbers in the order {ulx1, uly1, urx1, ury1, + lrx1, lry1, llx1, lly1, ulx2, uly2, urx2, ury2, lrx2, lry2, llx2, lly2} + or 8 numbers in the order {ulx2, uly2, urx2, ury2, lrx2, lry2, llx2, + lly2}. This function aborts on error. +-----------------------------------------------------------------------------*/ + + FILE * fp; + char * str; /* Entire file contents */ + int coords[16]; /* File as a list of up to 16 coordinates */ + char * c; + long int nread; + + /* Read the entire file. */ + fp = pm_openr(fname); + str = pm_read_unknown_size(fp, &nread); + REALLOCARRAY_NOFAIL(str, nread + 1); + str[nread] = '\0'; + pm_close(fp); + + /* Replace newlines and tabs with spaces to prettify error reporting. */ + for (c = str; *c != '\0'; ++c) + if (isspace(*c)) + *c = ' '; + + /* Read either {from, to} or just a {to} quadrilateral. */ + switch (parseCoords(str, coords)) { + case 16: + /* 16 integers: assign both the "from" and the "to" quadrilateral. */ + qmap->from.ul.x = coords[0]; + qmap->from.ul.y = coords[1]; + qmap->from.ur.x = coords[2]; + qmap->from.ur.y = coords[3]; + qmap->from.lr.x = coords[4]; + qmap->from.lr.y = coords[5]; + qmap->from.ll.x = coords[6]; + qmap->from.ll.y = coords[7]; + qmap->to.ul.x = coords[8]; + qmap->to.ul.y = coords[9]; + qmap->to.ur.x = coords[10]; + qmap->to.ur.y = coords[11]; + qmap->to.lr.x = coords[12]; + qmap->to.lr.y = coords[13]; + qmap->to.ll.x = coords[14]; + qmap->to.ll.y = coords[15]; + break; + case 8: + /* 8 integers: assign only the "to" quadrilateral. */ + memset((void *)&qmap->from, 0, sizeof(Quad)); + qmap->to.ul.x = coords[0]; + qmap->to.ul.y = coords[1]; + qmap->to.ur.x = coords[2]; + qmap->to.ur.y = coords[3]; + qmap->to.lr.x = coords[4]; + qmap->to.lr.y = coords[5]; + qmap->to.ll.x = coords[6]; + qmap->to.ll.y = coords[7]; + break; + default: + /* Any other number of integers: issue an error message. */ + pm_error("failed to parse \"%s\" as a list of either 8 or 16 integers", + str); + break; + } + + free(str); +} + + + +static void +parseCommandLine(int argc, + const char ** const argv, + struct CmdlineInfo * const cmdlineP ) { +/*---------------------------------------------------------------------------- + Parse program command line described in Unix standard form by argc + and argv. Return the information in the options as *cmdlineP. + + If command line is internally inconsistent (invalid options, etc.), + issue error message to stderr and abort program. +-----------------------------------------------------------------------------*/ + + optEntry *option_def; + /* Instructions to pm_optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + unsigned int mapFileSpec = 0, fromSpec = 0, toSpec = 0; + unsigned int viewSpec = 0, fillColorSpec = 0; + char *mapFile, *from, *to, *view; + + MALLOCARRAY_NOFAIL(option_def, 100); + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "mapfile", OPT_STRING, &mapFile, + &mapFileSpec, 0); + OPTENT3(0, "from", OPT_STRING, &from, + &fromSpec, 0); + OPTENT3(0, "to", OPT_STRING, &to, + &toSpec, 0); + OPTENT3(0, "view", OPT_STRING, &view, + &viewSpec, 0); + OPTENT3(0, "fill", OPT_STRING, &cmdlineP->fillColor, + &fillColorSpec, 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 */ + + pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and local variables. */ + + if (!fillColorSpec) + cmdlineP->fillColor = NULL; + + memset((void *)&cmdlineP->qmap, 0, sizeof(QuadMap)); + if (mapFileSpec) + readMapFile(mapFile, &cmdlineP->qmap); + if (fromSpec) + parseQuadString(from, &cmdlineP->qmap.from); + if (toSpec) + parseQuadString(to, &cmdlineP->qmap.to); + if (!mapFileSpec && !fromSpec && !toSpec && !viewSpec) + pm_error("You must specify at least one of " + "-mapfile, -qin, -qout, and -view"); + if (viewSpec) + parseViewString(view, &cmdlineP->bbox); + else + memset((void *)&cmdlineP->bbox, 0, sizeof(Quad)); + + if (argc < 2) + cmdlineP->inputFilespec = "-"; + else if (argc == 2) + cmdlineP->inputFilespec = argv[1]; + else + pm_error("Too many non-option arguments: %u. " + "Only argument is input file name", argc - 1); + + free((void *) option_def); +} + + + +static tuple +parseFillColor(const struct pam * const pamP, + const struct CmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Parse the fill color into the correct format for the given PAM metadata. +-----------------------------------------------------------------------------*/ + + tuple rgb; + tuple fillColor; + + if (!cmdlineP->fillColor) { + pnm_createBlackTuple(pamP, &fillColor); + return fillColor; + } + + rgb = pnm_parsecolor(cmdlineP->fillColor, pamP->maxval); + fillColor = pnm_allocpamtuple(pamP); + switch (pamP->depth) { + case 1: + /* Grayscale */ + fillColor[0] = (rgb[PAM_RED_PLANE]*299 + + rgb[PAM_GRN_PLANE]*587 + + rgb[PAM_BLU_PLANE]*114)/1000; + break; + case 2: + /* Grayscale + alpha */ + fillColor[0] = (rgb[PAM_RED_PLANE]*299 + + rgb[PAM_GRN_PLANE]*587 + + rgb[PAM_BLU_PLANE]*114)/1000; + fillColor[PAM_GRAY_TRN_PLANE] = pamP->maxval; + break; + case 3: + /* RGB */ + pnm_assigntuple(pamP, fillColor, rgb); + break; + case 4: + /* RGB + alpha */ + pnm_assigntuple(pamP, fillColor, rgb); + fillColor[PAM_TRN_PLANE] = pamP->maxval; + break; + default: + pm_error("unexpected image depth %d", pamP->depth); + break; + } + + return fillColor; +} + + + +static tuple ** +initOutputImage(const struct pam * const pamP, + const struct CmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Allocate and initialize the output image. +-----------------------------------------------------------------------------*/ + + tuple fillColor; /* Fill color to use for unused coordinates */ + tuple ** outImg; /* Output image */ + unsigned int row; + + outImg = pnm_allocpamarray(pamP); + + fillColor = parseFillColor(pamP, cmdlineP); + for (row = 0; row < pamP->height; ++row) { + unsigned int col; + + for (col = 0; col < pamP->width; ++col) { + pnm_assigntuple(pamP, outImg[row][col], fillColor); + } + } + + free((void *) fillColor); + return outImg; +} + + + +static void +computeSteps(const Quad * const qfrom, + const Quad * const qto, + double * const ustep, + double * const vstep) { +/*---------------------------------------------------------------------------- + Compute increments for u and v as these range from 0.0 to 1.0. +-----------------------------------------------------------------------------*/ + + double fx0, fx1, fxd; + double tx0, tx1, txd; + double fy0, fy1, fyd; + double ty0, ty1, tyd; + + /* Compute ustep as the inverse of the maximum possible x delta across + either the "from" or "to" quadrilateral. */ + fx0 = MIN4((double)qfrom->ur.x, + (double)qfrom->ul.x, + (double)qfrom->lr.x, + (double)qfrom->ll.x); + fx1 = MAX4((double)qfrom->ur.x, + (double)qfrom->ul.x, + (double)qfrom->lr.x, + (double)qfrom->ll.x); + fxd = fx1 - fx0; + tx0 = MIN4((double)qto->ur.x, + (double)qto->ul.x, + (double)qto->lr.x, + (double)qto->ll.x); + tx1 = MAX4((double)qto->ur.x, + (double)qto->ul.x, + (double)qto->lr.x, + (double)qto->ll.x); + txd = tx1 - tx0; + if (fxd == 0.0 && txd == 0.0) + *ustep = 1.0; /* Arbitrary nonzero step */ + *ustep = 0.5/MAX(fxd, txd); + /* Divide into 0.5 instead of 1.0 for additional smoothing. */ + + /* Compute vstep as the inverse of the maximum possible y delta across + either the "from" or "to" quadrilateral + . */ + fy0 = MIN4((double)qfrom->ur.y, + (double)qfrom->ul.y, + (double)qfrom->lr.y, + (double)qfrom->ll.y); + fy1 = MAX4((double)qfrom->ur.y, + (double)qfrom->ul.y, + (double)qfrom->lr.y, + (double)qfrom->ll.y); + fyd = fy1 - fy0; + ty0 = MIN4((double)qto->ur.y, + (double)qto->ul.y, + (double)qto->lr.y, + (double)qto->ll.y); + ty1 = MAX4((double)qto->ur.y, + (double)qto->ul.y, + (double)qto->lr.y, + (double)qto->ll.y); + tyd = ty1 - ty0; + if (fyd == 0.0 && tyd == 0.0) + *vstep = 1.0; /* Arbitrary nonzero step */ + *vstep = 0.5/MAX(fyd, tyd); + /* Divide into 0.5 instead of 1.0 for additional smoothing. */ +} + + + +static Quad * +prepareQuadrilateral(const struct pam * const pamP, + const Quad * const qdata) { +/*---------------------------------------------------------------------------- + If a quadrilateral has all zero points, replace it with a quadrilateral + of the full size of the image. The caller should free the result. +-----------------------------------------------------------------------------*/ + + Quad * qcopy; + + MALLOCVAR_NOFAIL(qcopy); + + if (qdata->ul.x == 0 && qdata->ul.y == 0 && + qdata->ur.x == 0 && qdata->ur.y == 0 && + qdata->ll.x == 0 && qdata->ll.y == 0 && + qdata->lr.x == 0 && qdata->lr.y == 0) { + /* Set the quadrilateral to the image's bounding box. */ + memset((void *)qcopy, 0, sizeof(Quad)); + qcopy->ur.x = pamP->width - 1; + qcopy->lr.x = pamP->width - 1; + qcopy->lr.y = pamP->height - 1; + qcopy->ll.y = pamP->height - 1; + } else { + /* Use the quadrilateral as specified. */ + memcpy(qcopy, qdata, sizeof(Quad)); + } + + return qcopy; +} + + +static void +coordsAtPercent(const Quad * const quad, + double const u, + double const v, + int * const x, + int * const y) { +/*---------------------------------------------------------------------------- + Return the (x, y) coordinates that lie at (u%, v%) from the upper left to + the lower right of a given quadrilateral. +-----------------------------------------------------------------------------*/ + + *x = (int) nearbyint((1.0 - u)*(1.0 - v)*quad->ul.x + + u*(1.0 - v)*quad->ur.x + + u*v*quad->lr.x + + (1.0 - u)*v*quad->ll.x); + *y = (int) nearbyint((1.0 - u)*(1.0 - v)*quad->ul.y + + u*(1.0 - v)*quad->ur.y + + u*v*quad->lr.y + + (1.0 - u)*v*quad->ll.y); +} + + + +static void +computeBoundingBox(const Quad * const q, + Quad * const bbox) { +/*---------------------------------------------------------------------------- + Compute the bounding box of a given quadrilateral. +-----------------------------------------------------------------------------*/ + + bbox->ul.x = bbox->ll.x = MIN4(q->ul.x, q->ur.x, q->lr.x, q->ll.x); + bbox->ul.y = bbox->ur.y = MIN4(q->ul.y, q->ur.y, q->lr.y, q->ll.y); + bbox->ur.x = bbox->lr.x = MAX4(q->ul.x, q->ur.x, q->lr.x, q->ll.x); + bbox->ll.y = bbox->lr.y = MAX4(q->ul.y, q->ur.y, q->lr.y, q->ll.y); +} + + + +static void +mapQuadrilaterals(const struct pam * const inPamP, + const struct pam * const outPamP, + const Quad * const qfrom, + const Quad * const qto, + tuple ** const inImg, + tuple ** const outImg, + int const xofs, + int const yofs) { +/*---------------------------------------------------------------------------- + Map the quadrilateral in the source image to the quadrilateral in the + target image. This is the function that implemens pamhomography's + primary functionality. +-----------------------------------------------------------------------------*/ + + sample ** channel; + /* Aggregated values for a single channel */ + unsigned long ** tally; + /* Number of values at each coordinate in the above */ + double ustep, vstep; + /* Steps to use when iterating from 0.0 to 1.0 */ + double u, v; + unsigned int plane, row, col; + + MALLOCARRAY2_NOFAIL(channel, outPamP->height, outPamP->width); + MALLOCARRAY2_NOFAIL(tally, outPamP->height, outPamP->width); + + computeSteps(qfrom, qto, &ustep, &vstep); + + for (plane = 0; plane < outPamP->depth; ++plane) { + /* Reset the channel colors and tally for each plane, */ + for (row = 0; row < outPamP->height; ++row) + for (col = 0; col < outPamP->width; ++col) { + channel[row][col] = 0; + tally[row][col] = 0; + } + + /* Iterate from 0% to 100% in the y dimension. */ + for (v = 0.0; v <= 1.0; v += vstep) { + /* Iterate from 0% to 100% in the x dimension. */ + for (u = 0.0; u <= 1.0; u += ustep) { + int x0, y0; /* "From" coordinate */ + int x1, y1; /* "To" coordinate */ + + /* Map (u%, v%) of one quadrilateral to (u%, v%) of the + other quadrilateral. */ + coordsAtPercent(qfrom, u, v, &x0, &y0); + coordsAtPercent(qto, u, v, &x1, &y1); + + /* Copy the source image's (x0, y0) to the destination + image's (x1, y1) in the current plane. */ + x1 += xofs; + y1 += yofs; + if (x0 >= 0 && y0 >= 0 && + x0 < inPamP->width && y0 < inPamP->height && + x1 >= 0 && y1 >= 0 && + x1 < outPamP->width && y1 < outPamP->height) { + channel[y1][x1] += inImg[y0][x0][plane]; + tally[y1][x1]++; + } + } + } + + /* Assign the current plane in the output image the average color + at each point. */ + for (row = 0; row < outPamP->height; ++row) + for (col = 0; col < outPamP->width; ++col) + if (tally[row][col] != 0) + outImg[row][col][plane] = + (channel[row][col] + tally[row][col]/2) / + tally[row][col]; + } + + pm_freearray2((void ** const)tally); + pm_freearray2((void ** const)channel); + free((void *)qto); + free((void *)qfrom); +} + + + +static void +processFile(FILE * const ifP, + const struct CmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Read the input image, create the output image, and map a quadrilateral in + the former to a quadrilateral in the latter. +-----------------------------------------------------------------------------*/ + + struct pam inPam; /* PAM metadata for the input file */ + struct pam outPam; /* PAM metadata for the output file */ + tuple ** inImg; /* Input image */ + tuple ** outImg; /* Output image */ + Quad *qfrom, *qto; /* Source and target quadrilaterals */ + Quad bbox; /* Bounding box around the transformed input image */ + + inImg = pnm_readpam(ifP, &inPam, PAM_STRUCT_SIZE(tuple_type)); + + /* Extract quadrilaterals and populate them with the image bounds + if necessary. */ + qfrom = prepareQuadrilateral(&inPam, &cmdlineP->qmap.from); + qto = prepareQuadrilateral(&inPam, &cmdlineP->qmap.to); + + /* Allocate storage for the target image. */ + if (cmdlineP->bbox.ul.x == 0 && cmdlineP->bbox.ul.y == 0 && + cmdlineP->bbox.lr.x == 0 && cmdlineP->bbox.lr.y == 0) + /* User did not specify a target bounding box. Compute optimal + dimensions. */ + computeBoundingBox(qto, &bbox); + else + /* User specified a target bounding box. Use it. */ + bbox = cmdlineP->bbox; + outPam = inPam; + outPam.file = stdout; + outPam.width = bbox.lr.x - bbox.ul.x + 1; + outPam.height = bbox.lr.y - bbox.ul.y + 1; + outImg = initOutputImage(&outPam, cmdlineP); + + mapQuadrilaterals(&inPam, &outPam, + qfrom, qto, + inImg, outImg, + -bbox.ul.x, -bbox.ul.y); + + pnm_writepam(&outPam, outImg); + + pnm_freepamarray(outImg, &outPam); + pnm_freepamarray(inImg, &inPam); +} + + + +int +main(int argc, const char *argv[]) { + + struct CmdlineInfo cmdline; /* Parsed command line */ + FILE * ifP; + + pm_proginit(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + processFile(ifP, &cmdline); + + pm_close(ifP); + + return 0; +} -- cgit 1.4.1