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