/* pamscale.c - rescale (resample) a PNM image This program evolved out of Jef Poskanzer's program Pnmscale from his Pbmplus package (which was derived from Poskanzer's 1989 Ppmscale). The resampling logic was taken from Michael Reinelt's program Pnmresample, somewhat recoded to follow Netpbm conventions. Michael submitted that for inclusion in Netpbm in December 2003. The frame of the program is by Bryan Henderson, and the old scaling algorithm is based on that in Jef Poskanzer's Pnmscale, but completely rewritten by Bryan Henderson ca. 2000. Plenty of other people contributed code changes over the years. Copyright (C) 2003 by Michael Reinelt Copyright (C) 1989, 1991 by Jef Poskanzer. 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. */ #define _XOPEN_SOURCE 500 /* get M_PI in math.h */ #include #include #include #include #include #include #include "pm_c_util.h" #include "mallocvar.h" #include "shhopt.h" #include "pam.h" /****************************/ /****************************/ /********* filters **********/ /****************************/ /****************************/ /* Most of the filters are FIR (finite impulse respone), but some ** (sinc, bessel) are IIR (infinite impulse respone). ** They should be windowed with hanning, hamming, blackman or ** kaiser window. ** For sinc and bessel the blackman window will be used by default. */ #define EPSILON 1e-7 /* x^2 and x^3 helper functions */ static __inline__ double pow2(double const x) { return x * x; } static __inline__ double pow3(double const x) { return x * x * x; } /* box, pulse, Fourier window, */ /* box function also known as rectangle function */ /* 1st order (constant) b-spline */ #define radius_point (0.0) #define radius_box (0.5) static double filter_box(double const x) { double const absx = x < 0.0 ? -x : x; return (absx <= 0.5) ? 1.0 : 0.0; } /* triangle, Bartlett window, */ /* triangle function also known as lambda function */ /* 2nd order (linear) b-spline */ #define radius_triangle (1.0) static double filter_triangle(double const x) { double const absx = x < 0.0 ? -x : x; return absx < 1.0 ? 1.0 - absx : 0.0; } /* 3rd order (quadratic) b-spline */ #define radius_quadratic (1.5) static double filter_quadratic(double const x) { double const absx = x < 0.0 ? -x : x; return absx < 0.5 ? 0.75 - pow2(absx) : absx < 1.5 ? 0.50 * pow2(absx - 1.5) : 0.0; } /* 4th order (cubic) b-spline */ #define radius_cubic (2.0) static double filter_cubic(double const x) { double const absx = x < 0.0 ? -x : x; return absx < 1.0 ? 0.5 * pow3(absx) - pow2(absx) + 2.0/3.0 : absx < 2.0 ? pow3(2.0-absx)/6.0 : 0.0; } /* Catmull-Rom spline, Overhauser spline */ #define radius_catrom (2.0) static double filter_catrom(double const x) { double const absx = x < 0.0 ? -x : x; return absx < 1.0 ? 1.5 * pow3(absx) - 2.5 * pow2(absx) + 1.0 : absx < 2.0 ? -0.5 * pow3(absx) + 2.5 * pow2(absx) - 4.0 * absx + 2.0 : 0.0; } /* Mitchell & Netravali's two-param cubic */ /* see Mitchell&Netravali, */ /* "Reconstruction Filters in Computer Graphics", SIGGRAPH 88 */ #define radius_mitchell (2.0) static double filter_mitchell(double x) { double const b = 1.0/3.0; double const c = 1.0/3.0; double const p0 = ( 6.0 - 2.0*b ) / 6.0; double const p2 = (-18.0 + 12.0*b + 6.0*c) / 6.0; double const p3 = ( 12.0 - 9.0*b - 6.0*c) / 6.0; double const q0 = ( 8.0*b + 24.0*c) / 6.0; double const q1 = ( - 12.0*b - 48.0*c) / 6.0; double const q2 = ( 6.0*b + 30.0*c) / 6.0; double const q3 = ( - b - 6.0*c) / 6.0; double const absx = x < 0.0 ? -x : x; return absx < 1.0 ? p3 * pow3(absx) + p2 * pow2(absx) + p0 : absx < 2.0 ? q3 * pow3(absx) + q2 * pow2(absx) + q1 * absx + q0 : 0.0; } /* Gaussian filter (infinite) */ #define radius_gauss (1.25) static double filter_gauss(double const x) { return exp(-2.0*pow2(x)) * sqrt(2.0/M_PI); } /* sinc, perfect lowpass filter (infinite) */ #define radius_sinc (4.0) static double filter_sinc(double const x) { /* Note: Some people say sinc(x) is sin(x)/x. Others say it's sin(PI*x)/(PI*x), a horizontal compression of the former which is zero at integer values. We use the latter, whose Fourier transform is a canonical rectangle function (edges at -1/2, +1/2, height 1). */ return x == 0.0 ? 1.0 : sin(M_PI*x)/(M_PI*x); } /* Bessel (for circularly symm. 2-d filt, infinite) */ /* See Pratt "Digital Image Processing" p. 97 for Bessel functions */ #define radius_bessel (3.2383) static double filter_bessel(double const x) { return x == 0.0 ? M_PI/4.0 : j1(M_PI * x) / (2.0 * x); } /* Hanning window (infinite) */ #define radius_hanning (1.0) static double filter_hanning(double const x) { return 0.5 * cos(M_PI * x) + 0.5; } /* Hamming window (infinite) */ #define radius_hamming (1.0) static double filter_hamming(double const x) { return 0.46 * cos(M_PI * x) + 0.54; } /* Blackman window (infinite) */ #define radius_blackman (1.0) static double filter_blackman(double const x) { return 0.5 * cos(M_PI * x) + 0.08 * cos(2.0 * M_PI * x) + 0.42; } /* parameterized Kaiser window (infinite) */ /* from Oppenheim & Schafer, Hamming */ #define radius_kaiser (1.0) /* modified zeroth order Bessel function of the first kind. */ static double bessel_i0(double const x) { int i; double sum, y, t; sum = 1.0; y = pow2(x)/4.0; t = y; for (i=2; t>EPSILON; ++i) { sum += t; t *= (double)y/pow2(i); } return sum; } static double filter_kaiser(double const x) { /* typically 4 < a < 9 */ /* param a trades off main lobe width (sharpness) */ /* for side lobe amplitude (ringing) */ double const a = 6.5; double const i0a = 1.0/bessel_i0(a); return i0a * bessel_i0(a * sqrt(1.0-pow2(x))); } /* normal distribution (infinite) */ /* Normal(x) = Gaussian(x/2)/2 */ #define radius_normal (1.0) static double filter_normal(double const x) { return exp(-pow2(x)/2.0) / sqrt(2.0*M_PI); } /* Hermite filter */ #define radius_hermite (1.0) static double filter_hermite(double const x) { /* f(x) = 2|x|^3 - 3|x|^2 + 1, -1 <= x <= 1 */ double const absx = x < 0.0 ? -x : x; return absx < 1.0 ? 2.0 * pow3(absx) - 3.0 * pow2(absx) + 1.0 : 0.0; } /* Lanczos filter */ #define radius_lanczos (3.0) static double filter_lanczos(double const x) { double const absx = x < 0.0 ? -x : x; return x < 3.0 ? filter_sinc(absx) * filter_sinc(absx/3.0) : 0.0; } typedef struct { const char *name; double (*function)(double); double radius; /* This is how far from the Y axis (on either side) the function has significant value. (You can use this to limit how much of your domain you bother to compute the function over). */ bool windowed; } filter; static filter Filters[] = { { "point", filter_box, radius_point, false }, { "box", filter_box, radius_box, false }, { "triangle", filter_triangle, radius_triangle, false }, { "quadratic", filter_quadratic, radius_quadratic, false }, { "cubic", filter_cubic, radius_cubic, false }, { "catrom", filter_catrom, radius_catrom, false }, { "mitchell", filter_mitchell, radius_mitchell, false }, { "gauss", filter_gauss, radius_gauss, false }, { "sinc", filter_sinc, radius_sinc, true }, { "bessel", filter_bessel, radius_bessel, true }, { "hanning", filter_hanning, radius_hanning, false }, { "hamming", filter_hamming, radius_hamming, false }, { "blackman", filter_blackman, radius_blackman, false }, { "kaiser", filter_kaiser, radius_kaiser, false }, { "normal", filter_normal, radius_normal, false }, { "hermite", filter_hermite, radius_hermite, false }, { "lanczos", filter_lanczos, radius_lanczos, false }, { NULL }, }; typedef double (*basicFunction_t)(double); /****************************/ /****************************/ /****** end of filters ******/ /****************************/ /****************************/ enum scaleType {SCALE_SEPARATE, SCALE_BOXFIT, SCALE_BOXFILL, SCALE_PIXELMAX}; /* This is a way of specifying the output dimensions. SCALE_SEPARATE means specify the horizontal and vertical scaling separately. One or both may be unspecified. SCALE_BOXFIT means specify height and width of a box, and the image must be scaled, preserving aspect ratio, to the largest size that will fit in the box. Some of the box may be empty. SCALE_BOXFILL means specify height and width of a box, and the image must be scaled, preserving aspect ratio, to the smallest size that completely fills the box. Some of the image may be outside the box. SCALE_PIXELMAX means specify the maximum number of pixels the result should have and scale preserving aspect ratio and maximizing image size. */ struct CmdlineInfo { /* All the information the user supplied in the command line, * in a form easy for the program to use. */ const char * inputFileName; /* Filespec of input file */ unsigned int reportonly; unsigned int nomix; basicFunction_t filterFunction; /* NULL if not using resample method */ basicFunction_t windowFunction; /* Meaningful only when filterFunction != NULL */ double filterRadius; /* Meaningful only when filterFunction != NULL */ enum scaleType scaleType; /* 'xsize' and 'ysize' are numbers of pixels. Their meaning depends upon 'scaleType'. for SCALE_BOXFIT and SCALE_BOXFILL, they are the box dimensions. For SCALE_SEPARATE, they are the separate dimensions, or zero to indicate unspecified. For SCALE_PIXELMAX, they are meaningless. */ unsigned int xsize; unsigned int ysize; /* 'xscale' and 'yscale' are meaningful only for scaleType == SCALE_SEPARATE and only where the corresponding xsize/ysize is unspecified. 0.0 means unspecified. */ float xscale; float yscale; /* 'pixels' is meaningful only for scaleType == SCALE_PIXELMAX */ unsigned int pixels; unsigned int linear; unsigned int verbose; }; static void lookupFilterByName(const char * const filtername, filter * const filterP) { unsigned int i; bool found; found = false; /* initial assumption */ for (i=0; Filters[i].name; ++i) { if (strcmp(filtername, Filters[i].name) == 0) { *filterP = Filters[i]; found = true; } } if (!found) { unsigned int i; char known_filters[1024]; strcpy(known_filters, ""); for (i = 0; Filters[i].name; ++i) { const char * const name = Filters[i].name; if (strlen(known_filters) + strlen(name) + 1 + 1 < sizeof(known_filters)) { strcat(known_filters, name); strcat(known_filters, " "); } } pm_error("No such filter as '%s'. Known filter names are: %s", filtername, known_filters); } } static void processFilterOptions(unsigned int const filterSpec, const char filterOpt[], unsigned int const windowSpec, const char windowOpt[], struct CmdlineInfo * const cmdlineP) { if (filterSpec) { filter baseFilter; lookupFilterByName(filterOpt, &baseFilter); cmdlineP->filterFunction = baseFilter.function; cmdlineP->filterRadius = baseFilter.radius; if (windowSpec) { filter windowFilter; lookupFilterByName(windowOpt, &windowFilter); if (cmdlineP->windowFunction == filter_box) cmdlineP->windowFunction = NULL; else cmdlineP->windowFunction = windowFilter.function; } else { /* Default for most filters is no window. Those that _require_ a window function get Blackman. */ if (baseFilter.windowed) cmdlineP->windowFunction = filter_blackman; else cmdlineP->windowFunction = NULL; } } else cmdlineP->filterFunction = NULL; } static void parseSizeParm(const char * const sizeString, const char * const description, unsigned int * const sizeP) { char * endptr; long int sizeLong; sizeLong = strtol(sizeString, &endptr, 10); if (strlen(sizeString) > 0 && *endptr != '\0') pm_error("%s size argument not an integer: '%s'", description, sizeString); else if (sizeLong > INT_MAX - 2) pm_error("%s size argument is too large " "for computations: %ld", description, sizeLong); else if (sizeLong <= 0) pm_error("%s size argument is not positive: %ld", description, sizeLong); else *sizeP = (unsigned int) sizeLong; } static void parseXyParms(int const argc, const char ** const argv, struct CmdlineInfo * const cmdlineP) { /* parameters are box width (columns), box height (rows), and optional filespec */ if (argc-1 < 2) pm_error("You must supply at least two parameters with " "-xyfit/xyfill/xysize: " "x and y dimensions of the bounding box."); else if (argc-1 > 3) pm_error("Too many arguments. With -xyfit/xyfill/xysize, " "you need 2 or 3 arguments."); else { parseSizeParm(argv[1], "horizontal", &cmdlineP->xsize); parseSizeParm(argv[2], "vertical", &cmdlineP->ysize); if (argc-1 < 3) cmdlineP->inputFileName = "-"; else cmdlineP->inputFileName = argv[3]; } } static void parseScaleParms(int const argc, const char ** const argv, struct CmdlineInfo * const cmdlineP) { /*---------------------------------------------------------------------------- Parse the parameters as a scale factor and optional filespec (e.g. 'pamscale .5' or 'pamscale .5 testimg.ppm'). -----------------------------------------------------------------------------*/ if (argc-1 < 1) pm_error("With no dimension options, you must supply at least " "one parameter: the scale factor."); else { cmdlineP->xscale = cmdlineP->yscale = atof(argv[1]); if (cmdlineP->xscale <= 0.0) pm_error("The scale parameter %s is not a positive number.", argv[1]); else { if (argc-1 < 2) cmdlineP->inputFileName = "-"; else { cmdlineP->inputFileName = argv[2]; if (argc-1 > 2) pm_error("Too many arguments. There are at most two " "arguments with this set of options: " "scale factor and input file name. " "You specified %u", argc-1); } } } } static void parseFilespecOnlyParms(int const argc, const char ** const argv, struct CmdlineInfo * const cmdlineP) { /* Only parameter allowed is optional filespec */ if (argc-1 < 1) cmdlineP->inputFileName = "-"; else cmdlineP->inputFileName = argv[1]; } static void parseCommandLine(int argc, const char ** 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. Note that the strings we return are stored in the storage that was passed to us as the argv array. We also trash *argv. --------------------------------------------------------------------------*/ optEntry * option_def; optStruct3 opt; /* Instructions to pm_optParseOptions3 on how to parse our options. */ unsigned int option_def_index; unsigned int xyfit, xyfill; int xsize, ysize, pixels; int reduce; float xscale, yscale; const char * filterOpt; const char * window; unsigned int filterSpec, windowSpec; unsigned int xscaleSpec, yscaleSpec, xsizeSpec, ysizeSpec; unsigned int pixelsSpec, reduceSpec; MALLOCARRAY_NOFAIL(option_def, 100); option_def_index = 0; /* incremented by OPTENT3 */ OPTENT3(0, "xsize", OPT_UINT, &xsize, &xsizeSpec, 0); OPTENT3(0, "width", OPT_UINT, &xsize, &xsizeSpec, 0); OPTENT3(0, "ysize", OPT_UINT, &ysize, &ysizeSpec, 0); OPTENT3(0, "height", OPT_UINT, &ysize, &ysizeSpec, 0); OPTENT3(0, "xscale", OPT_FLOAT, &xscale, &xscaleSpec, 0); OPTENT3(0, "yscale", OPT_FLOAT, &yscale, &yscaleSpec, 0); OPTENT3(0, "pixels", OPT_UINT, &pixels, &pixelsSpec, 0); OPTENT3(0, "reduce", OPT_UINT, &reduce, &reduceSpec, 0); OPTENT3(0, "xysize", OPT_FLAG, NULL, &xyfit, 0); OPTENT3(0, "xyfit", OPT_FLAG, NULL, &xyfit, 0); OPTENT3(0, "xyfill", OPT_FLAG, NULL, &xyfill, 0); OPTENT3(0, "verbose", OPT_FLAG, NULL, &cmdlineP->verbose, 0); OPTENT3(0, "filter", OPT_STRING, &filterOpt, &filterSpec, 0); OPTENT3(0, "window", OPT_STRING, &window, &windowSpec, 0); OPTENT3(0, "nomix", OPT_FLAG, NULL, &cmdlineP->nomix, 0); OPTENT3(0, "linear", OPT_FLAG, NULL, &cmdlineP->linear, 0); OPTENT3(0, "reportonly", OPT_FLAG, NULL, &cmdlineP->reportonly, 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 (cmdlineP->nomix && filterSpec) pm_error("You cannot specify both -nomix and -filter."); processFilterOptions(filterSpec, filterOpt, windowSpec, window, cmdlineP); if (xsizeSpec && xsize == 0) pm_error("-xsize/width must be greater than zero."); if (ysizeSpec && ysize == 0) pm_error("-ysize/height must be greater than zero."); if (xscaleSpec && xscale <= 0.0) pm_error("-xscale must be greater than zero."); if (yscaleSpec && yscale <= 0.0) pm_error("-yscale must be greater than zero."); if (reduceSpec && reduce <= 0) pm_error("-reduce must be greater than zero."); if (xsizeSpec && xscaleSpec) pm_error("Cannot specify both -xsize/width and -xscale."); if (ysizeSpec && yscaleSpec) pm_error("Cannot specify both -ysize/height and -yscale."); if ((xyfit || xyfill) && (xsizeSpec || xscaleSpec || ysizeSpec || yscaleSpec || reduceSpec || pixelsSpec) ) pm_error("Cannot specify -xyfit/xyfill/xysize with other " "dimension options."); if (xyfit && xyfill) pm_error("Cannot specify both -xyfit and -xyfill"); if (pixelsSpec && (xsizeSpec || xscaleSpec || ysizeSpec || yscaleSpec || reduceSpec) ) pm_error("Cannot specify -pixels with other dimension options."); if (reduceSpec && (xsizeSpec || xscaleSpec || ysizeSpec || yscaleSpec) ) pm_error("Cannot specify -reduce with other dimension options."); if (pixelsSpec && pixels == 0) pm_error("-pixels must be greater than zero"); /* Get the program parameters */ if (xyfit || xyfill) { cmdlineP->scaleType = xyfit ? SCALE_BOXFIT : SCALE_BOXFILL; parseXyParms(argc, argv, cmdlineP); } else if (reduceSpec) { cmdlineP->scaleType = SCALE_SEPARATE; parseFilespecOnlyParms(argc, argv, cmdlineP); cmdlineP->xsize = cmdlineP->ysize = 0; cmdlineP->xscale = cmdlineP->yscale = ((double) 1.0) / ((double) reduce); pm_message("reducing by %d gives scale factor of %f.", reduce, cmdlineP->xscale); } else if (pixelsSpec) { cmdlineP->scaleType = SCALE_PIXELMAX; parseFilespecOnlyParms(argc, argv, cmdlineP); cmdlineP->pixels = pixels; } else if (xsizeSpec || xscaleSpec || ysizeSpec || yscaleSpec) { cmdlineP->scaleType = SCALE_SEPARATE; parseFilespecOnlyParms(argc, argv, cmdlineP); cmdlineP->xsize = xsizeSpec ? xsize : 0; cmdlineP->ysize = ysizeSpec ? ysize : 0; cmdlineP->xscale = xscaleSpec ? xscale : 0.0; cmdlineP->yscale = yscaleSpec ? yscale : 0.0; } else { cmdlineP->scaleType = SCALE_SEPARATE; parseScaleParms(argc, argv, cmdlineP); cmdlineP->xsize = cmdlineP->ysize = 0; } } static void computeOutputDimensions(struct CmdlineInfo const cmdline, unsigned int const cols, unsigned int const rows, int * const newcolsP, int * const newrowsP) { double newcolsD, newrowsD; /* Intermediate calculation of the output dimensions, in double precision floating point to avoid arithmetic overflow. */ unsigned int newcols, newrows; /* The output dimensions we return */ switch(cmdline.scaleType) { case SCALE_PIXELMAX: { if (rows * cols <= cmdline.pixels) { newrowsD = rows; newcolsD = cols; } else { const double scale = sqrt( (float) cmdline.pixels / ((float) cols * (float) rows)); newrowsD = rows * scale; newcolsD = cols * scale; } } break; case SCALE_BOXFIT: case SCALE_BOXFILL: { double const aspect_ratio = (float) cols / (float) rows; double const box_aspect_ratio = (float) cmdline.xsize / (float) cmdline.ysize; if ((box_aspect_ratio > aspect_ratio && cmdline.scaleType == SCALE_BOXFIT) || (box_aspect_ratio < aspect_ratio && cmdline.scaleType == SCALE_BOXFILL)) { newrowsD = cmdline.ysize; newcolsD = newrowsD * aspect_ratio; } else { newcolsD = cmdline.xsize; newrowsD = newcolsD / aspect_ratio; } } break; case SCALE_SEPARATE: { if (cmdline.xsize) newcolsD = cmdline.xsize; else if (cmdline.xscale) newcolsD = cmdline.xscale * cols; else if (cmdline.ysize) newcolsD = cols * ((float) cmdline.ysize/rows); else newcolsD = cols; if (cmdline.ysize) newrowsD = cmdline.ysize; else if (cmdline.yscale) newrowsD = cmdline.yscale * rows; else if (cmdline.xsize) newrowsD = rows * ((float) cmdline.xsize/cols); else newrowsD = rows; } } /* If the rounding yields a zero dimension, we fudge it up to 1. We do this rather than considering it a specification error (and dying) because it's friendlier to automated processes that work on arbitrary input. It saves them having to check their numbers to avoid catastrophe. */ newcols = MAX(1, ROUNDU(newcolsD)); newrows = MAX(1, ROUNDU(newrowsD)); if (newcols > INT_MAX - 2) pm_error("output image width (%u) too large for computations", newcols); if (newrows > INT_MAX - 2) pm_error("output image height (%u) too large for computation", newrows); *newcolsP = newcols; *newrowsP = newrows; } /****************************/ /****************************/ /******* resampling *********/ /****************************/ /****************************/ /* The resample code was inspired by Paul Heckbert's zoom program. ** http://www.cs.cmu.edu/~ph/zoom */ struct filterFunction { /*---------------------------------------------------------------------------- A function to convolve with the samples. -----------------------------------------------------------------------------*/ basicFunction_t basicFunction; /* The basic shape of the function. Its horizontal scale is designed to filter out frequencies greater than 1. */ basicFunction_t windowFunction; /* A function to multiply by basicFunction(). NULL if none. */ double windowScaler; /* Factor by which to compress windowFunction() horizontally */ double horizontalScaler; /* Factor by which to compress basicFunction() * windowFunction horizontally. Note that compressing horizontally in the sample domain is equivalent to expanding horizontally (and shrinking vertically) in the frequency domain. I.e. values less than unity have the effect of chopping out high frequencies. */ double radius; /* A final filter. filterFunction(x) is zero for |x| > radius regardless of what the rest of the members say. Implementation note: This is important because windowFunction(), out of laziness, doesn't do the whole job of windowing. It is not zero beyond the cutoff points as it should be. If not for that, radius would only be a hint to describe what the other members already do, so the convolver knows where to stop. */ }; typedef struct { /* A term of the linear combination of input rows that makes up an output row. I.e. an input row and its weight. Alternatively, the analogous thing for a column. */ int position; /* Row/column number in the input image */ double weight; /* Weight to be given to that row/col. In [0, 1]. */ } WEIGHT; typedef struct { /* A description of the linear combination of input rows that generates a particular output row. An output row is a weighted average of some input rows. E.g. Row 2 of the output might be composed of 50% of Row 2 of the input and 50% of Row 3 of the input. Alternatively, the analogous thing for columns. */ unsigned int nWeight; /* Number of elements in 'Weight'. They're consecutive, starting at index 0. */ unsigned int allocWeight; /* Number of allocated frames in 'Weight' */ WEIGHT *Weight; /* The terms of the linear combination. Has 'nWeight' elements. The coefficients (weights) of each add up to unity. */ } WLIST; typedef struct { /* This identifies a row of the input image. */ int rowNumber; /* The row number in the input image of the row. -1 means no row. */ tuple *tuplerow; /* The tuples of the row. If rowNumber = -1, these are arbitrary, but allocated, tuples. */ } SCANLINE; typedef struct { /* A vertical window of a raster */ int width; /* Width of the window, in columns */ int height; /* Height of the window, in rows */ SCANLINE *line; /* An array of 'height' elements, malloced. This identifies the lines of the input image that compose the window. The index order is NOT the order of the rows in the image. E.g. line[0] isn't always the topmost row of the window. Rather, the rows are arranged in a cycle and you have to know independently where the topmost one is. E.g. the rows of a 5 line window with topmost row at index 3 might be: line[0] = Row 24 line[1] = Row 25 line[2] = Row 26 line[3] = Row 22 line[4] = Row 23 */ } SCAN; static int appendWeight(WLIST * const WList, int const index, double const weight) { /*---------------------------------------------------------------------------- Add a weighting of 'weight' for index 'index' to the weight list 'WList'. -----------------------------------------------------------------------------*/ if (weight == 0.0) { /* A weight of 0 in the list is redundant, so we don't add it. A weight entry says "Add W fraction of the pixel at index I," so where W is 0, it's the same as not having the entry at all. */ } else { unsigned int const n = WList->nWeight; assert(WList->allocWeight >= n+1); WList->Weight[n].position = index; WList->Weight[n].weight = weight; ++WList->nWeight; } return 0; } static sample unnormalize(double const normalized, sample const maxval) { /* Take care here, the conversion of any floating point value <= -1.0 to an unsigned type is _undefined_. See ISO 9899:1999 section 6.3.1.4. Not only is it undefined it also does the wrong thing in actual practice, EG on Darwin PowerPC (my iBook running OS X) negative values clamp to maxval. We get negative values because some of the filters (EG catrom) have negative weights. */ return MIN(maxval, (sample)(MAX(0.0, (normalized*maxval + 0.5)))); } static void initWeightList(WLIST * const weightListP, unsigned int const maxWeights) { weightListP->nWeight = 0; weightListP->allocWeight = maxWeights; MALLOCARRAY(weightListP->Weight, maxWeights); if (weightListP->Weight == NULL) pm_error("Out of memory allocating a %u-element weight list.", maxWeights); } static void createWeightList(unsigned int const targetPos, unsigned int const sourceSize, double const scale, struct filterFunction filter, WLIST * const weightListP) { /*---------------------------------------------------------------------------- Create a weight list for computing target pixel number 'targetPos' from a set of source pixels. These pixels are a line of pixels either horizontally or vertically. The weight list is a list of weights to give each source pixel in the set. The source pixel set is a window of source pixels centered on some point. The weights are defined by the function 'filter' of the position within the window, and normalized to add up to 1.0. Technically, the window is infinite, but we know that the filter function is zero beyond a certain distance from the center of the window. For example, assume 'targetPos' is 5. That means we're computing weights for either Column 5 or Row 5 of the target image. Assume it's Column 5. Assume 'radius' is 1. That means a window of two pixels' worth of a source row determines the color of the Column 5 pixel of a target row. Assume 'filter' is a triangle function -- 1 at 0, sloping down to 0 at -1 and 1. Now assume that the scale factor is 2 -- the target image will be twice the size of the source image. That means the two-pixel-wide window of the source row that affects Column 5 of the target row, which is centered at target position 5.5, is centered at source position 5.5/2 = 2.75. So it goes from source position 1.75 to 3.75. That means the window covers 1/4 of Column 1, all of Column 2, and 3/4 of Column 3 of the source row. We want to calculate 3 weights, one to be applied to each source pixel in computing the target pixel. Ideally, we would compute the average height of the filter function over each source pixel region. But that's too hard. So we approximate by assuming that the filter function is constant within each region, at the value the function has at the _center_ of the region. So for the Column 1 region, which goes from 1.75 to 2.00, centered -.875 from the center of the window, we assume a constant function value of triangle(-.875), which equals .125. For the 2.00-3.00 region, we get triangle(-.25) = .75. For the 3.00-3.75 region, we get triangle(.625) = .375. So the weights for the 3 regions, which we get by multiplying this constant function value by the width of the region and normalizing so they add up to 1 are: Source Column 1: .125*.25 / 1.0625 = .029 Source Column 2: .75*1.00 / 1.0625 = .706 Source Column 3: .375*.75 / 1.0625 = .265 These are the weights we return. Thus, if we assume that the source pixel 1 has value 10, source pixel 2 has value 20, and source pixel 3 has value 30, Caller would compute target pixel 5's value as 10*.029 + 20*.706 + 30*.265 = 22 -----------------------------------------------------------------------------*/ /* 'windowCenter', is the continuous position within the source of the center of the window that will influence target pixel 'targetPos'. 'left' and 'right' are the edges of the window. 'leftPixel' and 'rightPixel' are the pixel positions of the pixels at the edges of that window. Note that if we're doing vertical weights, "left" and "right" mean top and bottom. */ double const windowCenter = ((double)targetPos + 0.5) / scale; double left = MAX(0.0, windowCenter - filter.radius - EPSILON); unsigned int const leftPixel = floor(left); double right = MIN((double)sourceSize - EPSILON, windowCenter + filter.radius + EPSILON); unsigned int const rightPixel = floor(right); double norm; unsigned int j; initWeightList(weightListP, rightPixel - leftPixel + 1); /* calculate weights */ norm = 0.0; /* initial value */ for (j = leftPixel; j <= rightPixel; ++j) { /* Calculate the weight that source pixel 'j' will have in the value of target pixel 'targetPos'. */ double const regionLeft = MAX(left, (double)j); double const regionRight = MIN(right, (double)(j + 1)); double const regionWidth = regionRight - regionLeft; double const regionCenter = (regionRight + regionLeft) / 2; double const dist = regionCenter - windowCenter; double weight; weight = filter.basicFunction(filter.horizontalScaler * dist); if (filter.windowFunction) weight *= filter.windowFunction( filter.horizontalScaler * filter.windowScaler * dist); assert(regionWidth <= 1.0); weight *= regionWidth; norm += weight; appendWeight(weightListP, j, weight); } if (norm == 0.0) pm_error("INTERNAL ERROR: No source pixels contribute to target " "pixel %u", targetPos); /* normalize the weights so they add up to 1.0 */ if (norm != 1.0) { unsigned int n; for (n = 0; n < weightListP->nWeight; ++n) { weightListP->Weight[n].weight /= norm; } } } static void createWeightListSet(unsigned int const sourceSize, unsigned int const targetSize, struct filterFunction const filterFunction, WLIST ** const weightListSetP) { /*---------------------------------------------------------------------------- Create the set of weight lists that will effect the resample. This is where the actual work of resampling gets done. The weight list set is a bunch of factors one can multiply by the pixels in a region to effect a resampling. Multiplying by these factors effects all of the following transformations on the original pixels: 1) Filter out any frequencies that are artifacts of the original sampling. We assume a perfect sampling was done, which means the original continuous dataset had a maximum frequency of 1/2 of the original sample rate and anything above that is an artifact of the sampling. So we filter out anything above 1/2 of the original sample rate (sample rate == pixel resolution). 2) Filter out any frequencies that are too high to be captured by the new sampling -- i.e. frequencies above 1/2 the new sample rate. This is the information we must lose because of low sample rate. 3) Sample the result at the new sample rate. We do all three of these steps in a single convolution of the original pixels. Steps (1) and (2) can be combined into a single frequency domain rectangle function. A frequency domain rectangle function is a pixel domain sinc function, which is what we assume 'filterFunction' is. We get Step 3 by computing the convolution only at the new sample points. I don't know what any of this means when 'filterFunction' is not sinc. Maybe it just means some approximation or additional filtering steps are happening. -----------------------------------------------------------------------------*/ double const scale = (double)targetSize / sourceSize; WLIST *weightListSet; /* malloc'ed */ unsigned int targetPos; MALLOCARRAY_NOFAIL(weightListSet, targetSize); for (targetPos = 0; targetPos < targetSize; ++targetPos) createWeightList(targetPos, sourceSize, scale, filterFunction, &weightListSet[targetPos]); *weightListSetP = weightListSet; } static struct filterFunction makeFilterFunction(double const scale, basicFunction_t basicFunction, double const basicRadius, basicFunction_t windowFunction) { /*---------------------------------------------------------------------------- Create a function to convolve with the samples (so it isn't actually a filter function, but the Fourier transform of a filter function. A filter function is something you multiply by in the frequency domain) to create a function from which one can resample. Convolving with this function will achieve two goals: 1) filter out high frequencies that are artifacts of the original sampling (i.e. the turning of a continuous function into a staircase function); 2) filter out frequencies higher than half the resample rate, so that the resample will be a perfect sampling of it, and not have aliasing. To make the calculation even more efficient, we take advantage of the fact that the weight list doesn't depend on the particular old and new sample rates at all except -- all that's important is their ratio (which is 'scale'). So we assume the original sample rate is 1 and the new sample rate is 'scale'. -----------------------------------------------------------------------------*/ double const freqLimit = MIN(1.0, scale); /* We're going to cut out any frequencies above this, to accomplish Steps (1) and (2) above. */ struct filterFunction retval; retval.basicFunction = basicFunction; retval.windowFunction = windowFunction; retval.horizontalScaler = freqLimit; /* Our 'windowFunction' argument is a function normalized to the domain (-1, 1). We need to scale it horizontally to fit the basic filter function. We assume the radius of the filter function is the area to which the window should fit (i.e. zero beyond the radius, nonzero inside the radius). But that's really a misuse of radius, because radius is supposed to be just the distance beyond which we can assume for convenience that the filter function is zero, possibly giving up some precision. But note that 'windowFunction' isn't zero outside (-1, 1), even though the actual window function is supposed to be. Hence, scaling the window function exactly to the radius stops our calculations from noticing the wrong values outside (-1, 1) -- we'll never use them. */ retval.windowScaler = 1/basicRadius; retval.radius = basicRadius / retval.horizontalScaler; return retval; } static void destroyWeightListSet(WLIST * const weightListSet, unsigned int const size) { unsigned int i; for (i = 0; i < size; ++i) free(weightListSet[i].Weight); free(weightListSet); } static void createScanBuf(struct pam * const pamP, double const maxRowWeights, bool const verbose, SCAN * const scanbufP) { SCAN scanbuf; unsigned int lineNumber; scanbuf.width = pamP->width; scanbuf.height = maxRowWeights; MALLOCARRAY_NOFAIL(scanbuf.line, scanbuf.height); for (lineNumber = 0; lineNumber < scanbuf.height; ++lineNumber) { scanbuf.line[lineNumber].rowNumber = -1; scanbuf.line[lineNumber].tuplerow = pnm_allocpamrow(pamP); } if (verbose) pm_message("scanline buffer: %d lines of %d pixels", scanbuf.height, scanbuf.width); *scanbufP = scanbuf; } static void destroyScanbuf(SCAN const scanbuf) { unsigned int lineNumber; for (lineNumber = 0; lineNumber < scanbuf.height; ++lineNumber) pnm_freepamrow(scanbuf.line[lineNumber].tuplerow); free(scanbuf.line); } static void resampleDimensionMessage(struct pam * const inpamP, struct pam * const outpamP) { pm_message ("resampling from %d*%d to %d*%d (%f, %f)", inpamP->width, inpamP->height, outpamP->width, outpamP->height, (double)outpamP->width/inpamP->width, (double)outpamP->height/inpamP->height); } static void addInPixel(const struct pam * const pamP, tuple const tuple, float const weight, bool const haveOpacity, unsigned int const opacityPlane, double * const accum) { /*---------------------------------------------------------------------------- Add into accum[] the values from the tuple 'tuple', weighted by 'weight'. accum[P] is the accumulated normalized sample value for Plane P. Iff 'haveOpacity', Plane 'opacityPlane' of the tuple is an opacity (alpha, transparency) plane. -----------------------------------------------------------------------------*/ unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) { double const normalizedSample = (double)tuple[plane]/pamP->maxval; double opacityAdjustment; if (haveOpacity && plane != opacityPlane) opacityAdjustment = (double)tuple[opacityPlane]/pamP->maxval; else opacityAdjustment = 1; accum[plane] += opacityAdjustment * normalizedSample * weight; } } static void generateOutputTuple(const struct pam * const pamP, double const accum[], bool const haveOpacity, unsigned int const opacityPlane, tuple * const tupleP) { /*---------------------------------------------------------------------------- Convert the values accum[] accumulated for a pixel by outputOneResampledRow() to a bona fide PAM tuple as *tupleP, as described by *pamP. accum[P] is the pixel's plane P value normalized (i.e. in range 0..1). -----------------------------------------------------------------------------*/ unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) { float opacityAdjustedSample; /* normalized sample value */ if (haveOpacity && plane != opacityPlane) { if (accum[opacityPlane] < EPSILON) { opacityAdjustedSample = 0.0; } else opacityAdjustedSample = accum[plane] / accum[opacityPlane]; } else opacityAdjustedSample = accum[plane]; (*tupleP)[plane] = unnormalize(opacityAdjustedSample, pamP->maxval); } } static void outputOneResampledRow(const struct pam * const outpamP, SCAN const scanbuf, WLIST const YW, const WLIST * const XWeight, tuple * const line, double * const accum) { /*---------------------------------------------------------------------------- From the data in 'scanbuf' and weights in 'YW' and 'XWeight', generate one output row for the image described by *outpamP and output it. An output pixel is a weighted average of the pixels in a certain rectangle of the input. 'YW' and 'XWeight' describe those weights for each column of the row we are to output. 'line' and 'accum' are just working space that Caller provides us with to save us the time of allocating it. 'line' is at least big enough to hold an output row; 'accum' is at least outpamP->depth big. -----------------------------------------------------------------------------*/ unsigned int col; int haveOpacity; /* There is an opacity plane */ unsigned int opacityPlane; /* Plane number of opacity plane, if any */ pnm_getopacity(outpamP, &haveOpacity, &opacityPlane); /* We accumulate intensity in accum[], where accum[P] is the intensity for Plane P. These are normalized values (i.e. in the range 0..1 */ for (col = 0; col < outpamP->width; ++col) { WLIST const XW = XWeight[col]; unsigned int i; { unsigned int plane; for (plane = 0; plane < outpamP->depth; ++plane) accum[plane] = 0.0; } for (i = 0; i < YW.nWeight; ++i) { int const yp = YW.Weight[i].position; float const yw = YW.Weight[i].weight; int const slot = yp % scanbuf.height; unsigned int j; for (j = 0; j < XW.nWeight; ++j) { int const xp = XW.Weight[j].position; tuple const tuple = scanbuf.line[slot].tuplerow[xp]; addInPixel(outpamP, tuple, yw * XW.Weight[j].weight, haveOpacity, opacityPlane, accum); } } generateOutputTuple(outpamP, accum, haveOpacity, opacityPlane, &line[col]); } pnm_writepamrow(outpamP, line); } static bool scanbufContainsTheRows(SCAN const scanbuf, WLIST const rowWeights) { /*---------------------------------------------------------------------------- Return true iff scanbuf 'scanbuf' contains every row mentioned in 'rowWeights'. It might contain additional rows besides. -----------------------------------------------------------------------------*/ bool missingRow; unsigned int i; for (i = 0, missingRow = false; i < rowWeights.nWeight && !missingRow; ++i) { unsigned int const inputRow = rowWeights.Weight[i].position; unsigned int const slot = inputRow % scanbuf.height; /* This is the number of the slot in the scanbuf that would have the input row in question if the scanbuf has the row at all. */ if (scanbuf.line[slot].rowNumber != inputRow) { /* Nope, this slot has some other row or no row at all. So the row we're looking for isn't in the scanbuf. */ missingRow = true; } } return !missingRow; } static void createWeightLists(struct pam * const inpamP, struct pam * const outpamP, basicFunction_t const filterFunction, double const filterRadius, basicFunction_t const windowFunction, WLIST ** const horizWeightP, WLIST ** const vertWeightP, unsigned int * const maxRowWeightsP) { /*---------------------------------------------------------------------------- This is the function that actually does the resampling. Note that it does it without ever looking at the source or target pixels! It produces a simple set of numbers that Caller can blindly apply to the source pixels to get target pixels. -----------------------------------------------------------------------------*/ struct filterFunction horizFilter, vertFilter; horizFilter = makeFilterFunction( (double)outpamP->width/inpamP->width, filterFunction, filterRadius, windowFunction); createWeightListSet(inpamP->width, outpamP->width, horizFilter, horizWeightP); vertFilter = makeFilterFunction( (double)outpamP->height/inpamP->height, filterFunction, filterRadius, windowFunction); createWeightListSet(inpamP->height, outpamP->height, vertFilter, vertWeightP); *maxRowWeightsP = ceil(2.0*(vertFilter.radius+EPSILON) + 1 + EPSILON); } static void resample(struct pam * const inpamP, struct pam * const outpamP, basicFunction_t const filterFunction, double const filterRadius, basicFunction_t const windowFunction, bool const verbose, bool const linear) { /*--------------------------------------------------------------------------- Resample the image in the input file, described by *inpamP, so as to create the image in the output file, described by *outpamP. Input and output differ by height, width, and maxval only. Use the resampling filter function 'filterFunction', applied over radius 'filterRadius'. The input file is positioned past the header, to the beginning of the raster. The output file is too. ---------------------------------------------------------------------------*/ int inputRow, outputRow; WLIST * horizWeight; WLIST * vertWeight; SCAN scanbuf; unsigned int maxRowWeights; tuple * line; /* This is just work space for outputOneResampledRow() */ double * weight; /* This is just work space for outputOneResampledRow() */ if (linear) pm_error("You cannot use the resampling scaling method on " "linear input."); createWeightLists(inpamP, outpamP, filterFunction, filterRadius, windowFunction, &horizWeight, &vertWeight, &maxRowWeights); createScanBuf(inpamP, maxRowWeights, verbose, &scanbuf); if (verbose) resampleDimensionMessage(inpamP, outpamP); line = pnm_allocpamrow(outpamP); MALLOCARRAY_NOFAIL(weight, outpamP->depth); outputRow = 0; for (inputRow = 0; inputRow < inpamP->height; ++inputRow) { bool needMoreInput; /* We've output as much as we can using the rows that are in the scanbuf; it's time to move the window. Or fill it in the first place. */ unsigned int scanbufSlot; /* Read source row; add it to the scanbuf */ scanbufSlot = inputRow % scanbuf.height; scanbuf.line[scanbufSlot].rowNumber = inputRow; pnm_readpamrow(inpamP, scanbuf.line[scanbufSlot].tuplerow); /* Output all the rows we can make out of the current contents of the scanbuf. Might be none. */ needMoreInput = false; /* initial assumption */ while (outputRow < outpamP->height && !needMoreInput) { WLIST const rowWeights = vertWeight[outputRow]; /* The description of what makes up our current output row; i.e. what fractions of which input rows combine to create this output row. */ assert(rowWeights.nWeight <= scanbuf.height); if (scanbufContainsTheRows(scanbuf, rowWeights)) { outputOneResampledRow(outpamP, scanbuf, rowWeights, horizWeight, line, weight); ++outputRow; } else needMoreInput = true; } } if (outputRow != outpamP->height) pm_error("INTERNAL ERROR: assembled only %u of the required %u " "output rows.", outputRow, outpamP->height); pnm_freepamrow(line); destroyScanbuf(scanbuf); destroyWeightListSet(horizWeight, outpamP->width); destroyWeightListSet(vertWeight, outpamP->height); } /****************************/ /****************************/ /**** end of resampling *****/ /****************************/ /****************************/ static void zeroNewRow(struct pam * const pamP, tuplen * const tuplenrow) { unsigned int col; for (col = 0; col < pamP->width; ++col) { unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) tuplenrow[col][plane] = 0.0; } } static void accumOutputCol(struct pam * const pamP, tuplen const intuplen, float const fraction, tuplen const accumulator) { /*---------------------------------------------------------------------------- Add fraction 'fraction' of the pixel indicated by 'intuplen' to the pixel accumulator 'accumulator'. 'intuplen' and 'accumulator' are not a standard libnetpbm tuplen. It is proportional to light intensity. The foreground color component samples are proportional to light intensity, and have opacity factored in. -----------------------------------------------------------------------------*/ unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) accumulator[plane] += fraction * intuplen[plane]; } static void horizontalScale(tuplen * const inputtuplenrow, tuplen * const newtuplenrow, struct pam * const inpamP, struct pam * const outpamP, float const xscale, float * const stretchP) { /*---------------------------------------------------------------------------- Take the input row 'inputtuplenrow', described by *inpamP, and scale it by a factor of 'xscale', to create the output row 'newtuplenrow', described by *outpamP. Because of arithmetic imprecision, we may have to stretch slightly the contents of the last pixel of the output row to make a full pixel. Return as *stretchP the fraction of a pixel by which we had to stretch in this way. Assume maxval and depth of input and output are the same. -----------------------------------------------------------------------------*/ float fraccoltofill, fraccolleft; unsigned int col; unsigned int newcol; newcol = 0; fraccoltofill = 1.0; /* Output column is "empty" now */ zeroNewRow(outpamP, newtuplenrow); for (col = 0; col < inpamP->width; ++col) { /* Process one tuple from input ('inputtuplenrow') */ fraccolleft = xscale; /* Output all columns, if any, that can be filled using information from this input column, in addition to what's already in the output column. */ while (fraccolleft >= fraccoltofill) { /* Generate one output pixel in 'newtuplerow'. It will consist of anything accumulated from prior input pixels in accumulator[], plus a fraction of the current input pixel. */ assert(newcol < outpamP->width); accumOutputCol(inpamP, inputtuplenrow[col], fraccoltofill, newtuplenrow[newcol]); fraccolleft -= fraccoltofill; /* Set up to start filling next output column */ ++newcol; fraccoltofill = 1.0; } /* There's not enough left in the current input pixel to fill up a whole output column, so just accumulate the remainder of the pixel into the current output column. Because of rounding, we may have a tiny bit of pixel left and have run out of output pixels. In that case, we throw away what's left. */ if (fraccolleft > 0.0 && newcol < outpamP->width) { accumOutputCol(inpamP, inputtuplenrow[col], fraccolleft, newtuplenrow[newcol]); fraccoltofill -= fraccolleft; } } if (newcol < outpamP->width-1 || newcol > outpamP->width) pm_error("Internal error: last column filled is %d, but %d " "is the rightmost output column.", newcol, outpamP->width-1); if (newcol < outpamP->width) { /* We were still working on the last output column when we ran out of input columns. This would be because of rounding down, and we should be missing only a tiny fraction of that last output column. Just fill in the missing color with the color of the rightmost input pixel. */ accumOutputCol(inpamP, inputtuplenrow[inpamP->width-1], fraccoltofill, newtuplenrow[newcol]); *stretchP = fraccoltofill; } else *stretchP = 0.0; } static void zeroAccum(struct pam * const pamP, tuplen * const accumulator) { unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) { unsigned int col; for (col = 0; col < pamP->width; ++col) accumulator[col][plane] = 0.0; } } static void accumOutputRow(struct pam * const pamP, tuplen * const tuplenrow, float const fraction, tuplen * const accumulator) { /*---------------------------------------------------------------------------- Take 'fraction' times the samples in row 'tuplenrow' and add it to 'accumulator' in the same way as accumOutputCol(). 'fraction' is less than 1.0. -----------------------------------------------------------------------------*/ unsigned int plane; for (plane = 0; plane < pamP->depth; ++plane) { unsigned int col; for (col = 0; col < pamP->width; ++col) accumulator[col][plane] += fraction * tuplenrow[col][plane]; } } static void readARow(struct pam * const pamP, tuplen * const tuplenRow, const pnm_transformMap * const transform) { /*---------------------------------------------------------------------------- Read a row from the input file described by *pamP, as values (for the foreground color) proportional to light intensity, with opacity included. By contrast, a simple libnetpbm read would give the same numbers you find in a PAM: gamma-adjusted values for the foreground color component and scaled as if opaque. The latter means that full red would have a red intensity of 1.0 even if the pixel is only 75% opaque. We, on the other hand, would return red intensity of .75 in that case. The opacity plane we return is the same as a simple libnetpbm read would return. We ASSUME that the transform 'transform' is that necessary to effect the conversion to intensity-linear values and normalize. If it is NULL, we ASSUME that they already are intensity-proportional and just need to be normalized. -----------------------------------------------------------------------------*/ tuple * tupleRow; tupleRow = pnm_allocpamrow(pamP); pnm_readpamrow(pamP, tupleRow); pnm_normalizeRow(pamP, tupleRow, transform, tuplenRow); pnm_applyopacityrown(pamP, tuplenRow); pnm_freepamrow(tupleRow); } static void writeARow(struct pam * const pamP, tuplen * const tuplenRow, const pnm_transformMap * const transform) { /*---------------------------------------------------------------------------- Write a row to the output file described by *pamP, from values proportional to light intensity with opacity included (i.e. the same kind of number you would get form readARow()). We ASSUME that the transform 'transform' is that necessary to effect the conversion to brightness-linear unnormalized values. If it is NULL, we ASSUME that they already are brightness-proportional and just need to be unnormalized. We destroy *tuplenRow in the process. -----------------------------------------------------------------------------*/ tuple * tupleRow; tupleRow = pnm_allocpamrow(pamP); pnm_unapplyopacityrown(pamP, tuplenRow); pnm_unnormalizeRow(pamP, tuplenRow, transform, tupleRow); pnm_writepamrow(pamP, tupleRow); pnm_freepamrow(tupleRow); } static void issueStretchWarning(bool const verbose, double const fracrowtofill) { /* We need another input row to fill up this output row, but there aren't any more. That's because of rounding down on our scaling arithmetic. So we go ahead with the data from the last row we read, which amounts to stretching out the last output row. */ if (verbose) pm_message("%f of bottom row stretched because of " "arithmetic imprecision", fracrowtofill); } static void scaleHorizontallyAndOutputRow(struct pam * const inpamP, struct pam * const outpamP, tuplen * const rowAccumulator, const pnm_transformMap * const transform, tuplen * const newtuplenrow, float const xscale, unsigned int const row, bool const verbose) { /*---------------------------------------------------------------------------- Scale the row in 'rowAccumulator' horizontally by factor 'xscale' and output it. 'newtuplenrow' is work space Caller provides us. It is at least wide enough to hold one output row. -----------------------------------------------------------------------------*/ if (outpamP->width == inpamP->width) /* shortcut X scaling */ writeARow(outpamP, rowAccumulator, transform); /* This destroys 'rowAccumulator' */ else { float stretch; horizontalScale(rowAccumulator, newtuplenrow, inpamP, outpamP, xscale, &stretch); if (verbose && row == 0) pm_message("%f of right column stretched because of " "arithmetic imprecision", stretch); writeARow(outpamP, newtuplenrow, transform); /* This destroys 'newtuplenrow' */ } } static void createTransforms(struct pam * const inpamP, struct pam * const outpamP, bool const assumeLinear, const pnm_transformMap ** const inputTransformP, const pnm_transformMap ** const outputTransformP) { if (assumeLinear) { *inputTransformP = NULL; *outputTransformP = NULL; } else { *inputTransformP = pnm_createungammatransform(inpamP); *outputTransformP = pnm_creategammatransform(outpamP); } } static void destroyTransforms(const pnm_transformMap * const inputTransform, const pnm_transformMap * const outputTransform) { if (inputTransform) free((void*)inputTransform); if (outputTransform) free((void*)outputTransform); } static void scaleWithMixing(struct pam * const inpamP, struct pam * const outpamP, float const xscale, float const yscale, bool const assumeLinear, bool const verbose) { /*---------------------------------------------------------------------------- Scale the image described by *inpamP by xscale horizontally and yscale vertically and write the result as the image described by *outpamP. The input file is positioned past the header, to the beginning of the raster. The output file is too. Mix colors from input rows together in the output rows. 'assumeLinear' means to assume that the sample values in the input image vary from standard PAM in that they are proportional to intensity, (This makes the computation a lot faster, so you might use this even if the samples are actually standard PAM, to get approximate but fast results). -----------------------------------------------------------------------------*/ /* Here's how we think of the color mixing scaling operation: First, I'll describe scaling in one dimension. Assume we have a one row image. A raster row is ordinarily a sequence of discrete pixels which have no width and no distance between them -- only a sequence. Instead, think of the raster row as a bunch of pixels 1 unit wide adjacent to each other. For example, we are going to scale a 100 pixel row to a 150 pixel row. Imagine placing the input row right above the output row and stretching it so it is the same size as the output row. It still contains 100 pixels, but they are 1.5 units wide each. Our goal is to make the output row look as much as possible like the stretched input row, while observing that a pixel can be only one color. Output Pixel 0 is completely covered by Input Pixel 0, so we make Output Pixel 0 the same color as Input Pixel 0. Output Pixel 1 is covered half by Input Pixel 0 and half by Input Pixel 1. So we make Output Pixel 1 a 50/50 mix of Input Pixels 0 and 1. If you stand back far enough, input and output will look the same. This works for all scale factors, both scaling up and scaling down. For images with an opacity plane, imagine Input Pixel 0's foreground is fully opaque red (1,0,0,1), and Input Pixel 1 is fully transparent (foreground irrelevant) (0,0,0,0). We make Output Pixel 0's foreground fully opaque red as before. Output Pixel 1 is covered half by Input Pixel 0 and half by Input Pixel 1, so it is 50% opaque; but its foreground color is still red: (1,0,0,0.5). The output foreground color is the opacity and coverage weighted average of the input foreground colors, and the output opacity is the coverage weighted average of the input opacities. This program always stretches or squeezes the input row to be the same length as the output row; The output row's pixels are always 1 unit wide. The same thing works in the vertical direction. We think of rows as stacked strips of 1 unit height. We conceptually stretch the image vertically first (same process as above, but in place of a single-color pixels, we have a vector of colors). Then we take each row this vertical stretching generates and stretch it horizontally. */ tuplen * tuplenrow; /* An input row */ tuplen * newtuplenrow; /* Working space */ float rowsleft; /* The number of rows of output that need to be formed from the current input row (the one in tuplerow[]), less the number that have already been formed (either in accumulator[] or output to the file). This can be fractional because of the way we define rows as having height. */ float fracrowtofill; /* The fraction of the current output row (the one in vertScaledRow[]) that hasn't yet been filled in from an input row. */ tuplen * rowAccumulator; /* The red, green, and blue color intensities so far accumulated from input rows for the current output row. The ultimate value of this is an output row after vertical scaling, but before horizontal scaling. */ int rowsread; /* Number of rows of the input file that have been read */ int row; const pnm_transformMap * inputTransform; const pnm_transformMap * outputTransform; tuplenrow = pnm_allocpamrown(inpamP); rowAccumulator = pnm_allocpamrown(inpamP); rowsread = 0; rowsleft = 0.0; fracrowtofill = 1.0; newtuplenrow = pnm_allocpamrown(outpamP); createTransforms(inpamP, outpamP, assumeLinear, &inputTransform, &outputTransform); for (row = 0; row < outpamP->height; ++row) { /* First scale Y from tuplerow[] into rowAccumulator[]. */ zeroAccum(inpamP, rowAccumulator); if (outpamP->height == inpamP->height) { /* shortcut Y scaling */ readARow(inpamP, rowAccumulator, inputTransform); } else { while (fracrowtofill > 0) { if (rowsleft <= 0.0) { if (rowsread < inpamP->height) { readARow(inpamP, tuplenrow, inputTransform); ++rowsread; } else issueStretchWarning(verbose, fracrowtofill); rowsleft = yscale; } if (rowsleft < fracrowtofill) { accumOutputRow(inpamP, tuplenrow, rowsleft, rowAccumulator); fracrowtofill -= rowsleft; rowsleft = 0.0; } else { accumOutputRow(inpamP, tuplenrow, fracrowtofill, rowAccumulator); rowsleft = rowsleft - fracrowtofill; fracrowtofill = 0.0; } } fracrowtofill = 1.0; } /* 'rowAccumulator' now contains the contents of a single output row, but not yet horizontally scaled. Scale it now horizontally and write it out. */ scaleHorizontallyAndOutputRow(inpamP, outpamP, rowAccumulator, outputTransform, newtuplenrow, xscale, row, verbose); /* Destroys rowAccumulator */ } destroyTransforms(inputTransform, outputTransform); pnm_freepamrown(rowAccumulator); pnm_freepamrown(newtuplenrow); pnm_freepamrown(tuplenrow); } static void scaleWithoutMixing(const struct pam * const inpamP, const struct pam * const outpamP, float const xscale, float const yscale) { /*---------------------------------------------------------------------------- Scale the image described by *inpamP by xscale horizontally and yscale vertically and write the result as the image described by *outpamP. The input file is positioned past the header, to the beginning of the raster. The output file is too. Don't mix colors from different input pixels together in the output pixels. Each output pixel is an exact copy of some corresponding input pixel. -----------------------------------------------------------------------------*/ tuple * tuplerow; /* An input row */ tuple * newtuplerow; unsigned int row; /* The number of the next row to be output */ unsigned int rowInInput; /* The number of the next row to be read from the input */ assert(outpamP->maxval == inpamP->maxval); assert(outpamP->depth == inpamP->depth); tuplerow = pnm_allocpamrow(inpamP); rowInInput = 0; newtuplerow = pnm_allocpamrow(outpamP); for (row = 0; row < outpamP->height; ++row) { unsigned int col; unsigned int const inputRow = (int) (row / yscale); /* The number of the input row that we will use for this output row. */ for (; rowInInput <= inputRow; ++rowInInput) pnm_readpamrow(inpamP, tuplerow); for (col = 0; col < outpamP->width; ++col) { unsigned int const inputCol = (int) (col / xscale); pnm_assigntuple(inpamP, newtuplerow[col], tuplerow[inputCol]); } pnm_writepamrow(outpamP, newtuplerow); } /* Read off and discard rest of rows, because whatever is supplying the input stream may expect it to be consumed and because Caller will expect the stream to be positioned to the next image. */ for (; rowInInput < inpamP->height; ++rowInInput) pnm_readpamrow(inpamP, tuplerow); pnm_freepamrow(tuplerow); pnm_freepamrow(newtuplerow); } static void skipImage(struct pam * const pamP) { tuple * tuplerow; unsigned int row; tuplerow = pnm_allocpamrow(pamP); for (row = 0; row < pamP->height; ++row) pnm_readpamrow(pamP, tuplerow); pnm_freepamrow(tuplerow); } static void scale(FILE * const ifP, struct pam * const inpamP, struct pam * const outpamP, float const xscale, float const yscale, struct CmdlineInfo const cmdline) { pnm_writepaminit(outpamP); if (cmdline.nomix) { if (cmdline.verbose) pm_message("Using nomix method"); scaleWithoutMixing(inpamP, outpamP, xscale, yscale); } else if (!cmdline.filterFunction) { if (cmdline.verbose) pm_message("Using simple pixel mixing rescaling method"); scaleWithMixing(inpamP, outpamP, xscale, yscale, cmdline.linear, cmdline.verbose); } else { if (cmdline.verbose) pm_message("Using general filter method"); resample(inpamP, outpamP, cmdline.filterFunction, cmdline.filterRadius, cmdline.windowFunction, cmdline.verbose, cmdline.linear); } } static void pamscale(FILE * const ifP, FILE * const ofP, struct CmdlineInfo const cmdline) { struct pam inpam, outpam; float xscale, yscale; pnm_readpaminit(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type)); outpam = inpam; /* initial value */ outpam.file = ofP; if (PNM_FORMAT_TYPE(inpam.format) == PBM_TYPE && !cmdline.nomix) { outpam.format = PGM_TYPE; outpam.maxval = PGM_MAXMAXVAL; pm_message("promoting from PBM to PGM"); } else { outpam.format = inpam.format; outpam.maxval = inpam.maxval; } computeOutputDimensions(cmdline, inpam.width, inpam.height, &outpam.width, &outpam.height); xscale = (float) outpam.width / inpam.width; yscale = (float) outpam.height / inpam.height; if (cmdline.verbose) { pm_message("Scaling by %f horizontally to %d columns.", xscale, outpam.width); pm_message("Scaling by %f vertically to %d rows.", yscale, outpam.height); } if (xscale * inpam.width < outpam.width - 1 || yscale * inpam.height < outpam.height - 1) pm_error("Arithmetic precision of this program is inadequate to " "do the specified scaling. Use a smaller input image " "or a slightly different scale factor."); if (cmdline.reportonly) { printf("%d %d %f %f %d %d\n", inpam.width, inpam.height, xscale, yscale, outpam.width, outpam.height); skipImage(&inpam); } else scale(ifP, &inpam, &outpam, xscale, yscale, cmdline); } int main(int argc, const char **argv ) { struct CmdlineInfo cmdline; FILE * ifP; int eof; pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); ifP = pm_openr(cmdline.inputFileName); eof = false; while (!eof) { pamscale(ifP, stdout, cmdline); pnm_nextimage(ifP, &eof); } pm_close(ifP); pm_close(stdout); return 0; }