home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsm / netpbmsca / pnm / c / pnmnlfilt < prev    next >
Encoding:
Text File  |  1993-10-04  |  35.3 KB  |  854 lines
  1. /* pnmnlfilt.c - 4 in 1 (2 non-linear) filter
  2. **             - smooth an anyimage
  3. **             - do alpha trimmed mean filtering on an anyimage
  4. **             - do optimal estimation smoothing on an anyimage
  5. **             - do edge enhancement on an anyimage
  6. **
  7. ** Version 1.0
  8. **
  9. ** The implementation of an alpha-trimmed mean filter
  10. ** is based on the description in IEEE CG&A May 1990
  11. ** Page 23 by Mark E. Lee and Richard A. Redner.
  12. **
  13. ** The paper recommends using a hexagon sampling region around each
  14. ** pixel being processed, allowing an effective sub pixel radius to be
  15. ** specified. The hexagon values are sythesised by area sampling the
  16. ** rectangular pixels with a hexagon grid. The seven hexagon values
  17. ** obtained from the 3x3 pixel grid are used to compute the alpha
  18. ** trimmed mean. Note that an alpha value of 0.0 gives a conventional
  19. ** mean filter (where the radius controls the contribution of
  20. ** surrounding pixels), while a value of 0.5 gives a median filter.
  21. ** Although there are only seven values to trim from before finding
  22. ** the mean, the algorithm has been extended from that described in
  23. ** CG&A by using interpolation, to allow a continuous selection of
  24. ** alpha value between and including 0.0 to 0.5  The useful values
  25. ** for radius are between 0.3333333 (where the filter will have no
  26. ** effect because only one pixel is sampled), to 1.0, where all
  27. ** pixels in the 3x3 grid are sampled.
  28. **
  29. ** The optimal estimation filter is taken from an article "Converting Dithered
  30. ** Images Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November
  31. ** 1992, and this article references "Digital Image Enhancement andNoise Filtering by
  32. ** Use of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and
  33. ** Machine Intelligence, March 1980.
  34. **
  35. ** Also borrow the  technique used in pgmenhance(1) to allow edge
  36. ** enhancement if the alpha value is negative.
  37. **
  38. ** Author:
  39. **         Graeme W. Gill, 30th Jan 1993
  40. **         graeme@labtam.oz.au
  41. */
  42.  
  43. #include <math.h>
  44. #include "pnm.h"
  45.  
  46. double hex_area ARGS((double, double, double, double, double));
  47. int atfilt_setup ARGS((double, double, double));
  48. int atfilt0 ARGS((int *)), atfilt1 ARGS((int *)), atfilt2 ARGS((int *));
  49. int atfilt3 ARGS((int *)), atfilt4 ARGS((int *)), atfilt5 ARGS((int *));
  50. int (*atfuncs[6]) ARGS((int *)) = {atfilt0,atfilt1,atfilt2,atfilt3,atfilt4,atfilt5};
  51.  
  52. xelval omaxval; /* global so that pixel processing code can get at it quickly */
  53. int noisevariance;      /* global so that pixel processing code can get at it quickly */
  54.  
  55. int
  56. main( argc, argv )
  57. int argc;
  58. char* argv[];
  59.         {
  60.         double radius=0.0,alpha= -1.0;
  61.         FILE* ifp;
  62.         int rows, cols, format, oformat, row, col;
  63.         xelval maxval;
  64.         int (*atfunc) ARGS((int *));
  65.         char* usage = "alpha radius pnmfile\n\
  66.  0.0 <= alpha <= 0.5 for alpha trimmed mean -or- \n\
  67.  1.0 <= alpha <= 2.0 for optimal estimation -or- \n\
  68.  -0.1 >= alpha >= -0.9 for edge enhancement\n\
  69.  0.3333 <= radius <= 1.0 specify effective radius\n";
  70.  
  71.  
  72.         pnm_init( &argc, argv );
  73.  
  74.         if ( argc < 3 || argc > 4 )
  75.                 pm_usage( usage );
  76.  
  77.         if ( sscanf( argv[1], "%lf", &alpha ) != 1 )
  78.                 pm_usage( usage );
  79.         if ( sscanf( argv[2], "%lf", &radius ) != 1 )
  80.                 pm_usage( usage );
  81.  
  82.         if ((alpha > -0.1 && alpha < 0.0) || (alpha > 0.5 && alpha < 1.0))
  83.                 pm_error( "Alpha must be in range 0.0 <= alpha <= 0.5 for alpha trimmed mean" );
  84.         if (alpha > 2.0)
  85.                 pm_error( "Alpha must be in range 1.0 <= alpha <= 2.0 for optimal estimation" );
  86.         if (alpha < -0.9 || (alpha > -0.1 && alpha < 0.0))
  87.                 pm_error( "Alpha must be in range -0.9 <= alpha <= -0.1 for edge enhancement" );
  88.         if (radius < 0.333 || radius > 1.0)
  89.                 pm_error( "Radius must be in range 0.333333333 <= radius <= 1.0" );
  90.  
  91.         if ( argc == 4 )
  92.                 ifp = pm_openr( argv[3] );
  93.         else
  94.                 ifp = stdin;
  95.  
  96.         pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
  97.  
  98.         oformat = PNM_FORMAT_TYPE(format);
  99.         omaxval = PPM_MAXMAXVAL;        /* force output to max precision */
  100.  
  101.         atfunc = atfuncs[atfilt_setup(alpha,radius,(double)omaxval/(double)maxval)];
  102.  
  103.         if ( oformat < PGM_TYPE )
  104.                 {
  105.                 if (PGM_MAXMAXVAL > PPM_MAXMAXVAL)
  106.                         pm_error( "Can't handle pgm input file as maxval is too big" );
  107.                 oformat = RPGM_FORMAT;
  108.                 pm_message( "promoting file to PGM" );
  109.                 }
  110.         pnm_writepnminit( stdout, cols, rows, omaxval, oformat, 0 );
  111.  
  112.         if ( PNM_FORMAT_TYPE(oformat) == PPM_TYPE )
  113.                 {
  114.                 xel *irows[3];
  115.                 xel *irow0, *irow1, *irow2, *orow;
  116.                 int pr[9],pg[9],pb[9];          /* 3x3 neighbor pixel values */
  117.                 int r,g,b;
  118.  
  119.                 irows[0] = pnm_allocrow( cols );
  120.                 irows[1] = pnm_allocrow( cols );
  121.                 irows[2] = pnm_allocrow( cols );
  122.                 irow0 = irows[0];
  123.                 irow1 = irows[1];
  124.                 irow2 = irows[2];
  125.                 orow = pnm_allocrow( cols );
  126.  
  127.                 for ( row = 0; row < rows; row++ )
  128.                         {
  129.                         int po,no;              /* offsets for left and right colums in 3x3 */
  130.                         register xel *ip0, *ip1, *ip2, *op;
  131.  
  132.                         if (row == 0)
  133.                                 {
  134.                                 irow0 = irow1;
  135.                                 pnm_readpnmrow( ifp, irow1, cols, maxval, format );
  136.                                 }
  137.                         if (row == (rows-1))
  138.                                 irow2 = irow1;
  139.                         else
  140.                                 pnm_readpnmrow( ifp, irow2, cols, maxval, format );
  141.  
  142.                         for (col = cols-1,po= col>0?1:0,no=0,ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
  143.                              col >= 0;
  144.                              col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0)
  145.                                 {
  146.                                 /* grab 3x3 pixel values */
  147.                                 pr[0] = PPM_GETR( *ip1 );
  148.                                 pg[0] = PPM_GETG( *ip1 );
  149.                                 pb[0] = PPM_GETB( *ip1 );
  150.                                 pr[1] = PPM_GETR( *(ip1-no) );
  151.                                 pg[1] = PPM_GETG( *(ip1-no) );
  152.                                 pb[1] = PPM_GETB( *(ip1-no) );
  153.                                 pr[5] = PPM_GETR( *(ip1+po) );
  154.                                 pg[5] = PPM_GETG( *(ip1+po) );
  155.                                 pb[5] = PPM_GETB( *(ip1+po) );
  156.                                 pr[3] = PPM_GETR( *(ip2) );
  157.                                 pg[3] = PPM_GETG( *(ip2) );
  158.                                 pb[3] = PPM_GETB( *(ip2) );
  159.                                 pr[2] = PPM_GETR( *(ip2-no) );
  160.                                 pg[2] = PPM_GETG( *(ip2-no) );
  161.                                 pb[2] = PPM_GETB( *(ip2-no) );
  162.                                 pr[4] = PPM_GETR( *(ip2+po) );
  163.                                 pg[4] = PPM_GETG( *(ip2+po) );
  164.                                 pb[4] = PPM_GETB( *(ip2+po) );
  165.                                 pr[6] = PPM_GETR( *(ip0+po) );
  166.                                 pg[6] = PPM_GETG( *(ip0+po) );
  167.                                 pb[6] = PPM_GETB( *(ip0+po) );
  168.                                 pr[8] = PPM_GETR( *(ip0-no) );
  169.                                 pg[8] = PPM_GETG( *(ip0-no) );
  170.                                 pb[8] = PPM_GETB( *(ip0-no) );
  171.                                 pr[7] = PPM_GETR( *(ip0) );
  172.                                 pg[7] = PPM_GETG( *(ip0) );
  173.                                 pb[7] = PPM_GETB( *(ip0) );
  174.                                 r = (*atfunc)(pr);
  175.                                 g = (*atfunc)(pg);
  176.                                 b = (*atfunc)(pb);
  177.                                 PPM_ASSIGN( *op, r, g, b );
  178.                                 }
  179.                         pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
  180.                         if (irow1 == irows[2])
  181.                                 {
  182.                                 irow1 = irows[0];
  183.                                 irow2 = irows[1];
  184.                                 irow0 = irows[2];
  185.                                 }
  186.                         else if (irow1 == irows[1])
  187.                                 {
  188.                                 irow2 = irows[0];
  189.                                 irow0 = irows[1];
  190.                                 irow1 = irows[2];
  191.                                 }
  192.                         else    /* must be at irows[0] */
  193.                                 {
  194.                                 irow0 = irows[0];
  195.                                 irow1 = irows[1];
  196.                                 irow2 = irows[2];
  197.                                 }
  198.                         }
  199.                 }
  200.         else    /* Else must be PGM */
  201.                 {
  202.                 xel *irows[3];
  203.                 xel *irow0, *irow1, *irow2, *orow;
  204.                 int p[9];               /* 3x3 neighbor pixel values */
  205.                 int pv;
  206.                 int promote;
  207.  
  208.                 irows[0] = pnm_allocrow( cols );
  209.                 irows[1] = pnm_allocrow( cols );
  210.                 irows[2] = pnm_allocrow( cols );
  211.                 irow0 = irows[0];
  212.                 irow1 = irows[1];
  213.                 irow2 = irows[2];
  214.                 orow = pnm_allocrow( cols );
  215.                 /* we scale maxval to omaxval */
  216.                 promote = ( PNM_FORMAT_TYPE(format) != PNM_FORMAT_TYPE(oformat) );
  217.  
  218.                 for ( row = 0; row < rows; row++ )
  219.                         {
  220.                         int po,no;              /* offsets for left and right colums in 3x3 */
  221.                         register xel *ip0, *ip1, *ip2, *op;
  222.  
  223.                         if (row == 0)
  224.                                 {
  225.                                 irow0 = irow1;
  226.                                 pnm_readpnmrow( ifp, irow1, cols, maxval, format );
  227.                                 if ( promote )
  228.                                         pnm_promoteformatrow( irow1, cols, maxval, format, maxval, oformat );
  229.                                 }
  230.                         if (row == (rows-1))
  231.                                 irow2 = irow1;
  232.                         else
  233.                                 {
  234.                                 pnm_readpnmrow( ifp, irow2, cols, maxval, format );
  235.                                 if ( promote )
  236.                                         pnm_promoteformatrow( irow2, cols, maxval, format, maxval, oformat );
  237.                                 }
  238.  
  239.                         for (col = cols-1,po= col>0?1:0,no=0,ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
  240.                              col >= 0;
  241.                              col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0)
  242.                                 {
  243.                                 /* grab 3x3 pixel values */
  244.                                 p[0] = PNM_GET1( *ip1 );
  245.                                 p[1] = PNM_GET1( *(ip1-no) );
  246.                                 p[5] = PNM_GET1( *(ip1+po) );
  247.                                 p[3] = PNM_GET1( *(ip2) );
  248.                                 p[2] = PNM_GET1( *(ip2-no) );
  249.                                 p[4] = PNM_GET1( *(ip2+po) );
  250.                                 p[6] = PNM_GET1( *(ip0+po) );
  251.                                 p[8] = PNM_GET1( *(ip0-no) );
  252.                                 p[7] = PNM_GET1( *(ip0) );
  253.                                 pv = (*atfunc)(p);
  254.                                 PNM_ASSIGN1( *op, pv );
  255.                                 }
  256.                         pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
  257.                         if (irow1 == irows[2])
  258.                                 {
  259.                                 irow1 = irows[0];
  260.                                 irow2 = irows[1];
  261.                                 irow0 = irows[2];
  262.                                 }
  263.                         else if (irow1 == irows[1])
  264.                                 {
  265.                                 irow2 = irows[0];
  266.                                 irow0 = irows[1];
  267.                                 irow1 = irows[2];
  268.                                 }
  269.                         else    /* must be at irows[0] */
  270.                                 {
  271.                                 irow0 = irows[0];
  272.                                 irow1 = irows[1];
  273.                                 irow2 = irows[2];
  274.                                 }
  275.                         }
  276.                 }
  277.         pm_close( ifp );
  278.  
  279.         exit( 0 );
  280.         }
  281.  
  282. #define MXIVAL PPM_MAXMAXVAL    /* maximum input value */
  283. #define NOIVAL (MXIVAL + 1)             /* number of possible input values */
  284.  
  285. #define SCALEB 8                                /* scale bits */
  286. #define SCALE (1 << SCALEB)     /* scale factor */
  287. #define MXSVAL (MXIVAL * SCALE) /* maximum scaled values */
  288.  
  289. #define CSCALEB 2                               /* coarse scale bits */
  290. #define CSCALE (1 << CSCALEB)   /* coarse scale factor */
  291. #define MXCSVAL (MXIVAL * CSCALE)       /* maximum coarse scaled values */
  292. #define NOCSVAL (MXCSVAL + 1)   /* number of coarse scaled values */
  293. #define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB))  /* convert from scaled to coarse scaled */
  294. #define CSCTOSC(x) ((x) << (SCALEB - CSCALEB))  /* convert from course scaled to scaled */
  295.  
  296. #ifndef MAXINT
  297. # define MAXINT 0x7fffffff      /* assume this is a 32 bit machine */
  298. #endif
  299.  
  300. /* round and scale floating point to scaled integer */
  301. #define ROUND(x) ((int)(((x) * (double)SCALE) + 0.5))
  302. /* round and un-scale scaled integer value */
  303. #define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB)       /* rounded un-scale */
  304. #define UNSCALE(x) ((x) >> SCALEB)
  305.  
  306.  
  307. /* We restrict radius to the values: 0.333333 <= radius <= 1.0 */
  308. /* so that no fewer and no more than a 3x3 grid of pixels around */
  309. /* the pixel in question needs to be read. Given this, we only */
  310. /* need 3 or 4 weightings per hexagon, as follows: */
  311. /*                  _ _                         */
  312. /* Virtical hex:   |_|_|  1 2                   */
  313. /*                 |X|_|  0 3                   */
  314. /*                                       _      */
  315. /*              _                      _|_|   1 */
  316. /* Middle hex: |_| 1  Horizontal hex: |X|_| 0 2 */
  317. /*             |X| 0                    |_|   3 */
  318. /*             |_| 2                            */
  319.  
  320. /* all filters */
  321. int V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL];        /* vertical hex */
  322. int M0[NOIVAL],M1[NOIVAL],M2[NOIVAL];           /* middle hex */
  323. int H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL];        /* horizontal hex */
  324.  
  325. /* alpha trimmed and edge enhancement only */
  326. int ALFRAC[NOIVAL * 8];                 /* fractional alpha divider table */
  327.  
  328. /* optimal estimation only */
  329. int AVEDIV[7 * NOCSVAL];                /* divide by 7 to give average value */
  330. int SQUARE[2 * NOCSVAL];                /* scaled square lookup table */
  331.  
  332. /* Table initialisation function - return alpha range */
  333. int
  334. atfilt_setup(alpha,radius,maxscale)
  335. double alpha,radius,maxscale;   /* alpha, radius and amount to scale input pixel values */
  336.         {
  337.         /* other function variables */
  338.         int alpharange;                 /* alpha range value 0 - 3 */
  339.         double meanscale;               /* scale for finding mean */
  340.         double mmeanscale;              /* scale for finding mean - midle hex */
  341.         double alphafraction;   /* fraction of next largest/smallest to subtract from sum */
  342.  
  343.         /* do setup */
  344.  
  345.         if (alpha >= 0.0 && alpha < 1.0)        /* alpha trimmed mean */
  346.                 {
  347.                 double noinmean;
  348.                 /* number of elements (out of a possible 7) used in the mean */
  349.                 noinmean = ((0.5 - alpha) * 12.0) + 1.0;
  350.                 mmeanscale = meanscale = maxscale/noinmean;
  351.                 if (alpha == 0.0)                    /* mean filter */
  352.                         {
  353.                         alpharange = 0;
  354.                         alphafraction = 0.0;            /* not used */
  355.                         }
  356.                 else if (alpha < (1.0/6.0))     /* mean of 5 to 7 middle values */
  357.                         {
  358.                         alpharange = 1;
  359.                         alphafraction = (7.0 - noinmean)/2.0;
  360.                         }
  361.                 else if (alpha < (1.0/3.0))     /* mean of 3 to 5 middle values */
  362.                         {
  363.                         alpharange = 2;
  364.                         alphafraction = (5.0 - noinmean)/2.0;
  365.                         }
  366.                 else                            /* mean of 1 to 3 middle values */
  367.                         {                                                       /* alpha == 0.5 == median filter */
  368.                         alpharange = 3;
  369.                         alphafraction = (3.0 - noinmean)/2.0;
  370.                         }
  371.                 }
  372.         else if (alpha > 0.5)   /* optimal estimation - alpha controls noise variance threshold. */
  373.                 {
  374.                 int i;
  375.                 double noinmean = 7.0;
  376.                 alpharange = 5;                 /* edge enhancement function */
  377.                 alpha -= 1.0;                   /* normalise it to 0.0 -> 1.0 */
  378.                 mmeanscale = meanscale = maxscale;      /* compute scaled hex values */
  379.                 alphafraction = 1.0/noinmean;   /* Set up 1:1 division lookup - not used */
  380.                 noisevariance = alpha * (double)omaxval;
  381.                 noisevariance = noisevariance * noisevariance / 8.0;    /* estimate of noise variance */
  382.                 /* set yp optimal estimation specific stuff */
  383.                 for (i=0; i < (7 * NOCSVAL); i++)       /* divide scaled value by 7 lookup */
  384.                         {
  385.                         AVEDIV[i] = CSCTOSC(i)/7;       /* scaled divide by 7 */
  386.                         }
  387.                 for (i=0; i < (2 * NOCSVAL); i++)  /* compute square and rescale by (val >> (2 * SCALEB + 2)) table */
  388.                         {
  389.                         int val;
  390.                         val = CSCTOSC(i - NOCSVAL); /* NOCSVAL offset to cope with -ve input values */
  391.                         SQUARE[i] = (val * val) >> (2 * SCALEB + 2);
  392.                         }
  393.                 }
  394.         else    /* edge enhancement function */
  395.                 {
  396.                 alpharange = 4;                 /* edge enhancement function */
  397.                 alpha = -alpha;                 /* turn it the right way up */
  398.                 meanscale = maxscale * (-alpha/((1.0 - alpha) * 7.0)); /* mean of 7 and scaled by -alpha/(1-alpha) */
  399.                 mmeanscale = maxscale * (1.0/(1.0 - alpha) + meanscale);        /* middle pixel has 1/(1-alpha) as well */
  400.                 alphafraction = 0.0;    /* not used */
  401.                 }
  402.  
  403.                 /* Setup pixel weighting tables - note we pre-compute mean division here too. */
  404.                 {
  405.                 int i;
  406.                 double hexhoff,hexvoff;
  407.                 double tabscale,mtabscale;
  408.                 double v0,v1,v2,v3,m0,m1,m2,me0,me1,me2,h0,h1,h2,h3;
  409.  
  410.                 hexhoff = radius/2;                 /* horizontal offset of virtical hex centers */
  411.                 hexvoff = 3.0 * radius/sqrt(12.0);      /* virtical offset of virtical hex centers */
  412.                 /* scale tables to normalise by hexagon area, and number of hexes used in mean */
  413.                 tabscale = meanscale / (radius * hexvoff);
  414.                 mtabscale = mmeanscale / (radius * hexvoff);
  415.                 v0 = hex_area(0.0,0.0,hexhoff,hexvoff,radius) * tabscale;
  416.                 v1 = hex_area(0.0,1.0,hexhoff,hexvoff,radius) * tabscale;
  417.                 v2 = hex_area(1.0,1.0,hexhoff,hexvoff,radius) * tabscale;
  418.                 v3 = hex_area(1.0,0.0,hexhoff,hexvoff,radius) * tabscale;
  419.                 m0 = hex_area(0.0,0.0,0.0,0.0,radius) * mtabscale;
  420.                 m1 = hex_area(0.0,1.0,0.0,0.0,radius) * mtabscale;
  421.                 m2 = hex_area(0.0,-1.0,0.0,0.0,radius) * mtabscale;
  422.                 h0 = hex_area(0.0,0.0,radius,0.0,radius) * tabscale;
  423.                 h1 = hex_area(1.0,1.0,radius,0.0,radius) * tabscale;
  424.                 h2 = hex_area(1.0,0.0,radius,0.0,radius) * tabscale;
  425.                 h3 = hex_area(1.0,-1.0,radius,0.0,radius) * tabscale;
  426.  
  427.                 for (i=0; i <= MXIVAL; i++)
  428.                         {
  429.                         double fi;
  430.                         fi = (double)i;
  431.                         V0[i] = ROUND(fi * v0);
  432.                         V1[i] = ROUND(fi * v1);
  433.                         V2[i] = ROUND(fi * v2);
  434.                         V3[i] = ROUND(fi * v3);
  435.                         M0[i] = ROUND(fi * m0);
  436.                         M1[i] = ROUND(fi * m1);
  437.                         M2[i] = ROUND(fi * m2);
  438.                         H0[i] = ROUND(fi * h0);
  439.                         H1[i] = ROUND(fi * h1);
  440.                         H2[i] = ROUND(fi * h2);
  441.                         H3[i] = ROUND(fi * h3);
  442.                         }
  443.                 /* set up alpha fraction lookup table used on big/small */
  444.                 for (i=0; i < (NOIVAL * 8); i++)
  445.                         {
  446.                         ALFRAC[i] = ROUND((double)i * alphafraction);
  447.                         }
  448.                 }
  449.         return alpharange;
  450.         }
  451.  
  452.  
  453. /* Core pixel processing function - hand it 3x3 pixels and return result. */
  454. /* Mean filter */
  455. int
  456. atfilt0(p)
  457. int *p;         /* 9 pixel values from 3x3 neighbors */
  458.         {
  459.         int retv;
  460.         /* map to scaled hexagon values */
  461.         retv = M0[p[0]] + M1[p[3]] + M2[p[7]];
  462.         retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  463.         retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  464.         retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  465.         retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  466.         retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  467.         retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  468.         return UNSCALE(retv);
  469.         }
  470.  
  471. /* Mean of 5 - 7 middle values */
  472. int
  473. atfilt1(p)
  474. int *p;         /* 9 pixel values from 3x3 neighbors */
  475.         {
  476.         int h0,h1,h2,h3,h4,h5,h6;       /* hexagon values    2 3   */
  477.                                     /*                  1 0 4  */
  478.                                     /*                   6 5   */
  479.         int big,small;
  480.         /* map to scaled hexagon values */
  481.         h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  482.         h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  483.         h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  484.         h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  485.         h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  486.         h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  487.         h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  488.         /* sum values and also discover the largest and smallest */
  489.         big = small = h0;
  490. #define CHECK(xx) \
  491.         h0 += xx; \
  492.         if (xx > big) \
  493.                 big = xx; \
  494.         else if (xx < small) \
  495.                 small = xx;
  496.         CHECK(h1)
  497.         CHECK(h2)
  498.         CHECK(h3)
  499.         CHECK(h4)
  500.         CHECK(h5)
  501.         CHECK(h6)
  502. #undef CHECK
  503.         /* Compute mean of middle 5-7 values */
  504.         return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]);
  505.         }
  506.  
  507. /* Mean of 3 - 5 middle values */
  508. int
  509. atfilt2(p)
  510. int *p;         /* 9 pixel values from 3x3 neighbors */
  511.         {
  512.         int h0,h1,h2,h3,h4,h5,h6;       /* hexagon values    2 3   */
  513.                                     /*                  1 0 4  */
  514.                                     /*                   6 5   */
  515.         int big0,big1,small0,small1;
  516.         /* map to scaled hexagon values */
  517.         h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  518.         h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  519.         h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  520.         h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  521.         h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  522.         h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  523.         h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  524.         /* sum values and also discover the 2 largest and 2 smallest */
  525.         big0 = small0 = h0;
  526.         small1 = MAXINT;
  527.         big1 = 0;
  528. #define CHECK(xx) \
  529.         h0 += xx; \
  530.         if (xx > big1) \
  531.                 { \
  532.                 if (xx > big0) \
  533.                         { \
  534.                         big1 = big0; \
  535.                         big0 = xx; \
  536.                         } \
  537.                 else \
  538.                         big1 = xx; \
  539.                 } \
  540.         if (xx < small1) \
  541.                 { \
  542.                 if (xx < small0) \
  543.                         { \
  544.                         small1 = small0; \
  545.                         small0 = xx; \
  546.                         } \
  547.                 else \
  548.                         small1 = xx; \
  549.                 }
  550.         CHECK(h1)
  551.         CHECK(h2)
  552.         CHECK(h3)
  553.         CHECK(h4)
  554.         CHECK(h5)
  555.         CHECK(h6)
  556. #undef CHECK
  557.         /* Compute mean of middle 3-5 values */
  558.         return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]);
  559.         }
  560.  
  561. /* Mean of 1 - 3 middle values. If only 1 value, then this is a median filter. */
  562. int
  563. atfilt3(p)
  564. int *p;         /* 9 pixel values from 3x3 neighbors */
  565.         {
  566.         int h0,h1,h2,h3,h4,h5,h6;       /* hexagon values    2 3   */
  567.                                     /*                  1 0 4  */
  568.                                     /*                   6 5   */
  569.         int big0,big1,big2,small0,small1,small2;
  570.         /* map to scaled hexagon values */
  571.         h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  572.         h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  573.         h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  574.         h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  575.         h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  576.         h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  577.         h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  578.         /* sum values and also discover the 3 largest and 3 smallest */
  579.         big0 = small0 = h0;
  580.         small1 = small2 = MAXINT;
  581.         big1 = big2 = 0;
  582. #define CHECK(xx) \
  583.         h0 += xx; \
  584.         if (xx > big2) \
  585.                 { \
  586.                 if (xx > big1) \
  587.                         { \
  588.                         if (xx > big0) \
  589.                                 { \
  590.                                 big2 = big1; \
  591.                                 big1 = big0; \
  592.                                 big0 = xx; \
  593.                                 } \
  594.                         else \
  595.                                 { \
  596.                                 big2 = big1; \
  597.                                 big1 = xx; \
  598.                                 } \
  599.                         } \
  600.                 else \
  601.                         big2 = xx; \
  602.                 } \
  603.         if (xx < small2) \
  604.                 { \
  605.                 if (xx < small1) \
  606.                         { \
  607.                         if (xx < small0) \
  608.                                 { \
  609.                                 small2 = small1; \
  610.                                 small1 = small0; \
  611.                                 small0 = xx; \
  612.                                 } \
  613.                         else \
  614.                                 { \
  615.                                 small2 = small1; \
  616.                                 small1 = xx; \
  617.                                 } \
  618.                         } \
  619.                 else \
  620.                         small2 = xx; \
  621.                 }
  622.         CHECK(h1)
  623.         CHECK(h2)
  624.         CHECK(h3)
  625.         CHECK(h4)
  626.         CHECK(h5)
  627.         CHECK(h6)
  628. #undef CHECK
  629.         /* Compute mean of middle 1-3 values */
  630.         return  UNSCALE(h0 -big0 -big1 -small0 -small1 -ALFRAC[(big2 + small2)>>SCALEB]);
  631.         }
  632.  
  633. /* Edge enhancement */
  634. /* notice we use the global omaxval */
  635. int
  636. atfilt4(p)
  637. int *p;         /* 9 pixel values from 3x3 neighbors */
  638.         {
  639.         int hav;
  640.         /* map to scaled hexagon values and compute enhance value */
  641.         hav = M0[p[0]] + M1[p[3]] + M2[p[7]];
  642.         hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  643.         hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  644.         hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  645.         hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  646.         hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  647.         hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  648.         if (hav < 0)
  649.                 hav = 0;
  650.         hav = UNSCALE(hav);
  651.         if (hav > omaxval)
  652.                 hav = omaxval;
  653.         return hav;
  654.         }
  655.  
  656. /* Optimal estimation - do smoothing in inverse proportion */
  657. /* to the local variance. */
  658. /* notice we use the globals noisevariance and omaxval*/
  659. int
  660. atfilt5(p)
  661. int *p;         /* 9 pixel values from 3x3 neighbors */
  662.         {
  663.         int mean,variance,temp;
  664.         int h0,h1,h2,h3,h4,h5,h6;       /* hexagon values    2 3   */
  665.                                     /*                  1 0 4  */
  666.                                     /*                   6 5   */
  667.         /* map to scaled hexagon values */
  668.         h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  669.         h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  670.         h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  671.         h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  672.         h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  673.         h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  674.         h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  675.         mean = h0 + h1 + h2 + h3 + h4 + h5 + h6;
  676.         mean = AVEDIV[SCTOCSC(mean)];   /* compute scaled mean by dividing by 7 */
  677.         temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)];  /* compute scaled variance */
  678.         temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* and rescale to keep */
  679.         temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* within 32 bit limits */
  680.         temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  681.         temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  682.         temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  683.         temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];        /* (temp = h0 - mean) */
  684.         if (variance != 0)      /* avoid possible divide by 0 */
  685.                 temp = mean + (variance * temp) / (variance + noisevariance);   /* optimal estimate */
  686.         else temp = h0;
  687.         if (temp < 0)
  688.                 temp = 0;
  689.         temp = RUNSCALE(temp);
  690.         if (temp > omaxval)
  691.                 temp = omaxval;
  692.         return temp;
  693.         }
  694.  
  695. /* ************************************************** */
  696. /* Hexagon intersecting square area functions */
  697. /* Compute the area of the intersection of a triangle */
  698. /* and a rectangle */
  699.  
  700. double triang_area ARGS((double, double, double, double, double, double, double, double, int));
  701. double rectang_area ARGS((double, double, double, double, double, double, double, double));
  702.  
  703. /* Triangle orientation is per geometric axes (not graphical axies) */
  704.  
  705. #define NW 0    /* North west triangle /| */
  706. #define NE 1    /* North east triangle |\ */
  707. #define SW 2    /* South west triangle \| */
  708. #define SE 3    /* South east triangle |/ */
  709. #define STH 2
  710. #define EST 1
  711.  
  712. #define SWAPI(a,b) (t = a, a = -b, b = -t)
  713.  
  714. /* compute the area of overlap of a hexagon diameter d, */
  715. /* centered at hx,hy, with a unit square of center sx,sy. */
  716. double
  717. hex_area(sx,sy,hx,hy,d)
  718. double sx,sy;   /* square center */
  719. double hx,hy,d; /* hexagon center and diameter */
  720.         {
  721.         double hx0,hx1,hx2,hy0,hy1,hy2,hy3;
  722.         double sx0,sx1,sy0,sy1;
  723.  
  724.         /* compute square co-ordinates */
  725.         sx0 = sx - 0.5;
  726.         sy0 = sy - 0.5;
  727.         sx1 = sx + 0.5;
  728.         sy1 = sy + 0.5;
  729.         /* compute hexagon co-ordinates */
  730.         hx0 = hx - d/2.0;
  731.         hx1 = hx;
  732.         hx2 = hx + d/2.0;
  733.         hy0 = hy - 0.5773502692 * d;    /* d / sqrt(3) */
  734.         hy1 = hy - 0.2886751346 * d;    /* d / sqrt(12) */
  735.         hy2 = hy + 0.2886751346 * d;    /* d / sqrt(12) */
  736.         hy3 = hy + 0.5773502692 * d;    /* d / sqrt(3) */
  737.  
  738.         return
  739.                 triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) +
  740.                 triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) +
  741.                 rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) +
  742.                 triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) +
  743.                 triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE);
  744.         }
  745.  
  746. double
  747. triang_area(rx0,ry0,rx1,ry1,tx0,ty0,tx1,ty1,tt)
  748. double rx0,ry0,rx1,ry1;         /* rectangle boundaries */
  749. double tx0,ty0,tx1,ty1;         /* triangle boundaries */
  750. int tt;                                         /* triangle type */
  751.         {
  752.         double a,b,c,d;
  753.         double lx0,ly0,lx1,ly1;
  754.         /* Convert everything to a NW triangle */
  755.         if (tt & STH)
  756.                 {
  757.                 double t;
  758.         SWAPI(ry0,ry1);
  759.         SWAPI(ty0,ty1);
  760.                 }
  761.         if (tt & EST)
  762.                 {
  763.                 double t;
  764.         SWAPI(rx0,rx1);
  765.         SWAPI(tx0,tx1);
  766.                 }
  767.         /* Compute overlapping box */
  768.         if (tx0 > rx0)
  769.                 rx0 = tx0;
  770.         if (ty0 > ry0)
  771.                 ry0 = ty0;
  772.         if (tx1 < rx1)
  773.                 rx1 = tx1;
  774.         if (ty1 < ry1)
  775.                 ry1 = ty1;
  776.         if (rx1 <= rx0 || ry1 <= ry0)
  777.                 return 0.0;
  778.         /* Need to compute diagonal line intersection with the box */
  779.         /* First compute co-efficients to formulas x = a + by and y = c + dx */
  780.         b = (tx1 - tx0)/(ty1 - ty0);
  781.         a = tx0 - b * ty0;
  782.         d = (ty1 - ty0)/(tx1 - tx0);
  783.         c = ty0 - d * tx0;
  784.  
  785.         /* compute top or right intersection */
  786.         tt = 0;
  787.         ly1 = ry1;
  788.         lx1 = a + b * ly1;
  789.         if (lx1 <= rx0)
  790.                 return (rx1 - rx0) * (ry1 - ry0);
  791.         else if (lx1 > rx1)     /* could be right hand side */
  792.                 {
  793.                 lx1 = rx1;
  794.                 ly1 = c + d * lx1;
  795.                 if (ly1 <= ry0)
  796.                         return (rx1 - rx0) * (ry1 - ry0);
  797.                 tt = 1; /* right hand side intersection */
  798.                 }
  799.         /* compute left or bottom intersection */
  800.         lx0 = rx0;
  801.         ly0 = c + d * lx0;
  802.         if (ly0 >= ry1)
  803.                 return (rx1 - rx0) * (ry1 - ry0);
  804.         else if (ly0 < ry0)     /* could be right hand side */
  805.                 {
  806.                 ly0 = ry0;
  807.                 lx0 = a + b * ly0;
  808.                 if (lx0 >= rx1)
  809.                         return (rx1 - rx0) * (ry1 - ry0);
  810.                 tt |= 2;        /* bottom intersection */
  811.                 }
  812.  
  813.         if (tt == 0)    /* top and left intersection */
  814.                 {       /* rectangle minus triangle */
  815.                 return ((rx1 - rx0) * (ry1 - ry0))
  816.                      - (0.5 * (lx1 - rx0) * (ry1 - ly0));
  817.                 }
  818.         else if (tt == 1)       /* right and left intersection */
  819.                 {
  820.                 return ((rx1 - rx0) * (ly0 - ry0))
  821.                      + (0.5 * (rx1 - rx0) * (ly1 - ly0));
  822.                 }
  823.         else if (tt == 2)       /* top and bottom intersection */
  824.                 {
  825.                 return ((rx1 - lx1) * (ry1 - ry0))
  826.                      + (0.5 * (lx1 - lx0) * (ry1 - ry0));
  827.                 }
  828.         else /* tt == 3 */      /* right and bottom intersection */
  829.                 {       /* triangle */
  830.                 return (0.5 * (rx1 - lx0) * (ly1 - ry0));
  831.                 }
  832.         }
  833.  
  834. /* Compute rectangle area */
  835. double
  836. rectang_area(rx0,ry0,rx1,ry1,tx0,ty0,tx1,ty1)
  837. double rx0,ry0,rx1,ry1;         /* rectangle boundaries */
  838. double tx0,ty0,tx1,ty1;         /* rectangle boundaries */
  839.         {
  840.         /* Compute overlapping box */
  841.         if (tx0 > rx0)
  842.                 rx0 = tx0;
  843.         if (ty0 > ry0)
  844.                 ry0 = ty0;
  845.         if (tx1 < rx1)
  846.                 rx1 = tx1;
  847.         if (ty1 < ry1)
  848.                 ry1 = ty1;
  849.         if (rx1 <= rx0 || ry1 <= ry0)
  850.                 return 0.0;
  851.         return (rx1 - rx0) * (ry1 - ry0);
  852.         }
  853.  
  854.