diff options
Diffstat (limited to 'editor/pnmconvol.c')
| -rw-r--r-- | editor/pnmconvol.c | 1989 |
1 files changed, 1989 insertions, 0 deletions
diff --git a/editor/pnmconvol.c b/editor/pnmconvol.c new file mode 100644 index 00000000..0bf44ce3 --- /dev/null +++ b/editor/pnmconvol.c @@ -0,0 +1,1989 @@ +/* pnmconvol.c - general MxN convolution on a PNM image +** +** Version 2.0.1 January 30, 1995 +** +** Major rewriting by Mike Burns +** Copyright (C) 1994, 1995 by Mike Burns (burns@chem.psu.edu) +** +** 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. +*/ + +/* A change history is at the bottom */ + +#include "pnm.h" +#include "shhopt.h" +#include "mallocvar.h" + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* '-' if stdin */ + const char *kernelFilespec; + unsigned int nooffset; +}; + +static void +parseCommandLine(int argc, char ** argv, + struct cmdlineInfo *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; + /* Instructions to optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + MALLOCARRAY_NOFAIL(option_def, 100); + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "nooffset", OPT_FLAG, NULL, + &cmdlineP->nooffset, 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 */ + + optParseOptions3( &argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (argc-1 < 1) + pm_error("Need at least one argument: file specification of the " + "convolution kernel image."); + + cmdlineP->kernelFilespec = argv[1]; + + if (argc-1 >= 2) + cmdlineP->inputFilespec = argv[2]; + else + cmdlineP->inputFilespec = "-"; + + if (argc-1 > 2) + pm_error("Too many arguments. Only acceptable arguments are: " + "convolution file name and input file name"); +} + + +/* Macros to verify that r,g,b values are within proper range */ + +#define CHECK_GRAY \ + if ( tempgsum < 0L ) g = 0; \ + else if ( tempgsum > maxval ) g = maxval; \ + else g = tempgsum; + +#define CHECK_RED \ + if ( temprsum < 0L ) r = 0; \ + else if ( temprsum > maxval ) r = maxval; \ + else r = temprsum; + +#define CHECK_GREEN \ + if ( tempgsum < 0L ) g = 0; \ + else if ( tempgsum > maxval ) g = maxval; \ + else g = tempgsum; + +#define CHECK_BLUE \ + if ( tempbsum < 0L ) b = 0; \ + else if ( tempbsum > maxval ) b = maxval; \ + else b = tempbsum; + +struct convolveType { + void (*ppmConvolver)(const float ** const rweights, + const float ** const gweights, + const float ** const bweights); + void (*pgmConvolver)(const float ** const weights); +}; + +static FILE* ifp; +static int crows, ccols, ccolso2, crowso2; +static int cols, rows; +static xelval maxval; +static int format, newformat; + + + +static void +computeWeights(xel * const * const cxels, + int const ccols, + int const crows, + int const cformat, + xelval const cmaxval, + bool const offsetPgm, + float *** const rweightsP, + float *** const gweightsP, + float *** const bweightsP) { +/*---------------------------------------------------------------------------- + Compute the convolution matrix in normalized form from the PGM + form. Each element of the output matrix is the actual weight we give an + input pixel -- i.e. the thing by which we multiple a value from the + input image. + + 'offsetPgm' means the PGM convolution matrix is defined in offset form so + that it can represent negative values. E.g. with maxval 100, 50 means + 0, 100 means 50, and 0 means -50. If 'offsetPgm' is false, 0 means 0 + and there are no negative weights. +-----------------------------------------------------------------------------*/ + double const scale = (offsetPgm ? 2.0 : 1.0) / cmaxval; + double const offset = offsetPgm ? - 1.0 : 0.0; + + float** rweights; + float** gweights; + float** bweights; + + float rsum, gsum, bsum; + + unsigned int crow; + + /* Set up the normalized weights. */ + rweights = (float**) pm_allocarray(ccols, crows, sizeof(float)); + gweights = (float**) pm_allocarray(ccols, crows, sizeof(float)); + bweights = (float**) pm_allocarray(ccols, crows, sizeof(float)); + + rsum = gsum = bsum = 0.0; /* initial value */ + + for (crow = 0; crow < crows; ++crow) { + unsigned int ccol; + for (ccol = 0; ccol < ccols; ++ccol) { + switch (PNM_FORMAT_TYPE(cformat)) { + case PPM_TYPE: + rsum += rweights[crow][ccol] = + (PPM_GETR(cxels[crow][ccol]) * scale + offset); + gsum += gweights[crow][ccol] = + (PPM_GETG(cxels[crow][ccol]) * scale + offset); + bsum += bweights[crow][ccol] = + (PPM_GETB(cxels[crow][ccol]) * scale + offset); + break; + + default: + gsum += gweights[crow][ccol] = + (PNM_GET1(cxels[crow][ccol]) * scale + offset); + break; + } + } + } + *rweightsP = rweights; + *gweightsP = gweights; + *bweightsP = bweights; + + switch (PNM_FORMAT_TYPE(format)) { + case PPM_TYPE: + if (rsum < 0.9 || rsum > 1.1 || gsum < 0.9 || gsum > 1.1 || + bsum < 0.9 || bsum > 1.1) { + pm_message("WARNING - this convolution matrix is biased. " + "red, green, and blue average weights: %f, %f, %f " + "(unbiased would be 1).", + rsum, gsum, bsum); + + if (rsum < 0 && gsum < 0 && bsum < 0) + pm_message("Maybe you want the -nooffset option?"); + } + break; + + default: + if (gsum < 0.9 || gsum > 1.1) + pm_message("WARNING - this convolution matrix is biased. " + "average weight = %f (unbiased would be 1)", + gsum); + break; + } +} + + + +/* General PGM Convolution +** +** No useful redundancy in convolution matrix. +*/ + +static void +pgm_general_convolve(const float ** const weights) { + xel** xelbuf; + xel* outputrow; + xelval g; + int row; + xel **rowptr, *temprptr; + int toprow, temprow; + int i, irow; + long tempgsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + a row output buffer. + */ + xelbuf = pnm_allocarray(cols, crows); + outputrow = pnm_allocrow(cols); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray(1, crows); + + pnm_writepnminit(stdout, cols, rows, maxval, newformat, 0); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for (row = 0; row < crows - 1; ++row) { + pnm_readpnmrow(ifp, xelbuf[row], cols, maxval, format); + if (PNM_FORMAT_TYPE(format) != newformat) + pnm_promoteformatrow(xelbuf[row], cols, maxval, format, + maxval, newformat); + /* Write out just the part we're not going to convolve. */ + if (row < crowso2) + pnm_writepnmrow(stdout, xelbuf[row], cols, maxval, newformat, 0); + } + + /* Now the rest of the image - read in the row at the end of + xelbuf, and convolve and write out the row in the middle. + */ + for (; row < rows; ++row) { + int col; + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow(ifp, xelbuf[temprow], cols, maxval, format); + if (PNM_FORMAT_TYPE(format) != newformat) + pnm_promoteformatrow(xelbuf[temprow], cols, maxval, format, + maxval, newformat); + + /* Arrange rowptr to eliminate the use of mod function to determine + which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for (col = 0; col < cols; ++col) { + if (col < ccolso2 || col >= cols - ccolso2) + outputrow[col] = rowptr[crowso2][col]; + else { + int const leftcol = col - ccolso2; + int crow; + float gsum; + gsum = 0.0; + for (crow = 0; crow < crows; ++crow) { + int ccol; + temprptr = rowptr[crow] + leftcol; + for (ccol = 0; ccol < ccols; ++ccol) + gsum += PNM_GET1(*(temprptr + ccol)) + * weights[crow][ccol]; + } + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow(stdout, outputrow, cols, maxval, newformat, 0); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for (irow = crowso2 + 1; irow < crows; ++irow) + pnm_writepnmrow(stdout, rowptr[irow], cols, maxval, newformat, 0 ); +} + + + +/* PGM Mean Convolution +** +** This is the common case where you just want the target pixel replaced with +** the average value of its neighbors. This can work much faster than the +** general case because you can reduce the number of floating point operations +** that are required since all the weights are the same. You will only need +** to multiply by the weight once, not for every pixel in the convolution +** matrix. +** +** This algorithm works by creating sums for each column of crows height for +** the whole width of the image. Then add ccols column sums together to obtain +** the total sum of the neighbors and multiply that sum by the weight. As you +** move right to left to calculate the next pixel, take the total sum you just +** generated, add in the value of the next column and subtract the value of the +** leftmost column. Multiply that by the weight and that's it. As you move +** down a row, calculate new column sums by using previous sum for that column +** and adding in pixel on current row and subtracting pixel in top row. +** +*/ + + +static void +pgm_mean_convolve(const float ** const weights) { + float const gmeanweight = weights[0][0]; + + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval g; + int row, crow; + xel **rowptr, *temprptr; + int leftcol; + int i, irow; + int toprow, temprow; + int subrow, addrow; + int subcol, addcol; + long gisum; + int tempcol, crowsp1; + long tempgsum; + long *gcolumnsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. MEAN uses an extra row. */ + xelbuf = pnm_allocarray( cols, crows + 1 ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf. MEAN uses an extra row. */ + rowptr = (xel **) pnm_allocarray( 1, crows + 1); + + /* Allocate space for intermediate column sums */ + gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) ); + for ( col = 0; col < cols; ++col ) + gcolumnsum[col] = 0L; + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* Do first real row only */ + subrow = crows; + addrow = crows - 1; + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + gisum = 0L; + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + for ( ccol = 0; ccol < ccols; ++ccol ) + gcolumnsum[leftcol + ccol] += + PNM_GET1( *(temprptr + ccol) ); + } + for ( ccol = 0; ccol < ccols; ++ccol) + gisum += gcolumnsum[leftcol + ccol]; + tempgsum = (float) gisum * gmeanweight + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + /* Column numbers to subtract or add to isum */ + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + for ( crow = 0; crow < crows; ++crow ) + gcolumnsum[addcol] += PNM_GET1( rowptr[crow][addcol] ); + gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol]; + tempgsum = (float) gisum * gmeanweight + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + + ++row; + /* For all subsequent rows do it this way as the columnsums have been + ** generated. Now we can use them to reduce further calculations. + */ + crowsp1 = crows + 1; + for ( ; row < rows; ++row ) + { + toprow = row + 1; + temprow = row % (crows + 1); + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + /* This rearrangement using crows+1 rowptrs and xelbufs will cause + ** rowptr[0..crows-1] to always hold active xelbufs and for + ** rowptr[crows] to always hold the oldest (top most) xelbuf. + */ + temprow = (toprow + 1) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + gisum = 0L; + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + tempcol = leftcol + ccol; + gcolumnsum[tempcol] = gcolumnsum[tempcol] + - PNM_GET1( rowptr[subrow][ccol] ) + + PNM_GET1( rowptr[addrow][ccol] ); + gisum += gcolumnsum[tempcol]; + } + tempgsum = (float) gisum * gmeanweight + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + /* Column numbers to subtract or add to isum */ + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + gcolumnsum[addcol] = gcolumnsum[addcol] + - PNM_GET1( rowptr[subrow][addcol] ) + + PNM_GET1( rowptr[addrow][addcol] ); + gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol]; + tempgsum = (float) gisum * gmeanweight + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + +/* PGM Horizontal Convolution +** +** Similar idea to using columnsums of the Mean and Vertical convolution, +** but uses temporary sums of row values. Need to multiply by weights crows +** number of times. Each time a new line is started, must recalculate the +** initials rowsums for the newest row only. Uses queue to still access +** previous row sums. +** +*/ + +static void +pgm_horizontal_convolve(const float ** const weights) { + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval g; + int row, crow; + xel **rowptr, *temprptr; + int leftcol; + int i, irow; + int temprow; + int subcol, addcol; + float gsum; + int addrow, subrow; + long **growsum, **growsumptr; + int crowsp1; + long tempgsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. */ + xelbuf = pnm_allocarray( cols, crows + 1 ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray( 1, crows + 1); + + /* Allocate intermediate row sums. HORIZONTAL uses an extra row. */ + /* crows current rows and 1 extra for newest added row. */ + growsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) ); + growsumptr = (long **) pnm_allocarray( 1, crows + 1); + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* First row only */ + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = (row + 1) % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for ( crow = 0; crow < crows; ++crow ) + growsumptr[crow] = growsum[crow]; + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + gsum = 0.0; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + growsumptr[crow][leftcol] = 0L; + for ( ccol = 0; ccol < ccols; ++ccol ) + growsumptr[crow][leftcol] += + PNM_GET1( *(temprptr + ccol) ); + gsum += growsumptr[crow][leftcol] * weights[crow][0]; + } + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + gsum = 0.0; + leftcol = col - ccolso2; + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + growsumptr[crow][leftcol] = growsumptr[crow][subcol] + - PNM_GET1( rowptr[crow][subcol] ) + + PNM_GET1( rowptr[crow][addcol] ); + gsum += growsumptr[crow][leftcol] * weights[crow][0]; + } + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + + + /* For all subsequent rows */ + + subrow = crows; + addrow = crows - 1; + crowsp1 = crows + 1; + ++row; + for ( ; row < rows; ++row ) + { + temprow = row % crowsp1; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = (row + 2) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + { + rowptr[i] = xelbuf[irow]; + growsumptr[i] = growsum[irow]; + } + for (irow = 0; irow < temprow; ++irow, ++i) + { + rowptr[i] = xelbuf[irow]; + growsumptr[i] = growsum[irow]; + } + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + gsum = 0.0; + leftcol = col - ccolso2; + growsumptr[addrow][leftcol] = 0L; + for ( ccol = 0; ccol < ccols; ++ccol ) + growsumptr[addrow][leftcol] += + PNM_GET1( rowptr[addrow][leftcol + ccol] ); + for ( crow = 0; crow < crows; ++crow ) + gsum += growsumptr[crow][leftcol] * weights[crow][0]; + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + gsum = 0.0; + leftcol = col - ccolso2; + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + growsumptr[addrow][leftcol] = growsumptr[addrow][subcol] + - PNM_GET1( rowptr[addrow][subcol] ) + + PNM_GET1( rowptr[addrow][addcol] ); + for ( crow = 0; crow < crows; ++crow ) + gsum += growsumptr[crow][leftcol] * weights[crow][0]; + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + +/* PGM Vertical Convolution +** +** Uses column sums as in Mean Convolution. +** +*/ + + +static void +pgm_vertical_convolve(const float ** const weights) { + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval g; + int row, crow; + xel **rowptr, *temprptr; + int leftcol; + int i, irow; + int toprow, temprow; + int subrow, addrow; + int tempcol; + float gsum; + long *gcolumnsum; + int crowsp1; + int addcol; + long tempgsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. VERTICAL uses an extra row. */ + xelbuf = pnm_allocarray( cols, crows + 1 ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray( 1, crows + 1 ); + + /* Allocate space for intermediate column sums */ + gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) ); + for ( col = 0; col < cols; ++col ) + gcolumnsum[col] = 0L; + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* Now the rest of the image - read in the row at the end of + ** xelbuf, and convolve and write out the row in the middle. + */ + /* For first row only */ + + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + /* Arrange rowptr to eliminate the use of mod function to determine + ** which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + gsum = 0.0; + leftcol = col - ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + for ( ccol = 0; ccol < ccols; ++ccol ) + gcolumnsum[leftcol + ccol] += + PNM_GET1( *(temprptr + ccol) ); + } + for ( ccol = 0; ccol < ccols; ++ccol) + gsum += gcolumnsum[leftcol + ccol] * weights[0][ccol]; + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + gsum = 0.0; + leftcol = col - ccolso2; + addcol = col + ccolso2; + for ( crow = 0; crow < crows; ++crow ) + gcolumnsum[addcol] += PNM_GET1( rowptr[crow][addcol] ); + for ( ccol = 0; ccol < ccols; ++ccol ) + gsum += gcolumnsum[leftcol + ccol] * weights[0][ccol]; + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + + /* For all subsequent rows */ + subrow = crows; + addrow = crows - 1; + crowsp1 = crows + 1; + ++row; + for ( ; row < rows; ++row ) + { + toprow = row + 1; + temprow = row % (crows +1); + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + /* Arrange rowptr to eliminate the use of mod function to determine + ** which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = (toprow + 1) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + gsum = 0.0; + leftcol = col - ccolso2; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + tempcol = leftcol + ccol; + gcolumnsum[tempcol] = gcolumnsum[tempcol] + - PNM_GET1( rowptr[subrow][ccol] ) + + PNM_GET1( rowptr[addrow][ccol] ); + gsum = gsum + gcolumnsum[tempcol] * weights[0][ccol]; + } + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + else + { + gsum = 0.0; + leftcol = col - ccolso2; + addcol = col + ccolso2; + gcolumnsum[addcol] = gcolumnsum[addcol] + - PNM_GET1( rowptr[subrow][addcol] ) + + PNM_GET1( rowptr[addrow][addcol] ); + for ( ccol = 0; ccol < ccols; ++ccol ) + gsum += gcolumnsum[leftcol + ccol] * weights[0][ccol]; + tempgsum = gsum + 0.5; + CHECK_GRAY; + PNM_ASSIGN1( outputrow[col], g ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + + + +/* PPM General Convolution Algorithm +** +** No redundancy in convolution matrix. Just use brute force. +** See pgm_general_convolve() for more details. +*/ + +static void +ppm_general_convolve(const float ** const rweights, + const float ** const gweights, + const float ** const bweights) { + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval r, g, b; + int row, crow; + float rsum, gsum, bsum; + xel **rowptr, *temprptr; + int toprow, temprow; + int i, irow; + int leftcol; + long temprsum, tempgsum, tempbsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. */ + xelbuf = pnm_allocarray( cols, crows ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray( 1, crows ); + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* Now the rest of the image - read in the row at the end of + ** xelbuf, and convolve and write out the row in the middle. + */ + for ( ; row < rows; ++row ) + { + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + /* Arrange rowptr to eliminate the use of mod function to determine + ** which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else + { + leftcol = col - ccolso2; + rsum = gsum = bsum = 0.0; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + rsum += PPM_GETR( *(temprptr + ccol) ) + * rweights[crow][ccol]; + gsum += PPM_GETG( *(temprptr + ccol) ) + * gweights[crow][ccol]; + bsum += PPM_GETB( *(temprptr + ccol) ) + * bweights[crow][ccol]; + } + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + +/* PPM Mean Convolution +** +** Same as pgm_mean_convolve() but for PPM. +** +*/ + +static void +ppm_mean_convolve(const float ** const rweights, + const float ** const gweights, + const float ** const bweights) { + /* All weights of a single color are the same so just grab any one + of them. + */ + float const rmeanweight = rweights[0][0]; + float const gmeanweight = gweights[0][0]; + float const bmeanweight = bweights[0][0]; + + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval r, g, b; + int row, crow; + xel **rowptr, *temprptr; + int leftcol; + int i, irow; + int toprow, temprow; + int subrow, addrow; + int subcol, addcol; + long risum, gisum, bisum; + long temprsum, tempgsum, tempbsum; + int tempcol, crowsp1; + long *rcolumnsum, *gcolumnsum, *bcolumnsum; + + + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. MEAN uses an extra row. */ + xelbuf = pnm_allocarray( cols, crows + 1 ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf. MEAN uses an extra row. */ + rowptr = (xel **) pnm_allocarray( 1, crows + 1); + + /* Allocate space for intermediate column sums */ + rcolumnsum = (long *) pm_allocrow( cols, sizeof(long) ); + gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) ); + bcolumnsum = (long *) pm_allocrow( cols, sizeof(long) ); + for ( col = 0; col < cols; ++col ) + { + rcolumnsum[col] = 0L; + gcolumnsum[col] = 0L; + bcolumnsum[col] = 0L; + } + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* Do first real row only */ + subrow = crows; + addrow = crows - 1; + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + risum = 0L; + gisum = 0L; + bisum = 0L; + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + rcolumnsum[leftcol + ccol] += + PPM_GETR( *(temprptr + ccol) ); + gcolumnsum[leftcol + ccol] += + PPM_GETG( *(temprptr + ccol) ); + bcolumnsum[leftcol + ccol] += + PPM_GETB( *(temprptr + ccol) ); + } + } + for ( ccol = 0; ccol < ccols; ++ccol) + { + risum += rcolumnsum[leftcol + ccol]; + gisum += gcolumnsum[leftcol + ccol]; + bisum += bcolumnsum[leftcol + ccol]; + } + temprsum = (float) risum * rmeanweight + 0.5; + tempgsum = (float) gisum * gmeanweight + 0.5; + tempbsum = (float) bisum * bmeanweight + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + else + { + /* Column numbers to subtract or add to isum */ + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + rcolumnsum[addcol] += PPM_GETR( rowptr[crow][addcol] ); + gcolumnsum[addcol] += PPM_GETG( rowptr[crow][addcol] ); + bcolumnsum[addcol] += PPM_GETB( rowptr[crow][addcol] ); + } + risum = risum - rcolumnsum[subcol] + rcolumnsum[addcol]; + gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol]; + bisum = bisum - bcolumnsum[subcol] + bcolumnsum[addcol]; + temprsum = (float) risum * rmeanweight + 0.5; + tempgsum = (float) gisum * gmeanweight + 0.5; + tempbsum = (float) bisum * bmeanweight + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + + ++row; + /* For all subsequent rows do it this way as the columnsums have been + ** generated. Now we can use them to reduce further calculations. + */ + crowsp1 = crows + 1; + for ( ; row < rows; ++row ) + { + toprow = row + 1; + temprow = row % (crows + 1); + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + /* This rearrangement using crows+1 rowptrs and xelbufs will cause + ** rowptr[0..crows-1] to always hold active xelbufs and for + ** rowptr[crows] to always hold the oldest (top most) xelbuf. + */ + temprow = (toprow + 1) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + risum = 0L; + gisum = 0L; + bisum = 0L; + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + tempcol = leftcol + ccol; + rcolumnsum[tempcol] = rcolumnsum[tempcol] + - PPM_GETR( rowptr[subrow][ccol] ) + + PPM_GETR( rowptr[addrow][ccol] ); + risum += rcolumnsum[tempcol]; + gcolumnsum[tempcol] = gcolumnsum[tempcol] + - PPM_GETG( rowptr[subrow][ccol] ) + + PPM_GETG( rowptr[addrow][ccol] ); + gisum += gcolumnsum[tempcol]; + bcolumnsum[tempcol] = bcolumnsum[tempcol] + - PPM_GETB( rowptr[subrow][ccol] ) + + PPM_GETB( rowptr[addrow][ccol] ); + bisum += bcolumnsum[tempcol]; + } + temprsum = (float) risum * rmeanweight + 0.5; + tempgsum = (float) gisum * gmeanweight + 0.5; + tempbsum = (float) bisum * bmeanweight + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + else + { + /* Column numbers to subtract or add to isum */ + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + rcolumnsum[addcol] = rcolumnsum[addcol] + - PPM_GETR( rowptr[subrow][addcol] ) + + PPM_GETR( rowptr[addrow][addcol] ); + risum = risum - rcolumnsum[subcol] + rcolumnsum[addcol]; + gcolumnsum[addcol] = gcolumnsum[addcol] + - PPM_GETG( rowptr[subrow][addcol] ) + + PPM_GETG( rowptr[addrow][addcol] ); + gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol]; + bcolumnsum[addcol] = bcolumnsum[addcol] + - PPM_GETB( rowptr[subrow][addcol] ) + + PPM_GETB( rowptr[addrow][addcol] ); + bisum = bisum - bcolumnsum[subcol] + bcolumnsum[addcol]; + temprsum = (float) risum * rmeanweight + 0.5; + tempgsum = (float) gisum * gmeanweight + 0.5; + tempbsum = (float) bisum * bmeanweight + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + +/* PPM Horizontal Convolution +** +** Same as pgm_horizontal_convolve() +** +**/ + +static void +ppm_horizontal_convolve(const float ** const rweights, + const float ** const gweights, + const float ** const bweights) { + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval r, g, b; + int row, crow; + xel **rowptr, *temprptr; + int leftcol; + int i, irow; + int temprow; + int subcol, addcol; + float rsum, gsum, bsum; + int addrow, subrow; + long **rrowsum, **rrowsumptr; + long **growsum, **growsumptr; + long **browsum, **browsumptr; + int crowsp1; + long temprsum, tempgsum, tempbsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. */ + xelbuf = pnm_allocarray( cols, crows + 1 ); + outputrow = pnm_allocrow( cols ); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray( 1, crows + 1); + + /* Allocate intermediate row sums. HORIZONTAL uses an extra row */ + rrowsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) ); + rrowsumptr = (long **) pnm_allocarray( 1, crows + 1); + growsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) ); + growsumptr = (long **) pnm_allocarray( 1, crows + 1); + browsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) ); + browsumptr = (long **) pnm_allocarray( 1, crows + 1); + + pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 ); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for ( row = 0; row < crows - 1; ++row ) + { + pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[row], cols, maxval, format, maxval, newformat ); + /* Write out just the part we're not going to convolve. */ + if ( row < crowso2 ) + pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 ); + } + + /* First row only */ + temprow = row % crows; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = (row + 1) % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for ( crow = 0; crow < crows; ++crow ) + { + rrowsumptr[crow] = rrowsum[crow]; + growsumptr[crow] = growsum[crow]; + browsumptr[crow] = browsum[crow]; + } + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + leftcol = col - ccolso2; + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + for ( crow = 0; crow < crows; ++crow ) + { + temprptr = rowptr[crow] + leftcol; + rrowsumptr[crow][leftcol] = 0L; + growsumptr[crow][leftcol] = 0L; + browsumptr[crow][leftcol] = 0L; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + rrowsumptr[crow][leftcol] += + PPM_GETR( *(temprptr + ccol) ); + growsumptr[crow][leftcol] += + PPM_GETG( *(temprptr + ccol) ); + browsumptr[crow][leftcol] += + PPM_GETB( *(temprptr + ccol) ); + } + rsum += rrowsumptr[crow][leftcol] * rweights[crow][0]; + gsum += growsumptr[crow][leftcol] * gweights[crow][0]; + bsum += browsumptr[crow][leftcol] * bweights[crow][0]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + else + { + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + leftcol = col - ccolso2; + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + for ( crow = 0; crow < crows; ++crow ) + { + rrowsumptr[crow][leftcol] = rrowsumptr[crow][subcol] + - PPM_GETR( rowptr[crow][subcol] ) + + PPM_GETR( rowptr[crow][addcol] ); + rsum += rrowsumptr[crow][leftcol] * rweights[crow][0]; + growsumptr[crow][leftcol] = growsumptr[crow][subcol] + - PPM_GETG( rowptr[crow][subcol] ) + + PPM_GETG( rowptr[crow][addcol] ); + gsum += growsumptr[crow][leftcol] * gweights[crow][0]; + browsumptr[crow][leftcol] = browsumptr[crow][subcol] + - PPM_GETB( rowptr[crow][subcol] ) + + PPM_GETB( rowptr[crow][addcol] ); + bsum += browsumptr[crow][leftcol] * bweights[crow][0]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + + + /* For all subsequent rows */ + + subrow = crows; + addrow = crows - 1; + crowsp1 = crows + 1; + ++row; + for ( ; row < rows; ++row ) + { + temprow = row % crowsp1; + pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format ); + if ( PNM_FORMAT_TYPE(format) != newformat ) + pnm_promoteformatrow( + xelbuf[temprow], cols, maxval, format, maxval, newformat ); + + temprow = (row + 2) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + { + rowptr[i] = xelbuf[irow]; + rrowsumptr[i] = rrowsum[irow]; + growsumptr[i] = growsum[irow]; + browsumptr[i] = browsum[irow]; + } + for (irow = 0; irow < temprow; ++irow, ++i) + { + rowptr[i] = xelbuf[irow]; + rrowsumptr[i] = rrowsum[irow]; + growsumptr[i] = growsum[irow]; + browsumptr[i] = browsum[irow]; + } + + for ( col = 0; col < cols; ++col ) + { + if ( col < ccolso2 || col >= cols - ccolso2 ) + outputrow[col] = rowptr[crowso2][col]; + else if ( col == ccolso2 ) + { + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + leftcol = col - ccolso2; + rrowsumptr[addrow][leftcol] = 0L; + growsumptr[addrow][leftcol] = 0L; + browsumptr[addrow][leftcol] = 0L; + for ( ccol = 0; ccol < ccols; ++ccol ) + { + rrowsumptr[addrow][leftcol] += + PPM_GETR( rowptr[addrow][leftcol + ccol] ); + growsumptr[addrow][leftcol] += + PPM_GETG( rowptr[addrow][leftcol + ccol] ); + browsumptr[addrow][leftcol] += + PPM_GETB( rowptr[addrow][leftcol + ccol] ); + } + for ( crow = 0; crow < crows; ++crow ) + { + rsum += rrowsumptr[crow][leftcol] * rweights[crow][0]; + gsum += growsumptr[crow][leftcol] * gweights[crow][0]; + bsum += browsumptr[crow][leftcol] * bweights[crow][0]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + else + { + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + leftcol = col - ccolso2; + subcol = col - ccolso2 - 1; + addcol = col + ccolso2; + rrowsumptr[addrow][leftcol] = rrowsumptr[addrow][subcol] + - PPM_GETR( rowptr[addrow][subcol] ) + + PPM_GETR( rowptr[addrow][addcol] ); + growsumptr[addrow][leftcol] = growsumptr[addrow][subcol] + - PPM_GETG( rowptr[addrow][subcol] ) + + PPM_GETG( rowptr[addrow][addcol] ); + browsumptr[addrow][leftcol] = browsumptr[addrow][subcol] + - PPM_GETB( rowptr[addrow][subcol] ) + + PPM_GETB( rowptr[addrow][addcol] ); + for ( crow = 0; crow < crows; ++crow ) + { + rsum += rrowsumptr[crow][leftcol] * rweights[crow][0]; + gsum += growsumptr[crow][leftcol] * gweights[crow][0]; + bsum += browsumptr[crow][leftcol] * bweights[crow][0]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN( outputrow[col], r, g, b ); + } + } + pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 ); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for ( irow = crowso2 + 1; irow < crows; ++irow ) + pnm_writepnmrow( + stdout, rowptr[irow], cols, maxval, newformat, 0 ); + + } + + +/* PPM Vertical Convolution +** +** Same as pgm_vertical_convolve() +** +*/ + +static void +ppm_vertical_convolve(const float ** const rweights, + const float ** const gweights, + const float ** const bweights) { + int ccol, col; + xel** xelbuf; + xel* outputrow; + xelval r, g, b; + int row, crow; + xel **rowptr, *temprptr; + int i, irow; + int toprow, temprow; + int subrow, addrow; + int tempcol; + long *rcolumnsum, *gcolumnsum, *bcolumnsum; + int crowsp1; + int addcol; + long temprsum, tempgsum, tempbsum; + + /* Allocate space for one convolution-matrix's worth of rows, plus + ** a row output buffer. VERTICAL uses an extra row. */ + xelbuf = pnm_allocarray(cols, crows + 1); + outputrow = pnm_allocrow(cols); + + /* Allocate array of pointers to xelbuf */ + rowptr = (xel **) pnm_allocarray(1, crows + 1); + + /* Allocate space for intermediate column sums */ + MALLOCARRAY_NOFAIL(rcolumnsum, cols); + MALLOCARRAY_NOFAIL(gcolumnsum, cols); + MALLOCARRAY_NOFAIL(bcolumnsum, cols); + + for (col = 0; col < cols; ++col) { + rcolumnsum[col] = 0L; + gcolumnsum[col] = 0L; + bcolumnsum[col] = 0L; + } + + pnm_writepnminit(stdout, cols, rows, maxval, newformat, 0); + + /* Read in one convolution-matrix's worth of image, less one row. */ + for (row = 0; row < crows - 1; ++row) { + pnm_readpnmrow(ifp, xelbuf[row], cols, maxval, format); + if (PNM_FORMAT_TYPE(format) != newformat) + pnm_promoteformatrow(xelbuf[row], cols, maxval, format, + maxval, newformat); + /* Write out just the part we're not going to convolve. */ + if (row < crowso2) + pnm_writepnmrow(stdout, xelbuf[row], cols, maxval, newformat, 0); + } + + /* Now the rest of the image - read in the row at the end of + ** xelbuf, and convolve and write out the row in the middle. + */ + /* For first row only */ + + toprow = row + 1; + temprow = row % crows; + pnm_readpnmrow(ifp, xelbuf[temprow], cols, maxval, format); + if (PNM_FORMAT_TYPE(format) != newformat) + pnm_promoteformatrow(xelbuf[temprow], cols, maxval, format, maxval, + newformat); + + /* Arrange rowptr to eliminate the use of mod function to determine + ** which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = toprow % crows; + i = 0; + for (irow = temprow; irow < crows; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for (col = 0; col < cols; ++col) { + if (col < ccolso2 || col >= cols - ccolso2) + outputrow[col] = rowptr[crowso2][col]; + else if (col == ccolso2) { + int const leftcol = col - ccolso2; + float rsum, gsum, bsum; + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + for (crow = 0; crow < crows; ++crow) { + temprptr = rowptr[crow] + leftcol; + for (ccol = 0; ccol < ccols; ++ccol) { + rcolumnsum[leftcol + ccol] += + PPM_GETR(*(temprptr + ccol)); + gcolumnsum[leftcol + ccol] += + PPM_GETG(*(temprptr + ccol)); + bcolumnsum[leftcol + ccol] += + PPM_GETB(*(temprptr + ccol)); + } + } + for (ccol = 0; ccol < ccols; ++ccol) { + rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol]; + gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol]; + bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN(outputrow[col], r, g, b); + } else { + int const leftcol = col - ccolso2; + float rsum, gsum, bsum; + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + addcol = col + ccolso2; + for (crow = 0; crow < crows; ++crow) { + rcolumnsum[addcol] += PPM_GETR( rowptr[crow][addcol]); + gcolumnsum[addcol] += PPM_GETG( rowptr[crow][addcol]); + bcolumnsum[addcol] += PPM_GETB( rowptr[crow][addcol]); + } + for (ccol = 0; ccol < ccols; ++ccol) { + rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol]; + gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol]; + bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN(outputrow[col], r, g, b); + } + } + pnm_writepnmrow(stdout, outputrow, cols, maxval, newformat, 0); + + /* For all subsequent rows */ + subrow = crows; + addrow = crows - 1; + crowsp1 = crows + 1; + ++row; + for (; row < rows; ++row) { + toprow = row + 1; + temprow = row % (crows +1); + pnm_readpnmrow(ifp, xelbuf[temprow], cols, maxval, format); + if (PNM_FORMAT_TYPE(format) != newformat) + pnm_promoteformatrow(xelbuf[temprow], cols, maxval, format, + maxval, newformat); + + /* Arrange rowptr to eliminate the use of mod function to determine + ** which row of xelbuf is 0...crows. Mod function can be very costly. + */ + temprow = (toprow + 1) % crowsp1; + i = 0; + for (irow = temprow; irow < crowsp1; ++i, ++irow) + rowptr[i] = xelbuf[irow]; + for (irow = 0; irow < temprow; ++irow, ++i) + rowptr[i] = xelbuf[irow]; + + for (col = 0; col < cols; ++col) { + if (col < ccolso2 || col >= cols - ccolso2) + outputrow[col] = rowptr[crowso2][col]; + else if (col == ccolso2) { + int const leftcol = col - ccolso2; + float rsum, gsum, bsum; + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + + for (ccol = 0; ccol < ccols; ++ccol) { + tempcol = leftcol + ccol; + rcolumnsum[tempcol] = rcolumnsum[tempcol] + - PPM_GETR(rowptr[subrow][ccol]) + + PPM_GETR(rowptr[addrow][ccol]); + rsum = rsum + rcolumnsum[tempcol] * rweights[0][ccol]; + gcolumnsum[tempcol] = gcolumnsum[tempcol] + - PPM_GETG(rowptr[subrow][ccol]) + + PPM_GETG(rowptr[addrow][ccol]); + gsum = gsum + gcolumnsum[tempcol] * gweights[0][ccol]; + bcolumnsum[tempcol] = bcolumnsum[tempcol] + - PPM_GETB(rowptr[subrow][ccol]) + + PPM_GETB(rowptr[addrow][ccol]); + bsum = bsum + bcolumnsum[tempcol] * bweights[0][ccol]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN(outputrow[col], r, g, b); + } else { + int const leftcol = col - ccolso2; + float rsum, gsum, bsum; + rsum = 0.0; + gsum = 0.0; + bsum = 0.0; + addcol = col + ccolso2; + rcolumnsum[addcol] = rcolumnsum[addcol] + - PPM_GETR(rowptr[subrow][addcol]) + + PPM_GETR(rowptr[addrow][addcol]); + gcolumnsum[addcol] = gcolumnsum[addcol] + - PPM_GETG(rowptr[subrow][addcol]) + + PPM_GETG(rowptr[addrow][addcol]); + bcolumnsum[addcol] = bcolumnsum[addcol] + - PPM_GETB(rowptr[subrow][addcol]) + + PPM_GETB(rowptr[addrow][addcol]); + for (ccol = 0; ccol < ccols; ++ccol) { + rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol]; + gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol]; + bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol]; + } + temprsum = rsum + 0.5; + tempgsum = gsum + 0.5; + tempbsum = bsum + 0.5; + CHECK_RED; + CHECK_GREEN; + CHECK_BLUE; + PPM_ASSIGN(outputrow[col], r, g, b); + } + } + pnm_writepnmrow(stdout, outputrow, cols, maxval, newformat, 0); + } + + /* Now write out the remaining unconvolved rows in xelbuf. */ + for (irow = crowso2 + 1; irow < crows; ++irow) + pnm_writepnmrow(stdout, rowptr[irow], cols, maxval, newformat, 0); + +} + + + +static void +determineConvolveType(xel * const * const cxels, + struct convolveType * const typeP) { +/*---------------------------------------------------------------------------- + Determine which form of convolution is best. The general form always + works, but with some special case convolution matrices, faster forms + of convolution are possible. + + We don't check for the case that one of the PPM colors can have + differing types. We handle only cases where all PPMs are of the same + special case. +-----------------------------------------------------------------------------*/ + int horizontal, vertical; + int tempcxel, rtempcxel, gtempcxel, btempcxel; + int crow, ccol; + + switch (PNM_FORMAT_TYPE(format)) { + case PPM_TYPE: + horizontal = TRUE; /* initial assumption */ + crow = 0; + while (horizontal && (crow < crows)) { + ccol = 1; + rtempcxel = PPM_GETR(cxels[crow][0]); + gtempcxel = PPM_GETG(cxels[crow][0]); + btempcxel = PPM_GETB(cxels[crow][0]); + while (horizontal && (ccol < ccols)) { + if ((PPM_GETR(cxels[crow][ccol]) != rtempcxel) || + (PPM_GETG(cxels[crow][ccol]) != gtempcxel) || + (PPM_GETB(cxels[crow][ccol]) != btempcxel)) + horizontal = FALSE; + ++ccol; + } + ++crow; + } + + vertical = TRUE; /* initial assumption */ + ccol = 0; + while (vertical && (ccol < ccols)) { + crow = 1; + rtempcxel = PPM_GETR(cxels[0][ccol]); + gtempcxel = PPM_GETG(cxels[0][ccol]); + btempcxel = PPM_GETB(cxels[0][ccol]); + while (vertical && (crow < crows)) { + if ((PPM_GETR(cxels[crow][ccol]) != rtempcxel) | + (PPM_GETG(cxels[crow][ccol]) != gtempcxel) | + (PPM_GETB(cxels[crow][ccol]) != btempcxel)) + vertical = FALSE; + ++crow; + } + ++ccol; + } + break; + + default: + horizontal = TRUE; /* initial assumption */ + crow = 0; + while (horizontal && (crow < crows)) { + ccol = 1; + tempcxel = PNM_GET1(cxels[crow][0]); + while (horizontal && (ccol < ccols)) { + if (PNM_GET1(cxels[crow][ccol]) != tempcxel) + horizontal = FALSE; + ++ccol; + } + ++crow; + } + + vertical = TRUE; /* initial assumption */ + ccol = 0; + while (vertical && (ccol < ccols)) { + crow = 1; + tempcxel = PNM_GET1(cxels[0][ccol]); + while (vertical && (crow < crows)) { + if (PNM_GET1(cxels[crow][ccol]) != tempcxel) + vertical = FALSE; + ++crow; + } + ++ccol; + } + break; + } + + /* Which type do we have? */ + if (horizontal && vertical) { + typeP->ppmConvolver = ppm_mean_convolve; + typeP->pgmConvolver = pgm_mean_convolve; + } else if (horizontal) { + typeP->ppmConvolver = ppm_horizontal_convolve; + typeP->pgmConvolver = pgm_horizontal_convolve; + } else if (vertical) { + typeP->ppmConvolver = ppm_vertical_convolve; + typeP->pgmConvolver = pgm_vertical_convolve; + } else { + typeP->ppmConvolver = ppm_general_convolve; + typeP->pgmConvolver = pgm_general_convolve; + } +} + + + +static void +convolveIt(int const format, + struct convolveType const convolveType, + const float** const rweights, + const float** const gweights, + const float** const bweights) { + + switch (PNM_FORMAT_TYPE(format)) { + case PPM_TYPE: + convolveType.ppmConvolver(rweights, gweights, bweights); + break; + + default: + convolveType.pgmConvolver(gweights); + } +} + + + +int +main(int argc, char * argv[]) { + + struct cmdlineInfo cmdline; + FILE* cifp; + xel** cxels; + int cformat; + xelval cmaxval; + struct convolveType convolveType; + float ** rweights; + float ** gweights; + float ** bweights; + + pnm_init(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + cifp = pm_openr(cmdline.kernelFilespec); + + /* Read in the convolution matrix. */ + cxels = pnm_readpnm(cifp, &ccols, &crows, &cmaxval, &cformat); + pm_close(cifp); + + if (ccols % 2 != 1 || crows % 2 != 1) + pm_error("the convolution matrix must have an odd number of " + "rows and columns" ); + + ccolso2 = ccols / 2; + crowso2 = crows / 2; + + ifp = pm_openr(cmdline.inputFilespec); + + pnm_readpnminit(ifp, &cols, &rows, &maxval, &format); + if (cols < ccols || rows < crows) + pm_error("the image is smaller than the convolution matrix" ); + + newformat = MAX(PNM_FORMAT_TYPE(cformat), PNM_FORMAT_TYPE(format)); + if (PNM_FORMAT_TYPE(cformat) != newformat) + pnm_promoteformat(cxels, ccols, crows, cmaxval, cformat, + cmaxval, newformat); + if (PNM_FORMAT_TYPE(format) != newformat) { + switch (PNM_FORMAT_TYPE(newformat)) { + case PPM_TYPE: + if (PNM_FORMAT_TYPE(format) != newformat) + pm_message("promoting to PPM"); + break; + case PGM_TYPE: + if (PNM_FORMAT_TYPE(format) != newformat) + pm_message("promoting to PGM"); + break; + } + } + + computeWeights(cxels, ccols, crows, newformat, cmaxval, !cmdline.nooffset, + &rweights, &gweights, &bweights); + + /* Handle certain special cases when runtime can be improved. */ + + determineConvolveType(cxels, &convolveType); + + convolveIt(format, convolveType, + (const float **)rweights, + (const float **)gweights, + (const float **)bweights); + + pm_close(stdout); + pm_close(ifp); + return 0; +} + + + +/****************************************************************************** + SOME CHANGE HISTORY +******************************************************************************* + + Version 2.0.1 Changes + --------------------- + Fixed four lines that were improperly allocated as sizeof( float ) when they + should have been sizeof( long ). + + Version 2.0 Changes + ------------------- + + Version 2.0 was written by Mike Burns (derived from Jef Poskanzer's + original) in January 1995. + + Reduce run time by general optimizations and handling special cases of + convolution matrices. Program automatically determines if convolution + matrix is one of the types it can make use of so no extra command line + arguments are necessary. + + Examples of convolution matrices for the special cases are + + Mean Horizontal Vertical + x x x x x x x y z + x x x y y y x y z + x x x z z z x y z + + I don't know if the horizontal and vertical ones are of much use, but + after working on the mean convolution, it gave me ideas for the other two. + + Some other compiler dependent optimizations + ------------------------------------------- + Created separate functions as code was getting too large to put keep both + PGM and PPM cases in same function and also because SWITCH statement in + inner loop can take progressively more time the larger the size of the + convolution matrix. GCC is affected this way. + + Removed use of MOD (%) operator from innermost loop by modifying manner in + which the current xelbuf[] is chosen. + + This is from the file pnmconvol.README, dated August 1995, extracted in + April 2000, which was in the March 1994 Netpbm release: + + ----------------------------------------------------------------------------- + This is a faster version of the pnmconvol.c program that comes with netpbm. + There are no changes to the command line arguments, so this program can be + dropped in without affecting the way you currently run it. An updated man + page is also included. + + My original intention was to improve the running time of applying a + neighborhood averaging convolution matrix to an image by using a different + algorithm, but I also improved the run time of performing the general + convolution by optimizing that code. The general convolution runs in 1/4 to + 1/2 of the original time and neighborhood averaging runs in near constant + time for the convolution masks I tested (3x3, 5x5, 7x7, 9x9). + + Sample times for two computers are below. Times are in seconds as reported + by /bin/time for a 512x512 pgm image. + + Matrix IBM RS6000 SUN IPC + Size & Type 220 + + 3x3 + original pnmconvol 6.3 18.4 + new general case 3.1 6.0 + new average case 1.8 2.6 + + 5x5 + original pnmconvol 11.9 44.4 + new general case 5.6 11.9 + new average case 1.8 2.6 + + 7x7 + original pnmconvol 20.3 82.9 + new general case 9.4 20.7 + new average case 1.8 2.6 + + 9x9 + original pnmconvol 30.9 132.4 + new general case 14.4 31.8 + new average case 1.8 2.6 + + + Send all questions/comments/bugs to me at burns@chem.psu.edu. + + - Mike + + ---------------------------------------------------------------------------- + Mike Burns System Administrator + burns@chem.psu.edu Department of Chemistry + (814) 863-2123 The Pennsylvania State University + ---------------------------------------------------------------------------- + +*/ |