home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / NETWORK / netpbm_src.lzh / NETPBM / PNM / pnmconvol.c < prev    next >
C/C++ Source or Header  |  1996-11-18  |  55KB  |  1,905 lines
  1. /* pnmconvol.c - general MxN convolution on a portable anymap
  2. **
  3. ** Version 2.0.1 January 30, 1995
  4. **
  5. ** Major rewriting by Mike Burns
  6. ** Copyright (C) 1994, 1995 by Mike Burns (burns@chem.psu.edu)
  7. **
  8. ** Copyright (C) 1989, 1991 by Jef Poskanzer.
  9. **
  10. ** Permission to use, copy, modify, and distribute this software and its
  11. ** documentation for any purpose and without fee is hereby granted, provided
  12. ** that the above copyright notice appear in all copies and that both that
  13. ** copyright notice and this permission notice appear in supporting
  14. ** documentation.  This software is provided "as is" without express or
  15. ** implied warranty.
  16. */
  17.  
  18. /* Version 2.0.1 Changes
  19. ** ---------------------
  20. ** Fixed four lines that were improperly allocated as sizeof( float ) when they
  21. ** should have been sizeof( long ).
  22. **
  23. ** Version 2.0 Changes
  24. ** -------------------
  25. ** Reduce run time by general optimizations and handling special cases of
  26. ** convolution matrices.  Program automatically determines if convolution 
  27. ** matrix is one of the types it can make use of so no extra command line
  28. ** arguments are necessary.
  29. **
  30. ** Examples of convolution matrices for the special cases are
  31. **
  32. **    Mean       Horizontal    Vertical
  33. **    x x x        x x x        x y z
  34. **    x x x        y y y        x y z
  35. **    x x x        z z z        x y z
  36. **
  37. ** I don't know if the horizontal and vertical ones are of much use, but
  38. ** after working on the mean convolution, it gave me ideas for the other two.
  39. **
  40. ** Some other compiler dependent optimizations
  41. ** -------------------------------------------
  42. ** Created separate functions as code was getting too large to put keep both
  43. ** PGM and PPM cases in same function and also because SWITCH statement in
  44. ** inner loop can take progressively more time the larger the size of the 
  45. ** convolution matrix.  GCC is affected this way.
  46. **
  47. ** Removed use of MOD (%) operator from innermost loop by modifying manner in
  48. ** which the current xelbuf[] is chosen.
  49. ** 
  50. */
  51.  
  52. #include "pnm.h"
  53.  
  54. /* Macros to verify that r,g,b values are within proper range */
  55.  
  56. #define CHECK_GRAY \
  57.     if ( tempgsum < 0L ) g = 0; \
  58.     else if ( tempgsum > maxval ) g = maxval; \
  59.     else g = tempgsum;
  60.  
  61. #define CHECK_RED \
  62.     if ( temprsum < 0L ) r = 0; \
  63.     else if ( temprsum > maxval ) r = maxval; \
  64.     else r = temprsum;
  65.  
  66. #define CHECK_GREEN \
  67.     if ( tempgsum < 0L ) g = 0; \
  68.     else if ( tempgsum > maxval ) g = maxval; \
  69.     else g = tempgsum;
  70.  
  71. #define CHECK_BLUE \
  72.     if ( tempbsum < 0L ) b = 0; \
  73.     else if ( tempbsum > maxval ) b = maxval; \
  74.     else b = tempbsum;
  75.  
  76.  
  77. static int check_convolve_type ARGS((xel **cxels));
  78. static void pgm_general_convolve ARGS((void));
  79. static void ppm_general_convolve ARGS((void));
  80. static void pgm_mean_convolve ARGS((void));
  81. static void ppm_mean_convolve ARGS((void));
  82. static void pgm_vertical_convolve ARGS((void));
  83. static void ppm_vertical_convolve ARGS((void));
  84. static void pgm_horizontal_convolve ARGS((void));
  85. static void ppm_horizontal_convolve ARGS((void));
  86.  
  87. #define TRUE    1
  88. #define FALSE    0
  89.  
  90. #define GENERAL_CONVOLVE    0
  91. #define MEAN_CONVOLVE        1
  92. #define HORIZONTAL_CONVOLVE    2
  93. #define VERTICAL_CONVOLVE    3
  94.  
  95. static FILE* ifp;
  96. static float** rweights;
  97. static float** gweights;
  98. static float** bweights;
  99. static int crows, ccols, ccolso2, crowso2;
  100. static int cols, rows;
  101. static xelval maxval;
  102. static int format, newformat;
  103. static float rmeanweight, gmeanweight, bmeanweight;
  104.  
  105. int
  106. main( argc, argv )
  107.     int argc;
  108.     char* argv[];
  109.     {
  110.     FILE* cifp;
  111.     xel** cxels;
  112.     int argn, cformat;
  113.     int crow, row;
  114.     register int ccol, col;
  115.     xelval cmaxval;
  116.     xelval g;
  117.     float gsum;
  118.     float rsum, bsum;
  119.     char* usage = "<convolutionfile> [pnmfile]";
  120.     int convolve_type;
  121.  
  122.     pnm_init( &argc, argv );
  123.  
  124.     argn = 1;
  125.  
  126.     if ( argn == argc )
  127.     pm_usage( usage );
  128.     cifp = pm_openr( argv[argn] );
  129.     ++argn;
  130.  
  131.     if ( argn != argc )
  132.     {
  133.     ifp = pm_openr( argv[argn] );
  134.     ++argn;
  135.     }
  136.     else
  137.     ifp = stdin;
  138.  
  139.     if ( argn != argc )
  140.     pm_usage( usage );
  141.  
  142.     pnm_pbmmaxval = PNM_MAXMAXVAL;  /* use larger value for better results */
  143.  
  144.     /* Read in the convolution matrix. */
  145.     cxels = pnm_readpnm( cifp, &ccols, &crows, &cmaxval, &cformat );
  146.     pm_close( cifp );
  147.     if ( ccols % 2 != 1 || crows % 2 != 1 )
  148.     pm_error(
  149.      "the convolution matrix must have an odd number of rows and columns" );
  150.     ccolso2 = ccols / 2;
  151.     crowso2 = crows / 2;
  152.  
  153.     pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
  154.     if ( cols < ccols || rows < crows )
  155.     pm_error(
  156.         "the image is smaller than the convolution matrix" );
  157.  
  158.     newformat = max( PNM_FORMAT_TYPE(cformat), PNM_FORMAT_TYPE(format) );
  159.     if ( PNM_FORMAT_TYPE(cformat) != newformat )
  160.     pnm_promoteformat( cxels, ccols, crows, cmaxval, cformat, cmaxval, newformat );
  161.     if ( PNM_FORMAT_TYPE(format) != newformat )
  162.         {
  163.         switch ( PNM_FORMAT_TYPE(newformat) )
  164.             {
  165.             case PPM_TYPE:
  166.             if ( PNM_FORMAT_TYPE(format) != newformat )
  167.                 pm_message( "promoting to PPM" );
  168.             break;
  169.             case PGM_TYPE:
  170.             if ( PNM_FORMAT_TYPE(format) != newformat )
  171.                 pm_message( "promoting to PGM" );
  172.             break;
  173.             }
  174.         }
  175.  
  176.  
  177.     /* Set up the normalized weights. */
  178.     rweights = (float**) pm_allocarray( ccols, crows, sizeof(float) );
  179.     gweights = (float**) pm_allocarray( ccols, crows, sizeof(float) );
  180.     bweights = (float**) pm_allocarray( ccols, crows, sizeof(float) );
  181.     rsum = gsum = bsum = 0;
  182.     for ( crow = 0; crow < crows; ++crow )
  183.     for ( ccol = 0; ccol < ccols; ++ccol )
  184.         {
  185.         switch ( PNM_FORMAT_TYPE(format) )
  186.         {
  187.         case PPM_TYPE:
  188.         rsum += rweights[crow][ccol] =
  189.             ( PPM_GETR(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
  190.         gsum += gweights[crow][ccol] =
  191.             ( PPM_GETG(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
  192.         bsum += bweights[crow][ccol] =
  193.             ( PPM_GETB(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
  194.         break;
  195.  
  196.         default:
  197.         gsum += gweights[crow][ccol] =
  198.             ( PNM_GET1(cxels[crow][ccol]) * 2.0 / cmaxval - 1.0 );
  199.         break;
  200.         }
  201.         }
  202.  
  203.     /* For mean_convolve routines.  All weights of a single color are the same
  204.     ** so just grab any one of them.
  205.     */
  206.     rmeanweight = rweights[0][0];
  207.     gmeanweight = gweights[0][0];
  208.     bmeanweight = bweights[0][0];
  209.  
  210.     switch ( PNM_FORMAT_TYPE(format) )
  211.     {
  212.     case PPM_TYPE:
  213.     if ( rsum < 0.9 || rsum > 1.1 || gsum < 0.9 || gsum > 1.1 ||
  214.          bsum < 0.9 || bsum > 1.1 )
  215.         pm_message(
  216.         "WARNING - this convolution matrix is biased" );
  217.     break;
  218.  
  219.     default:
  220.     if ( gsum < 0.9 || gsum > 1.1 )
  221.         pm_message(
  222.          "WARNING - this convolution matrix is biased" );
  223.     break;
  224.     }
  225.  
  226.     /* Handle certain special cases when runtime can be improved. */
  227.  
  228.     convolve_type = check_convolve_type(cxels);
  229.  
  230.     switch ( PNM_FORMAT_TYPE(format) )
  231.     {
  232.     case PPM_TYPE:
  233.     switch (convolve_type)
  234.         {
  235.         case MEAN_CONVOLVE:
  236.         ppm_mean_convolve();
  237.         break;
  238.  
  239.         case HORIZONTAL_CONVOLVE:
  240.         ppm_horizontal_convolve();
  241.         break;
  242.  
  243.         case VERTICAL_CONVOLVE:
  244.         ppm_vertical_convolve();
  245.         break;
  246.  
  247.         default: /* GENERAL_CONVOLVE */
  248.         ppm_general_convolve();
  249.         break;  
  250.         }
  251.     break;
  252.  
  253.     default:
  254.     switch (convolve_type)
  255.         {
  256.         case MEAN_CONVOLVE:
  257.         pgm_mean_convolve();
  258.         break;
  259.  
  260.         case HORIZONTAL_CONVOLVE:
  261.         pgm_horizontal_convolve();
  262.         break;
  263.  
  264.         case VERTICAL_CONVOLVE:
  265.         pgm_vertical_convolve();
  266.         break;
  267.  
  268.         default: /* GENERAL_CONVOLVE */
  269.         pgm_general_convolve();
  270.         break;  
  271.         }
  272.     break;
  273.     } /* PNM_TYPE */
  274.  
  275.     pm_close( ifp );
  276.     exit( 0 );
  277.     }
  278.  
  279.  
  280.  
  281. /* check_convolve_type
  282. **
  283. ** Determine if the convolution matrix is one of the special cases that
  284. ** can be processed faster than the general form.
  285. **
  286. ** Does not check for the case where one of the PPM colors can have differing 
  287. ** types.  Only handles cases where all PPM's are of the same special case.
  288. */
  289.  
  290. static int
  291. check_convolve_type (cxels)
  292.     xel **cxels;
  293.     {
  294.     int convolve_type;
  295.     int horizontal, vertical;
  296.     int tempcxel, rtempcxel, gtempcxel, btempcxel;
  297.     int crow, ccol;
  298.  
  299.     switch ( PNM_FORMAT_TYPE(format) )
  300.     {
  301.     case PPM_TYPE:
  302.     horizontal = TRUE;
  303.     crow = 0;
  304.     while ( horizontal && (crow < crows) )
  305.         {
  306.         ccol = 1;
  307.         rtempcxel = PPM_GETR(cxels[crow][0]);
  308.         gtempcxel = PPM_GETG(cxels[crow][0]);
  309.         btempcxel = PPM_GETB(cxels[crow][0]);
  310.         while ( horizontal && (ccol < ccols) )
  311.         {
  312.         if (( PPM_GETR(cxels[crow][ccol]) != rtempcxel ) |
  313.             ( PPM_GETG(cxels[crow][ccol]) != gtempcxel ) |
  314.             ( PPM_GETB(cxels[crow][ccol]) != btempcxel )) 
  315.             horizontal = FALSE;
  316.         ++ccol;
  317.         }
  318.         ++crow;
  319.         }
  320.  
  321.     vertical = TRUE;
  322.     ccol = 0;
  323.     while ( vertical && (ccol < ccols) )
  324.         {
  325.         crow = 1;
  326.         rtempcxel = PPM_GETR(cxels[0][ccol]);
  327.         gtempcxel = PPM_GETG(cxels[0][ccol]);
  328.         btempcxel = PPM_GETB(cxels[0][ccol]);
  329.         while ( vertical && (crow < crows) )
  330.         {
  331.         if (( PPM_GETR(cxels[crow][ccol]) != rtempcxel ) |
  332.             ( PPM_GETG(cxels[crow][ccol]) != gtempcxel ) |
  333.             ( PPM_GETB(cxels[crow][ccol]) != btempcxel ))
  334.             vertical = FALSE;
  335.         ++crow;
  336.         }
  337.         ++ccol;
  338.         }
  339.     break;
  340.  
  341.     default:
  342.     horizontal = TRUE;
  343.     crow = 0;
  344.     while ( horizontal && (crow < crows) )
  345.         {
  346.         ccol = 1;
  347.         tempcxel = PNM_GET1(cxels[crow][0]);
  348.         while ( horizontal && (ccol < ccols) )
  349.         {
  350.         if ( PNM_GET1(cxels[crow][ccol]) != tempcxel )
  351.             horizontal = FALSE;
  352.         ++ccol;
  353.         }
  354.         ++crow;
  355.         }
  356.  
  357.     vertical = TRUE;
  358.     ccol = 0;
  359.     while ( vertical && (ccol < ccols) )
  360.         {
  361.         crow = 1;
  362.         tempcxel = PNM_GET1(cxels[0][ccol]);
  363.         while ( vertical && (crow < crows) )
  364.         {
  365.         if ( PNM_GET1(cxels[crow][ccol]) != tempcxel )
  366.             vertical = FALSE;
  367.         ++crow;
  368.         }
  369.         ++ccol;
  370.         }
  371.     break;
  372.     }
  373.  
  374.     /* Which type do we have? */
  375.     if ( horizontal && vertical )
  376.     convolve_type = MEAN_CONVOLVE;
  377.     else if ( horizontal )
  378.     convolve_type = HORIZONTAL_CONVOLVE;
  379.     else if ( vertical )
  380.     convolve_type = VERTICAL_CONVOLVE;
  381.     else
  382.     convolve_type = GENERAL_CONVOLVE;
  383.  
  384.     return (convolve_type);
  385.     }
  386.  
  387.  
  388.  
  389.  
  390. /* General PGM Convolution
  391. **
  392. ** No useful redundancy in convolution matrix.
  393. */
  394.  
  395. static void
  396. pgm_general_convolve()
  397.     {
  398.     register int ccol, col;
  399.     xel** xelbuf;
  400.     xel* outputrow;
  401.     xel x, y;
  402.     xelval g;
  403.     int row, crow;
  404.     float gsum;
  405.     xel **rowptr, *temprptr;
  406.     int toprow, temprow;
  407.     int i, irow;
  408.     int leftcol;
  409.     long tempgsum;
  410.  
  411.     /* Allocate space for one convolution-matrix's worth of rows, plus
  412.     ** a row output buffer. */
  413.     xelbuf = pnm_allocarray( cols, crows );
  414.     outputrow = pnm_allocrow( cols );
  415.  
  416.     /* Allocate array of pointers to xelbuf */
  417.     rowptr = (xel **) pnm_allocarray( 1, crows );
  418.  
  419.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  420.  
  421.     /* Read in one convolution-matrix's worth of image, less one row. */
  422.     for ( row = 0; row < crows - 1; ++row )
  423.     {
  424.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  425.     if ( PNM_FORMAT_TYPE(format) != newformat )
  426.         pnm_promoteformatrow(
  427.         xelbuf[row], cols, maxval, format, maxval, newformat );
  428.     /* Write out just the part we're not going to convolve. */
  429.     if ( row < crowso2 )
  430.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  431.     }
  432.  
  433.     /* Now the rest of the image - read in the row at the end of
  434.     ** xelbuf, and convolve and write out the row in the middle.
  435.     */
  436.     for ( ; row < rows; ++row )
  437.     {
  438.     toprow = row + 1;
  439.     temprow = row % crows;
  440.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  441.     if ( PNM_FORMAT_TYPE(format) != newformat )
  442.         pnm_promoteformatrow(
  443.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  444.  
  445.     /* Arrange rowptr to eliminate the use of mod function to determine
  446.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  447.     */
  448.     temprow = toprow % crows;
  449.     i = 0;
  450.     for (irow = temprow; irow < crows; ++i, ++irow)
  451.         rowptr[i] = xelbuf[irow];
  452.     for (irow = 0; irow < temprow; ++irow, ++i)
  453.         rowptr[i] = xelbuf[irow];
  454.  
  455.     for ( col = 0; col < cols; ++col )
  456.         {
  457.         if ( col < ccolso2 || col >= cols - ccolso2 )
  458.         outputrow[col] = rowptr[crowso2][col];
  459.         else
  460.         {
  461.         leftcol = col - ccolso2;
  462.         gsum = 0.0;
  463.         for ( crow = 0; crow < crows; ++crow )
  464.             {
  465.             temprptr = rowptr[crow] + leftcol;
  466.             for ( ccol = 0; ccol < ccols; ++ccol )
  467.             gsum += PNM_GET1( *(temprptr + ccol) )
  468.                 * gweights[crow][ccol];
  469.             }
  470.             tempgsum = gsum + 0.5;
  471.             CHECK_GRAY;
  472.             PNM_ASSIGN1( outputrow[col], g );
  473.         }
  474.         }
  475.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  476.     }
  477.  
  478.     /* Now write out the remaining unconvolved rows in xelbuf. */
  479.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  480.     pnm_writepnmrow(
  481.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  482.  
  483.     pm_close( stdout );
  484.     }
  485.  
  486.  
  487. /* PGM Mean Convolution
  488. **
  489. ** This is the common case where you just want the target pixel replaced with
  490. ** the average value of its neighbors.  This can work much faster than the
  491. ** general case because you can reduce the number of floating point operations
  492. ** that are required since all the weights are the same.  You will only need
  493. ** to multiply by the weight once, not for every pixel in the convolution
  494. ** matrix.
  495. **
  496. ** This algorithm works by creating sums for each column of crows height for
  497. ** the whole width of the image.  Then add ccols column sums together to obtain
  498. ** the total sum of the neighbors and multiply that sum by the weight.  As you
  499. ** move right to left to calculate the next pixel, take the total sum you just
  500. ** generated, add in the value of the next column and subtract the value of the
  501. ** leftmost column.  Multiply that by the weight and that's it.  As you move
  502. ** down a row, calculate new column sums by using previous sum for that column
  503. ** and adding in pixel on current row and subtracting pixel in top row.
  504. **
  505. */
  506.  
  507.  
  508. static void
  509. pgm_mean_convolve()
  510.     {
  511.     register int ccol, col;
  512.     xel** xelbuf;
  513.     xel* outputrow;
  514.     xelval g;
  515.     int row, crow;
  516.     xel **rowptr, *temprptr;
  517.     int leftcol;
  518.     int i, irow;
  519.     int toprow, temprow;
  520.     int subrow, addrow;
  521.     int subcol, addcol;
  522.     long gisum;
  523.     int tempcol, crowsp1;
  524.     long tempgsum;
  525.     long *gcolumnsum;
  526.  
  527.     /* Allocate space for one convolution-matrix's worth of rows, plus
  528.     ** a row output buffer.  MEAN uses an extra row. */
  529.     xelbuf = pnm_allocarray( cols, crows + 1 );
  530.     outputrow = pnm_allocrow( cols );
  531.  
  532.     /* Allocate array of pointers to xelbuf. MEAN uses an extra row. */
  533.     rowptr = (xel **) pnm_allocarray( 1, crows + 1);
  534.  
  535.     /* Allocate space for intermediate column sums */
  536.     gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  537.     for ( col = 0; col < cols; ++col )
  538.     gcolumnsum[col] = 0L;
  539.  
  540.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  541.  
  542.     /* Read in one convolution-matrix's worth of image, less one row. */
  543.     for ( row = 0; row < crows - 1; ++row )
  544.     {
  545.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  546.     if ( PNM_FORMAT_TYPE(format) != newformat )
  547.         pnm_promoteformatrow(
  548.         xelbuf[row], cols, maxval, format, maxval, newformat );
  549.     /* Write out just the part we're not going to convolve. */
  550.     if ( row < crowso2 )
  551.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  552.     }
  553.  
  554.     /* Do first real row only */
  555.     subrow = crows;
  556.     addrow = crows - 1;
  557.     toprow = row + 1;
  558.     temprow = row % crows;
  559.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  560.     if ( PNM_FORMAT_TYPE(format) != newformat )
  561.     pnm_promoteformatrow(
  562.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  563.  
  564.     temprow = toprow % crows;
  565.     i = 0;
  566.     for (irow = temprow; irow < crows; ++i, ++irow)
  567.     rowptr[i] = xelbuf[irow];
  568.     for (irow = 0; irow < temprow; ++irow, ++i)
  569.     rowptr[i] = xelbuf[irow];
  570.  
  571.     gisum = 0L;
  572.     for ( col = 0; col < cols; ++col )
  573.     {
  574.     if ( col < ccolso2 || col >= cols - ccolso2 )
  575.         outputrow[col] = rowptr[crowso2][col];
  576.     else if ( col == ccolso2 )
  577.         {
  578.         leftcol = col - ccolso2;
  579.         for ( crow = 0; crow < crows; ++crow )
  580.         {
  581.         temprptr = rowptr[crow] + leftcol;
  582.         for ( ccol = 0; ccol < ccols; ++ccol )
  583.             gcolumnsum[leftcol + ccol] += 
  584.             PNM_GET1( *(temprptr + ccol) );
  585.         }
  586.         for ( ccol = 0; ccol < ccols; ++ccol)
  587.         gisum += gcolumnsum[leftcol + ccol];
  588.         tempgsum = (float) gisum * gmeanweight + 0.5;
  589.         CHECK_GRAY;
  590.         PNM_ASSIGN1( outputrow[col], g );
  591.         }
  592.     else
  593.         {
  594.         /* Column numbers to subtract or add to isum */
  595.         subcol = col - ccolso2 - 1;
  596.         addcol = col + ccolso2;  
  597.         for ( crow = 0; crow < crows; ++crow )
  598.         gcolumnsum[addcol] += PNM_GET1( rowptr[crow][addcol] );
  599.         gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol];
  600.         tempgsum = (float) gisum * gmeanweight + 0.5;
  601.         CHECK_GRAY;
  602.         PNM_ASSIGN1( outputrow[col], g );
  603.         }
  604.     }
  605.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  606.  
  607.     ++row;
  608.     /* For all subsequent rows do it this way as the columnsums have been
  609.     ** generated.  Now we can use them to reduce further calculations.
  610.     */
  611.     crowsp1 = crows + 1;
  612.     for ( ; row < rows; ++row )
  613.     {
  614.     toprow = row + 1;
  615.     temprow = row % (crows + 1);
  616.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  617.     if ( PNM_FORMAT_TYPE(format) != newformat )
  618.         pnm_promoteformatrow(
  619.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  620.  
  621.     /* This rearrangement using crows+1 rowptrs and xelbufs will cause
  622.     ** rowptr[0..crows-1] to always hold active xelbufs and for 
  623.     ** rowptr[crows] to always hold the oldest (top most) xelbuf.
  624.     */
  625.     temprow = (toprow + 1) % crowsp1;
  626.     i = 0;
  627.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  628.         rowptr[i] = xelbuf[irow];
  629.     for (irow = 0; irow < temprow; ++irow, ++i)
  630.         rowptr[i] = xelbuf[irow];
  631.  
  632.     gisum = 0L;
  633.     for ( col = 0; col < cols; ++col )
  634.         {
  635.         if ( col < ccolso2 || col >= cols - ccolso2 )
  636.         outputrow[col] = rowptr[crowso2][col];
  637.         else if ( col == ccolso2 )
  638.         {
  639.         leftcol = col - ccolso2;
  640.         for ( ccol = 0; ccol < ccols; ++ccol )
  641.             {
  642.             tempcol = leftcol + ccol;
  643.             gcolumnsum[tempcol] = gcolumnsum[tempcol]
  644.             - PNM_GET1( rowptr[subrow][ccol] )
  645.             + PNM_GET1( rowptr[addrow][ccol] );
  646.             gisum += gcolumnsum[tempcol];
  647.             }
  648.         tempgsum = (float) gisum * gmeanweight + 0.5;
  649.         CHECK_GRAY;
  650.         PNM_ASSIGN1( outputrow[col], g );
  651.         }
  652.         else
  653.         {
  654.         /* Column numbers to subtract or add to isum */
  655.         subcol = col - ccolso2 - 1;
  656.         addcol = col + ccolso2;  
  657.         gcolumnsum[addcol] = gcolumnsum[addcol]
  658.             - PNM_GET1( rowptr[subrow][addcol] )
  659.             + PNM_GET1( rowptr[addrow][addcol] );
  660.         gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol];
  661.         tempgsum = (float) gisum * gmeanweight + 0.5;
  662.         CHECK_GRAY;
  663.         PNM_ASSIGN1( outputrow[col], g );
  664.         }
  665.         }
  666.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  667.     }
  668.  
  669.     /* Now write out the remaining unconvolved rows in xelbuf. */
  670.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  671.     pnm_writepnmrow(
  672.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  673.  
  674.     pm_close( stdout );
  675.     }
  676.  
  677.  
  678. /* PGM Horizontal Convolution
  679. **
  680. ** Similar idea to using columnsums of the Mean and Vertical convolution,
  681. ** but uses temporary sums of row values.  Need to multiply by weights crows
  682. ** number of times.  Each time a new line is started, must recalculate the
  683. ** initials rowsums for the newest row only.  Uses queue to still access
  684. ** previous row sums.
  685. **
  686. */
  687.  
  688. static void
  689. pgm_horizontal_convolve()
  690.     {
  691.     register int ccol, col;
  692.     xel** xelbuf;
  693.     xel* outputrow;
  694.     xel x;
  695.     xelval g;
  696.     int row, crow;
  697.     xel **rowptr, *temprptr;
  698.     int leftcol;
  699.     int i, irow;
  700.     int temprow;
  701.     int subcol, addcol;
  702.     float gsum;
  703.     int addrow, subrow;
  704.     long **growsum, **growsumptr;
  705.     int crowsp1;
  706.     long tempgsum;
  707.  
  708.     /* Allocate space for one convolution-matrix's worth of rows, plus
  709.     ** a row output buffer. */
  710.     xelbuf = pnm_allocarray( cols, crows + 1 );
  711.     outputrow = pnm_allocrow( cols );
  712.  
  713.     /* Allocate array of pointers to xelbuf */
  714.     rowptr = (xel **) pnm_allocarray( 1, crows + 1);
  715.  
  716.     /* Allocate intermediate row sums.  HORIZONTAL uses an extra row. */
  717.     /* crows current rows and 1 extra for newest added row.           */
  718.     growsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) );
  719.     growsumptr = (long **) pnm_allocarray( 1, crows + 1);
  720.  
  721.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  722.  
  723.     /* Read in one convolution-matrix's worth of image, less one row. */
  724.     for ( row = 0; row < crows - 1; ++row )
  725.     {
  726.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  727.     if ( PNM_FORMAT_TYPE(format) != newformat )
  728.         pnm_promoteformatrow(
  729.         xelbuf[row], cols, maxval, format, maxval, newformat );
  730.     /* Write out just the part we're not going to convolve. */
  731.     if ( row < crowso2 )
  732.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  733.     }
  734.  
  735.     /* First row only */
  736.     temprow = row % crows;
  737.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  738.     if ( PNM_FORMAT_TYPE(format) != newformat )
  739.     pnm_promoteformatrow(
  740.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  741.  
  742.     temprow = (row + 1) % crows;
  743.     i = 0;
  744.     for (irow = temprow; irow < crows; ++i, ++irow)
  745.     rowptr[i] = xelbuf[irow];
  746.     for (irow = 0; irow < temprow; ++irow, ++i)
  747.     rowptr[i] = xelbuf[irow];
  748.  
  749.     for ( crow = 0; crow < crows; ++crow )
  750.     growsumptr[crow] = growsum[crow];
  751.  
  752.     for ( col = 0; col < cols; ++col )
  753.     {
  754.     if ( col < ccolso2 || col >= cols - ccolso2 )
  755.         outputrow[col] = rowptr[crowso2][col];
  756.     else if ( col == ccolso2 )
  757.         {
  758.         leftcol = col - ccolso2;
  759.         gsum = 0.0;
  760.         for ( crow = 0; crow < crows; ++crow )
  761.         {
  762.         temprptr = rowptr[crow] + leftcol;
  763.         growsumptr[crow][leftcol] = 0L;
  764.         for ( ccol = 0; ccol < ccols; ++ccol )
  765.             growsumptr[crow][leftcol] += 
  766.                 PNM_GET1( *(temprptr + ccol) );
  767.         gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  768.         }
  769.         tempgsum = gsum + 0.5;
  770.         CHECK_GRAY;
  771.         PNM_ASSIGN1( outputrow[col], g );
  772.         }
  773.     else
  774.         {
  775.         gsum = 0.0;
  776.         leftcol = col - ccolso2;
  777.         subcol = col - ccolso2 - 1;
  778.         addcol = col + ccolso2;
  779.         for ( crow = 0; crow < crows; ++crow )
  780.         {
  781.         growsumptr[crow][leftcol] = growsumptr[crow][subcol]
  782.             - PNM_GET1( rowptr[crow][subcol] )
  783.             + PNM_GET1( rowptr[crow][addcol] );
  784.         gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  785.         }
  786.         tempgsum = gsum + 0.5;
  787.         CHECK_GRAY;
  788.         PNM_ASSIGN1( outputrow[col], g );
  789.         }
  790.         }
  791.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  792.  
  793.  
  794.     /* For all subsequent rows */
  795.  
  796.     subrow = crows;
  797.     addrow = crows - 1;
  798.     crowsp1 = crows + 1;
  799.     ++row;
  800.     for ( ; row < rows; ++row )
  801.     {
  802.     temprow = row % crowsp1;
  803.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  804.     if ( PNM_FORMAT_TYPE(format) != newformat )
  805.         pnm_promoteformatrow(
  806.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  807.  
  808.     temprow = (row + 2) % crowsp1;
  809.     i = 0;
  810.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  811.         {
  812.         rowptr[i] = xelbuf[irow];
  813.         growsumptr[i] = growsum[irow];
  814.         }
  815.     for (irow = 0; irow < temprow; ++irow, ++i)
  816.         {
  817.         rowptr[i] = xelbuf[irow];
  818.         growsumptr[i] = growsum[irow];
  819.         }
  820.  
  821.     for ( col = 0; col < cols; ++col )
  822.         {
  823.         if ( col < ccolso2 || col >= cols - ccolso2 )
  824.         outputrow[col] = rowptr[crowso2][col];
  825.         else if ( col == ccolso2 )
  826.         {
  827.         gsum = 0.0;
  828.         leftcol = col - ccolso2;
  829.         growsumptr[addrow][leftcol] = 0L;
  830.         for ( ccol = 0; ccol < ccols; ++ccol )
  831.             growsumptr[addrow][leftcol] += 
  832.             PNM_GET1( rowptr[addrow][leftcol + ccol] );
  833.         for ( crow = 0; crow < crows; ++crow )
  834.             gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  835.         tempgsum = gsum + 0.5;
  836.         CHECK_GRAY;
  837.         PNM_ASSIGN1( outputrow[col], g );
  838.         }
  839.         else
  840.         {
  841.         gsum = 0.0;
  842.         leftcol = col - ccolso2;
  843.         subcol = col - ccolso2 - 1;
  844.         addcol = col + ccolso2;  
  845.         growsumptr[addrow][leftcol] = growsumptr[addrow][subcol]
  846.             - PNM_GET1( rowptr[addrow][subcol] )
  847.             + PNM_GET1( rowptr[addrow][addcol] );
  848.         for ( crow = 0; crow < crows; ++crow )
  849.             gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  850.         tempgsum = gsum + 0.5;
  851.         CHECK_GRAY;
  852.         PNM_ASSIGN1( outputrow[col], g );
  853.         }
  854.         }
  855.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  856.     }
  857.  
  858.     /* Now write out the remaining unconvolved rows in xelbuf. */
  859.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  860.     pnm_writepnmrow(
  861.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  862.  
  863.     pm_close( stdout );
  864.     }
  865.  
  866.  
  867. /* PGM Vertical Convolution
  868. **
  869. ** Uses column sums as in Mean Convolution.
  870. **
  871. */
  872.  
  873.  
  874. static void
  875. pgm_vertical_convolve()
  876.     {
  877.     register int ccol, col;
  878.     xel** xelbuf;
  879.     xel* outputrow;
  880.     xelval g;
  881.     int row, crow;
  882.     xel **rowptr, *temprptr;
  883.     int leftcol;
  884.     int i, irow;
  885.     int toprow, temprow;
  886.     int subrow, addrow;
  887.     int tempcol;
  888.     float gsum;
  889.     long *gcolumnsum;
  890.     int crowsp1;
  891.     int addcol;
  892.     long tempgsum;
  893.  
  894.     /* Allocate space for one convolution-matrix's worth of rows, plus
  895.     ** a row output buffer. VERTICAL uses an extra row. */
  896.     xelbuf = pnm_allocarray( cols, crows + 1 );
  897.     outputrow = pnm_allocrow( cols );
  898.  
  899.     /* Allocate array of pointers to xelbuf */
  900.     rowptr = (xel **) pnm_allocarray( 1, crows + 1 );
  901.  
  902.     /* Allocate space for intermediate column sums */
  903.     gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  904.     for ( col = 0; col < cols; ++col )
  905.     gcolumnsum[col] = 0L;
  906.  
  907.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  908.  
  909.     /* Read in one convolution-matrix's worth of image, less one row. */
  910.     for ( row = 0; row < crows - 1; ++row )
  911.     {
  912.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  913.     if ( PNM_FORMAT_TYPE(format) != newformat )
  914.         pnm_promoteformatrow(
  915.         xelbuf[row], cols, maxval, format, maxval, newformat );
  916.     /* Write out just the part we're not going to convolve. */
  917.     if ( row < crowso2 )
  918.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  919.     }
  920.  
  921.     /* Now the rest of the image - read in the row at the end of
  922.     ** xelbuf, and convolve and write out the row in the middle.
  923.     */
  924.     /* For first row only */
  925.  
  926.     toprow = row + 1;
  927.     temprow = row % crows;
  928.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  929.     if ( PNM_FORMAT_TYPE(format) != newformat )
  930.     pnm_promoteformatrow(
  931.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  932.  
  933.     /* Arrange rowptr to eliminate the use of mod function to determine
  934.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  935.     */
  936.     temprow = toprow % crows;
  937.     i = 0;
  938.     for (irow = temprow; irow < crows; ++i, ++irow)
  939.     rowptr[i] = xelbuf[irow];
  940.     for (irow = 0; irow < temprow; ++irow, ++i)
  941.     rowptr[i] = xelbuf[irow];
  942.  
  943.     for ( col = 0; col < cols; ++col )
  944.     {
  945.     if ( col < ccolso2 || col >= cols - ccolso2 )
  946.         outputrow[col] = rowptr[crowso2][col];
  947.     else if ( col == ccolso2 )
  948.         {
  949.         gsum = 0.0;
  950.         leftcol = col - ccolso2;
  951.         for ( crow = 0; crow < crows; ++crow )
  952.         {
  953.         temprptr = rowptr[crow] + leftcol;
  954.         for ( ccol = 0; ccol < ccols; ++ccol )
  955.             gcolumnsum[leftcol + ccol] += 
  956.             PNM_GET1( *(temprptr + ccol) );
  957.         }
  958.         for ( ccol = 0; ccol < ccols; ++ccol)
  959.         gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  960.         tempgsum = gsum + 0.5;
  961.         CHECK_GRAY;
  962.         PNM_ASSIGN1( outputrow[col], g );
  963.         }
  964.     else
  965.         {
  966.         gsum = 0.0;
  967.         leftcol = col - ccolso2;
  968.         addcol = col + ccolso2;  
  969.         for ( crow = 0; crow < crows; ++crow )
  970.         gcolumnsum[addcol] += PNM_GET1( rowptr[crow][addcol] );
  971.         for ( ccol = 0; ccol < ccols; ++ccol )
  972.         gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  973.         tempgsum = gsum + 0.5;
  974.         CHECK_GRAY;
  975.         PNM_ASSIGN1( outputrow[col], g );
  976.         }
  977.     }
  978.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  979.  
  980.     /* For all subsequent rows */
  981.     subrow = crows;
  982.     addrow = crows - 1;
  983.     crowsp1 = crows + 1;
  984.     ++row;
  985.     for ( ; row < rows; ++row )
  986.     {
  987.     toprow = row + 1;
  988.     temprow = row % (crows +1);
  989.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  990.     if ( PNM_FORMAT_TYPE(format) != newformat )
  991.         pnm_promoteformatrow(
  992.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  993.  
  994.     /* Arrange rowptr to eliminate the use of mod function to determine
  995.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  996.     */
  997.     temprow = (toprow + 1) % crowsp1;
  998.     i = 0;
  999.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  1000.         rowptr[i] = xelbuf[irow];
  1001.     for (irow = 0; irow < temprow; ++irow, ++i)
  1002.         rowptr[i] = xelbuf[irow];
  1003.  
  1004.     for ( col = 0; col < cols; ++col )
  1005.         {
  1006.         if ( col < ccolso2 || col >= cols - ccolso2 )
  1007.         outputrow[col] = rowptr[crowso2][col];
  1008.         else if ( col == ccolso2 )
  1009.         {
  1010.         gsum = 0.0;
  1011.         leftcol = col - ccolso2;
  1012.         for ( ccol = 0; ccol < ccols; ++ccol )
  1013.             {
  1014.             tempcol = leftcol + ccol;
  1015.             gcolumnsum[tempcol] = gcolumnsum[tempcol] 
  1016.             - PNM_GET1( rowptr[subrow][ccol] )
  1017.             + PNM_GET1( rowptr[addrow][ccol] );
  1018.             gsum = gsum + gcolumnsum[tempcol] * gweights[0][ccol];
  1019.             }
  1020.         tempgsum = gsum + 0.5;
  1021.         CHECK_GRAY;
  1022.         PNM_ASSIGN1( outputrow[col], g );
  1023.         }
  1024.         else
  1025.         {
  1026.         gsum = 0.0;
  1027.         leftcol = col - ccolso2;
  1028.         addcol = col + ccolso2;
  1029.         gcolumnsum[addcol] = gcolumnsum[addcol]
  1030.             - PNM_GET1( rowptr[subrow][addcol] )
  1031.             + PNM_GET1( rowptr[addrow][addcol] );
  1032.         for ( ccol = 0; ccol < ccols; ++ccol )
  1033.             gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  1034.         tempgsum = gsum + 0.5;
  1035.         CHECK_GRAY;
  1036.         PNM_ASSIGN1( outputrow[col], g );
  1037.         }
  1038.         }
  1039.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1040.     }
  1041.  
  1042.     /* Now write out the remaining unconvolved rows in xelbuf. */
  1043.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  1044.     pnm_writepnmrow(
  1045.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  1046.  
  1047.     pm_close( stdout );
  1048.     }
  1049.  
  1050.  
  1051.  
  1052.  
  1053. /* PPM General Convolution Algorithm
  1054. **
  1055. ** No redundancy in convolution matrix.  Just use brute force.
  1056. ** See pgm_general_convolve() for more details.
  1057. */
  1058.  
  1059. static void
  1060. ppm_general_convolve()
  1061.     {
  1062.     register int ccol, col;
  1063.     xel** xelbuf;
  1064.     xel* outputrow;
  1065.     xel x, y;
  1066.     xelval r, g, b;
  1067.     int row, crow;
  1068.     float rsum, gsum, bsum;
  1069.     xel **rowptr, *temprptr;
  1070.     int toprow, temprow;
  1071.     int i, irow;
  1072.     int leftcol;
  1073.     long temprsum, tempgsum, tempbsum;
  1074.  
  1075.     /* Allocate space for one convolution-matrix's worth of rows, plus
  1076.     ** a row output buffer. */
  1077.     xelbuf = pnm_allocarray( cols, crows );
  1078.     outputrow = pnm_allocrow( cols );
  1079.  
  1080.     /* Allocate array of pointers to xelbuf */
  1081.     rowptr = (xel **) pnm_allocarray( 1, crows );
  1082.  
  1083.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  1084.  
  1085.     /* Read in one convolution-matrix's worth of image, less one row. */
  1086.     for ( row = 0; row < crows - 1; ++row )
  1087.     {
  1088.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  1089.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1090.         pnm_promoteformatrow(
  1091.         xelbuf[row], cols, maxval, format, maxval, newformat );
  1092.     /* Write out just the part we're not going to convolve. */
  1093.     if ( row < crowso2 )
  1094.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  1095.     }
  1096.  
  1097.     /* Now the rest of the image - read in the row at the end of
  1098.     ** xelbuf, and convolve and write out the row in the middle.
  1099.     */
  1100.     for ( ; row < rows; ++row )
  1101.     {
  1102.     toprow = row + 1;
  1103.     temprow = row % crows;
  1104.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1105.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1106.         pnm_promoteformatrow(
  1107.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1108.  
  1109.     /* Arrange rowptr to eliminate the use of mod function to determine
  1110.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  1111.     */
  1112.     temprow = toprow % crows;
  1113.     i = 0;
  1114.     for (irow = temprow; irow < crows; ++i, ++irow)
  1115.         rowptr[i] = xelbuf[irow];
  1116.     for (irow = 0; irow < temprow; ++irow, ++i)
  1117.         rowptr[i] = xelbuf[irow];
  1118.  
  1119.     for ( col = 0; col < cols; ++col )
  1120.         {
  1121.         if ( col < ccolso2 || col >= cols - ccolso2 )
  1122.         outputrow[col] = rowptr[crowso2][col];
  1123.         else
  1124.         {
  1125.         leftcol = col - ccolso2;
  1126.         rsum = gsum = bsum = 0.0;
  1127.         for ( crow = 0; crow < crows; ++crow )
  1128.             {
  1129.             temprptr = rowptr[crow] + leftcol;
  1130.             for ( ccol = 0; ccol < ccols; ++ccol )
  1131.             {
  1132.             rsum += PPM_GETR( *(temprptr + ccol) )
  1133.                 * rweights[crow][ccol];
  1134.             gsum += PPM_GETG( *(temprptr + ccol) )
  1135.                 * gweights[crow][ccol];
  1136.             bsum += PPM_GETB( *(temprptr + ccol) )
  1137.                 * bweights[crow][ccol];
  1138.             }
  1139.             }
  1140.             temprsum = rsum + 0.5;
  1141.             tempgsum = gsum + 0.5;
  1142.             tempbsum = bsum + 0.5;
  1143.             CHECK_RED;
  1144.             CHECK_GREEN;
  1145.             CHECK_BLUE;
  1146.             PPM_ASSIGN( outputrow[col], r, g, b );
  1147.         }
  1148.         }
  1149.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1150.     }
  1151.  
  1152.     /* Now write out the remaining unconvolved rows in xelbuf. */
  1153.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  1154.     pnm_writepnmrow(
  1155.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  1156.  
  1157.     pm_close( stdout );
  1158.     }
  1159.  
  1160.  
  1161. /* PPM Mean Convolution
  1162. **
  1163. ** Same as pgm_mean_convolve() but for PPM.
  1164. **
  1165. */
  1166.  
  1167. static void
  1168. ppm_mean_convolve()
  1169.     {
  1170.     register int ccol, col;
  1171.     xel** xelbuf;
  1172.     xel* outputrow;
  1173.     xelval r, g, b;
  1174.     int row, crow;
  1175.     xel **rowptr, *temprptr;
  1176.     int leftcol;
  1177.     int i, irow;
  1178.     int toprow, temprow;
  1179.     int subrow, addrow;
  1180.     int subcol, addcol;
  1181.     long risum, gisum, bisum;
  1182.     float rsum, gsum, bsum;
  1183.     long temprsum, tempgsum, tempbsum;
  1184.     int tempcol, crowsp1;
  1185.     long *rcolumnsum, *gcolumnsum, *bcolumnsum;
  1186.  
  1187.     /* Allocate space for one convolution-matrix's worth of rows, plus
  1188.     ** a row output buffer.  MEAN uses an extra row. */
  1189.     xelbuf = pnm_allocarray( cols, crows + 1 );
  1190.     outputrow = pnm_allocrow( cols );
  1191.  
  1192.     /* Allocate array of pointers to xelbuf. MEAN uses an extra row. */
  1193.     rowptr = (xel **) pnm_allocarray( 1, crows + 1);
  1194.  
  1195.     /* Allocate space for intermediate column sums */
  1196.     rcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1197.     gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1198.     bcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1199.     for ( col = 0; col < cols; ++col )
  1200.     {
  1201.     rcolumnsum[col] = 0L;
  1202.     gcolumnsum[col] = 0L;
  1203.     bcolumnsum[col] = 0L;
  1204.     }
  1205.  
  1206.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  1207.  
  1208.     /* Read in one convolution-matrix's worth of image, less one row. */
  1209.     for ( row = 0; row < crows - 1; ++row )
  1210.     {
  1211.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  1212.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1213.         pnm_promoteformatrow(
  1214.         xelbuf[row], cols, maxval, format, maxval, newformat );
  1215.     /* Write out just the part we're not going to convolve. */
  1216.     if ( row < crowso2 )
  1217.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  1218.     }
  1219.  
  1220.     /* Do first real row only */
  1221.     subrow = crows;
  1222.     addrow = crows - 1;
  1223.     toprow = row + 1;
  1224.     temprow = row % crows;
  1225.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1226.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1227.     pnm_promoteformatrow(
  1228.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1229.  
  1230.     temprow = toprow % crows;
  1231.     i = 0;
  1232.     for (irow = temprow; irow < crows; ++i, ++irow)
  1233.     rowptr[i] = xelbuf[irow];
  1234.     for (irow = 0; irow < temprow; ++irow, ++i)
  1235.     rowptr[i] = xelbuf[irow];
  1236.  
  1237.     risum = 0L;
  1238.     gisum = 0L;
  1239.     bisum = 0L;
  1240.     for ( col = 0; col < cols; ++col )
  1241.     {
  1242.     if ( col < ccolso2 || col >= cols - ccolso2 )
  1243.         outputrow[col] = rowptr[crowso2][col];
  1244.     else if ( col == ccolso2 )
  1245.         {
  1246.         leftcol = col - ccolso2;
  1247.         for ( crow = 0; crow < crows; ++crow )
  1248.         {
  1249.         temprptr = rowptr[crow] + leftcol;
  1250.         for ( ccol = 0; ccol < ccols; ++ccol )
  1251.             {
  1252.             rcolumnsum[leftcol + ccol] += 
  1253.             PPM_GETR( *(temprptr + ccol) );
  1254.             gcolumnsum[leftcol + ccol] += 
  1255.             PPM_GETG( *(temprptr + ccol) );
  1256.             bcolumnsum[leftcol + ccol] += 
  1257.             PPM_GETB( *(temprptr + ccol) );
  1258.             }
  1259.         }
  1260.         for ( ccol = 0; ccol < ccols; ++ccol)
  1261.         {
  1262.         risum += rcolumnsum[leftcol + ccol];
  1263.         gisum += gcolumnsum[leftcol + ccol];
  1264.         bisum += bcolumnsum[leftcol + ccol];
  1265.         }
  1266.         temprsum = (float) risum * rmeanweight + 0.5;
  1267.         tempgsum = (float) gisum * gmeanweight + 0.5;
  1268.         tempbsum = (float) bisum * bmeanweight + 0.5;
  1269.         CHECK_RED;
  1270.         CHECK_GREEN;
  1271.         CHECK_BLUE;
  1272.         PPM_ASSIGN( outputrow[col], r, g, b );
  1273.         }
  1274.     else
  1275.         {
  1276.         /* Column numbers to subtract or add to isum */
  1277.         subcol = col - ccolso2 - 1;
  1278.         addcol = col + ccolso2;  
  1279.         for ( crow = 0; crow < crows; ++crow )
  1280.         {
  1281.         rcolumnsum[addcol] += PPM_GETR( rowptr[crow][addcol] );
  1282.         gcolumnsum[addcol] += PPM_GETG( rowptr[crow][addcol] );
  1283.         bcolumnsum[addcol] += PPM_GETB( rowptr[crow][addcol] );
  1284.         }
  1285.         risum = risum - rcolumnsum[subcol] + rcolumnsum[addcol];
  1286.         gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol];
  1287.         bisum = bisum - bcolumnsum[subcol] + bcolumnsum[addcol];
  1288.         temprsum = (float) risum * rmeanweight + 0.5;
  1289.         tempgsum = (float) gisum * gmeanweight + 0.5;
  1290.         tempbsum = (float) bisum * bmeanweight + 0.5;
  1291.         CHECK_RED;
  1292.         CHECK_GREEN;
  1293.         CHECK_BLUE;
  1294.         PPM_ASSIGN( outputrow[col], r, g, b );
  1295.         }
  1296.     }
  1297.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1298.  
  1299.     ++row;
  1300.     /* For all subsequent rows do it this way as the columnsums have been
  1301.     ** generated.  Now we can use them to reduce further calculations.
  1302.     */
  1303.     crowsp1 = crows + 1;
  1304.     for ( ; row < rows; ++row )
  1305.     {
  1306.     toprow = row + 1;
  1307.     temprow = row % (crows + 1);
  1308.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1309.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1310.         pnm_promoteformatrow(
  1311.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1312.  
  1313.     /* This rearrangement using crows+1 rowptrs and xelbufs will cause
  1314.     ** rowptr[0..crows-1] to always hold active xelbufs and for 
  1315.     ** rowptr[crows] to always hold the oldest (top most) xelbuf.
  1316.     */
  1317.     temprow = (toprow + 1) % crowsp1;
  1318.     i = 0;
  1319.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  1320.         rowptr[i] = xelbuf[irow];
  1321.     for (irow = 0; irow < temprow; ++irow, ++i)
  1322.         rowptr[i] = xelbuf[irow];
  1323.  
  1324.     risum = 0L;
  1325.     gisum = 0L;
  1326.     bisum = 0L;
  1327.     for ( col = 0; col < cols; ++col )
  1328.         {
  1329.         if ( col < ccolso2 || col >= cols - ccolso2 )
  1330.         outputrow[col] = rowptr[crowso2][col];
  1331.         else if ( col == ccolso2 )
  1332.         {
  1333.         leftcol = col - ccolso2;
  1334.         for ( ccol = 0; ccol < ccols; ++ccol )
  1335.             {
  1336.             tempcol = leftcol + ccol;
  1337.             rcolumnsum[tempcol] = rcolumnsum[tempcol]
  1338.             - PPM_GETR( rowptr[subrow][ccol] )
  1339.             + PPM_GETR( rowptr[addrow][ccol] );
  1340.             risum += rcolumnsum[tempcol];
  1341.             gcolumnsum[tempcol] = gcolumnsum[tempcol]
  1342.             - PPM_GETG( rowptr[subrow][ccol] )
  1343.             + PPM_GETG( rowptr[addrow][ccol] );
  1344.             gisum += gcolumnsum[tempcol];
  1345.             bcolumnsum[tempcol] = bcolumnsum[tempcol]
  1346.             - PPM_GETB( rowptr[subrow][ccol] )
  1347.             + PPM_GETB( rowptr[addrow][ccol] );
  1348.             bisum += bcolumnsum[tempcol];
  1349.             }
  1350.         temprsum = (float) risum * rmeanweight + 0.5;
  1351.         tempgsum = (float) gisum * gmeanweight + 0.5;
  1352.         tempbsum = (float) bisum * bmeanweight + 0.5;
  1353.         CHECK_RED;
  1354.         CHECK_GREEN;
  1355.         CHECK_BLUE;
  1356.         PPM_ASSIGN( outputrow[col], r, g, b );
  1357.         }
  1358.         else
  1359.         {
  1360.         /* Column numbers to subtract or add to isum */
  1361.         subcol = col - ccolso2 - 1;
  1362.         addcol = col + ccolso2;  
  1363.         rcolumnsum[addcol] = rcolumnsum[addcol]
  1364.             - PPM_GETR( rowptr[subrow][addcol] )
  1365.             + PPM_GETR( rowptr[addrow][addcol] );
  1366.         risum = risum - rcolumnsum[subcol] + rcolumnsum[addcol];
  1367.         gcolumnsum[addcol] = gcolumnsum[addcol]
  1368.             - PPM_GETG( rowptr[subrow][addcol] )
  1369.             + PPM_GETG( rowptr[addrow][addcol] );
  1370.         gisum = gisum - gcolumnsum[subcol] + gcolumnsum[addcol];
  1371.         bcolumnsum[addcol] = bcolumnsum[addcol]
  1372.             - PPM_GETB( rowptr[subrow][addcol] )
  1373.             + PPM_GETB( rowptr[addrow][addcol] );
  1374.         bisum = bisum - bcolumnsum[subcol] + bcolumnsum[addcol];
  1375.         temprsum = (float) risum * rmeanweight + 0.5;
  1376.         tempgsum = (float) gisum * gmeanweight + 0.5;
  1377.         tempbsum = (float) bisum * bmeanweight + 0.5;
  1378.         CHECK_RED;
  1379.         CHECK_GREEN;
  1380.         CHECK_BLUE;
  1381.         PPM_ASSIGN( outputrow[col], r, g, b );
  1382.         }
  1383.         }
  1384.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1385.     }
  1386.  
  1387.     /* Now write out the remaining unconvolved rows in xelbuf. */
  1388.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  1389.     pnm_writepnmrow(
  1390.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  1391.  
  1392.     pm_close( stdout );
  1393.     }
  1394.  
  1395.  
  1396. /* PPM Horizontal Convolution
  1397. **
  1398. ** Same as pgm_horizontal_convolve()
  1399. **
  1400. **/
  1401.  
  1402. static void
  1403. ppm_horizontal_convolve()
  1404.     {
  1405.     register int ccol, col;
  1406.     xel** xelbuf;
  1407.     xel* outputrow;
  1408.     xel x;
  1409.     xelval r, g, b;
  1410.     int row, crow;
  1411.     xel **rowptr, *temprptr;
  1412.     int leftcol;
  1413.     int i, irow;
  1414.     int temprow;
  1415.     int subcol, addcol;
  1416.     float rsum, gsum, bsum;
  1417.     int addrow, subrow;
  1418.     long **rrowsum, **rrowsumptr;
  1419.     long **growsum, **growsumptr;
  1420.     long **browsum, **browsumptr;
  1421.     int crowsp1;
  1422.     long temprsum, tempgsum, tempbsum;
  1423.  
  1424.     /* Allocate space for one convolution-matrix's worth of rows, plus
  1425.     ** a row output buffer. */
  1426.     xelbuf = pnm_allocarray( cols, crows + 1 );
  1427.     outputrow = pnm_allocrow( cols );
  1428.  
  1429.     /* Allocate array of pointers to xelbuf */
  1430.     rowptr = (xel **) pnm_allocarray( 1, crows + 1);
  1431.  
  1432.     /* Allocate intermediate row sums.  HORIZONTAL uses an extra row */
  1433.     rrowsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) );
  1434.     rrowsumptr = (long **) pnm_allocarray( 1, crows + 1);
  1435.     growsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) );
  1436.     growsumptr = (long **) pnm_allocarray( 1, crows + 1);
  1437.     browsum = (long **) pm_allocarray( cols, crows + 1, sizeof(long) );
  1438.     browsumptr = (long **) pnm_allocarray( 1, crows + 1);
  1439.  
  1440.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  1441.  
  1442.     /* Read in one convolution-matrix's worth of image, less one row. */
  1443.     for ( row = 0; row < crows - 1; ++row )
  1444.     {
  1445.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  1446.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1447.         pnm_promoteformatrow(
  1448.         xelbuf[row], cols, maxval, format, maxval, newformat );
  1449.     /* Write out just the part we're not going to convolve. */
  1450.     if ( row < crowso2 )
  1451.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  1452.     }
  1453.  
  1454.     /* First row only */
  1455.     temprow = row % crows;
  1456.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1457.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1458.     pnm_promoteformatrow(
  1459.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1460.  
  1461.     temprow = (row + 1) % crows;
  1462.     i = 0;
  1463.     for (irow = temprow; irow < crows; ++i, ++irow)
  1464.     rowptr[i] = xelbuf[irow];
  1465.     for (irow = 0; irow < temprow; ++irow, ++i)
  1466.     rowptr[i] = xelbuf[irow];
  1467.  
  1468.     for ( crow = 0; crow < crows; ++crow )
  1469.     {
  1470.     rrowsumptr[crow] = rrowsum[crow];
  1471.     growsumptr[crow] = growsum[crow];
  1472.     browsumptr[crow] = browsum[crow];
  1473.     }
  1474.  
  1475.     for ( col = 0; col < cols; ++col )
  1476.     {
  1477.     if ( col < ccolso2 || col >= cols - ccolso2 )
  1478.         outputrow[col] = rowptr[crowso2][col];
  1479.     else if ( col == ccolso2 )
  1480.         {
  1481.         leftcol = col - ccolso2;
  1482.         rsum = 0.0;
  1483.         gsum = 0.0;
  1484.         bsum = 0.0;
  1485.         for ( crow = 0; crow < crows; ++crow )
  1486.         {
  1487.         temprptr = rowptr[crow] + leftcol;
  1488.         rrowsumptr[crow][leftcol] = 0L;
  1489.         growsumptr[crow][leftcol] = 0L;
  1490.         browsumptr[crow][leftcol] = 0L;
  1491.         for ( ccol = 0; ccol < ccols; ++ccol )
  1492.             {
  1493.             rrowsumptr[crow][leftcol] += 
  1494.                 PPM_GETR( *(temprptr + ccol) );
  1495.             growsumptr[crow][leftcol] += 
  1496.                 PPM_GETG( *(temprptr + ccol) );
  1497.             browsumptr[crow][leftcol] += 
  1498.                 PPM_GETB( *(temprptr + ccol) );
  1499.             }
  1500.         rsum += rrowsumptr[crow][leftcol] * rweights[crow][0];
  1501.         gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  1502.         bsum += browsumptr[crow][leftcol] * bweights[crow][0];
  1503.         }
  1504.         temprsum = rsum + 0.5;
  1505.         tempgsum = gsum + 0.5;
  1506.         tempbsum = bsum + 0.5;
  1507.         CHECK_RED;
  1508.         CHECK_GREEN;
  1509.         CHECK_BLUE;
  1510.         PPM_ASSIGN( outputrow[col], r, g, b );
  1511.         }
  1512.     else
  1513.         {
  1514.         rsum = 0.0;
  1515.         gsum = 0.0;
  1516.         bsum = 0.0;
  1517.         leftcol = col - ccolso2;
  1518.         subcol = col - ccolso2 - 1;
  1519.         addcol = col + ccolso2;
  1520.         for ( crow = 0; crow < crows; ++crow )
  1521.         {
  1522.         rrowsumptr[crow][leftcol] = rrowsumptr[crow][subcol]
  1523.             - PPM_GETR( rowptr[crow][subcol] )
  1524.             + PPM_GETR( rowptr[crow][addcol] );
  1525.         rsum += rrowsumptr[crow][leftcol] * rweights[crow][0];
  1526.         growsumptr[crow][leftcol] = growsumptr[crow][subcol]
  1527.             - PPM_GETG( rowptr[crow][subcol] )
  1528.             + PPM_GETG( rowptr[crow][addcol] );
  1529.         gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  1530.         browsumptr[crow][leftcol] = browsumptr[crow][subcol]
  1531.             - PPM_GETB( rowptr[crow][subcol] )
  1532.             + PPM_GETB( rowptr[crow][addcol] );
  1533.         bsum += browsumptr[crow][leftcol] * bweights[crow][0];
  1534.         }
  1535.         temprsum = rsum + 0.5;
  1536.         tempgsum = gsum + 0.5;
  1537.         tempbsum = bsum + 0.5;
  1538.         CHECK_RED;
  1539.         CHECK_GREEN;
  1540.         CHECK_BLUE;
  1541.         PPM_ASSIGN( outputrow[col], r, g, b );
  1542.         }
  1543.         }
  1544.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1545.  
  1546.  
  1547.     /* For all subsequent rows */
  1548.  
  1549.     subrow = crows;
  1550.     addrow = crows - 1;
  1551.     crowsp1 = crows + 1;
  1552.     ++row;
  1553.     for ( ; row < rows; ++row )
  1554.     {
  1555.     temprow = row % crowsp1;
  1556.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1557.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1558.         pnm_promoteformatrow(
  1559.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1560.  
  1561.     temprow = (row + 2) % crowsp1;
  1562.     i = 0;
  1563.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  1564.         {
  1565.         rowptr[i] = xelbuf[irow];
  1566.         rrowsumptr[i] = rrowsum[irow];
  1567.         growsumptr[i] = growsum[irow];
  1568.         browsumptr[i] = browsum[irow];
  1569.         }
  1570.     for (irow = 0; irow < temprow; ++irow, ++i)
  1571.         {
  1572.         rowptr[i] = xelbuf[irow];
  1573.         rrowsumptr[i] = rrowsum[irow];
  1574.         growsumptr[i] = growsum[irow];
  1575.         browsumptr[i] = browsum[irow];
  1576.         }
  1577.  
  1578.     for ( col = 0; col < cols; ++col )
  1579.         {
  1580.         if ( col < ccolso2 || col >= cols - ccolso2 )
  1581.         outputrow[col] = rowptr[crowso2][col];
  1582.         else if ( col == ccolso2 )
  1583.         {
  1584.         rsum = 0.0;
  1585.         gsum = 0.0;
  1586.         bsum = 0.0;
  1587.         leftcol = col - ccolso2;
  1588.         rrowsumptr[addrow][leftcol] = 0L;
  1589.         growsumptr[addrow][leftcol] = 0L;
  1590.         browsumptr[addrow][leftcol] = 0L;
  1591.         for ( ccol = 0; ccol < ccols; ++ccol )
  1592.             {
  1593.             rrowsumptr[addrow][leftcol] += 
  1594.             PPM_GETR( rowptr[addrow][leftcol + ccol] );
  1595.             growsumptr[addrow][leftcol] += 
  1596.             PPM_GETG( rowptr[addrow][leftcol + ccol] );
  1597.             browsumptr[addrow][leftcol] += 
  1598.             PPM_GETB( rowptr[addrow][leftcol + ccol] );
  1599.             }
  1600.         for ( crow = 0; crow < crows; ++crow )
  1601.             {
  1602.             rsum += rrowsumptr[crow][leftcol] * rweights[crow][0];
  1603.             gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  1604.             bsum += browsumptr[crow][leftcol] * bweights[crow][0];
  1605.             }
  1606.         temprsum = rsum + 0.5;
  1607.         tempgsum = gsum + 0.5;
  1608.         tempbsum = bsum + 0.5;
  1609.         CHECK_RED;
  1610.         CHECK_GREEN;
  1611.         CHECK_BLUE;
  1612.         PPM_ASSIGN( outputrow[col], r, g, b );
  1613.         }
  1614.         else
  1615.         {
  1616.         rsum = 0.0;
  1617.         gsum = 0.0;
  1618.         bsum = 0.0;
  1619.         leftcol = col - ccolso2;
  1620.         subcol = col - ccolso2 - 1;
  1621.         addcol = col + ccolso2;  
  1622.         rrowsumptr[addrow][leftcol] = rrowsumptr[addrow][subcol]
  1623.             - PPM_GETR( rowptr[addrow][subcol] )
  1624.             + PPM_GETR( rowptr[addrow][addcol] );
  1625.         growsumptr[addrow][leftcol] = growsumptr[addrow][subcol]
  1626.             - PPM_GETG( rowptr[addrow][subcol] )
  1627.             + PPM_GETG( rowptr[addrow][addcol] );
  1628.         browsumptr[addrow][leftcol] = browsumptr[addrow][subcol]
  1629.             - PPM_GETB( rowptr[addrow][subcol] )
  1630.             + PPM_GETB( rowptr[addrow][addcol] );
  1631.         for ( crow = 0; crow < crows; ++crow )
  1632.             {
  1633.             rsum += rrowsumptr[crow][leftcol] * rweights[crow][0];
  1634.             gsum += growsumptr[crow][leftcol] * gweights[crow][0];
  1635.             bsum += browsumptr[crow][leftcol] * bweights[crow][0];
  1636.             }
  1637.         temprsum = rsum + 0.5;
  1638.         tempgsum = gsum + 0.5;
  1639.         tempbsum = bsum + 0.5;
  1640.         CHECK_RED;
  1641.         CHECK_GREEN;
  1642.         CHECK_BLUE;
  1643.         PPM_ASSIGN( outputrow[col], r, g, b );
  1644.         }
  1645.         }
  1646.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1647.     }
  1648.  
  1649.     /* Now write out the remaining unconvolved rows in xelbuf. */
  1650.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  1651.     pnm_writepnmrow(
  1652.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  1653.  
  1654.     pm_close( stdout );
  1655.     }
  1656.  
  1657.  
  1658. /* PPM Vertical Convolution
  1659. **
  1660. ** Same as pgm_vertical_convolve()
  1661. **
  1662. */
  1663.  
  1664. static void
  1665. ppm_vertical_convolve()
  1666.     {
  1667.     register int ccol, col;
  1668.     xel** xelbuf;
  1669.     xel* outputrow;
  1670.     xelval r, g, b;
  1671.     int row, crow;
  1672.     xel **rowptr, *temprptr;
  1673.     int leftcol;
  1674.     int i, irow;
  1675.     int toprow, temprow;
  1676.     int subrow, addrow;
  1677.     int tempcol;
  1678.     float rsum, gsum, bsum;
  1679.     long *rcolumnsum, *gcolumnsum, *bcolumnsum;
  1680.     int crowsp1;
  1681.     int addcol;
  1682.     long temprsum, tempgsum, tempbsum;
  1683.  
  1684.     /* Allocate space for one convolution-matrix's worth of rows, plus
  1685.     ** a row output buffer. VERTICAL uses an extra row. */
  1686.     xelbuf = pnm_allocarray( cols, crows + 1 );
  1687.     outputrow = pnm_allocrow( cols );
  1688.  
  1689.     /* Allocate array of pointers to xelbuf */
  1690.     rowptr = (xel **) pnm_allocarray( 1, crows + 1 );
  1691.  
  1692.     /* Allocate space for intermediate column sums */
  1693.     rcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1694.     gcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1695.     bcolumnsum = (long *) pm_allocrow( cols, sizeof(long) );
  1696.     for ( col = 0; col < cols; ++col )
  1697.     {
  1698.     rcolumnsum[col] = 0L;
  1699.     gcolumnsum[col] = 0L;
  1700.     bcolumnsum[col] = 0L;
  1701.     }
  1702.  
  1703.     pnm_writepnminit( stdout, cols, rows, maxval, newformat, 0 );
  1704.  
  1705.     /* Read in one convolution-matrix's worth of image, less one row. */
  1706.     for ( row = 0; row < crows - 1; ++row )
  1707.     {
  1708.     pnm_readpnmrow( ifp, xelbuf[row], cols, maxval, format );
  1709.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1710.         pnm_promoteformatrow(
  1711.         xelbuf[row], cols, maxval, format, maxval, newformat );
  1712.     /* Write out just the part we're not going to convolve. */
  1713.     if ( row < crowso2 )
  1714.         pnm_writepnmrow( stdout, xelbuf[row], cols, maxval, newformat, 0 );
  1715.     }
  1716.  
  1717.     /* Now the rest of the image - read in the row at the end of
  1718.     ** xelbuf, and convolve and write out the row in the middle.
  1719.     */
  1720.     /* For first row only */
  1721.  
  1722.     toprow = row + 1;
  1723.     temprow = row % crows;
  1724.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1725.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1726.     pnm_promoteformatrow(
  1727.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1728.  
  1729.     /* Arrange rowptr to eliminate the use of mod function to determine
  1730.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  1731.     */
  1732.     temprow = toprow % crows;
  1733.     i = 0;
  1734.     for (irow = temprow; irow < crows; ++i, ++irow)
  1735.     rowptr[i] = xelbuf[irow];
  1736.     for (irow = 0; irow < temprow; ++irow, ++i)
  1737.     rowptr[i] = xelbuf[irow];
  1738.  
  1739.     for ( col = 0; col < cols; ++col )
  1740.     {
  1741.     if ( col < ccolso2 || col >= cols - ccolso2 )
  1742.         outputrow[col] = rowptr[crowso2][col];
  1743.     else if ( col == ccolso2 )
  1744.         {
  1745.         rsum = 0.0;
  1746.         gsum = 0.0;
  1747.         bsum = 0.0;
  1748.         leftcol = col - ccolso2;
  1749.         for ( crow = 0; crow < crows; ++crow )
  1750.         {
  1751.         temprptr = rowptr[crow] + leftcol;
  1752.         for ( ccol = 0; ccol < ccols; ++ccol )
  1753.             {
  1754.             rcolumnsum[leftcol + ccol] += 
  1755.             PPM_GETR( *(temprptr + ccol) );
  1756.             gcolumnsum[leftcol + ccol] += 
  1757.             PPM_GETG( *(temprptr + ccol) );
  1758.             bcolumnsum[leftcol + ccol] += 
  1759.             PPM_GETB( *(temprptr + ccol) );
  1760.             }
  1761.         }
  1762.         for ( ccol = 0; ccol < ccols; ++ccol)
  1763.         {
  1764.         rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol];
  1765.         gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  1766.         bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol];
  1767.         }
  1768.         temprsum = rsum + 0.5;
  1769.         tempgsum = gsum + 0.5;
  1770.         tempbsum = bsum + 0.5;
  1771.         CHECK_RED;
  1772.         CHECK_GREEN;
  1773.         CHECK_BLUE;
  1774.         PPM_ASSIGN( outputrow[col], r, g, b );
  1775.         }
  1776.     else
  1777.         {
  1778.         rsum = 0.0;
  1779.         gsum = 0.0;
  1780.         bsum = 0.0;
  1781.         leftcol = col - ccolso2;
  1782.         addcol = col + ccolso2;  
  1783.         for ( crow = 0; crow < crows; ++crow )
  1784.         {
  1785.         rcolumnsum[addcol] += PPM_GETR( rowptr[crow][addcol] );
  1786.         gcolumnsum[addcol] += PPM_GETG( rowptr[crow][addcol] );
  1787.         bcolumnsum[addcol] += PPM_GETB( rowptr[crow][addcol] );
  1788.         }
  1789.         for ( ccol = 0; ccol < ccols; ++ccol )
  1790.         {
  1791.         rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol];
  1792.         gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  1793.         bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol];
  1794.         }
  1795.         temprsum = rsum + 0.5;
  1796.         tempgsum = gsum + 0.5;
  1797.         tempbsum = bsum + 0.5;
  1798.         CHECK_RED;
  1799.         CHECK_GREEN;
  1800.         CHECK_BLUE;
  1801.         PPM_ASSIGN( outputrow[col], r, g, b );
  1802.         }
  1803.     }
  1804.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1805.  
  1806.     /* For all subsequent rows */
  1807.     subrow = crows;
  1808.     addrow = crows - 1;
  1809.     crowsp1 = crows + 1;
  1810.     ++row;
  1811.     for ( ; row < rows; ++row )
  1812.     {
  1813.     toprow = row + 1;
  1814.     temprow = row % (crows +1);
  1815.     pnm_readpnmrow( ifp, xelbuf[temprow], cols, maxval, format );
  1816.     if ( PNM_FORMAT_TYPE(format) != newformat )
  1817.         pnm_promoteformatrow(
  1818.         xelbuf[temprow], cols, maxval, format, maxval, newformat );
  1819.  
  1820.     /* Arrange rowptr to eliminate the use of mod function to determine
  1821.     ** which row of xelbuf is 0...crows.  Mod function can be very costly.
  1822.     */
  1823.     temprow = (toprow + 1) % crowsp1;
  1824.     i = 0;
  1825.     for (irow = temprow; irow < crowsp1; ++i, ++irow)
  1826.         rowptr[i] = xelbuf[irow];
  1827.     for (irow = 0; irow < temprow; ++irow, ++i)
  1828.         rowptr[i] = xelbuf[irow];
  1829.  
  1830.     for ( col = 0; col < cols; ++col )
  1831.         {
  1832.         if ( col < ccolso2 || col >= cols - ccolso2 )
  1833.         outputrow[col] = rowptr[crowso2][col];
  1834.         else if ( col == ccolso2 )
  1835.         {
  1836.         rsum = 0.0;
  1837.         gsum = 0.0;
  1838.         bsum = 0.0;
  1839.         leftcol = col - ccolso2;
  1840.         for ( ccol = 0; ccol < ccols; ++ccol )
  1841.             {
  1842.             tempcol = leftcol + ccol;
  1843.             rcolumnsum[tempcol] = rcolumnsum[tempcol] 
  1844.             - PPM_GETR( rowptr[subrow][ccol] )
  1845.             + PPM_GETR( rowptr[addrow][ccol] );
  1846.             rsum = rsum + rcolumnsum[tempcol] * rweights[0][ccol];
  1847.             gcolumnsum[tempcol] = gcolumnsum[tempcol] 
  1848.             - PPM_GETG( rowptr[subrow][ccol] )
  1849.             + PPM_GETG( rowptr[addrow][ccol] );
  1850.             gsum = gsum + gcolumnsum[tempcol] * gweights[0][ccol];
  1851.             bcolumnsum[tempcol] = bcolumnsum[tempcol] 
  1852.             - PPM_GETB( rowptr[subrow][ccol] )
  1853.             + PPM_GETB( rowptr[addrow][ccol] );
  1854.             bsum = bsum + bcolumnsum[tempcol] * bweights[0][ccol];
  1855.             }
  1856.         temprsum = rsum + 0.5;
  1857.         tempgsum = gsum + 0.5;
  1858.         tempbsum = bsum + 0.5;
  1859.         CHECK_RED;
  1860.         CHECK_GREEN;
  1861.         CHECK_BLUE;
  1862.         PPM_ASSIGN( outputrow[col], r, g, b );
  1863.         }
  1864.         else
  1865.         {
  1866.         rsum = 0.0;
  1867.         gsum = 0.0;
  1868.         bsum = 0.0;
  1869.         leftcol = col - ccolso2;
  1870.         addcol = col + ccolso2;
  1871.         rcolumnsum[addcol] = rcolumnsum[addcol]
  1872.             - PPM_GETR( rowptr[subrow][addcol] )
  1873.             + PPM_GETR( rowptr[addrow][addcol] );
  1874.         gcolumnsum[addcol] = gcolumnsum[addcol]
  1875.             - PPM_GETG( rowptr[subrow][addcol] )
  1876.             + PPM_GETG( rowptr[addrow][addcol] );
  1877.         bcolumnsum[addcol] = bcolumnsum[addcol]
  1878.             - PPM_GETB( rowptr[subrow][addcol] )
  1879.             + PPM_GETB( rowptr[addrow][addcol] );
  1880.         for ( ccol = 0; ccol < ccols; ++ccol )
  1881.             {
  1882.             rsum += rcolumnsum[leftcol + ccol] * rweights[0][ccol];
  1883.             gsum += gcolumnsum[leftcol + ccol] * gweights[0][ccol];
  1884.             bsum += bcolumnsum[leftcol + ccol] * bweights[0][ccol];
  1885.             }
  1886.         temprsum = rsum + 0.5;
  1887.         tempgsum = gsum + 0.5;
  1888.         tempbsum = bsum + 0.5;
  1889.         CHECK_RED;
  1890.         CHECK_GREEN;
  1891.         CHECK_BLUE;
  1892.         PPM_ASSIGN( outputrow[col], r, g, b );
  1893.         }
  1894.         }
  1895.     pnm_writepnmrow( stdout, outputrow, cols, maxval, newformat, 0 );
  1896.     }
  1897.  
  1898.     /* Now write out the remaining unconvolved rows in xelbuf. */
  1899.     for ( irow = crowso2 + 1; irow < crows; ++irow )
  1900.     pnm_writepnmrow(
  1901.             stdout, rowptr[irow], cols, maxval, newformat, 0 );
  1902.  
  1903.     pm_close( stdout );
  1904.     }
  1905.