/* pnmconvol.c - general MxN convolution on a portable anymap
**
** 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"

/* 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;


static int check_convolve_type ARGS((xel **cxels));
static void pgm_general_convolve ARGS((void));
static void ppm_general_convolve ARGS((void));
static void pgm_mean_convolve ARGS((void));
static void ppm_mean_convolve ARGS((void));
static void pgm_vertical_convolve ARGS((void));
static void ppm_vertical_convolve ARGS((void));
static void pgm_horizontal_convolve ARGS((void));
static void ppm_horizontal_convolve ARGS((void));

#define TRUE    1
#define FALSE   0

#define GENERAL_CONVOLVE    0
#define MEAN_CONVOLVE       1
#define HORIZONTAL_CONVOLVE 2
#define VERTICAL_CONVOLVE   3

static FILE* ifp;
static float** rweights;
static float** gweights;
static float** bweights;
static int crows, ccols, ccolso2, crowso2;
static int cols, rows;
static xelval maxval;
static int format, newformat;
static float rmeanweight, gmeanweight, bmeanweight;



static void
computeWeights(xel ** const cxels, int const ccols, int const crows,
               int const cformat, xelval const cmaxval,
               float *** const rweightsP,
               float *** const gweightsP,
               float *** const bweightsP ) {

    float** rweights;
    float** gweights;
    float** bweights;

    float rsum, gsum, bsum;

    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;

    for ( crow = 0; crow < crows; ++crow ) {
        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]) * 2.0 / cmaxval - 1.0 );
                gsum += gweights[crow][ccol] =
                    ( PPM_GETG(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
                bsum += bweights[crow][ccol] =
                    ( PPM_GETB(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
                break;
                
            default:
                gsum += gweights[crow][ccol] =
                    ( PNM_GET1(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
                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);
        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;
    }
}


int
main(int argc, char * argv[] )
    {
    FILE* cifp;
    xel** cxels;
    int argn, cformat;
    xelval cmaxval;
    char* usage = "<convolutionfile> [pnmfile]";
    int convolve_type;

    pnm_init( &argc, argv );

    argn = 1;

    if ( argn == argc )
    pm_usage( usage );
    cifp = pm_openr( argv[argn] );
    ++argn;

    if ( argn != argc )
    {
    ifp = pm_openr( argv[argn] );
    ++argn;
    }
    else
    ifp = stdin;

    if ( argn != argc )
    pm_usage( usage );

    /* 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;

    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,
                   &rweights, &gweights, &bweights);
    /* For mean_convolve routines.  All weights of a single color are the same
    ** so just grab any one of them.
    */
    rmeanweight = rweights[0][0];
    gmeanweight = gweights[0][0];
    bmeanweight = bweights[0][0];

    /* Handle certain special cases when runtime can be improved. */

    convolve_type = check_convolve_type(cxels);

    switch ( PNM_FORMAT_TYPE(format) )
    {
    case PPM_TYPE:
    switch (convolve_type)
        {
        case MEAN_CONVOLVE:
        ppm_mean_convolve();
        break;

        case HORIZONTAL_CONVOLVE:
        ppm_horizontal_convolve();
        break;

        case VERTICAL_CONVOLVE:
        ppm_vertical_convolve();
        break;

        default: /* GENERAL_CONVOLVE */
        ppm_general_convolve();
        break;  
        }
    break;

    default:
    switch (convolve_type)
        {
        case MEAN_CONVOLVE:
        pgm_mean_convolve();
        break;

        case HORIZONTAL_CONVOLVE:
        pgm_horizontal_convolve();
        break;

        case VERTICAL_CONVOLVE:
        pgm_vertical_convolve();
        break;

        default: /* GENERAL_CONVOLVE */
        pgm_general_convolve();
        break;  
        }
    break;
    } /* PNM_TYPE */

    pm_close( ifp );
    exit( 0 );
    }



/* check_convolve_type
**
** Determine if the convolution matrix is one of the special cases that
** can be processed faster than the general form.
**
** Does not check for the case where one of the PPM colors can have differing 
** types.  Only handles cases where all PPM's are of the same special case.
*/

static int
check_convolve_type (cxels)
    xel **cxels;
    {
    int convolve_type;
    int horizontal, vertical;
    int tempcxel, rtempcxel, gtempcxel, btempcxel;
    int crow, ccol;

    switch ( PNM_FORMAT_TYPE(format) )
    {
    case PPM_TYPE:
    horizontal = TRUE;
    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;
    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;
    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;
    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 )
    convolve_type = MEAN_CONVOLVE;
    else if ( horizontal )
    convolve_type = HORIZONTAL_CONVOLVE;
    else if ( vertical )
    convolve_type = VERTICAL_CONVOLVE;
    else
    convolve_type = GENERAL_CONVOLVE;

    return (convolve_type);
    }




/* General PGM Convolution
**
** No useful redundancy in convolution matrix.
*/

static void
pgm_general_convolve()
    {
    register int ccol, col;
    xel** xelbuf;
    xel* outputrow;
    xelval g;
    int row, crow;
    float gsum;
    xel **rowptr, *temprptr;
    int toprow, temprow;
    int i, irow;
    int leftcol;
    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 )
    {
    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;
        gsum = 0.0;
        for ( crow = 0; crow < crows; ++crow )
            {
            temprptr = rowptr[crow] + leftcol;
            for ( ccol = 0; ccol < ccols; ++ccol )
            gsum += PNM_GET1( *(temprptr + ccol) )
                * gweights[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 );

    pm_close( stdout );
    }


/* 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()
    {
    register 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 );

    pm_close( stdout );
    }


/* 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()
    {
    register 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] * gweights[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] * gweights[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] * gweights[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] * gweights[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 );

    pm_close( stdout );
    }


/* PGM Vertical Convolution
**
** Uses column sums as in Mean Convolution.
**
*/


static void
pgm_vertical_convolve()
    {
    register 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] * gweights[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] * gweights[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] * gweights[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] * gweights[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 );

    pm_close( stdout );
    }




/* 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()
    {
    register 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 );

    pm_close( stdout );
    }


/* PPM Mean Convolution
**
** Same as pgm_mean_convolve() but for PPM.
**
*/

static void
ppm_mean_convolve()
    {
    register 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 );

    pm_close( stdout );
    }


/* PPM Horizontal Convolution
**
** Same as pgm_horizontal_convolve()
**
**/

static void
ppm_horizontal_convolve()
    {
    register 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 );

    pm_close( stdout );
    }


/* PPM Vertical Convolution
**
** Same as pgm_vertical_convolve()
**
*/

static void
ppm_vertical_convolve()
    {
    register 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 tempcol;
    float rsum, gsum, bsum;
    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 */
    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 );
    }

    /* 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 )
        {
        rsum = 0.0;
        gsum = 0.0;
        bsum = 0.0;
        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)
        {
        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
        {
        rsum = 0.0;
        gsum = 0.0;
        bsum = 0.0;
        leftcol = col - ccolso2;
        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 )
        {
        rsum = 0.0;
        gsum = 0.0;
        bsum = 0.0;
        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] );
            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
        {
        rsum = 0.0;
        gsum = 0.0;
        bsum = 0.0;
        leftcol = col - ccolso2;
        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 );

    pm_close( stdout );
    }

/******************************************************************************
                            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
 ----------------------------------------------------------------------------

*/








