/* ** Version 1.0 September 28, 1996 ** ** Copyright (C) 1996 by Mike Burns ** ** Adapted to Netpbm 2005.08.10 by Bryan Henderson ** ** 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. */ /* References ** ---------- ** The select k'th value implementation is based on Algorithm 489 by ** Robert W. Floyd from the "Collected Algorithms from ACM" Volume II. ** ** The histogram sort is based is described in the paper "A Fast Two- ** Dimensional Median Filtering Algorithm" in "IEEE Transactions on ** Acoustics, Speech, and Signal Processing" Vol. ASSP-27, No. 1, February ** 1979. The algorithm I more closely followed is found in "Digital ** Image Processing Algorithms" by Ioannis Pitas. */ #include #include "pm_c_util.h" #include "pgm.h" #include "shhopt.h" #include "mallocvar.h" #include "nstring.h" enum MedianMethod {MEDIAN_UNSPECIFIED, SELECT_MEDIAN, HISTOGRAM_SORT_MEDIAN}; #define MAX_MEDIAN_TYPES 2 struct CmdlineInfo { /* All the information the user supplied in the command line, in a form easy for the program to use. */ const char * inputFileName; unsigned int width; unsigned int height; unsigned int cutoff; enum MedianMethod type; }; static int const forceplain = 0; /* Global variables common to each median sort routine. */ static gray ** grays; /* The convolution buffer. This is a circular buffer that contains the rows of the input image that are being convolved into the current output row. */ static gray * grayrow; /* A buffer for building the current output row */ static void parseCommandLine(int argc, const char ** 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 pm_optParseOptions3 on how to parse our options. */ optStruct3 opt; unsigned int option_def_index; unsigned int widthSpec, heightSpec, cutoffSpec, typeSpec; const char * type; MALLOCARRAY_NOFAIL(option_def, 100); option_def_index = 0; /* incremented by OPTENT3 */ OPTENT3(0, "width", OPT_UINT, &cmdlineP->width, &widthSpec, 0); OPTENT3(0, "height", OPT_UINT, &cmdlineP->height, &heightSpec, 0); OPTENT3(0, "cutoff", OPT_UINT, &cmdlineP->cutoff, &cutoffSpec, 0); OPTENT3(0, "type", OPT_STRING, &type, &typeSpec, 0); opt.opt_table = option_def; opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ opt.allowNegNum = FALSE; /* We may have 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 (widthSpec) { if (cmdlineP->width < 1) pm_error("-width must be at least 1"); } else cmdlineP->width = 3; if (heightSpec) { if (cmdlineP->height < 1) pm_error("-height must be at least 1"); } else cmdlineP->height = 3; if (!cutoffSpec) cmdlineP->cutoff = 250; if (typeSpec) { if (streq(type, "histogram_sort")) cmdlineP->type = HISTOGRAM_SORT_MEDIAN; else if (streq(type, "select")) cmdlineP->type = SELECT_MEDIAN; else pm_error("Invalid value '%s' for -type. Valid values are " "'histogram_sort' and 'select'", type); } else cmdlineP->type = MEDIAN_UNSPECIFIED; if (argc-1 < 1) cmdlineP->inputFileName = "-"; else { cmdlineP->inputFileName = argv[1]; if (argc-1 > 1) pm_error ("Too many arguments. The only argument is " "the optional input file name"); } } static void setWindow(gray ** const convBuffer, unsigned int const crows, gray ** const cgrayrow, unsigned int const lastRow) { /*---------------------------------------------------------------------------- Set 'cgrayrow' so it points into the circular buffer 'convBuffer' such that cgrayrow[0] is the topmost row in the buffer, given that the bottommost row in the buffer is row number 'lastRow'. -----------------------------------------------------------------------------*/ unsigned int const windowTopRow = (lastRow + 1) % crows; unsigned int bufferRow; unsigned int wrow; wrow = 0; for (bufferRow = windowTopRow; bufferRow < crows; ++wrow, ++bufferRow) cgrayrow[wrow] = grays[bufferRow]; for (bufferRow = 0; bufferRow < windowTopRow; ++wrow, ++bufferRow) cgrayrow[wrow] = grays[bufferRow]; } static void select489(gray * const a, int * const parray, int const n, int const k) { gray t; int i, j, l, r; int ptmp; l = 0; r = n - 1; while ( r > l ) { t = a[parray[k]]; i = l; j = r; ptmp = parray[l]; parray[l] = parray[k]; parray[k] = ptmp; if ( a[parray[r]] > t ) { ptmp = parray[r]; parray[r] = parray[l]; parray[l] = ptmp; } while ( i < j ) { ptmp = parray[i]; parray[i] = parray[j]; parray[j] = ptmp; ++i; --j; while ( a[parray[i]] < t ) ++i; while ( a[parray[j]] > t ) --j; } if ( a[parray[l]] == t ) { ptmp = parray[l]; parray[l] = parray[j]; parray[j] = ptmp; } else { ++j; ptmp = parray[j]; parray[j] = parray[r]; parray[r] = ptmp; } if ( j <= k ) l = j + 1; if ( k <= j ) r = j - 1; } } static void selectMedian(FILE * const ifP, unsigned int const ccols, unsigned int const crows, unsigned int const cols, unsigned int const rows, int const format, gray const maxval, unsigned int const median, unsigned int const firstRow) { unsigned int const ccolso2 = ccols / 2; unsigned int const crowso2 = crows / 2; unsigned int const numValues = crows * ccols; unsigned int col; gray * garray; /* Array of the currently gray values */ int * parray; int * subcol; gray ** cgrayrow; unsigned int row; garray = pgm_allocrow(numValues); MALLOCARRAY(cgrayrow, crows); MALLOCARRAY(parray, numValues); MALLOCARRAY(subcol, cols); if (cgrayrow == NULL || parray == NULL || subcol == NULL) pm_error("Unable to allocate memory"); for (col = 0; col < cols; ++col) subcol[col] = (col - (ccolso2 + 1)) % ccols; /* Apply median to main part of image. */ for (row = firstRow; row < rows; ++row) { unsigned int crow; unsigned int col; pgm_readpgmrow(ifP, grays[row % crows], cols, maxval, format); setWindow(grays, crows, cgrayrow, row); for (col = 0; col < cols; ++col) { if (col < ccolso2 || col >= cols - ccolso2) { grayrow[col] = cgrayrow[crowso2][col]; } else if (col == ccolso2) { unsigned int const leftcol = col - ccolso2; unsigned int i; i = 0; for (crow = 0; crow < crows; ++crow) { gray * const temprptr = cgrayrow[crow] + leftcol; unsigned int ccol; for (ccol = 0; ccol < ccols; ++ccol) { garray[i] = *(temprptr + ccol); parray[i] = i; ++i; } } select489(garray, parray, numValues, median); grayrow[col] = garray[parray[median]]; } else { unsigned int const addcol = col + ccolso2; unsigned int crow; unsigned int tsum; for (crow = 0, tsum = 0; crow < crows; ++crow, tsum += ccols) garray[tsum + subcol[col]] = *(cgrayrow[crow] + addcol ); select489( garray, parray, numValues, median ); grayrow[col] = garray[parray[median]]; } } pgm_writepgmrow( stdout, grayrow, cols, maxval, forceplain ); } free(subcol); free(parray); free(cgrayrow); pgm_freerow(garray); } static void histogramSortMedian(FILE * const ifP, unsigned int const ccols, unsigned int const crows, unsigned int const cols, unsigned int const rows, int const format, gray const maxval, unsigned int const median, unsigned int const firstRow) { unsigned int const ccolso2 = ccols / 2; unsigned int const crowso2 = crows / 2; unsigned int const histmax = maxval + 1; unsigned int * hist; unsigned int mdn, ltmdn; gray * leftCol; gray * rghtCol; gray ** cgrayrow; /* The window of the image currently being convolved, with cgrayrow[0] being the top row of the window. Pointers into grays[] */ unsigned int row; /* Row number in input -- bottommost row in the window we're currently convolving */ MALLOCARRAY(cgrayrow, crows); MALLOCARRAY(hist, histmax); if (cgrayrow == NULL || hist == NULL) pm_error("Unable to allocate memory"); leftCol = pgm_allocrow(crows); rghtCol = pgm_allocrow(crows); /* Apply median to main part of image. */ for (row = firstRow; row < rows; ++row) { unsigned int col; unsigned int i; /* initialize hist[] */ for (i = 0; i < histmax; ++i) hist[i] = 0; pgm_readpgmrow(ifP, grays[row % crows], cols, maxval, format); setWindow(grays, crows, cgrayrow, row); for (col = 0; col < cols; ++col) { if (col < ccolso2 || col >= cols - ccolso2) grayrow[col] = cgrayrow[crowso2][col]; else if (col == ccolso2) { unsigned int const leftcol = col - ccolso2; unsigned int crow; i = 0; for (crow = 0; crow < crows; ++crow) { unsigned int ccol; gray * const temprptr = cgrayrow[crow] + leftcol; for (ccol = 0; ccol < ccols; ++ccol) { gray const g = *(temprptr + ccol); ++hist[g]; ++i; } } ltmdn = 0; for (mdn = 0; ltmdn <= median; ++mdn) ltmdn += hist[mdn]; --mdn; if (ltmdn > median) ltmdn -= hist[mdn]; grayrow[col] = mdn; } else { unsigned int const subcol = col - (ccolso2 + 1); unsigned int const addcol = col + ccolso2; unsigned int crow; for (crow = 0; crow < crows; ++crow) { leftCol[crow] = *(cgrayrow[crow] + subcol); rghtCol[crow] = *(cgrayrow[crow] + addcol); } for (crow = 0; crow < crows; ++crow) { { gray const g = leftCol[crow]; --hist[(unsigned int) g]; if ((unsigned int) g < mdn) --ltmdn; } { gray const g = rghtCol[crow]; ++hist[(unsigned int) g]; if ((unsigned int) g < mdn) ++ltmdn; } } if (ltmdn > median) do { --mdn; ltmdn -= hist[mdn]; } while (ltmdn > median); else { /* This one change from Pitas algorithm can reduce run ** time by up to 10%. */ while (ltmdn <= median) { ltmdn += hist[mdn]; ++mdn; } --mdn; if (ltmdn > median) ltmdn -= hist[mdn]; } grayrow[col] = mdn; } } pgm_writepgmrow(stdout, grayrow, cols, maxval, forceplain); } pgm_freerow(leftCol); pgm_freerow(rghtCol); free(hist); free(cgrayrow); } static void convolve(FILE * const ifP, unsigned int const cols, unsigned int const rows, gray const maxval, int const format, unsigned int const ccols, unsigned int const crows, enum MedianMethod const medianMethod, unsigned int const median) { unsigned int const crowso2 = crows / 2; unsigned int row; /* An even-size convolution window is biased toward the top and left. So if it is 8 rows, the window covers 4 rows above the target row and 3 rows below it, plus the target row itself. 'crowso2' is the number of the target row within the window. There are always 'crowso2' rows above it and either crowso2 or crowso2-1 rows below it. */ /* Allocate space for number of rows in mask size. */ grays = pgm_allocarray(cols, crows); grayrow = pgm_allocrow(cols); /* Prime the convolution window -- fill it except the last row */ for (row = 0; row < crows - 1; ++row) pgm_readpgmrow(ifP, grays[row], cols, maxval, format); /* Copy the top half out verbatim, since convolution kernel for these rows runs off the top of the image. */ for (row = 0; row < crowso2; ++row) pgm_writepgmrow(stdout, grays[row], cols, maxval, forceplain); switch (medianMethod) { case SELECT_MEDIAN: selectMedian(ifP, ccols, crows, cols, rows, format, maxval, median, crows-1); break; case HISTOGRAM_SORT_MEDIAN: histogramSortMedian(ifP, ccols, crows, cols, rows, format, maxval, median, crows-1); break; case MEDIAN_UNSPECIFIED: pm_error("INTERNAL ERROR: median unspecified"); } /* Copy the bottom half of the remaining convolution window verbatim, since convolution kernel for these rows runs off the bottom of the image. */ assert(crows >= crowso2 + 1); for (row = rows - (crows-crowso2-1); row < rows; ++row) pgm_writepgmrow(stdout, grays[row % crows], cols, maxval, forceplain); pgm_freearray(grays, crows); pgm_freerow(grayrow); } int main(int argc, const char * argv[]) { struct CmdlineInfo cmdline; FILE * ifP; int cols, rows; int format; gray maxval; unsigned int ccols, crows; unsigned int median; enum MedianMethod medianMethod; pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); ifP = pm_openr(cmdline.inputFileName); assert(cmdline.height > 0 && cmdline.width > 0); pgm_readpgminit(ifP, &cols, &rows, &maxval, &format); ccols = MIN(cmdline.width, cols+1); crows = MIN(cmdline.height, rows+1); pgm_writepgminit(stdout, cols, rows, maxval, forceplain); median = (crows * ccols) / 2; /* Choose which sort to run. */ if (cmdline.type == MEDIAN_UNSPECIFIED) { if ((maxval / ((ccols * crows) - 1)) < cmdline.cutoff) medianMethod = HISTOGRAM_SORT_MEDIAN; else medianMethod = SELECT_MEDIAN; } else medianMethod = cmdline.type; convolve(ifP, cols, rows, maxval, format, ccols, crows, medianMethod, median); pm_close(ifP); pm_close(stdout); return 0; }