cleanup

git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@1012 9d0c8265-081b-0410-96cb-a4ca84ce46f8

2 files changed, 908 insertions, 719 deletions

diff --git a/editor/pamflip.c b/editor/pamflip.c
index 1c479e54..61c6efce 100644
--- a/editor/pamflip.c
+++ b/editor/pamflip.c
@@ -745,7 +745,7 @@ typedef struct {
 /*----------------------------------------------------------------------------
    A description of the quilt of cells that make up the output image.
 -----------------------------------------------------------------------------*/
-    unsigned int ranks, files;
+    unsigned int rankCt, fileCt;
         /* Dimensions of the output image in cells */
     struct pam ** pam;
         /* pam[y][x] is the pam structure for the cell at rank y, file x
@@ -769,8 +769,12 @@ initOutCell(struct pam *     const outCellPamP,
             unsigned int     const inCellHeight,
             struct pam *     const inpamP,
             struct xformCore const xformCore) {
-
-    unsigned int outCellFiles, outCellRanks;
+/*----------------------------------------------------------------------------
+   Set up an output cell.  Create and open a temporary file to hold its
+   raster.  Figure out the dimensions of the cell.  Return a PAM structure
+   that describes the cell (including identifying that temporary file).
+-----------------------------------------------------------------------------*/
+    unsigned int outCellFileCt, outCellRankCt;
 
     *outCellPamP = *inpamP;  /* initial value */
 
@@ -779,10 +783,10 @@ initOutCell(struct pam *     const outCellPamP,
     outCellPamP->file = pm_tmpfile();
 
     xformDimensions(xformCore, inCellWidth, inCellHeight,
-                    &outCellFiles, &outCellRanks);
+                    &outCellFileCt, &outCellRankCt);
 
-    outCellPamP->width = outCellFiles;
-    outCellPamP->height = outCellRanks;
+    outCellPamP->width = outCellFileCt;
+    outCellPamP->height = outCellRankCt;
 }
 
 
@@ -792,9 +796,26 @@ createOutputMap(struct pam *       const inpamP,
                 unsigned int       const maxRows,
                 struct xformMatrix const cellXform,
                 struct xformCore   const xformCore) {
-
-    unsigned int const inCellFiles = 1;
-    unsigned int const inCellRanks = (inpamP->height + maxRows - 1) / maxRows;
+/*----------------------------------------------------------------------------
+   Create and return the output map.  That's a map of all the output cells
+   (from which the output image can be assembled, once those cells are filled
+   in).
+
+   The map contains the dimensions of each output cell as well as file
+   stream handles for the temporary files containing the pixels of each
+   cell.  We create and open those files.
+
+   Note that a complexity of determining cell dimensions (which we handle)
+   is that the input image isn't necessarily a whole multiple of the input
+   cell size, so the last cell may be short.
+
+   The map does not contain the mapping from input cells to output cells
+   (e.g. in a top-bottom transformation of a 10-cell image, input cell
+   0 maps to output cell 9); caller can compute that when needed from
+   'cellXform'.
+-----------------------------------------------------------------------------*/
+    unsigned int const inCellFileCt = 1;
+    unsigned int const inCellRankCt = (inpamP->height + maxRows - 1) / maxRows;
 
     outputMap * mapP;
     unsigned int outCellRank;
@@ -802,27 +823,27 @@ createOutputMap(struct pam *       const inpamP,
 
     MALLOCVAR_NOFAIL(mapP);
 
-    xformDimensions(xformCore, inCellFiles, inCellRanks,
-                    &mapP->files, &mapP->ranks);
+    xformDimensions(xformCore, inCellFileCt, inCellRankCt,
+                    &mapP->fileCt, &mapP->rankCt);
 
-    MALLOCARRAY(mapP->pam, mapP->ranks);
+    MALLOCARRAY(mapP->pam, mapP->rankCt);
     if (mapP->pam == NULL)
         pm_error("Could not allocate a cell array for %u ranks of cells.",
-                 mapP->ranks);
+                 mapP->rankCt);
 
-    for (outCellRank = 0; outCellRank < mapP->ranks; ++outCellRank) {
+    for (outCellRank = 0; outCellRank < mapP->rankCt; ++outCellRank) {
 
-        MALLOCARRAY(mapP->pam[outCellRank], mapP->files);
+        MALLOCARRAY(mapP->pam[outCellRank], mapP->fileCt);
 
         if (mapP->pam[outCellRank] == NULL)
             pm_error("Failed to allocate rank %u of the cell array, "
-                     "%u cells wide", outCellRank, mapP->files);
+                     "%u cells wide", outCellRank, mapP->fileCt);
     }
 
-    for (inCellRank = 0; inCellRank < inCellRanks; ++inCellRank) {
+    for (inCellRank = 0; inCellRank < inCellRankCt; ++inCellRank) {
         unsigned int const inCellFile = 0;
         unsigned int const inCellStartRow = inCellRank * maxRows;
-        unsigned int const inCellRows =
+        unsigned int const inCellRowCt =
             MIN(inpamP->height - inCellStartRow, maxRows);
 
         unsigned int outCellFile, outCellRank;
@@ -830,7 +851,7 @@ createOutputMap(struct pam *       const inpamP,
                        &outCellFile, &outCellRank);
     
         initOutCell(&mapP->pam[outCellRank][outCellFile],
-                    inpamP->width, inCellRows,
+                    inpamP->width, inCellRowCt,
                     inpamP, xformCore);
     }
     return mapP;
@@ -843,7 +864,7 @@ destroyOutputMap(outputMap * const mapP) {
 
     unsigned int outCellRank;
 
-    for (outCellRank = 0; outCellRank < mapP->ranks; ++outCellRank)
+    for (outCellRank = 0; outCellRank < mapP->rankCt; ++outCellRank)
         free(mapP->pam[outCellRank]);
 
     free(mapP->pam);
@@ -858,9 +879,9 @@ rewindCellFiles(outputMap * const outputMapP) {
 
     unsigned int outCellRank;
 
-    for (outCellRank = 0; outCellRank < outputMapP->ranks; ++outCellRank) {
+    for (outCellRank = 0; outCellRank < outputMapP->rankCt; ++outCellRank) {
         unsigned int outCellFile;
-        for (outCellFile = 0; outCellFile < outputMapP->files; ++outCellFile)
+        for (outCellFile = 0; outCellFile < outputMapP->fileCt; ++outCellFile)
             pm_seek(outputMapP->pam[outCellRank][outCellFile].file, 0);
     }
 }
@@ -872,9 +893,9 @@ closeCellFiles(outputMap * const outputMapP) {
 
     unsigned int outCellRank;
 
-    for (outCellRank = 0; outCellRank < outputMapP->ranks; ++outCellRank) {
+    for (outCellRank = 0; outCellRank < outputMapP->rankCt; ++outCellRank) {
         unsigned int outCellFile;
-        for (outCellFile = 0; outCellFile < outputMapP->files; ++outCellFile)
+        for (outCellFile = 0; outCellFile < outputMapP->fileCt; ++outCellFile)
             pm_close(outputMapP->pam[outCellRank][outCellFile].file);
     }
 }
@@ -932,7 +953,7 @@ stitchCellsToOutput(outputMap *  const outputMapP,
 
     tupleRow = pnm_allocpamrow(outpamP);
 
-    for (outRank = 0; outRank < outputMapP->ranks; ++outRank) {
+    for (outRank = 0; outRank < outputMapP->rankCt; ++outRank) {
         unsigned int const cellRows = outputMapP->pam[outRank][0].height;
             /* Number of rows in every cell in this rank */
 
@@ -944,7 +965,7 @@ stitchCellsToOutput(outputMap *  const outputMapP,
 
             outCol = 0;
 
-            for (outFile = 0; outFile < outputMapP->files; ++outFile) {
+            for (outFile = 0; outFile < outputMapP->fileCt; ++outFile) {
                 struct pam * const outCellPamP = 
                     &outputMapP->pam[outRank][outFile];
 
@@ -973,7 +994,7 @@ transformNonPbmChunk(struct pam *     const inpamP,
                      unsigned int     const maxRows,
                      bool             const verbose) {
 /*----------------------------------------------------------------------------
-  Same as transformNonPbmChunk(), except we read 'maxRows' rows of the
+  Same as transformNonPbmWhole(), except we read 'maxRows' rows of the
   input into memory at a time, storing intermediate results in temporary
   files, to limit our use of virtual and real memory.
 
@@ -981,28 +1002,33 @@ transformNonPbmChunk(struct pam *     const inpamP,
   header).
 
   We call the strip of 'maxRows' rows that we read a source cell.  We
-  transform that cell according to 'xformCore' to create to create a
+  transform that cell according to 'xformCore' to create a
   target cell.  We store all the target cells in temporary files.
   We consider the target cells to be arranged in a column matrix the
   same as the source cells within the source image; we transform that
   matrix according to 'xformCore'.  The resulting cell matrix is the
   target image.
 -----------------------------------------------------------------------------*/
-    unsigned int const inCellFiles = 1;
-    unsigned int const inCellRanks = (inpamP->height + maxRows - 1) / maxRows;
+    /* The cells of the source image ("inCell") are in a 1-column matrix.
+       "rank" is the vertical position of a cell in that matrix; "file" is
+       the horizontal position (always 0, of course).
+    */
+    unsigned int const inCellFileCt = 1;
+    unsigned int const inCellRankCt = (inpamP->height + maxRows - 1) / maxRows;
 
     struct xformMatrix cellXform;
     unsigned int inCellRank;
     outputMap * outputMapP;
 
     if (verbose)
-        pm_message("Transforming in %u chunks, using temp files", inCellRanks);
+        pm_message("Transforming in %u chunks, using temp files",
+                   inCellRankCt);
 
-    computeXformMatrix(&cellXform, inCellFiles, inCellRanks, xformCore);
+    computeXformMatrix(&cellXform, inCellFileCt, inCellRankCt, xformCore);
 
     outputMapP = createOutputMap(inpamP, maxRows, cellXform, xformCore);
 
-    for (inCellRank = 0; inCellRank < inCellRanks; ++inCellRank) {
+    for (inCellRank = 0; inCellRank < inCellRankCt; ++inCellRank) {
         unsigned int const inCellFile = 0;
         unsigned int const inCellStartRow = inCellRank * maxRows;
         unsigned int const inCellRows =
diff --git a/editor/pnmconvol.c b/editor/pnmconvol.c
index 031feb66..9d3a3ce9 100644
--- a/editor/pnmconvol.c
+++ b/editor/pnmconvol.c
@@ -1,6 +1,4 @@
-/* pnmconvol.c - general MxN convolution on a PNM image
-**
-** Version 2.0.1 January 30, 1995
+/* pnmconvol.c - general MxN convolution on a Netpbm image
 **
 ** Major rewriting by Mike Burns
 ** Copyright (C) 1994, 1995 by Mike Burns (burns@chem.psu.edu)
@@ -22,6 +20,7 @@
 #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.
@@ -81,53 +80,35 @@ parseCommandLine(int argc, char ** argv,
 }
 
 
-/* Macros to verify that r,g,b values are within proper range */
+struct convKernel {
+    unsigned int const cols;
+    unsigned int const rows;
+    float ** weight[3];
+};
 
-#define CHECK_GRAY \
-    if (tempgsum < 0L) g = 0; \
-    else if (tempgsum > maxval) g = maxval; \
-    else g = tempgsum;
+convKernelDestroy(struct convKernel * const convKernelP) {
 
-#define CHECK_RED \
-    if (temprsum < 0L) r = 0; \
-    else if (temprsum > maxval) r = maxval; \
-    else r = temprsum;
+    unsigned int plane;
 
-#define CHECK_GREEN \
-    if (tempgsum < 0L) g = 0; \
-    else if (tempgsum > maxval) g = maxval; \
-    else g = tempgsum;
+    for (plane = 0; plane < 3; ++plane)
+        pm_freearray(convKernelP->weight[plane]);
 
-#define CHECK_BLUE \
-    if (tempbsum < 0L) b = 0; \
-    else if (tempbsum > maxval) b = maxval; \
-    else b = tempbsum;
+    free(convKernelP);
+}
 
 struct convolveType {
-    void (*ppmConvolver)(const float ** const rweights,
-                         const float ** const gweights,
-                         const float ** const bweights);
-    void (*pgmConvolver)(const float ** const weights);
+    void (*convolver)(struct pam *       const inpamP,
+                      struct pam *       const outpamP,
+                      const convKernel * const convKernelP);
 };
 
-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) {
+computeKernel(struct pam *         const cpamP
+              tuple * const *      const ctuples, 
+              bool                 const offsetPgm,
+              struct convKernel ** const convKernelP) {
 /*----------------------------------------------------------------------------
    Compute the convolution matrix in normalized form from the PGM
    form.  Each element of the output matrix is the actual weight we give an
@@ -139,728 +120,906 @@ computeWeights(xel * const *   const cxels,
    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 scale = (offsetPgm ? 2.0 : 1.0) / cpamP->maxval;
     double const offset = offsetPgm ? - 1.0 : 0.0;
 
-    float** rweights;
-    float** gweights;
-    float** bweights;
+    struct convKernel * convKernelP;
+    float sum[3];
+    unsigned int plane;
+    unsigned int row;
 
-    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 */
+    MALLOCVAR_NOFAIL(convKernelP);
     
-    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;
-            }
+    convKernelP->weight[PAM_RED_PLANE] =
+        pm_allocarray(cpamP->width, cpamP->height,
+                      sizeof(convKernelP->weight[PAM_RED_PLANE][0][0]));
+    
+    for (plane = 0; plane < 3; ++plane)
+        sum[plane] = 0.0; /* initial value */
+    
+    for (row = 0; row < cpamP->height; ++row) {
+        unsigned int col;
+        for (col = 0; col < cpamP->width; ++col) {
+            unsigned int plane;
+            for (plane = 0; plane < cpamP->depth; ++plane) {
+                sum += convKernelP->weight[plane][row][col] =
+                    ctuples[row][col][plane] * scale + offset);
         }
     }
-    *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) {
+    switch (PNM_FORMAT_TYPE(cpamP->format)) {
+    case PPM_TYPE: {
+        unsigned int plane;
+        bool biased, negative;
+        for (plane = 0, biased = false, negative = false;
+             plane < cpamP->depth;
+             ++plane) {
+            if (sum[plane] < 0.9 || sum[plane] > 1.1)
+                biased = true;
+            if (sum[plane] < 0.0)
+                negative = true;
+        }
+    
+        if (biased) {
             pm_message("WARNING - this convolution matrix is biased.  " 
                        "red, green, and blue average weights: %f, %f, %f "
                        "(unbiased would be 1).",
-                       rsum, gsum, bsum);
+                       sum[PAM_RED_PLANE],
+                       sum[PAM_GRN_PLANE],
+                       sum[PAM_BLU_PLANE]);
 
-            if (rsum < 0 && gsum < 0 && bsum < 0)
+            if (negative)
                 pm_message("Maybe you want the -nooffset option?");
         }
-        break;
+    } break;
 
     default:
-        if (gsum < 0.9 || gsum > 1.1)
+        if (sum[0] < 0.9 || sum[0] > 1.1)
             pm_message("WARNING - this convolution matrix is biased.  "
                        "average weight = %f (unbiased would be 1)",
-                       gsum);
+                       sum[0]);
         break;
     }
 }
 
 
 
-/* General PGM Convolution
-**
-** No useful redundancy in convolution matrix.
-*/
+static tuple **
+allocRowbuf(struct pam * const pamP,
+            unsigned int const height) {
+
+    tuple ** rowbuf;
+
+    MALLOCARRAY(rowbuf, height);
+
+    if (rowbuf == NULL)
+        pm_error("Failed to allocate %u-row buffer");
+    else {
+        unsigned int row;
+    
+        for (row = 0; row < height; ++row)
+            rowbuf[i] = pnm_allocpamrow(pamP);
+    }
+
+    return rowbuf;
+}
+
+
 
 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;
+freeRowbuf(tuple **     const rowbuf,
+           unsigned int const height) {
 
-    /* Allocate space for one convolution-matrix's worth of rows, plus
-       a row output buffer.
-    */
-    xelbuf = pnm_allocarray(cols, crows);
-    outputrow = pnm_allocrow(cols);
+    unsigned int row;
 
-    /* Allocate array of pointers to xelbuf */
-    rowptr = (xel **) pnm_allocarray(1, crows);
+    for (row = 0; row < height; ++row)
+        pnm_freepamrow(rowbuf[row]);
 
-    pnm_writepnminit(stdout, cols, rows, maxval, newformat, 0);
+    free(rowbuf);
+}
 
-    /* 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];
+static void
+readAndScaleRow(struct pam * const inpamP,
+                tuple *      const inrow,
+                sample       const newMaxval,
+                unsigned int const newDepth) {
 
-        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);
-    }
+    pnm_readpamrow(inpamP, inrow);
 
-    /* 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 );
+    if (newMaxval != inpamP->maxval)
+        pnm_scaletuplerow(inpamP, inrow, inrow, newMaxval);
+
+    if (newDepth == 3 && inpamP->depth == 1)
+        pnm_makerowrgb(inpamP, inrow);
 }
 
 
 
-/* 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
+readInitialRowbuf(struct pam *        const inpamP,
+                  struct convKernel * const convKernelP,
+                  tuple **            const rowbuf,
+                  sample              const outputMaxval,
+                  unsigned int        const outputDepth) {
+/*----------------------------------------------------------------------------
+  Read in one convolution kernel's worth of image, less one row,
+  into the row buffer rowbuf[], starting at the beginning of the buffer.
+  
+  Scale the contents to maxval 'outputMaxval'.
+-----------------------------------------------------------------------------*/
+    unsigned int row;
+
+    for (row = 0; row < convKernelP->rows - 1; ++row)
+        readAndScaleRow(inpamP, rowbuf[row], outputMaxval, outputDepth);
+}
+
 
 
 static void
-pgm_mean_convolve(const float ** const weights) {
-    float const gmeanweight = weights[0][0];
+writeUnconvolvedTop(struct pam *        const outpamP,
+                    struct convKernel * const convKernelP,
+                    tuple **            const rowbuf) {
+/*----------------------------------------------------------------------------
+   Write out the top part that we can't convolve because the convolution
+   kernel runs off the top of the image.
 
-    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;
+   Assume those rows are in the window rowbuf[], with the top row of the
+   image as the first row in rowbuf[].
+-----------------------------------------------------------------------------*/
+    unsigned int row;
 
-    /* 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 );
+    for (row = 0; row < convKernelP->rows/2; ++row)
+        pnm_writepamrow(outpamP, rowbuf[row]);
+}
 
-    /* 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 );
+static void
+writeUnconvolvedBottom(struct pam *        const outpamP,
+                       struct convKernel * const convKernelP,
+                       tuple **            const circMap) {
+/*----------------------------------------------------------------------------
+  Write out the bottom part that we can't convolve because the convolution
+  kernel runs off the bottom of the image.
 
-    /* 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 );
-    }
+  Assume the end of the image is in the row buffer, mapped by 'circMap'
+  such that the top of the window is circMap[0].
+-----------------------------------------------------------------------------*/
+    unsigned int row;
 
-    /* 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 );
+    for (row = convKernelP->rows / 2 + 1; row < convKernelP->rows; ++row)
+        pnm_writepamrow(outpamP, circMap[row]);
+}
+
+
+
+static void
+setupCircMap(tuple **     const circMap,
+             tuple **     const rowbuf,
+             unsigned int const topRowbufRow) {
+/*----------------------------------------------------------------------------
+  Set up circMap[] to reflect the case that index 'topRowbufRow' of rowbuf[]
+  is for the topmost row in the window.
+-----------------------------------------------------------------------------*/
+    unsigned int row;
+    unsigned int i;
 
-    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 );
+    for (row = topRowbufRow; row < convKernelP->rows; ++i, ++row)
+        circMap[i] = rowbuf[row];
+
+    for (row = 0; row < topRowbufRow; ++row, ++i)
+        circMap[i] = rowbuf[row];
+}
+
+
+
+static void
+convolveGeneralRowPlane(struct pam *        const pamP,
+                        tuple **            const circMap,
+                        struct convKernel * const convKernelP,
+                        unsigned int        const plane,
+                        tuple *             const outputrow) {
+
+    unsigned int const crowso2 = convKernelP->rows / 2;
+    unsigned int const ccolso2 = convKernelP->cols / 2;
+
+    unsigned int col;
+    
+    for (col = 0; col < pamP->width; ++col) {
+        if (col < ccolso2 || col >= pamP->width - ccolso2)
+            outputrow[col][plane] = circMap[crowso2][col][plane];
+        else {
+            unsigned int const leftcol = col - ccolso2;
+            unsigned int crow;
+            float sum;
+            sum = 0.0;
+            for (crow = 0; crow < crows; ++crow) {
+                tuple ** const leftrptr = &circMap[crow][leftcol];
+                unsigned int ccol;
+                for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                    sum += leftrptr[ccol][plane] *
+                        convKernelP->weight[plane][crow][ccol];
+            }
+            outputrow[col][plane] = MIN(pamP->maxval, MAX(0, sum + 0.5));
         }
     }
-    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.
+
+static void
+convolveGeneral(struct pam *        const inpamP,
+                struct pam *        const outpamP,
+                struct convKernel * const convKernelP) {
+/*----------------------------------------------------------------------------
+   Do the convolution without taking advantage of any useful redundancy in the
+   convolution matrix.
+-----------------------------------------------------------------------------*/
+    tuple ** rowbuf;
+        /* A vertical window of the input image.  It holds as many rows as the
+           convolution kernel covers -- the rows we're currently using to
+           create output rows.  It is a circular buffer.
+        */
+    tuple ** circMap;
+        /* A map from image row number within window to element of rowbuf[].
+           E.g. if rowbuf[] if 5 rows high and rowbuf[2] contains the
+           topmost row, then circMap[0] == 2, circMap[1] = 3,
+           circMap[4] = 1.  You could calculate the same thing with a mod
+           function, but that is sometimes more expensive.
+        */
+    tuple * outputrow;
+        /* The convolved row to be output */
+    unsigned int row;
+
+    rowbuf = allocRowbuf(pamP, convKernelP->rows);
+    MALLOCARRAY_NOFAIL(circMap, convKernelP->rows);
+    outputrow = pnm_allocpamrow(outpamP);
+
+    pnm_writepaminit(outpamP);
+
+    assert(convKernelP->rows > 0);
+
+    readInitialRowbuf(inpamP, convKernelP, rowbuf, outpamP->maxval);
+
+    writeUnconvolvedTop(outpamP, convKernelP, rowbuf);
+
+    /* Now the rest of the image - read in the row at the bottom of the
+       window, then convolve and write out the row in the middle of the
+       window.
     */
-    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 (row = convKernelP->rows - 1; row < inpamP->height; ++row) {
+        unsigned int const rowbufRow = row % convKernelP->rows;
 
-    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 );
+        setupCircMap(circMap, rowbuf, (row + 1) % convKernelP->rows);
+
+        readAndScaleRow(inpamP, rowbuf[rowbufRow],
+                        outpamP->maxval, outpamP->depth);
+
+        for (plane = 0; plane < outpamP->depth; ++plane)
+            convolveGeneralRowPlane(outpamP, circMap, convKernelP, plane,
+                                    outputrow);
+
+        pnm_writepamrow(outpamP, outputrow);
     }
+    writeUnconvolvedBottom(outpamP, convKernelP, circMap);
 
-    /* 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 );
+    freeRowbuf(rowbuf, convKernelP->rows);
+}
+
+
+
+static sample **
+allocSum(unsigned int const depth,
+         unsigned int const size) {
 
+    sample ** sum;
+
+    MALLOCARRAY(sum, depth);
+
+    if (!sum)
+        pm_error("Could not allocate memory for %u planes of sums", depth);
+    else {
+        unsigned int plane;
+
+        for (plane = 0; plane < depth; ++plane) {
+            MALLOCARRAY(sum[plane], size);
+            
+            if (!sum[plane])
+                pm_error("Could not allocate memory for %u sums", size);
+        }
     }
+    return sum;
+}
 
 
-/* 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;
+freeSum(sample **    const sum,
+        unsigned int const depth) {
 
-    /* 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 );
+    unsigned int plane;
 
-    /* Allocate array of pointers to xelbuf */
-    rowptr = (xel **) pnm_allocarray( 1, crows + 1);
+    for (plane = 0; plane < depth; ++plane)
+        free(sum[plane]);
 
-    /* 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);
+    free(sum);
+}
 
-    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 );
+
+static void
+convolveAndComputeColumnSums(struct pam *        const inpamP,
+                             struct pam *        const outpamP,
+                             tuple **            const rows,
+                             struct convKernel * const convKernelP,
+                             unsigned int        conts plane,
+                             tuple *             const outputrow,
+                             sample **           const convColumnSum) {
+
+    unsigned int const ccolso2 = convKernelP->cols / 2;
+    float        const weight  = convKernelP->weight[plane][0][0];
+
+    unsigned int col;
+    sample matrixSum;
+
+    for (col = 0; col < inpamP->width; ++col)
+        convColumnSum[col] = 0;      /* Initial value */
+
+    for (col = 0, matrixSum = 0; col < inpamP->width; ++col) {
+        if (col < ccolso2 || col >= inpamP->width - ccolso2)
+            outputrow[col] = circMap[crowso2][col];
+        else if (col == ccolso2) {
+            unsigned int const leftcol = col - ccolso2;
+            unsigned int crow;
+            unsigned int ccol;
+            for (crow = 0; crow < convKernelP->rows; ++crow) {
+                tuple ** const temprptr = &circMap[crow][leftcol];
+                unsigned int ccol;
+                for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                    convColumnSum[leftcol + ccol] += temprptr[ccol][plane];
+            }
+            for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                matrixSum += convColumnSum[leftcol + ccol];
+            outputrow[col][plane] =
+                MIN(outpamP->maxval, MAX(0, matrixSum * weight + 0.5));
+        } else {
+            /* Column numbers to subtract or add to isum */
+            unsigned int const subcol = col - ccolso2 - 1;
+            unsigned int const addcol = col + ccolso2;  
+            unsigned int crow;
+            for (crow = 0; crow < convKernelP->rows; ++crow)
+                convColumnSum[addcol] += circMap[crow][addcol][plane];
+            matrixSum =
+                matrixSum - convColumnSum[subcol] + convColumnSum[addcol];
+            outputrow[col][plane] =
+                MIN(outpamP->maxval, MAX(0, matrixSum * weight + 0.5));
+        }
     }
+}
 
-    /* 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 );
+static void
+convolveAndComputeColumnSums(struct pam *        const inpamP,
+                             struct pam *        const outpamP,
+                             tuple **            const rows,
+                             struct convKernel * const convKernelP,
+                             tuple *             const outputrow,
+                             sample **           const convColumnSum) {
+/*----------------------------------------------------------------------------
+  Convolve the rows in rows[] -- one convolution kernel's worth, where
+  rows[0] is the top.  Put the result in outputrow[].
 
+  Along the way, add up the sum of each column and return that as
+  convColumnSum[].
+-----------------------------------------------------------------------------*/
+    unsigned int plane;
 
-    /* For all subsequent rows */
+    for (plane = 0; plane < outpamP->depth; ++plane) 
+        convolveAndComputeColumnSumsPlane(inpamP, outpamP,
+                                          rows, convKernelP, plane,
+                                          outputrow, convColumnSum[plane]);
+}
 
-    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 );
-        }
+static void
+convolveMeanRowPlane(struct pam *        const pamP,
+                     tuple **            const circMap,
+                     struct convKernel * const convKernelP,
+                     unsigned int        const plane,
+                     tuple *             const outputrow,
+                     sample *            const convColumnSum) {
+
+    unsigned int const ccolso2 = convKernelP->cols / 2;
+    float const weight = convKernelP->weight[plane][0][0];
+
+    unsigned int col;
+    sample gisum;
+
+    gisum = 0;
+    for (col = 0; col < cols; ++col) {
+        if (col < ccolso2 || col >= cols - ccolso2)
+            outputrow[col] = circMap[crowso2][col];
+        else if (col == ccolso2) {
+            unsigned int const leftcol = col - ccolso2;
+
+            unsigned int ccol;
+
+            for (ccol = 0; ccol < convKernelP->cols; ++ccol) {
+                sample * const thisColumnSumP =
+                    &convColumnSum[leftcol + ccol];
+                *thisColumnSumP = *thisColumnSumP
+                    - circMap[subrow][ccol][plane]
+                    + circMap[addrow][ccol][plane];
+                gisum += *thisColumnSumP;
+            }
+            outputrow[col][plane] =
+                MIN(outpamP->maxval, MAX(0, gisum * weight + 0.5));
+        } else {
+            /* Column numbers to subtract or add to isum */
+            unsigned int const subcol = col - ccolso2 - 1;
+            unsigned int const addcol = col + ccolso2;  
+
+            convColumnSum[addcol] = convColumnSum[addcol]
+                - circMap[subrow][addcol][plane]
+                + circMap[addrow][addcol][plane];
+            gisum = gisum - convColumnSum[subcol] + convColumnSum[addcol];
+
+            outputrow[col][plane] =
+                MIN(outpamP->maxval, MAX(0, gisum * weight + 0.5));
         }
-    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
+convolveMean(struct pam *        const inpamP,
+             struct pam *        const outpamP,
+             struct convKernel * const convKernelP) {
+/*----------------------------------------------------------------------------
+  Mean Convolution
+
+  This is for 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 as follows: At a certain vertical position in the
+  image, create sums for each column fragment of the convolution height all
+  the way across the image.  Then add those sums across the convolution width
+  to obtain the total sum over the convolution area and multiply that sum by
+  the weight.  As you move left to right, to calculate the next output 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.
+
+  We assume the convolution kernel is uniform -- same weights everywhere.
+
+  We assume the output is PGM and the input is PGM or PBM.
+-----------------------------------------------------------------------------*/
+    unsigned int const windowHeight = convKernelP->rows + 1;
+        /* The height of the window we keep in the row buffer.  The buffer
+           contains the rows covered by the convolution kernel, plus the row
+           immediately above that.  The latter is there because to compute
+           the sliding mean, we need to subtract off the row that the
+           convolution kernel just slid past.
+        */
+
+    tuple ** rowbuf;
+        /* Same as in convolvePgmGeneral */
+    tuple ** circMap;
+        /* Same as in convolvePgmGeneral */
+    tuple * outputrow;
+        /* Same as in convolvePgmGeneral */
+    unsigned int row;
+        /* Row number of next row to read in from the file */
+    sample matrixSum;
+        /* Sum of all pixels in current convolution window */
+    sample ** convColumnSum;  /* Malloc'd */
+        /* convColumnSum[plane][col] is the sum of Plane 'plane' of all the
+           pixels in the Column 'col' of the image within the current vertical
+           convolution window.  I.e. if our convolution kernel is 5 rows high
+           and we're now looking at Rows 10-15, convColumn[0][3] is the sum of
+           Plane 0 of Column 3, Rows 10-15.
+        */
+    unsigned int col;
+
+    rowbuf = allocRowbuf(pamP, windowHeight);
+    MALLOCARRAY_NOFAIL(circMap, windowHeight);
+    outputrow = pnm_allocpamrow(outpamP);
+
+    convColumnSum = allocSum(outpamP->depth, outpamP->width);
+
+    pnm_writepaminit(outpamP);
+
+    readInitialRowbuf(inpamP, convKernelP, rowbuf, outpamP->maxval);
+
+    writeUnconvolvedTop(outpamP, convKernelP, rowbuf);
+
+    /* Add a row to the window to have enough to convolve */
+    readAndScaleRow(inpamP, rowbuf[convKernelP->rows - 1],
+                    outpamP->maxval, outpamP->depth);
+
+    setupCircMap(circMap, rowbuf, 0);
+
+    /* Convolve the first window the long way */
+    convolveAndComputeColumnSums(inpamP, outpamP, circMap, convKernelP,
+                                 outputrow, convColumnSum);
+
+    /* Write that first convolved row */
+    pnm_writepamrow(outpamP, outputrow);
+
+    /* For all subsequent rows do it this way as the columnsums have been
+       generated.  Now we can use them to reduce further calculations.
+    */
+    for (row = convKernelP->rows; row < inpamP->height; ++row) {
+        unsigned int const subrow = 0;
+            /* Row just above convolution window -- what we subtract from
+               running sum
+            */
+        unsigned int const addrow = 1 + (convKernelP->rows - 1);
+            /* Bottom row of convolution window: What we add to running sum */
+
+        readAndScaleRow(inpamP, rowbuf[row % windowHeight],
+                        outpamP->maxval, outpamP->depth);
+
+        /* Remember the window is one row higher than the convolution
+           kernel.  The top row in the window is not part of this convolution.
+        */
+
+        setupCircMap(circMap, rowbuf, (row + 1) % windowHeight);
+
+        for (plane = 0; plane < outpamP->depth; ++plane)
+            convolveMeanRowPlane(outpamP, circMap, convKernelP, plane,
+                                 outputrow, convColumnSum[plane]);
+
+        pnm_writepamrow(outpamP, outputrow);
     }
+    writeUnconvolvedBottom(outpamP, convKernelP, circMap);
 
+    freeColumnSum(convColumnSum, outpamP->depth);
+    freeRowbuf(rowbuf, windowHeight);
+}
 
-/* 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;
+convolveHorizontalRowPlane(struct pam *        const outpamP,
+                           tuple **            const circMap,
+                           struct convKernel * const convKernelP,
+                           unsigned int        const plane,
+                           tuple *             const outputrow,
+                           sample *            const sumCircMap) {
 
-    /* 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 );
+    unsigned int col;
 
-    /* Allocate array of pointers to xelbuf */
-    rowptr = (xel **) pnm_allocarray( 1, crows + 1 );
+    for (col = 0; col < cols; ++col) {
+        if (col < ccolso2 || col >= cols - ccolso2)
+            outputrow[col][0] = circMap[crowso2][col][0];
+        else if (col == ccolso2) {
+            unsigned int const leftcol = col - ccolso2;
+            
+            float matrixSum;
+            unsigned int crow;
 
-    /* Allocate space for intermediate column sums */
-    gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
-    for ( col = 0; col < cols; ++col )
-    gcolumnsum[col] = 0L;
+            for (crow = 0, matrixSum = 0.0; crow < crows; ++crow) {
+                tuple ** const temprptr = circMap[crow] + leftcol;
 
-    pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
+                unsigned int ccol;
+                
+                sumCircMap[crow][leftcol] = 0L;
+                for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                    sumCircMap[crow][leftcol] += temprptr[ccol][0];
+                matrixSum += sumCircMap[crow][leftcol] *
+                    convKernelP->weight[crow][0][0];
+            }
+            outputrow[col][0] =
+                MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+        } else {
+            unsigned int const leftcol = col - ccolso2;
+            unsigned int const subcol  = col - ccolso2 - 1;
+            usnigned int const addcol  = col + ccolso2;
 
-    /* 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 );
+            float matrixSum;
+            unsigned int crow;
+
+            matrixSum = 0.0;
+
+            for (crow = 0; crow < crows; ++crow) {
+                sumCircMap[crow][leftcol] = sumCircMap[crow][subcol]
+                    - circMap[crow][subcol][0]
+                    + circMap[crow][addcol][0];
+                matrixSum += sumCircMap[crow][leftcol] *
+                    convKernelP->weight[crow][0][0];
+            }
+            outputrow[col][0] =
+                MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+        }
     }
+}
 
-    /* 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];
+static void
+convolveHorizontal(struct pam *        const inpamP,
+                   struct pam *        const outpamP,
+                   struct convKernel * const convKernelP) {
+/*----------------------------------------------------------------------------
+  Horizontal Convolution
 
-    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 );
-        }
+  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.
+-----------------------------------------------------------------------------*/
+    unsigned int const ccolso2 = convKernelP->cols / 2;
+    unsigned int const crowso2 = convKernelP->rows / 2;
+    unsigned int const windowHeight = convKernelP->rows + 1;
+        /* Same as in convolvePgmMean */
+
+    tuple ** rowbuf;
+        /* Same as in convolvePgmGeneral */
+    tuple ** circMap;
+        /* Same as in convolvePgmGeneral */
+    tuple * outputrow;
+        /* Same as in convolvePgmGeneral */
+    unsigned int row;
+        /* Row number of next row to read in from the file */
+    unsigned int plane;
+    sample ** convRowSum;  /* Malloc'd */
+    sample * sumCircMap;  /* Malloc'd */
+    unsigned int plane;
+
+    rowbuf = allocRowbuf(pamP, windowHeight);
+    MALLOCARRAY_NOFAIL(circMap, windowHeight);
+    outputrow = pnm_allocpamrow(outpamP);
+
+    convRowSum = allocSum(outpamP->depth, windowHeight);
+    MALLOCARRAY_NOFAIL(sumCircMap, windowHeight);
+
+    pnm_writepaminit(outpamP);
+
+    readInitialRowbuf(inpamP, convKernelP, rowbuf, outpamP->maxval);
+
+    writeUnconvolvedTop(outpamP, convKernelP, rowbuf);
+
+    /* Add a row to the window to have enough to convolve */
+    readAndScaleRow(inpamP, rowbuf[convKernelP->rows - 1],
+                    outpamP->maxval, outpamP->depth);
+
+    setupCircMap(circMap, rowbuf, 0);
+
+    for (plane = 0; plane < outpamP->depth; ++plane) {
+        unsigned int crow;
+
+        for (crow = 0; crow < crows; ++crow)
+            sumCircMap[crow] = convRowSum[plane][crow];
+ 
+        convolveHorizontalRowPlane(outpamP, circMap, convKernelP, plane,
+                                   outputrow, sumCircMap);
     }
-    pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
+    pnm_writepamrow(outpamP, outputrow);
 
-    /* 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 );
+        /* For all subsequent rows */
 
-    /* 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];
+        unsigned int const addrow = crows - 1;
+        unsigned int const windowHeight = crows + 1;
 
-    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];
+        for (row = convKernelP->rows ; row < rows; ++row) {
+            unsigned int const toprow = (row + 2) % windowHeight;
+
+            unsigned int col;
+
+            readAndScaleRow(inpamP, rowbuf[row % windowHeight],
+                            outpamP->maxval, outpamP->depth);
+            
+            i = 0;
+            for (irow = toprow; irow < windowHeight; ++i, ++irow) {
+                circMap[i] = rowbuf[irow];
+                sumCircMap[i] = convRowSum[plane][irow];
             }
-        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 );
-        }
+            for (irow = 0; irow < toprow; ++irow, ++i) {
+                circMap[i] = rowbuf[irow];
+                sumCircMap[i] = convRowSum[plane][irow];
+            }
+            
+            for (col = 0; col < cols; ++col) {
+                if (col < ccolso2 || col >= cols - ccolso2)
+                    outputrow[col] = circMap[crowso2][col];
+                else if (col == ccolso2) {
+                    unsigned int const leftcol = col - ccolso2;
+
+                    float matrixSum;
+
+                    {
+                        unsigned int ccol;
+                        sumCircMap[addrow][leftcol] = 0L;
+                        for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                            sumCircMap[addrow][leftcol] += 
+                                circMap[addrow][leftcol + ccol][0];
+                    }
+                    {
+                        unsigned int crow;
+                        for (crow = 0, matrixSum = 0.0; crow < crows; ++crow)
+                            matrixSum += sumCircMap[crow][leftcol] *
+                                convKernelP->weight[crow][0][0];
+                    }
+                    outputrow[col][0] =
+                        MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+                } else {
+                    unsigned int const leftcol = col - ccolso2;
+                    unsigned int const subcol  = col - ccolso2 - 1;
+                    unsigned int const addcol  = col + ccolso2;  
+
+                    float matrixSum;
+                    unsigned int crow;
+
+                    sumCircMap[addrow][leftcol] = sumCircMap[addrow][subcol]
+                        - circMap[addrow][subcol][0]
+                        + circMap[addrow][addcol][0];
+
+                    for (crow = 0, matrixSum = 0.0; crow < crows; ++crow)
+                        matrixSum += sumCircMap[crow][leftcol] *
+                            convKernelP->weight[crow][0][0];
+
+                    outputrow[col][0] =
+                        MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+                }
+            }
+            pnm_writepamrow(outpamP, outputrow);
         }
-    pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
     }
+    writeUnconvolvedBottom(outpamP, convKernelP, circMap);
+
+    freeRowbuf(rowbuf, windowHeight);
+}
 
-    /* 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
+convolvePgmVertical(struct pam *        const inpamP,
+                    struct pam *        const outpamP,
+                    struct convKernel * const convKernelP) {
+
+    /* Uses column sums as in mean convolution, above */
+
+    unsigned int const ccolso2 = convKernelP->cols / 2;
+    float const weight = convKernelP->weight[0][0][0];
+    unsigned int const windowHeight = convKernelP->rows + 1;
+        /* The height of the window we keep in the row buffer.  The buffer
+           contains the rows covered by the convolution kernel, plus the row
+           immediately above that.  The latter is there because to compute
+           the sliding mean, we need to subtract off the row that the
+           convolution kernel just slid past.
+        */
+
+    tuple ** rowbuf;
+        /* Same as in convolvePgmGeneral */
+    tuple ** circMap;
+        /* Same as in convolvePgmGeneral */
+    tuple * outputrow;
+        /* Same as in convolvePgmGeneral */
+    unsigned int row;
+        /* Row number of next row to read in from the file */
+    sample * convColumnSum;  /* Malloc'd */
+        /* convColumnSum[col] is the sum of all the pixels in the Column 'col'
+           of the image within the current vertical convolution window.
+           I.e. if our convolution kernel is 5 rows high and we're now looking
+           at Rows 10-15, convColumn[3] is the sum of Column 3, Rows 10-15.
+        */
+    sample matrixSum;
+    unsigned int col;
+
+    rowbuf = allocRowbuf(pamP, windowHeight);
+    MALLOCARRAY_NOFAIL(circMap, windowHeight);
+    outputrow = pnm_allocpamrow(outpamP);
+
+    MALLOCARRAY_NOFAIL(convColumnSum, inpamP->width);
+
+    pnm_writepaminit(outpamP);
+
+    readInitialRowbuf(inpamP, convKernelP, rowbuf, outpamP->maxval);
+
+    writeUnconvolvedTop(outpamP, convKernelP, rowbuf);
+
+    setupCircMap(circMap, rowbuf, 0);
+
+    /* Convolve the first window the long way */
+    convolveAndComputeColumnSums(inpamP, outpamP, circMap, convKernelP,
+                                 outputrow, convColumnSum);
+
+    /* Write that first convolved row */
+    pnm_writepamrow(outpamP, outputrow);
+
+    for (row = convKernelP->rows; row < inpamP->height; ++row) {
+        unsigned int const subrow = 0;
+            /* Row just above convolution window -- what we subtract from
+               running sum
+            */
+        unsigned int const addrow = 1 + (convKernelP->rows - 1);
+            /* Bottom row of convolution window: What we add to running sum */
+
+        unsigned int col;
+
+        readAndScaleRow(inpamP, rowbuf[row % windowHeight],
+                        outpamP->maxval, outpamP->depth);
+
+        /* Remember the window is one row higher than the convolution
+           kernel.  The top row in the window is not part of this convolution.
+        */
+
+        setupCircMap(circMap, rowbuf, (row + 1) % windowHeight);
+
+        for (col = 0; col < cols; ++col) {
+            if (col < ccolso2 || col >= cols - ccolso2)
+                outputrow[col] = circMap[crowso2][col];
+            else if (col == ccolso2) {
+                unsigned int const leftcol = col - ccolso2;
+
+                float matrixSum;
+                unsigned int crow;
+                unsigned int ccol;
+
+                for (crow = 0; crow < convKernelP->rows; ++crow) {
+                    unsigned int ccol;
+
+                    for (ccol = 0; ccol < convKernelP->cols; ++ccol)
+                        convColumnSum[leftcol + ccol] +=
+                            circMap[crow][ccol][0];
+                }
+                for (ccol = 0, matrixSum = 0.0;
+                     ccol < convKernelP->cols;
+                     ++ccol) {
+                    matrixSum += convColumnSum[leftcol + ccol] *
+                        convKernelP->weight[0][ccol][0];
+                }
+                outputrow[col][0] =
+                    MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+            } else {
+                unsigned int const leftcol = col - ccolso2;
+                unsigned int const addcol  = col + ccolso2;
+
+                float matrixSum;
+                unsigned int ccol;
+
+                convColumnSum[addcol] = convColumnSum[addcol]
+                    - circMap[subrow][addcol][0]
+                    + circMap[addrow][addcol][0];
+
+                for (ccol = 0, matrixSum = 0.0;
+                     ccol < convKernelP->cols;
+                     ++ccol) {
+                    matrixSum += convColumnSum[leftcol + ccol] *
+                        convKernelP->weight[0][ccol][0];
+                }
+                outputrow[col][0] =
+                    MIN(outpamP->maxval, MAX(0, matrixSum + 0.5));
+            }
+        }
+        pnm_writepamrow(outpamP, outputrow);
     }
+    writeUnconvolvedBottom(outpamP, convKernelP, circMap);
 
+    freeRowbuf(rowbuf, windowHeight);
+}
 
 
 
 /* PPM General Convolution Algorithm
 **
 ** No redundancy in convolution matrix.  Just use brute force.
-** See pgm_general_convolve() for more details.
+** See convolvePgmGeneral() for more details.
 */
 
 static void
@@ -1787,36 +1946,35 @@ determineConvolveType(xel * const *         const cxels,
     
     /* Which type do we have? */
     if (horizontal && vertical) {
-        typeP->ppmConvolver = ppm_mean_convolve;
-        typeP->pgmConvolver = pgm_mean_convolve;
+        typeP->ppmConvolver = convolvePpmMean;
+        typeP->pgmConvolver = convolvePgmMean;
     } else if (horizontal) {
-        typeP->ppmConvolver = ppm_horizontal_convolve;
-        typeP->pgmConvolver = pgm_horizontal_convolve;
+        typeP->ppmConvolver = convolvePpmHorizontal;
+        typeP->pgmConvolver = convolvePgmHorizontal;
     } else if (vertical) {
-        typeP->ppmConvolver = ppm_vertical_convolve;
-        typeP->pgmConvolver = pgm_vertical_convolve;
+        typeP->ppmConvolver = convolvePpmVertical;
+        typeP->pgmConvolver = convolvePgmVertical;
     } else {
-        typeP->ppmConvolver = ppm_general_convolve;
-        typeP->pgmConvolver = pgm_general_convolve;
+        typeP->ppmConvolver = convolvePpmGeneral;
+        typeP->pgmConvolver = convolvePgmGeneral;
     }
 }
 
 
 
 static void
-convolveIt(int                 const format,
+convolveIt(struct pam *        const inpamP,
+           struct pam *        const outpamP,
            struct convolveType const convolveType,
-           const float**       const rweights,
-           const float**       const gweights,
-           const float**       const bweights) {
+           struct convKernel * const convKernelP) {
 
-    switch (PNM_FORMAT_TYPE(format)) {
+    switch (PNM_FORMAT_TYPE(inpamP->format)) {
     case PPM_TYPE:
-        convolveType.ppmConvolver(rweights, gweights, bweights);
+        convolveType.ppmConvolver(inpamP, outpamP, convKernelP);
         break;
 
     default:
-        convolveType.pgmConvolver(gweights);
+        convolveType.pgmConvolver(inpamP, outpamP, convKernelP);
     }
 }
 
@@ -1826,24 +1984,25 @@ int
 main(int argc, char * argv[]) {
 
     struct cmdlineInfo cmdline;
-    FILE* cifp;
-    xel** cxels;
+    FILE * ifP;
+    FILE * cifP;
+    tuple ** ctuples;
     int cformat;
     xelval cmaxval;
     struct convolveType convolveType;
-    float ** rweights;
-    float ** gweights;
-    float ** bweights;
+    struct convKernel * convKernelP;
+    struct pam inpam;
+    sturct pam cpam;
 
     pnm_init(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    cifp = pm_openr(cmdline.kernelFilespec);
+    cifP = pm_openr(cmdline.kernelFilespec);
 
     /* Read in the convolution matrix. */
-    cxels = pnm_readpnm(cifp, &ccols, &crows, &cmaxval, &cformat);
-    pm_close(cifp);
+    ctuples = pnm_readpam(cifP, &cpam, PAM_STRUCT_SIZE(tuple_type));
+    pm_close(cifP);
 
     if (ccols % 2 != 1 || crows % 2 != 1)
         pm_error("the convolution matrix must have an odd number of "
@@ -1852,43 +2011,47 @@ main(int argc, char * argv[]) {
     ccolso2 = ccols / 2;
     crowso2 = crows / 2;
 
-    ifp = pm_openr(cmdline.inputFilespec);
+    ifP = pm_openr(cmdline.inputFilespec);
 
-    pnm_readpnminit(ifp, &cols, &rows, &maxval, &format);
-    if (cols < ccols || rows < crows)
+    pnm_readpaminit(ifP, &inpam, PAM_STRUCT_SIZE(allocation_depth));
+    if (inpam.cols < cpam.cols || inpam.rows < cpam.rows)
         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)
+    outpam = inpam;
+
+    outpam.format = MAX(PNM_FORMAT_TYPE(cpam.format),
+                        PNM_FORMAT_TYPE(inpam.format));
+    if (PNM_FORMAT_TYPE(cpam.format) != outpam.format)
         pnm_promoteformat(cxels, ccols, crows, cmaxval, cformat, 
-                          cmaxval, newformat);
-    if (PNM_FORMAT_TYPE(format) != newformat) {
-        switch (PNM_FORMAT_TYPE(newformat)) {
+                          cmaxval, outpam.format);
+    if (PNM_FORMAT_TYPE(inpam.format) != outpam.format) {
+        switch (PNM_FORMAT_TYPE(outpam.format)) {
         case PPM_TYPE:
-            if (PNM_FORMAT_TYPE(format) != newformat)
+            if (PNM_FORMAT_TYPE(inpam.format) != outpam.format)
                 pm_message("promoting to PPM");
             break;
         case PGM_TYPE:
-            if (PNM_FORMAT_TYPE(format) != newformat)
+            if (PNM_FORMAT_TYPE(inpam.format) != outpam.format)
                 pm_message("promoting to PGM");
             break;
         }
     }
 
-    computeWeights(cxels, ccols, crows, newformat, cmaxval, !cmdline.nooffset,
-                   &rweights, &gweights, &bweights);
+    pnm_setminallocationdepth(&inpam, MAX(inpam.depth, outpam.depth));
+
+    computeKernel(cpam, ctuples, outpam.format, !cmdline.nooffset,
+                  &convKernelP);
 
     /* Handle certain special cases when runtime can be improved. */
 
-    determineConvolveType(cxels, &convolveType);
+    determineConvolveType(ctuples, &convolveType);
 
-    convolveIt(format, convolveType, 
-               (const float **)rweights, 
-               (const float **)gweights, 
-               (const float **)bweights);
+    convolveIt(&inpam, &outpam, convolveType, convKernelP);
 
+    convKernelDestroy(convKernelP);
     pm_close(stdout);
-    pm_close(ifp);
+    pm_close(ifP);
+
     return 0;
 }