Article 13644 of alt.sources:

Archive-name: pgmmedian.c_1.0
Submitted-by: burns@cac.psu.edu

Description
-----------
pgmmedian is a pbmplus/netpbm program that applies a median filter to
an image to reduce noise such as impulse noise.  The following shar
file includes the source code, man page and README file for pgmmedian.

- Mike

----------------------------------------------------------------------------
Mike Burns                                                UNIX Systems Group 
burns@cac.psu.edu                              Center for Academic Computing 
(814) 863-5606                             The Pennsylvania State University

---- Cut Here and unpack ----
#!/bin/sh
# shar:	Shell Archiver  (v1.22)
#
#	Run the following text with /bin/sh to create:
#	  pgmmedian.README
#	  pgmmedian.1
#	  pgmmedian.c
#
sed 's/^X//' << 'SHAR_EOF' > pgmmedian.README &&
XREADME for pgmmedian 1.0
X
X
XDescription
X-----------
Xpgmmedian is a pbmplus/netpbm program that applies a median filter to
Xan image to reduce noise such as impulse noise.  It does not try to
Xcharacterize the noise to perform a better noise reduction.
X
X
XInstallation
X------------
XCopy C source and man page to your {pbmplus|netpbm}/pgm directory and then
Xmodify the Imakefile to include pgmmedian in MAN1, SRCS, OBJS, BINS and
Xadd the target line
X
XNormalPbmplusProgramTarget(pgmmedian)
X
XThen run xmkmf and make.
X
X
XInformation about the algorithms
X--------------------------------
XWhile implementing this median filter, I tried many types of sorts and select
Xk'th value routines to find the median value.  I eventually chose to 
Ximplement it with a histogram sort and a select as the histogram sort is 
Xfastest method in many cases especially with larger mask sizes and when
Xthe histogram sort is not the fastest, select was nearly the fastest and
Xeasier to implement than the other method I created.  References for both
Xof these are given in the man page and source code.
X
XAll regular sorting methods (quicksort, insertion, shell, merge) were too
Xslow.  I tried the select k'th from Sedgewick's book and the CACM 489 one
Xand found 489 to be faster.  I also tried the algorithm "Median Finding
Xon a 3 x 3 Grid" from Graphics Gems.  The histogram sort still beat it and
Xbeat the 5x5 version given in the appendix.  For combinations of small mask 
Xsizes and small maxvals the histogram sort is always fastest, but there is
Xa point where it becomes much slower as the maxval gets larger because
Xyou have to compare bins that are farther and farther apart on average.  
XThat is why I use a cutoff value to decide between using a histogram sort
Xor select.  See the timing tables below to see how poorly a histogram sort
Xcan perform for a 3x3 mask with maxval of 64k.  
X
XI modified the histogram sort to them use a linked list of bins rather than
Xallocating a column of 64k bins.  This way you only have to go to the next
Xbin in the linked list, rather than possibly go through 10,000 bins to find
Xthe median.  The added overhead of having to insert and delete bins from the
Xlinked list was a little steep and only saves time over the standard 
Xhistogram sort when the maxval is 64k.  And when the mask size is small, 
Xselect still beats it, so I decided not to also implement it in the final 
Xversion.
X
XI also designed my own method of finding the median that sorts each column,
Xthen merge sorts the columns together to find the median.  Then as the mask
Xmoves horizontally across the image, you already have the columes sorted.
XYou have to remove the proper elements from the left most column and merge
Xin the sorted values from the right most column.  It used a shell sort for 
Xthe sorting and binary searches to find and replace values from an array 
Xof pointers, etc.  It's nearly always 5-20% faster than select, but
Xquite a bit more messy to debug when something goes wrong, so I decided to
Xstick with select for those cases where it can beat the histogram sort.
XThis method and the histogram sort were the only methods I used that reused
Xthe data from the previous mask.
X
XThe following timings are from a Sun IPC and Sparc 20.  The row of numbers
Xacross the top is the result of (maxval/((cols X rows) -1) ).  The second
Xrow of numbers is the maxval of the image used for that run.
X
Xselect - select from CACM #489
Xmysort - the sort I devised described above
Xhist   - histogram sort
Xhist2  - histogram sort using linked lists
X
X
XImage size of 512x512 
X
XSun Sparc IPC
X-------------
X
XMask size of 3x3
X----------------
X   8               128    256    512    1024   2048   4096    8192
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e    13.5    13.8   15.4   16.2    16.7   17.9   19.1    19.7 
XVer19m    11.7    12.0   13.5   14.3    14.9   16.1   17.3    18.1
XVer20j    14.2    14.8   16.4   18.3    20.0   24.0   31.0    42.8
X
X
XMask size of 5x5
X----------------
X  24                      85     170     341    682   1365    2730
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e    21.3    21.8   23.2   24.1    24.7   25.7   26.8    27.6 
XVer19m    18.1    18.5   20.0   20.9    21.4   22.7   23.9    24.4
XVer20b    14.0    15.5   18.9   23.9    32.8   51.4   89.6   171.0
XVer20j    32.4    33.0   34.9   36.5    38.4   42.6   49.2    60.7 
X
X
XMask size of 9x9
X----------------
X  80                      25     51      102    204    409     819
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e    45.2    45.6   47.2   48.1    48.6   49.4   50.8    51.4 
XVer19m    41.6    42.0   43.4   44.3    44.8   46.1   47.4    47.8
XVer20b    14.0    15.5   18.9   23.9    32.8   51.4   89.6   171.0
XVer20j    32.4    33.0   34.9   36.5    38.4   42.6   49.2    60.7
X
X
XMask size of 13x13
X------------------
X 168                      12     24      48     97     195     390
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e    81.4    81.5   83.3   84.3    84.5   85.2   86.8    95.2
XVer19m    78.7    79.1   80.4   81.4    82.0   83.4   84.2    84.5
XVer20b    16.9    18.2   21.0   25.8    33.4   49.6   82.4   150.5
XVer20j    44.4    46.1   47.2   48.5    50.4   54.3   60.8    72.0 
X
X
XMask size of 17x17
X------------------
X 288                       7     14      28     56     113     227
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e   129.7   129.9  131.8  132.5   133.0  133.2  134.6   146.7
XVer19m   133.4   133.8  135.1  136.1   136.5  137.9  139.0   138.9
XVer20b    19.8    20.9   23.8   28.1    35.0   49.3   78.8   138.9
XVer20j    55.7    56.1   57.1   58.0    59.6   66.7   72.9    84.0 
X
X
XMask size of 25x25
X------------------
X 624                       3      6      13     26     52      105
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e   269.0   269.9  270.7  272.8   273.2  271.9  273.5   297.6
XVer19m   271.1   271.5  273.0  274.1   274.6  275.7  276.8   276.2
XVer20b    27.3    28.3   30.9   34.6    40.2   52.7   77.7   127.5
XVer20j    65.7    66.4   68.1   69.5    71.3   75.9   82.0    92.5 
X
X
XMask size of 35x35
X------------------
X1224                              3       6     13     26      53
X           511     1k     2k     4k      8k     16k    32k     64k
XVer14e   511.7   512.0  514.1  512.3   514.9  514.1  514.8   559.2
XVer19m   515.7   516.1  517.7  518.5   519.0  520.3  521.5   522.4
XVer20b    34.7    35.6   38.2   41.0    46.1   56.9   77.2   117.9
XVer20j    72.0    72.6   74.2   75.7    77.9   82.3   88.2    98.7 
X
X
XSun Sparc 20 model 71
X---------------------
X
XMask size of 3x3
X----------------
X   8               128    256    512   1024   2048   4096   8192
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect      2.8    2.9    3.3    3.6    3.6    3.9    4.2    4.5
Xmysort      2.4    2.4    2.7    2.9    2.9    3.2    3.4    3.7
Xhist        1.7    2.1    2.8    4.0    6.2   10.7   19.7   37.7
Xhist2       2.5    2.6    3.0    3.2    3.4    4.0    4.8    6.0
X
X
XMask size of 5x5
X----------------
X  24                      85     170    341    682   1365   2730
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect      4.6    4.7    5.0    5.4    5.4    5.7    6.0    6.1
Xmysort      3.7    3.8    4.1    4.4    4.5    4.7    4.9    5.1
Xhist        1.9    2.2    2.8    3.9    5.7    9.3   16.6   31.0
Xhist2       3.7    3.7    4.1    4.3    4.5    5.0    5.7    7.0
X
X
XMask size of 9x9
X----------------
X  80                      25     51     102    204    409    819
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect     10.4   10.5   11.0   11.1   11.3   11.6   11.9   12.1
Xmysort      8.9    8.9    9.3    9.5    9.5    9.8   10.0   10.2
Xhist        2.2    2.5    3.1    3.9    5.4    8.4   14.0   25.4
Xhist2       6.0    6.1    6.4    6.6    6.9    7.4    8.1    9.5
X
X
X
XMask size of 13x13
X------------------
X 168                      12     24     48     97     195    390
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect     20.0   19.8   20.2   20.8   20.5   20.7   21.2   21.5
Xmysort     17.1   17.1   17.5   17.7   17.7   18.1   18.2   18.3
Xhist        2.7    2.9    3.5    4.2    5.5    8.2   13.0   22.4
Xhist2       8.2    8.3    8.7    9.0    9.1    9.7   10.6   11.8
X
X
XMask size of 17x17
X------------------
X 288                       7     14     28     56     113    227
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect     32.6   32.7   33.2   33.7   33.4   33.7   34.1   34.4
Xmysort     28.5   28.5   30.0   29.0   29.1   29.3   29.6   29.6
Xhist        3.4    3.6    4.1    4.7    5.8    8.2   12.8   20.9
Xhist2      10.3   10.4   12.2   11.0   11.3   11.7   12.5   13.6
X
X
XMask size of 25x25
X------------------
X 624                       3      6     13     26     52     105
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect     69.6   69.7   70.5   70.4   70.8   70.8   70.9   71.2
Xmysort     60.2   60.2   60.8   60.7   60.7   61.0   61.3   61.4
Xhist        5.7    5.5    6.0    6.8    7.5    9.6   13.8   20.1
Xhist2      13.4   13.5   13.7   14.6   14.3   14.8   15.6   16.7
X
X
XMask size of 35x35
X------------------
X1224                              3      6     13     26      53
X            511    1k     2k     4k     8k     16k    32k    64k
Xselect    137.7  136.6  138.6  137.3  137.9  137.9  138.2  138.5
Xmysort    230.0  114.9  115.6  115.4  115.7  115.8  116.0  116.1
Xhist        7.3    7.5    7.9    9.8    9.2   10.7   14.0   19.3
Xhist2      15.1   15.3   17.5   15.8   19.1   16.4   17.2   18.6
X
XSend all questions/comments/bugs to me at burns@cac.psu.edu.
X
X- Mike
X
X----------------------------------------------------------------------------
XMike Burns                                                UNIX Systems Group 
Xburns@cac.psu.edu                              Center for Academic Computing 
X(814) 863-5606                             The Pennsylvania State University
X
X
SHAR_EOF
chmod 0644 pgmmedian.README || echo "restore of pgmmedian.README fails"
sed 's/^X//' << 'SHAR_EOF' > pgmmedian.1 &&
X.TH pgmmedian 1 "28 September 1996"
X.IX pgmmedian
X.SH NAME
Xpgmmedian - apply a median filter to a PGM file
X.SH SYNOPSIS
X.B pgmmedian
X.RB [ -cutoff
X.IR int ]
X.RB [ -size
X.IR width
X.IR height ]
X.RB [ -type
X.IR median_type ]
X.RI [ pnmfile ]
X.br
X.SH DESCRIPTION
XApplies a median filter to a PGM image.  The program automatically selects
Xwhether to use a histogram sort or select k'th value method to determine
Xmedian.  This automatic selection is configurable via the -cutoff option and
Xwhich method to use can be forced with the -type option.
X.SH OPTIONS
X.TP
X.B -cutoff int
XThe cutoff value to decide between using the histogram sort or select k'th 
Xvalue to find the median.  It is compared to (maxval / ((width * height) - 1)).
XThe default value is
X.IR 250.
XIf expression is less than the cutoff the histogram sort is used, otherwise
Xselect k'th value method is used.
X.TP
X.B -size width height
XWidth and height of median mask to apply.
X.TP
X.B -type median_type
XSelects which method to use to find median irregardless of cutoff value.  
XChoices currently are histogram_sort and select.
XOnly enough letters to be unique are required for the type option.
X.SH "BUGS"
XNone.
X.SH "REFERENCES"
X"Collected Algorithms from ACM" Volume II, Algorithm 489 by Robert W. Floyd
X.TP
X"A Fast Two\-Dimensional Median Filtering Algorithm" in "IEEE Transactions
Xon Acoustics, Speech, and Signal Processing" Vol. ASSP-27, No. 1, February
X1979
X.TP
X"Digital Image Processing Algorithms" by Ioannis Pitas,
XPrentice Hall, 1993  ISBN 0-13-145814-0
X.SH "SEE ALSO"
Xpgm(5), pgmnoise(1), pnmaddnoise(1), pnmsmooth(1)
X.SH AUTHOR
XCopyright (C) 1996 by Mike Burns <burns@cac.psu.edu>
X.\" Permission to use, copy, modify, and distribute this software and its
X.\" documentation for any purpose and without fee is hereby granted, provided
X.\" that the above copyright notice appear in all copies and that both that
X.\" copyright notice and this permission notice appear in supporting
X.\" documentation.  This software is provided "as is" without express or
X.\" implied warranty.
SHAR_EOF
chmod 0644 pgmmedian.1 || echo "restore of pgmmedian.1 fails"
sed 's/^X//' << 'SHAR_EOF' > pgmmedian.c &&
X/* pgmmedian.c - MxN median filter for PGM files.
X**
X** Version 1.0  September 28, 1996
X**
X** Copyright (C) 1996 by Mike Burns <burns@cac.psu.edu>
X**
X** Permission to use, copy, modify, and distribute this software and its
X** documentation for any purpose and without fee is hereby granted, provided
X** that the above copyright notice appear in all copies and that both that
X** copyright notice and this permission notice appear in supporting
X** documentation.  This software is provided "as is" without express or
X** implied warranty.
X*/
X
X/* References
X** ----------
X** The select k'th value implementation is based on Algorithm 489 by 
X** Robert W. Floyd from the "Collected Algorithms from ACM" Volume II.
X**
X** The histogram sort is based is described in the paper "A Fast Two-
X** Dimensional Median Filtering Algorithm" in "IEEE Transactions on 
X** Acoustics, Speech, and Signal Processing" Vol. ASSP-27, No. 1, February
X** 1979.  The algorithm I more closely followed is found in "Digital
X** Image Processing Algorithms" by Ioannis Pitas.
X*/
X
X
X#include "pgm.h"
X
X#define MAX_MEDIAN_TYPES	2
X#define SELECT_MEDIAN		0
X#define HISTOGRAM_SORT_MEDIAN	1
X
Xstatic void select_median ARGS(( FILE *ifp, int ccols, int crows,
X    int cols, int rows, int median ));
Xstatic void select_489 ARGS(( gray *a, int *parray, int n, int k ));
Xstatic void histogram_sort_median ARGS(( FILE *ifp, int ccols, int crows,
X    int cols, int rows, int median ));
X
X/* Global variables common to each median sort routine. */
Xint format, forceplain;
Xgray maxval;
Xgray **grays;
Xgray *grayrow;
Xgray **rowptr;
Xint ccolso2, crowso2;
Xint row;
X
Xint
Xmain( argc, argv )
X    int argc;
X    char *argv[];
X    {
X    FILE *ifp;
X    int argn;
X    int ccols, crows;
X    int median_type;
X    int i;
X    int median_cutoff;	
X    int cols, rows;
X    int median;
X    char *usage = "[-size width height] [-type median_type] [-cutoff value] [pgmfile]";
X
X    int median_type_selected;
X    char *median_name[] =
X	{
X	"select",
X	"histogram_sort",
X	};
X    int median_id[] = 
X	{
X	SELECT_MEDIAN,
X	HISTOGRAM_SORT_MEDIAN
X	};
X    int median_compare[] =
X	{
X	1,
X	2
X	};
X
X    pgm_init( &argc, argv );
X
X    ccols = 3;
X    crows = 3;
X    /* Cutoff to choose between select and histogram median. */
X    median_cutoff = 250;
X    median_type = SELECT_MEDIAN;
X    median_type_selected = 0;
X
X    argn = 1;
X    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
X	{
X	if ( pm_keymatch( argv[argn], "-cutoff", 2 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( "incorrect number of arguments for -cutoff option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -cutoff option: %s",
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    if ( (median_cutoff = atoi(argv[argn])) == 0 )
X		pm_error( "invalid cutoff specification: %s", argv[argn] );
X	    }
X	else if ( pm_keymatch( argv[argn], "-size", 2 ) )
X	    {
X	    ++argn;
X	    if ( argn+1 >= argc )
X		{
X		pm_message( "incorrect number of arguments for -size option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' || argv[argn+1][0] == '-' )
X		{
X		pm_message( "invalid arguments to -size option: %s %s",
X		    argv[argn], argv[argn+1] );
X		pm_usage( usage );
X		}
X	    if ( (ccols = atoi(argv[argn])) == 0 )
X		pm_error( "invalid width size specification: %s", argv[argn] );
X	    ++argn;
X	    if ( (crows = atoi(argv[argn])) == 0 )
X		pm_error( "invalid height size specification: %s",argv[argn] );
X	    if ( ccols % 2 != 1 || crows % 2 != 1 )
X		pm_error( "the convolution matrix must have an odd number of rows and columns" );
X	    }
X	else if ( pm_keymatch( argv[argn], "-type", 2 ) )
X	    {
X	    ++argn;
X	    if ( argn >= argc )
X		{
X		pm_message( "incorrect number of arguments for -type option" );
X		pm_usage( usage );
X		}
X	    else if ( argv[argn][0] == '-' )
X		{
X		pm_message( "invalid argument to -type option: %s",
X		    argv[argn] );
X		pm_usage( usage );
X		}
X	    i = 0;
X	    while ( ( i < MAX_MEDIAN_TYPES ) && !pm_keymatch( argv[argn],
X		median_name[i], median_compare[i] ) )
X		++i;
X	    if ( i == MAX_MEDIAN_TYPES )
X		{
X		pm_message( "invalid argument to -type option: %s",
X		argv[argn] );
X		pm_usage( usage );
X		}
X	    median_type = median_id[i];
X	    median_type_selected = 1;
X	    }
X	else
X	    pm_usage( usage );
X	++argn;
X	}
X
X    if ( argn < argc )
X	{
X	ifp = pm_openr( argv[argn] );
X	++argn;
X	}
X    else
X	ifp = stdin;
X
X    if ( argn != argc )
X	pm_usage( usage );
X
X    /* Code common to all sorts methods */
X    forceplain = 0;
X
X    ccolso2 = ccols / 2;
X    crowso2 = crows / 2;
X
X    pgm_readpgminit( ifp, &cols, &rows, &maxval, &format );
X    pgm_writepgminit( stdout, cols, rows, maxval, forceplain );
X
X    /* Allocate space for number of rows in mask size. */
X    grays = pgm_allocarray( cols, crows );
X    grayrow = pgm_allocrow( cols );
X
X    /* Allocate pointers to mask row buffer. */
X    rowptr = pgm_allocarray( 1, crows );
X
X    /* Read in and write out initial rows that won't get changed. */
X    for ( row = 0; row < crows - 1; ++row )
X	{
X	pgm_readpgmrow( ifp, grays[row], cols, maxval, format );
X	/* Write out the unchanged row. */
X        if ( row < crowso2 )
X	    pgm_writepgmrow( stdout, grays[row], cols, maxval, forceplain );
X	}
X
X    median = ( crows * ccols ) / 2;
X
X    /* Choose which sort to run. */
X    if ( ! median_type_selected ) {
X	if ( ( maxval / ( (ccols * crows) - 1 ) ) < median_cutoff )
X	    median_type = HISTOGRAM_SORT_MEDIAN;
X	else
X	    median_type = SELECT_MEDIAN;
X    }
X
X    switch ( median_type )
X	{
X	case SELECT_MEDIAN:
X	select_median( ifp, ccols, crows, cols, rows, median );
X	break;
X
X	case HISTOGRAM_SORT_MEDIAN:
X	histogram_sort_median( ifp, ccols, crows, cols, rows, median );
X	break;
X
X	default:
X	pm_error( "no median type" );
X	break;
X	}
X	
X    pm_close( ifp );
X    pm_close( stdout );
X
X    pgm_freearray( grays, crows );
X    pgm_freerow( grayrow );
X    pgm_freearray( rowptr, crows );
X
X    exit( 0 );
X    }
X
X
Xstatic void
Xselect_489( gray *a, int *parray, int n, int k )
X    {
X    gray t, v;
X    int i, j, l, r;
X    int ptmp, ttmp;
X
X    l = 0;
X    r = n - 1;
X    while ( r > l )
X	{
X	t = a[parray[k]];
X	ttmp = parray[k];
X	i = l;
X	j = r;
X	ptmp = parray[l];
X	parray[l] = parray[k];
X	parray[k] = ptmp;
X	if ( a[parray[r]] > t )
X	    {
X	    ptmp = parray[r];
X	    parray[r] = parray[l];
X	    parray[l] = ptmp;
X	    }
X	while ( i < j )
X	    {
X	    ptmp = parray[i];
X	    parray[i] = parray[j];
X	    parray[j] = ptmp;
X	    ++i;
X	    --j;
X	    while ( a[parray[i]] < t )
X		++i;
X	    while ( a[parray[j]] > t )
X		--j;
X	    }
X	if ( a[parray[l]] == t )
X	    {
X	    ptmp = parray[l];
X	    parray[l] = parray[j];
X	    parray[j] = ptmp;
X	    }
X	else
X	    {
X	    ++j;
X	    ptmp = parray[j];
X	    parray[j] = parray[r];
X	    parray[r] = ptmp;
X	    }
X	if ( j <= k )
X	    l = j + 1;
X	if ( k <= j )
X	    r = j - 1;
X	}
X    }
X
X
X
Xstatic void
Xselect_median( FILE *ifp, int ccols, int crows, int cols, int rows, int median )
X    {
X    register int ccol, col;
X    int crow;
X    gray g;
X    int rownum, irow, temprow;
X    gray *temprptr;
X    int i, leftcol;
X    int num_values;
X    gray *garray;
X
X    int *parray;
X    int addcol;
X    int *subcol;
X    int tsum;
X
X    /* Allocate storage for array of the current gray values. */
X    garray = pgm_allocrow( crows * ccols );
X
X    num_values = crows * ccols;
X
X    parray = (int *) pm_allocrow( crows * ccols, sizeof(int) );
X    subcol = (int *) pm_allocrow( cols, sizeof(int) );
X
X    for ( i = 0; i < cols; ++i )
X	subcol[i] = ( i - (ccolso2 + 1) ) % ccols;
X
X    /* Apply median to main part of image. */
X    for ( ; row < rows; ++row )
X	{
X	temprow = row % crows;
X	pgm_readpgmrow( ifp, grays[temprow], cols, maxval, format );
X
X        /* Rotate pointers to rows, so rows can be accessed in order. */
X        temprow = ( row + 1 ) % crows;
X	rownum = 0;
X	for ( irow = temprow; irow < crows; ++rownum, ++irow )
X	    rowptr[rownum] = grays[irow];
X	for ( irow = 0; irow < temprow; ++rownum, ++irow )
X	    rowptr[rownum] = grays[irow];
X
X	for ( col = 0; col < cols; ++col )
X	    {
X	    if ( col < ccolso2 || col >= cols - ccolso2 )
X		{
X		grayrow[col] = rowptr[crowso2][col];
X		}
X	    else if ( col == ccolso2 )
X		{
X		leftcol = col - ccolso2;
X		i = 0;
X		for ( crow = 0; crow < crows; ++crow )
X		    {
X		    temprptr = rowptr[crow] + leftcol;
X		    for ( ccol = 0; ccol < ccols; ++ccol )
X			{
X			garray[i] = *( temprptr + ccol );
X			parray[i] = i;
X			++i;
X			}
X		    }
X		select_489( garray, parray, num_values, median );
X		grayrow[col] = garray[parray[median]];
X		}
X	    else
X		{
X		addcol = col + ccolso2;
X		for ( crow = 0, tsum = 0; crow < crows; ++crow, tsum += ccols )
X		    {
X		    garray[tsum + subcol[col]] =
X			*(rowptr[crow] + addcol );
X		    }
X		select_489( garray, parray, num_values, median );
X		grayrow[col] = garray[parray[median]];
X		}
X	    }
X	pgm_writepgmrow( stdout, grayrow, cols, maxval, forceplain );
X	}
X
X    /* Write out remaining unchanged rows. */
X    for ( irow = crowso2 + 1; irow < crows; ++irow )
X	pgm_writepgmrow( stdout, rowptr[irow], cols, maxval, forceplain );
X
X    pgm_freerow( garray );
X    pm_freerow( (char *) parray );
X    pm_freerow( (char *) subcol );
X    }
X
X
X
Xstatic void
Xhistogram_sort_median( FILE *ifp, int ccols, int crows, int cols, int rows,
X    int median )
X    {
X    register int ccol, col;
X    int crow;
X    gray g;
X    int rownum, irow, temprow;
X    gray *temprptr;
X    int i, leftcol;
X
X    int *hist, histmax;
X    int mdn, ltmdn;
X    gray *left_col, *right_col;
X
X    int addcol, subcol;
X
X    histmax = maxval + 1;
X    hist = (int *) pm_allocrow( histmax, sizeof( int ) );
X    left_col = pgm_allocrow( crows );
X    right_col = pgm_allocrow( crows );
X
X    /* Apply median to main part of image. */
X    for ( ; row < rows; ++row )
X	{
X	/* initialize hist[] */
X	for ( i = 0; i < histmax; ++i )
X	    hist[i] = 0;
X
X	temprow = row % crows;
X	pgm_readpgmrow( ifp, grays[temprow], cols, maxval, format );
X
X        /* Rotate pointers to rows, so rows can be accessed in order. */
X        temprow = ( row + 1 ) % crows;
X	rownum = 0;
X	for ( irow = temprow; irow < crows; ++rownum, ++irow )
X	    rowptr[rownum] = grays[irow];
X	for ( irow = 0; irow < temprow; ++rownum, ++irow )
X	    rowptr[rownum] = grays[irow];
X
X	for ( col = 0; col < cols; ++col )
X	    {
X	    if ( col < ccolso2 || col >= cols - ccolso2 )
X		{
X		grayrow[col] = rowptr[crowso2][col];
X		}
X	    else if ( col == ccolso2 )
X		{
X		leftcol = col - ccolso2;
X		i = 0;
X		for ( crow = 0; crow < crows; ++crow )
X		    {
X		    temprptr = rowptr[crow] + leftcol;
X		    for ( ccol = 0; ccol < ccols; ++ccol )
X			{
X			g = *( temprptr + ccol );
X			++hist[g];
X			++i;
X			}
X		    }
X		ltmdn = 0;
X		for ( mdn = 0; ltmdn <= median; ++mdn )
X		    ltmdn += hist[mdn];
X		mdn--;
X		if ( ltmdn > median ) 
X		    ltmdn -= hist[mdn];
X
X		grayrow[col] = mdn;
X		}
X	    else 
X		{
X		subcol = col - ( ccolso2 + 1 );
X		addcol = col + ccolso2;
X		for ( crow = 0; crow < crows; ++crow )
X		    {
X		    left_col[crow] = *( rowptr[crow] + subcol );
X		    right_col[crow] = *( rowptr[crow] + addcol );
X		    }
X		for ( crow = 0; crow < crows; ++crow )
X		    {
X		    g = left_col[crow];
X		    hist[(int) g]--;
X		    if ( (int) g < mdn )
X			ltmdn--;
X		    g = right_col[crow];
X		    hist[(int) g]++;
X		    if ( (int) g < mdn )
X			ltmdn++;
X		    }
X
X		if ( ltmdn > median )
X		    do
X			{
X			mdn--;
X			ltmdn -= hist[mdn];
X			} while ( ltmdn > median );
X		else
X		    {
X		    /* This one change from Pitas algorithm can reduce run
X		    ** time by up to 10%.
X		    */
X		    while ( ltmdn <= median )
X			{
X			ltmdn += hist[mdn];
X			mdn++;
X			}
X		    mdn--;
X		    if ( ltmdn > median ) 
X			ltmdn -= hist[mdn];
X		    }
X		grayrow[col] = mdn;
X	
X		}
X	    }
X	pgm_writepgmrow( stdout, grayrow, cols, maxval, forceplain );
X	}
X
X    /* Write out remaining unchanged rows. */
X    for ( irow = crowso2 + 1; irow < crows; ++irow )
X	pgm_writepgmrow( stdout, rowptr[irow], cols, maxval, forceplain );
X
X    pm_freerow( (char *) hist );
X    pgm_freerow( left_col );
X    pgm_freerow( right_col );
X
X    }
SHAR_EOF
chmod 0644 pgmmedian.c || echo "restore of pgmmedian.c fails"
exit 0


