#include #include #include #include "pm_c_util.h" #include "mallocvar.h" #include "shhopt.h" #include "pam.h" enum function {FN_ADD, FN_SUBTRACT, FN_MULTIPLY, FN_DIVIDE, FN_DIFFERENCE, FN_MINIMUM, FN_MAXIMUM, FN_MEAN, FN_COMPARE, FN_AND, FN_OR, FN_NAND, FN_NOR, FN_XOR, FN_SHIFTLEFT, FN_SHIFTRIGHT }; static bool isDyadic(enum function const function) { bool retval; switch (function) { case FN_ADD: case FN_MEAN: case FN_AND: case FN_OR: case FN_XOR: retval = FALSE; break; case FN_SUBTRACT: case FN_DIFFERENCE: case FN_MINIMUM: case FN_MAXIMUM: case FN_COMPARE: case FN_MULTIPLY: case FN_DIVIDE: case FN_NAND: case FN_NOR: case FN_SHIFTLEFT: case FN_SHIFTRIGHT: retval = TRUE; break; } return retval; } struct CmdlineInfo { /* All the information the user supplied in the command line, in a form easy for the program to use. */ enum function function; unsigned int operandCt; const char ** operandFileNames; }; static void parseCommandLine(int argc, const char ** const argv, struct CmdlineInfo * const cmdlineP) { /*---------------------------------------------------------------------------- Note that the file spec array we return is stored in the storage that was passed to us as the argv array. -----------------------------------------------------------------------------*/ optEntry * option_def; /* Instructions to OptParseOptions2 on how to parse our options. */ optStruct3 opt; unsigned int option_def_index; unsigned int addSpec, subtractSpec, multiplySpec, divideSpec, differenceSpec, minimumSpec, maximumSpec, meanSpec, compareSpec, andSpec, orSpec, nandSpec, norSpec, xorSpec, shiftleftSpec, shiftrightSpec; MALLOCARRAY_NOFAIL(option_def, 100); option_def_index = 0; /* incremented by OPTENTRY */ OPTENT3(0, "add", OPT_FLAG, NULL, &addSpec, 0); OPTENT3(0, "subtract", OPT_FLAG, NULL, &subtractSpec, 0); OPTENT3(0, "multiply", OPT_FLAG, NULL, &multiplySpec, 0); OPTENT3(0, "divide", OPT_FLAG, NULL, ÷Spec, 0); OPTENT3(0, "difference", OPT_FLAG, NULL, &differenceSpec, 0); OPTENT3(0, "minimum", OPT_FLAG, NULL, &minimumSpec, 0); OPTENT3(0, "maximum", OPT_FLAG, NULL, &maximumSpec, 0); OPTENT3(0, "mean", OPT_FLAG, NULL, &meanSpec, 0); OPTENT3(0, "compare", OPT_FLAG, NULL, &compareSpec, 0); OPTENT3(0, "and", OPT_FLAG, NULL, &andSpec, 0); OPTENT3(0, "or", OPT_FLAG, NULL, &orSpec, 0); OPTENT3(0, "nand", OPT_FLAG, NULL, &nandSpec, 0); OPTENT3(0, "nor", OPT_FLAG, NULL, &norSpec, 0); OPTENT3(0, "xor", OPT_FLAG, NULL, &xorSpec, 0); OPTENT3(0, "shiftleft", OPT_FLAG, NULL, &shiftleftSpec, 0); OPTENT3(0, "shiftright", OPT_FLAG, NULL, &shiftrightSpec, 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 some of *cmdlineP and others. */ if (addSpec + subtractSpec + multiplySpec + divideSpec + differenceSpec + minimumSpec + maximumSpec + meanSpec + compareSpec + andSpec + orSpec + nandSpec + norSpec + xorSpec + shiftleftSpec + shiftrightSpec > 1) pm_error("You may specify only one function"); if (addSpec) cmdlineP->function = FN_ADD; else if (subtractSpec) cmdlineP->function = FN_SUBTRACT; else if (multiplySpec) cmdlineP->function = FN_MULTIPLY; else if (divideSpec) cmdlineP->function = FN_DIVIDE; else if (differenceSpec) cmdlineP->function = FN_DIFFERENCE; else if (minimumSpec) cmdlineP->function = FN_MINIMUM; else if (maximumSpec) cmdlineP->function = FN_MAXIMUM; else if (meanSpec) cmdlineP->function = FN_MEAN; else if (compareSpec) cmdlineP->function = FN_COMPARE; else if (andSpec) cmdlineP->function = FN_AND; else if (orSpec) cmdlineP->function = FN_OR; else if (nandSpec) cmdlineP->function = FN_NAND; else if (norSpec) cmdlineP->function = FN_NOR; else if (xorSpec) cmdlineP->function = FN_XOR; else if (shiftleftSpec) cmdlineP->function = FN_SHIFTLEFT; else if (shiftrightSpec) cmdlineP->function = FN_SHIFTRIGHT; else pm_error("You must specify a function (e.g. '-add')"); if (argc-1 < 2) pm_error("You must specify at least two arguments: the files which " "are the operands of the function. You specified %u", argc-1); else { if (isDyadic(cmdlineP->function) && argc-1 > 2) pm_error("You specified %u arguments, but a dyadic function. " "For a dyadic function, you must specify 2 arguments: " "the operands of the function", argc-1); else { cmdlineP->operandCt = argc-1; cmdlineP->operandFileNames = &argv[1]; } } } enum category { CATEGORY_FRACTIONAL_ARITH, /* Arithmetic in which each sample represents a the fraction sample/maxval. */ CATEGORY_BITSTRING, /* And, Or, etc. Maxval isn't a scale factor at all; it's a mask. */ CATEGORY_SHIFT /* Left argument is a bit string, but right argument is a whole number (left maxval is a mask; right maxval is meaningless). */ }; static enum category functionCategory(enum function const function) { enum category retval; switch (function) { case FN_ADD: case FN_SUBTRACT: case FN_DIFFERENCE: case FN_MINIMUM: case FN_MAXIMUM: case FN_MEAN: case FN_COMPARE: case FN_MULTIPLY: case FN_DIVIDE: retval = CATEGORY_FRACTIONAL_ARITH; break; case FN_AND: case FN_OR: case FN_NAND: case FN_NOR: case FN_XOR: retval = CATEGORY_BITSTRING; break; case FN_SHIFTLEFT: case FN_SHIFTRIGHT: retval = CATEGORY_SHIFT; break; } return retval; } static void computeOutputType(struct pam * const outpamP, struct pam const inpam1, struct pam const inpam2, enum function const function) { outpamP->size = sizeof(struct pam); outpamP->len = PAM_STRUCT_SIZE(tuple_type); outpamP->file = stdout; outpamP->format = MAX(inpam1.format, inpam2.format); outpamP->plainformat = FALSE; outpamP->height = inpam1.height; outpamP->width = inpam1.width; outpamP->depth = MAX(inpam1.depth, inpam2.depth); switch (functionCategory(function)) { case CATEGORY_FRACTIONAL_ARITH: if (function == FN_COMPARE) outpamP->maxval = 2; else outpamP->maxval = MAX(inpam1.maxval, inpam2.maxval); break; case CATEGORY_BITSTRING: if (inpam2.maxval != inpam1.maxval) pm_error("For a bit string operation, the maxvals of the " "left and right image must be the same. You have " "left=%u and right=%u", (unsigned)inpam1.maxval, (unsigned)inpam2.maxval); if (pm_bitstomaxval(pm_maxvaltobits(inpam1.maxval)) != inpam1.maxval) pm_error("For a bit string operation, the maxvals of the inputs " "must be a full binary count, i.e. a power of two " "minus one such as 0xff. You have 0x%x", (unsigned)inpam1.maxval); outpamP->maxval = inpam1.maxval; break; case CATEGORY_SHIFT: if (pm_bitstomaxval(pm_maxvaltobits(inpam1.maxval)) != inpam1.maxval) pm_error("For a bit shift operation, the maxval of the left " "input image " "must be a full binary count, i.e. a power of two " "minus one such as 0xff. You have 0x%x", (unsigned)inpam1.maxval); outpamP->maxval = inpam1.maxval; } outpamP->bytes_per_sample = (pm_maxvaltobits(outpamP->maxval)+7)/8; strcpy(outpamP->tuple_type, inpam1.tuple_type); } static samplen samplenSum(samplen const operands[], unsigned int const operandCt) { unsigned int i; samplen total; for (i = 1, total = operands[0]; i < operandCt; ++i) { total += operands[1]; if (total > 1.) total = 1.; } return total; } static samplen samplenMin(samplen const operands[], unsigned int const operandCt) { unsigned int i; samplen min; for (i = 1, min = operands[0]; i < operandCt; ++i) { if (operands[i] < min) min = operands[i]; } return min; } static samplen samplenMax(samplen const operands[], unsigned int const operandCt) { unsigned int i; samplen max; for (i = 1, max = operands[0]; i < operandCt; ++i) { if (operands[i] > max) max = operands[i]; } return max; } static samplen samplenMean(samplen const operands[], unsigned int const operandCt) { unsigned int i; double sum; for (i = 0, sum = 0.; i < operandCt; ++i) sum += operands[i]; return sum / operandCt; } static samplen samplenProduct(samplen const operands[], unsigned int const operandCt) { unsigned int i; double product; for (i = 1, product = operands[0]; i < operandCt; ++i) product *= operands[i]; return product; } static samplen applyNormalizedFunction(enum function const function, samplen const operands[], unsigned int const operandCt) { samplen result; switch (function) { case FN_ADD: result = samplenSum(operands, operandCt); break; case FN_SUBTRACT: result = MAX(0., operands[0] - operands[1]); break; case FN_MULTIPLY: result = samplenProduct(operands, operandCt); break; case FN_DIVIDE: result = (operands[1] > operands[0]) ? operands[0] / operands[1] : 1.; break; case FN_DIFFERENCE: result = operands[0] > operands[1] ? operands[0] - operands[1] : operands[1] - operands[0]; break; case FN_MINIMUM: result = samplenMin(operands, operandCt); break; case FN_MAXIMUM: result = samplenMax(operands, operandCt); break; case FN_MEAN: result = samplenMean(operands, operandCt); break; case FN_COMPARE: result = operands[0] > operands[1] ? 1. : operands[0] < operands[1] ? 0. : .5; break; default: pm_error("Internal error. applyNormalizedFunction() called " "for a function it doesn't know how to do: %u", function); } return result; } static void doNormalizedArith(struct pam * const inpam1P, struct pam * const inpam2P, struct pam * const outpamP, enum function const function) { /* Some of the logic in this subroutine is designed for future expansion into non-dyadic computations. But for now, all computations have exactly two operands. */ unsigned int const operandCt = 2; tuplen ** tuplerown; /* tuplerown[0] is the current row in the first operand image */ tuplen * tuplerownOut; unsigned int row; samplen * operands; /* operand[0] is the first operand in the current one-sample computation */ unsigned int * plane; /* plane[0] is the plane number in the first operand image for the current one-sample computation. plane[1] is the plane number in the second operand image, etc. */ MALLOCARRAY_NOFAIL(operands, operandCt); MALLOCARRAY_NOFAIL(plane, operandCt); MALLOCARRAY_NOFAIL(tuplerown, operandCt); tuplerown[0] = pnm_allocpamrown(inpam1P); tuplerown[1] = pnm_allocpamrown(inpam2P); tuplerownOut = pnm_allocpamrown(outpamP); for (row = 0; row < outpamP->height; ++row) { unsigned int col; pnm_readpamrown(inpam1P, tuplerown[0]); pnm_readpamrown(inpam2P, tuplerown[1]); for (col = 0; col < outpamP->width; ++col) { unsigned int outplane; for (outplane = 0; outplane < outpamP->depth; ++outplane) { unsigned int op; plane[0] = MIN(outplane, inpam1P->depth-1); plane[1] = MIN(outplane, inpam2P->depth-1); for (op = 0; op < operandCt; ++op) operands[op] = tuplerown[op][col][plane[op]]; tuplerownOut[col][outplane] = applyNormalizedFunction(function, operands, operandCt); assert(tuplerownOut[col][outplane] >= 0.); assert(tuplerownOut[col][outplane] <= 1.); } } pnm_writepamrown(outpamP, tuplerownOut); } pnm_freepamrown(tuplerown[0]); pnm_freepamrown(tuplerown[1]); free(tuplerown); pnm_freepamrown(tuplerownOut); free(plane); free(operands); } static sample sampleSum(sample const operands[], unsigned int const operandCt, sample const maxval) { unsigned int i; sample total; for (i = 1, total = operands[0]; i < operandCt; ++i) { total += operands[i]; if (total > maxval) total = maxval; } return total; } static sample sampleMin(sample const operands[], unsigned int const operandCt) { unsigned int i; sample min; for (i = 1, min = operands[0]; i < operandCt; ++i) { if (operands[i] < min) min = operands[i]; } return min; } static sample sampleMax(sample const operands[], unsigned int const operandCt) { unsigned int i; sample max; for (i = 1, max = operands[0]; i < operandCt; ++i) { if (operands[i] > max) max = operands[i]; } return max; } static sample sampleMean(sample const operands[], unsigned int const operandCt) { unsigned int i; unsigned int sum; for (i = 0, sum = 0; i < operandCt; ++i) { sum += operands[i]; if (UINT_MAX - operands[i] < sum) pm_error("Arithmetic overflow adding samples for mean"); } return ROUNDDIV(sum, operandCt); } static sample sampleProduct(sample const operands[], unsigned int const operandCt, sample const maxval) { unsigned int i; double product; for (i = 0, product = 1.0; i < operandCt; ++i) { product *= ((double)operands[i]/maxval); } return (sample)(product * maxval + 0.5); } static sample sampleAnd(sample const operands[], unsigned int const operandCt) { unsigned int i; sample accum; for (i = 1, accum = operands[0]; i < operandCt; ++i) { accum &= operands[i]; } return accum; } static sample sampleOr(sample const operands[], unsigned int const operandCt) { unsigned int i; sample accum; for (i = 1, accum = operands[0]; i < operandCt; ++i) { accum |= operands[i]; } return accum; } static sample sampleNand(sample const operands[], unsigned int const operandCt, sample const maxval) { unsigned int i; sample accum; for (i = 1, accum = operands[0]; i < operandCt; ++i) { accum &= operands[i]; } return ~accum & maxval; } static sample sampleNor(sample const operands[], unsigned int const operandCt, sample const maxval) { unsigned int i; sample accum; for (i = 1, accum = operands[0]; i < operandCt; ++i) { accum |= operands[i]; } return ~accum & maxval; } static sample sampleXor(sample const operands[], unsigned int const operandCt) { unsigned int i; sample accum; for (i = 1, accum = operands[0]; i < operandCt; ++i) { accum ^= operands[i]; } return accum; } static sample applyUnNormalizedFunction(enum function const function, sample const operands[], unsigned int const operandCt, sample const maxval) { /*---------------------------------------------------------------------------- Apply function 'function' to the arguments operands[] (there are 'operandCt' of them), assuming both are based on the same maxval 'maxval'. Return a value which is also a fraction of 'maxval'. Exception: for the shift operations, operands[1] is not based on any maxval. It is an absolute bit count. We assume 'operandCount' is sensible for 'function'. E.g. if 'function' is FN_DIFFERENCE, 'operandCt' is 2. For a function that has a concept of left and right argument, operands[0] is the left and operands[1] is the right. -----------------------------------------------------------------------------*/ sample result; switch (function) { case FN_ADD: result = sampleSum(operands, operandCt, maxval); break; case FN_SUBTRACT: result = MAX(0, (int)operands[0] - (int)operands[1]); break; case FN_DIFFERENCE: result = operands[0] > operands[1] ? operands[0] - operands[1] : operands[1] - operands[0]; break; case FN_MINIMUM: result = sampleMin(operands, operandCt); break; case FN_MAXIMUM: result = sampleMax(operands, operandCt); break; case FN_MEAN: result = sampleMean(operands, operandCt); break; case FN_COMPARE: result = operands[0] > operands[1] ? 2 : operands[0] < operands[1] ? 0 : 1; break; case FN_MULTIPLY: result = sampleProduct(operands, operandCt, maxval); break; case FN_DIVIDE: result = (operands[1] > operands[0]) ? ROUNDDIV(operands[0] * maxval, operands[1]) : maxval; break; case FN_AND: result = sampleAnd(operands, operandCt); break; case FN_OR: result = sampleOr(operands, operandCt); break; case FN_NAND: result = sampleNand(operands, operandCt, maxval); break; case FN_NOR: result = sampleNor(operands, operandCt, maxval); break; case FN_XOR: result = sampleXor(operands, operandCt); break; case FN_SHIFTLEFT: result = (operands[0] << operands[1]) & maxval; break; case FN_SHIFTRIGHT: result = operands[0] >> operands[1]; break; default: pm_error("Internal error. applyUnNormalizedFunction() called " "for a function it doesn't know how to do: %u", function); } return result; } static void doUnNormalizedArith(struct pam * const inpam1P, struct pam * const inpam2P, struct pam * const outpamP, enum function const function) { /*---------------------------------------------------------------------------- Take advantage of the fact that both inputs and the output use the same maxval to do the computation without time-consuming normalization of sample values. -----------------------------------------------------------------------------*/ /* Some of the logic in this subroutine is designed for future expansion into non-dyadic computations. But for now, all computations have exactly two operands. */ unsigned int const operandCt = 2; sample const maxval = outpamP->maxval; tuple ** tuplerow; /* tuplerow[0] is the current row in the first operand image */ tuple * tuplerowOut; unsigned int row; sample * operands; /* operand[0] is the first operand in the current one-sample computation */ unsigned int * plane; /* plane[0] is the plane number in the first operand image for the current one-sample computation. plane[1] is the plane number in the second operand image, etc. */ /* Input conditions: */ assert(inpam1P->maxval == maxval); assert(inpam2P->maxval == maxval); assert(outpamP->maxval == maxval); MALLOCARRAY_NOFAIL(operands, operandCt); MALLOCARRAY_NOFAIL(plane, operandCt); MALLOCARRAY_NOFAIL(tuplerow, operandCt); tuplerow[0] = pnm_allocpamrow(inpam1P); tuplerow[1] = pnm_allocpamrow(inpam2P); tuplerowOut = pnm_allocpamrow(outpamP); for (row = 0; row < outpamP->height; ++row) { unsigned int col; pnm_readpamrow(inpam1P, tuplerow[0]); pnm_readpamrow(inpam2P, tuplerow[1]); for (col = 0; col < outpamP->width; ++col) { unsigned int outplane; for (outplane = 0; outplane < outpamP->depth; ++outplane) { unsigned int op; plane[0] = MIN(outplane, inpam1P->depth-1); plane[1] = MIN(outplane, inpam2P->depth-1); for (op = 0; op < operandCt; ++op) operands[op] = tuplerow[op][col][plane[op]]; tuplerowOut[col][outplane] = applyUnNormalizedFunction(function, operands, operandCt, maxval); assert(tuplerowOut[col][outplane] >= 0); assert(tuplerowOut[col][outplane] <= outpamP->maxval); } } pnm_writepamrow(outpamP, tuplerowOut); } pnm_freepamrow(tuplerow[0]); pnm_freepamrow(tuplerow[1]); free(tuplerow); pnm_freepamrow(tuplerowOut); free(plane); free(operands); } int main(int argc, const char *argv[]) { struct CmdlineInfo cmdline; struct pam inpam1; struct pam inpam2; struct pam outpam; FILE * if1P; FILE * if2P; pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); if (cmdline.operandCt != 2) /* Code for > 2 operands not written yet */ pm_error("You specified %u operands. We understand only 2.", cmdline.operandCt); if1P = pm_openr(cmdline.operandFileNames[0]); if2P = pm_openr(cmdline.operandFileNames[1]); pnm_readpaminit(if1P, &inpam1, PAM_STRUCT_SIZE(tuple_type)); pnm_readpaminit(if2P, &inpam2, PAM_STRUCT_SIZE(tuple_type)); if (inpam1.width != inpam2.width || inpam1.height != inpam2.height) pm_error("The two images must be the same width and height. " "The first is %ux%ux%u, but the second is %ux%ux%u", inpam1.width, inpam1.height, inpam1.depth, inpam2.width, inpam2.height, inpam2.depth); if (inpam1.depth != 1 && inpam2.depth != 1 && inpam1.depth != inpam2.depth) pm_error("The two images must have the same depth or one of them " "must have depth 1. The first has depth %u and the second " "has depth %u", inpam1.depth, inpam2.depth); computeOutputType(&outpam, inpam1, inpam2, cmdline.function); pnm_writepaminit(&outpam); switch (functionCategory(cmdline.function)) { case CATEGORY_FRACTIONAL_ARITH: if (inpam1.maxval == inpam2.maxval) doUnNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function); else doNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function); break; case CATEGORY_BITSTRING: case CATEGORY_SHIFT: doUnNormalizedArith(&inpam1, &inpam2, &outpam, cmdline.function); break; } pm_close(if1P); pm_close(if2P); return 0; }