home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1994 March / Source_Code_CD-ROM_Walnut_Creek_March_1994.iso / compsrcs / misc / volume35 / pnmnlflt / part01 < prev    next >
Encoding:
Text File  |  1993-03-03  |  35.1 KB  |  1,061 lines
  1. Newsgroups: comp.sources.misc
  2. From: graeme@labtam.labtam.oz.au (Graeme Gill)
  3. Subject: v35i116:  pnmnlfilt - A nonlinear area filter for the pbm package, Part01/01
  4. Message-ID: <1993Mar4.190827.6812@sparky.imd.sterling.com>
  5. X-Md4-Signature: 7241a5f87217d6750c6d5ba0e9d5420b
  6. Date: Thu, 4 Mar 1993 19:08:27 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: graeme@labtam.labtam.oz.au (Graeme Gill)
  10. Posting-number: Volume 35, Issue 116
  11. Archive-name: pnmnlfilt/part01
  12. Environment: PBM
  13.  
  14. WHAT IS pnmnlfilt ?
  15.  
  16.     It is another utility for the pbm library.
  17.  
  18.     pnmnlfilt - non-linear filters: smooth, alpha trim mean, optimal
  19.     estimation smoothing, edge enhancement.
  20.  
  21.     This is something of a swiss army knife filter. It has 3 distinct
  22.     operating modes.
  23.  
  24.     The optimal estimation filter and the edge enhancement filter in
  25.     particular are good for turning quantized (ie. 8 bit color) images
  26.     back into 24 bit color images suitable for jpeg encoding.
  27.     It can do this with less blurring that ordinary smoothing
  28.     filters.  The filter may also be useful for smoothing out noise
  29.     artifacts of distributed ray tracing images. The various modes
  30.     of filtering work well in combination (eg. edge enhancement
  31.     after optimal estimation filtering etc.).
  32.  
  33.  
  34. INSTALLATION:
  35.  
  36.     Unpack the files.
  37.  
  38.     You will need to have already installed and compiled the pbm package.
  39.     pnmnlfilt.c and pnmnlfilt.1 should be copied to the pnm directory.
  40.  
  41.     The Makefile (or Imakefile if you are using them) file in pnm should be
  42.     modified to add "pnmnlfilt" to the "MATHBINARIES" = define.
  43.     Add "pnmnlfilt.1" to the "MANUALS1 =" define.
  44.  
  45.     Re-make the pnm tools in the usual way. The manual entry gives
  46.     details of the options and suggested usage of pnmnlfilt.
  47.  
  48.     Graeme W. Gill
  49.         graeme@labtam.oz.au
  50. -----
  51. #! /bin/sh
  52. # This is a shell archive.  Remove anything before this line, then feed it
  53. # into a shell via "sh file" or similar.  To overwrite existing files,
  54. # type "sh file -c".
  55. # Contents:  pnmnlfilt.1 pnmnlfilt.c
  56. # Wrapped by kent@sparky on Thu Mar  4 13:05:22 1993
  57. PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin ; export PATH
  58. echo If this archive is complete, you will see the following message:
  59. echo '          "shar: End of archive 1 (of 1)."'
  60. if test -f 'pnmnlfilt.1' -a "${1}" != "-c" ; then 
  61.   echo shar: Will not clobber existing file \"'pnmnlfilt.1'\"
  62. else
  63.   echo shar: Extracting \"'pnmnlfilt.1'\" \(5647 characters\)
  64.   sed "s/^X//" >'pnmnlfilt.1' <<'END_OF_FILE'
  65. X.TH pnmnlfilt 1 "5 February 1993"
  66. X.IX pnmnlfilt
  67. X.SH NAME
  68. Xpnmnlfilt - non-linear filters: smooth, alpha trim mean, optimal
  69. Xestimation smoothing, edge enhancement.
  70. X.SH SYNOPSIS
  71. X.B pnmnlfilt
  72. X.RI alpha
  73. X.RI radius
  74. X.RI [ pnmfile ]
  75. X.SH DESCRIPTION
  76. X.IX smoothing
  77. X.IX dithering
  78. X.IX alpha trim
  79. X.IX mean filter
  80. X.IX median filter
  81. X.IX optimal estimation
  82. XThis is something of a swiss army knife filter. It has 3 distinct operating
  83. Xmodes. In all of the modes each pixel in the image is examined and processed
  84. Xaccording to it and its surrounding pixels values. Rather than using the
  85. X9 pixels in a 3x3 block, 7 hexagonal area samples are taken, the size of
  86. Xthe hexagons being controlled by the radius parameter. A radius value of
  87. X0.3333 means that the 7 hexagons exactly fit into the center pixel (ie.
  88. Xthere will be no filtering effect). A radius value of 1.0 means that
  89. Xthe 7 hexagons exactly fit a 3x3 pixel array.
  90. X.SH Alpha trimmed mean filter.    (0.0 <= alpha <= 0.5)
  91. X.PP
  92. XThe value of the center pixel will be
  93. Xreplaced by the mean of the 7 hexagon values, but the 7 values are
  94. Xsorted by size and the top and bottom alpha portion of the 7 are
  95. Xexcluded from the mean.  This implies that an alpha value of 0.0 gives
  96. Xthe same sort of output as a normal convolution (ie. averaging or
  97. Xsmoothing filter), where radius will determine the "strength" of the
  98. Xfilter. A good value to start from for subtle filtering is alpha = 0.0, radius = 0.55
  99. XFor a more blatant effect, try alpha 0.0 and radius 1.0
  100. X.PP
  101. XAn alpha value of 0.5 will cause the median value of the
  102. X7 hexagons to be used to replace the center pixel value. This sort
  103. Xof filter is good for eliminating "pop" or single pixel noise from
  104. Xan image without spreading the noise out or smudging features on
  105. Xthe image. Judicious use of the radius parameter will fine tune the
  106. Xfiltering. Intermediate values of alpha give effects somewhere
  107. Xbetween smoothing and "pop" noise reduction. For subtle filtering
  108. Xtry starting with values of alpha = 0.4, radius = 0.6  For a more blatant
  109. Xeffect try alpha = 0.5, radius = 1.0
  110. X.SH Optimal estimation smoothing. (1.0 <= alpha <= 2.0)
  111. X.PP
  112. XThis type of filter applies a smoothing filter adaptively over the image.
  113. XFor each pixel the variance of the surrounding hexagon values is calculated,
  114. Xand the amount of smoothing is made inversely proportional to it. The idea
  115. Xis that if the variance is small then it is due to noise in the image, while
  116. Xif the variance is large, it is because of "wanted" image features. As usual
  117. Xthe radius parameter controls the effective radius, but it probably advisable to
  118. Xleave the radius between 0.8 and 1.0 for the variance calculation to be meaningful.
  119. XThe alpha parameter sets the noise threshold, over which less smoothing will be done.
  120. XThis means that small values of alpha will give the most subtle filtering effect,
  121. Xwhile large values will tend to smooth all parts of the image. You could start
  122. Xwith values like alpha = 1.2, radius = 1.0 and try increasing or decreasing the
  123. Xalpha parameter to get the desired effect. This type of filter is best for
  124. Xfiltering out dithering noise in both bitmap and color images.
  125. X.SH Edge enhancement. (-0.1 >= alpha >= -0.9)
  126. X.PP
  127. XThis is the opposite type of filter to the smoothing filter. It enhances
  128. Xedges. The alpha parameter controls the amount of edge enhancement, from
  129. Xsubtle (-0.1) to blatant (-0.9). The radius parameter controls the effective
  130. Xradius as usual, but useful values are between 0.5 and 0.9. Try starting
  131. Xwith values of alpha = 0.3, radius = 0.8
  132. X.SH Combination use.
  133. X.PP
  134. XThe various modes of 
  135. X.B pnmnlfilt
  136. Xcan be used one after the other to get the desired result. For instance to
  137. Xturn a monochrome dithered image into a grayscale image you could try
  138. Xone or two passes of the smoothing filter, followed by a pass of the optimal estimation
  139. Xfilter, then some subtle edge enhancement. Note that using edge enhancement is
  140. Xonly likely to be useful after one of the non-linear filters (alpha trimmed mean
  141. Xor optimal estimation filter), as edge enhancement is the direct opposite of
  142. Xsmoothing.
  143. X.PP
  144. XFor reducing color quantization noise in images (ie. turning .gif files back into
  145. X24 bit files) you could try a pass of the optimal estimation filter
  146. X(alpha 1.2, radius 1.0), a pass of the median filter (alpha 0.5, radius 0.55),
  147. Xand possibly a pass of the edge enhancement filter.
  148. XSeveral passes of the optimal estimation filter with declining alpha
  149. Xvalues are more effective than a single pass with a large alpha value.
  150. XAs usual, there is a tradeoff between filtering effectiveness and loosing
  151. Xdetail. Experimentation is encouraged.
  152. X.SH References:
  153. X.PP
  154. XThe alpha-trimmed mean filter is 
  155. Xbased on the description in IEEE CG&A May 1990 
  156. XPage 23 by Mark E. Lee and Richard A. Redner,
  157. Xand has been enhanced to allow continuous alpha adjustment.
  158. X.PP
  159. XThe optimal estimation filter is taken from an article "Converting Dithered
  160. XImages Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November
  161. X1992, and this article references "Digital Image Enhancement and Noise Filtering by
  162. XUse of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and
  163. XMachine Intelligence, March 1980.
  164. X.PP
  165. XThe edge enhancement details are from pgmenhance(1),
  166. Xwhich is taken from Philip R. Thompson's "xim"
  167. Xprogram, which in turn took it from section 6 of "Digital Halftones by
  168. XDot Diffusion", D. E. Knuth, ACM Transaction on Graphics Vol. 6, No. 4,
  169. XOctober 1987, which in turn got it from two 1976 papers by J. F. Jarvis
  170. Xet. al.
  171. X.SH "SEE ALSO"
  172. Xpgmenhance(1), pnmconvol(1), pnm(5)
  173. X.SH BUGS
  174. XIntegers and tables may overflow if PPM_MAXMAXVAL is greater than 255.
  175. X.SH AUTHOR
  176. XGraeme W. Gill    graeme@labtam.oz.au
  177. END_OF_FILE
  178.   if test 5647 -ne `wc -c <'pnmnlfilt.1'`; then
  179.     echo shar: \"'pnmnlfilt.1'\" unpacked with wrong size!
  180.   fi
  181.   # end of 'pnmnlfilt.1'
  182. fi
  183. if test -f 'pnmnlfilt.c' -a "${1}" != "-c" ; then 
  184.   echo shar: Will not clobber existing file \"'pnmnlfilt.c'\"
  185. else
  186.   echo shar: Extracting \"'pnmnlfilt.c'\" \(25912 characters\)
  187.   sed "s/^X//" >'pnmnlfilt.c' <<'END_OF_FILE'
  188. X/* pnmnlfilt.c - 4 in 1 (2 non-linear) filter
  189. X**             - smooth an anyimage
  190. X**             - do alpha trimmed mean filtering on an anyimage
  191. X**             - do optimal estimation smoothing on an anyimage
  192. X**             - do edge enhancement on an anyimage
  193. X**
  194. X** Version 1.0
  195. X**
  196. X** The implementation of an alpha-trimmed mean filter
  197. X** is based on the description in IEEE CG&A May 1990 
  198. X** Page 23 by Mark E. Lee and Richard A. Redner. 
  199. X**
  200. X** The paper recommends using a hexagon sampling region around each 
  201. X** pixel being processed, allowing an effective sub pixel radius to be 
  202. X** specified. The hexagon values are sythesised by area sampling the 
  203. X** rectangular pixels with a hexagon grid. The seven hexagon values 
  204. X** obtained from the 3x3 pixel grid are used to compute the alpha 
  205. X** trimmed mean. Note that an alpha value of 0.0 gives a conventional 
  206. X** mean filter (where the radius controls the contribution of 
  207. X** surrounding pixels), while a value of 0.5 gives a median filter. 
  208. X** Although there are only seven values to trim from before finding 
  209. X** the mean, the algorithm has been extended from that described in 
  210. X** CG&A by using interpolation, to allow a continuous selection of 
  211. X** alpha value between and including 0.0 to 0.5  The useful values 
  212. X** for radius are between 0.3333333 (where the filter will have no 
  213. X** effect because only one pixel is sampled), to 1.0, where all 
  214. X** pixels in the 3x3 grid are sampled. 
  215. X**
  216. X** The optimal estimation filter is taken from an article "Converting Dithered
  217. X** Images Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November
  218. X** 1992, and this article references "Digital Image Enhancement andNoise Filtering by
  219. X** Use of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and
  220. X** Machine Intelligence, March 1980.
  221. X**
  222. X** Also borrow the  technique used in pgmenhance(1) to allow edge
  223. X** enhancement if the alpha value is negative.
  224. X**
  225. X** Author:
  226. X**         Graeme W. Gill, 30th Jan 1993 
  227. X**         graeme@labtam.oz.au
  228. X*/
  229. X
  230. X#include "pnm.h"
  231. X#include <math.h>
  232. Xdouble hex_area();
  233. Xint atfilt_setup();
  234. Xint atfilt0(),atfilt1(),atfilt2(),atfilt3(),atfilt4(),atfilt5();
  235. Xint (*atfuncs[6])() = {atfilt0,atfilt1,atfilt2,atfilt3,atfilt4,atfilt5};
  236. X
  237. Xxelval omaxval;    /* global so that pixel processing code can get at it quickly */
  238. Xint noisevariance;    /* global so that pixel processing code can get at it quickly */
  239. X
  240. Xvoid
  241. Xmain( argc, argv )
  242. Xint argc;
  243. Xchar* argv[];
  244. X    {
  245. X    double radius=0.0,alpha= -1.0;
  246. X    FILE* ifp;
  247. X    int rows, cols, format, oformat, row, col;
  248. X    xelval maxval;
  249. X    int (*atfunc)();
  250. X    char* usage = "alpha radius pnmfile\n\
  251. X 0.0 <= alpha <= 0.5 for alpha trimmed mean -or- \n\
  252. X 1.0 <= alpha <= 2.0 for optimal estimation -or- \n\
  253. X -0.1 >= alpha >= -0.9 for edge enhancement\n\
  254. X 0.3333 <= radius <= 1.0 specify effective radius\n";
  255. X
  256. X    pnm_init( &argc, argv );
  257. X
  258. X    if ( argc < 3 || argc > 4 )
  259. X        pm_usage( usage );
  260. X
  261. X    if ( sscanf( argv[1], "%lf", &alpha ) != 1 )
  262. X        pm_usage( usage );
  263. X    if ( sscanf( argv[2], "%lf", &radius ) != 1 )
  264. X        pm_usage( usage );
  265. X
  266. X    if ((alpha > -0.1 && alpha < 0.0) || (alpha > 0.5 && alpha < 1.0))
  267. X        pm_error( "Alpha must be in range 0.0 <= alpha <= 0.5 for alpha trimmed mean" );
  268. X    if (alpha > 2.0)
  269. X        pm_error( "Alpha must be in range 1.0 <= alpha <= 2.0 for optimal estimation" );
  270. X    if (alpha < -0.9 || (alpha > -0.1 && alpha < 0.0))
  271. X        pm_error( "Alpha must be in range -0.9 <= alpha <= -0.1 for edge enhancement" );
  272. X    if (radius < 0.333 || radius > 1.0)
  273. X        pm_error( "Radius must be in range 0.333333333 <= radius <= 1.0" );
  274. X
  275. X    if ( argc == 4 )
  276. X        ifp = pm_openr( argv[3] );
  277. X    else
  278. X        ifp = stdin;
  279. X
  280. X    pnm_readpnminit( ifp, &cols, &rows, &maxval, &format );
  281. X
  282. X    oformat = PNM_FORMAT_TYPE(format);
  283. X    omaxval = PPM_MAXMAXVAL;    /* force output to max precision */
  284. X
  285. X    atfunc = atfuncs[atfilt_setup(alpha,radius,(double)omaxval/(double)maxval)];
  286. X
  287. X    if ( oformat < PGM_TYPE )
  288. X        {
  289. X        if (PGM_MAXMAXVAL > PPM_MAXMAXVAL)
  290. X            pm_error( "Can't handle pgm input file as maxval is too big" );
  291. X        oformat = RPGM_FORMAT;
  292. X        pm_message( "promoting file to PGM" );
  293. X        }
  294. X    pnm_writepnminit( stdout, cols, rows, omaxval, oformat, 0 );
  295. X
  296. X    if ( PNM_FORMAT_TYPE(oformat) == PPM_TYPE )
  297. X        {
  298. X        xel *irows[3];
  299. X        xel *irow0, *irow1, *irow2, *orow;
  300. X        int pr[9],pg[9],pb[9];        /* 3x3 neighbor pixel values */
  301. X        int r,g,b;
  302. X
  303. X        irows[0] = pnm_allocrow( cols );
  304. X        irows[1] = pnm_allocrow( cols );
  305. X        irows[2] = pnm_allocrow( cols );
  306. X        irow0 = irows[0];
  307. X        irow1 = irows[1];
  308. X        irow2 = irows[2];
  309. X        orow = pnm_allocrow( cols );
  310. X
  311. X        for ( row = 0; row < rows; row++ )
  312. X            {
  313. X            int po,no;        /* offsets for left and right colums in 3x3 */
  314. X            register xel *ip0, *ip1, *ip2, *op;
  315. X
  316. X            if (row == 0)
  317. X                {
  318. X                irow0 = irow1;
  319. X                pnm_readpnmrow( ifp, irow1, cols, maxval, format );
  320. X                }
  321. X            if (row == (rows-1))
  322. X                irow2 = irow1;
  323. X            else
  324. X                pnm_readpnmrow( ifp, irow2, cols, maxval, format );
  325. X
  326. X            for (col = cols-1,po= col>0?1:0,no=0,ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
  327. X                 col >= 0;
  328. X                 col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0)
  329. X                {
  330. X                /* grab 3x3 pixel values */
  331. X                pr[0] = PPM_GETR( *ip1 );
  332. X                pg[0] = PPM_GETG( *ip1 );
  333. X                pb[0] = PPM_GETB( *ip1 );
  334. X                pr[1] = PPM_GETR( *(ip1-no) );
  335. X                pg[1] = PPM_GETG( *(ip1-no) );
  336. X                pb[1] = PPM_GETB( *(ip1-no) );
  337. X                pr[5] = PPM_GETR( *(ip1+po) );
  338. X                pg[5] = PPM_GETG( *(ip1+po) );
  339. X                pb[5] = PPM_GETB( *(ip1+po) );
  340. X                pr[3] = PPM_GETR( *(ip2) );
  341. X                pg[3] = PPM_GETG( *(ip2) );
  342. X                pb[3] = PPM_GETB( *(ip2) );
  343. X                pr[2] = PPM_GETR( *(ip2-no) );
  344. X                pg[2] = PPM_GETG( *(ip2-no) );
  345. X                pb[2] = PPM_GETB( *(ip2-no) );
  346. X                pr[4] = PPM_GETR( *(ip2+po) );
  347. X                pg[4] = PPM_GETG( *(ip2+po) );
  348. X                pb[4] = PPM_GETB( *(ip2+po) );
  349. X                pr[6] = PPM_GETR( *(ip0+po) );
  350. X                pg[6] = PPM_GETG( *(ip0+po) );
  351. X                pb[6] = PPM_GETB( *(ip0+po) );
  352. X                pr[8] = PPM_GETR( *(ip0-no) );
  353. X                pg[8] = PPM_GETG( *(ip0-no) );
  354. X                pb[8] = PPM_GETB( *(ip0-no) );
  355. X                pr[7] = PPM_GETR( *(ip0) );
  356. X                pg[7] = PPM_GETG( *(ip0) );
  357. X                pb[7] = PPM_GETB( *(ip0) );
  358. X                r = (*atfunc)(pr);
  359. X                g = (*atfunc)(pg);
  360. X                b = (*atfunc)(pb);
  361. X                PPM_ASSIGN( *op, r, g, b );
  362. X                }
  363. X            pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
  364. X            if (irow1 == irows[2])
  365. X                {
  366. X                irow1 = irows[0];
  367. X                irow2 = irows[1];
  368. X                irow0 = irows[2];
  369. X                }
  370. X            else if (irow1 == irows[1])
  371. X                {
  372. X                irow2 = irows[0];
  373. X                irow0 = irows[1];
  374. X                irow1 = irows[2];
  375. X                }
  376. X            else    /* must be at irows[0] */
  377. X                {
  378. X                irow0 = irows[0];
  379. X                irow1 = irows[1];
  380. X                irow2 = irows[2];
  381. X                }
  382. X            }
  383. X        }
  384. X    else    /* Else must be PGM */
  385. X        {
  386. X        xel *irows[3];
  387. X        xel *irow0, *irow1, *irow2, *orow;
  388. X        int p[9];        /* 3x3 neighbor pixel values */
  389. X        int pv;
  390. X        int promote;
  391. X
  392. X        irows[0] = pnm_allocrow( cols );
  393. X        irows[1] = pnm_allocrow( cols );
  394. X        irows[2] = pnm_allocrow( cols );
  395. X        irow0 = irows[0];
  396. X        irow1 = irows[1];
  397. X        irow2 = irows[2];
  398. X        orow = pnm_allocrow( cols );
  399. X        /* we scale maxval to omaxval */
  400. X        promote = ( PNM_FORMAT_TYPE(format) != PNM_FORMAT_TYPE(oformat) );
  401. X
  402. X        for ( row = 0; row < rows; row++ )
  403. X            {
  404. X            int po,no;        /* offsets for left and right colums in 3x3 */
  405. X            register xel *ip0, *ip1, *ip2, *op;
  406. X
  407. X            if (row == 0)
  408. X                {
  409. X                irow0 = irow1;
  410. X                pnm_readpnmrow( ifp, irow1, cols, maxval, format );
  411. X                if ( promote )
  412. X                    pnm_promoteformatrow( irow1, cols, maxval, format, maxval, oformat );
  413. X                }
  414. X            if (row == (rows-1))
  415. X                irow2 = irow1;
  416. X            else
  417. X                {
  418. X                pnm_readpnmrow( ifp, irow2, cols, maxval, format );
  419. X                if ( promote )
  420. X                    pnm_promoteformatrow( irow2, cols, maxval, format, maxval, oformat );
  421. X                }
  422. X
  423. X            for (col = cols-1,po= col>0?1:0,no=0,ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
  424. X                 col >= 0;
  425. X                 col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0)
  426. X                {
  427. X                /* grab 3x3 pixel values */
  428. X                p[0] = PNM_GET1( *ip1 );
  429. X                p[1] = PNM_GET1( *(ip1-no) );
  430. X                p[5] = PNM_GET1( *(ip1+po) );
  431. X                p[3] = PNM_GET1( *(ip2) );
  432. X                p[2] = PNM_GET1( *(ip2-no) );
  433. X                p[4] = PNM_GET1( *(ip2+po) );
  434. X                p[6] = PNM_GET1( *(ip0+po) );
  435. X                p[8] = PNM_GET1( *(ip0-no) );
  436. X                p[7] = PNM_GET1( *(ip0) );
  437. X                pv = (*atfunc)(p);
  438. X                PNM_ASSIGN1( *op, pv );
  439. X                }
  440. X            pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
  441. X            if (irow1 == irows[2])
  442. X                {
  443. X                irow1 = irows[0];
  444. X                irow2 = irows[1];
  445. X                irow0 = irows[2];
  446. X                }
  447. X            else if (irow1 == irows[1])
  448. X                {
  449. X                irow2 = irows[0];
  450. X                irow0 = irows[1];
  451. X                irow1 = irows[2];
  452. X                }
  453. X            else    /* must be at irows[0] */
  454. X                {
  455. X                irow0 = irows[0];
  456. X                irow1 = irows[1];
  457. X                irow2 = irows[2];
  458. X                }
  459. X            }
  460. X        }
  461. X    pm_close( ifp );
  462. X
  463. X    exit( 0 );
  464. X    }
  465. X
  466. X#define MXIVAL PPM_MAXMAXVAL    /* maximum input value */
  467. X#define NOIVAL (MXIVAL + 1)        /* number of possible input values */
  468. X
  469. X#define SCALEB 8                /* scale bits */
  470. X#define SCALE (1 << SCALEB)    /* scale factor */
  471. X#define MXSVAL (MXIVAL * SCALE)    /* maximum scaled values */
  472. X
  473. X#define CSCALEB 2                /* coarse scale bits */
  474. X#define CSCALE (1 << CSCALEB)    /* coarse scale factor */
  475. X#define MXCSVAL (MXIVAL * CSCALE)    /* maximum coarse scaled values */
  476. X#define NOCSVAL (MXCSVAL + 1)    /* number of coarse scaled values */
  477. X#define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB))    /* convert from scaled to coarse scaled */
  478. X#define CSCTOSC(x) ((x) << (SCALEB - CSCALEB))    /* convert from course scaled to scaled */
  479. X
  480. X#ifndef MAXINT
  481. X# define MAXINT    0x7fffffff    /* assume this is a 32 bit machine */
  482. X#endif
  483. X
  484. X/* round and scale floating point to scaled integer */
  485. X#define ROUND(x) ((int)(((x) * (double)SCALE) + 0.5))
  486. X/* round and un-scale scaled integer value */
  487. X#define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB)    /* rounded un-scale */
  488. X#define UNSCALE(x) ((x) >> SCALEB)
  489. X
  490. X
  491. X/* We restrict radius to the values: 0.333333 <= radius <= 1.0 */
  492. X/* so that no fewer and no more than a 3x3 grid of pixels around */
  493. X/* the pixel in question needs to be read. Given this, we only */
  494. X/* need 3 or 4 weightings per hexagon, as follows: */
  495. X/*                  _ _                         */
  496. X/* Virtical hex:   |_|_|  1 2                   */
  497. X/*                 |X|_|  0 3                   */
  498. X/*                                       _      */
  499. X/*              _                      _|_|   1 */
  500. X/* Middle hex: |_| 1  Horizontal hex: |X|_| 0 2 */
  501. X/*             |X| 0                    |_|   3 */
  502. X/*             |_| 2                            */
  503. X
  504. X/* all filters */
  505. Xint V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL];    /* vertical hex */
  506. Xint M0[NOIVAL],M1[NOIVAL],M2[NOIVAL];            /* middle hex */
  507. Xint H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL];    /* horizontal hex */
  508. X
  509. X/* alpha trimmed and edge enhancement only */
  510. Xint ALFRAC[NOIVAL * 8];            /* fractional alpha divider table */
  511. X
  512. X/* optimal estimation only */
  513. Xint AVEDIV[7 * NOCSVAL];        /* divide by 7 to give average value */
  514. Xint SQUARE[2 * NOCSVAL];        /* scaled square lookup table */
  515. X
  516. X/* Table initialisation function - return alpha range */
  517. Xint 
  518. Xatfilt_setup(alpha,radius,maxscale)
  519. Xdouble alpha,radius,maxscale;    /* alpha, radius and amount to scale input pixel values */
  520. X    {
  521. X    /* other function variables */
  522. X    int alpharange;            /* alpha range value 0 - 3 */
  523. X    double meanscale;        /* scale for finding mean */
  524. X    double mmeanscale;        /* scale for finding mean - midle hex */
  525. X    double alphafraction;    /* fraction of next largest/smallest to subtract from sum */
  526. X
  527. X    /* do setup */
  528. X
  529. X    if (alpha >= 0.0 && alpha < 1.0)    /* alpha trimmed mean */
  530. X        {
  531. X        double noinmean;
  532. X        /* number of elements (out of a possible 7) used in the mean */
  533. X        noinmean = ((0.5 - alpha) * 12.0) + 1.0;
  534. X        mmeanscale = meanscale = maxscale/noinmean;
  535. X        if (alpha == 0.0)                 /* mean filter */
  536. X            {
  537. X            alpharange = 0;
  538. X            alphafraction = 0.0;        /* not used */
  539. X            }
  540. X        else if (alpha < (1.0/6.0))     /* mean of 5 to 7 middle values */
  541. X            {
  542. X            alpharange = 1;
  543. X            alphafraction = (7.0 - noinmean)/2.0;
  544. X            }
  545. X        else if (alpha < (1.0/3.0))     /* mean of 3 to 5 middle values */
  546. X            {
  547. X            alpharange = 2;
  548. X            alphafraction = (5.0 - noinmean)/2.0;
  549. X            }
  550. X        else                            /* mean of 1 to 3 middle values */
  551. X            {                            /* alpha == 0.5 == median filter */
  552. X            alpharange = 3;
  553. X            alphafraction = (3.0 - noinmean)/2.0;
  554. X            }
  555. X        }
  556. X    else if (alpha > 0.5)    /* optimal estimation - alpha controls noise variance threshold. */
  557. X        {
  558. X        int i;
  559. X        double noinmean = 7.0;
  560. X        alpharange = 5;            /* edge enhancement function */
  561. X        alpha -= 1.0;            /* normalise it to 0.0 -> 1.0 */
  562. X        mmeanscale = meanscale = maxscale;    /* compute scaled hex values */
  563. X        alphafraction = 1.0/noinmean;    /* Set up 1:1 division lookup - not used */
  564. X        noisevariance = alpha * (double)omaxval;
  565. X        noisevariance = noisevariance * noisevariance / 8.0;    /* estimate of noise variance */
  566. X        /* set yp optimal estimation specific stuff */
  567. X        for (i=0; i < (7 * NOCSVAL); i++)    /* divide scaled value by 7 lookup */
  568. X            {
  569. X            AVEDIV[i] = CSCTOSC(i)/7;    /* scaled divide by 7 */
  570. X            }
  571. X        for (i=0; i < (2 * NOCSVAL); i++)  /* compute square and rescale by (val >> (2 * SCALEB + 2)) table */
  572. X            {
  573. X            int val;
  574. X            val = CSCTOSC(i - NOCSVAL); /* NOCSVAL offset to cope with -ve input values */
  575. X            SQUARE[i] = (val * val) >> (2 * SCALEB + 2);
  576. X            }
  577. X        }
  578. X    else    /* edge enhancement function */
  579. X        {
  580. X        alpharange = 4;            /* edge enhancement function */
  581. X        alpha = -alpha;            /* turn it the right way up */
  582. X        meanscale = maxscale * (-alpha/((1.0 - alpha) * 7.0)); /* mean of 7 and scaled by -alpha/(1-alpha) */
  583. X        mmeanscale = maxscale * (1.0/(1.0 - alpha) + meanscale);    /* middle pixel has 1/(1-alpha) as well */
  584. X        alphafraction = 0.0;    /* not used */
  585. X        }
  586. X    
  587. X        /* Setup pixel weighting tables - note we pre-compute mean division here too. */
  588. X        {
  589. X        int i;
  590. X        double hexhoff,hexvoff;
  591. X        double tabscale,mtabscale;
  592. X        double v0,v1,v2,v3,m0,m1,m2,me0,me1,me2,h0,h1,h2,h3;
  593. X
  594. X        hexhoff = radius/2;                 /* horizontal offset of virtical hex centers */
  595. X        hexvoff = 3.0 * radius/sqrt(12.0);    /* virtical offset of virtical hex centers */
  596. X        /* scale tables to normalise by hexagon area, and number of hexes used in mean */
  597. X        tabscale = meanscale / (radius * hexvoff);
  598. X        mtabscale = mmeanscale / (radius * hexvoff);
  599. X        v0 = hex_area(0.0,0.0,hexhoff,hexvoff,radius) * tabscale;
  600. X        v1 = hex_area(0.0,1.0,hexhoff,hexvoff,radius) * tabscale;
  601. X        v2 = hex_area(1.0,1.0,hexhoff,hexvoff,radius) * tabscale;
  602. X        v3 = hex_area(1.0,0.0,hexhoff,hexvoff,radius) * tabscale;
  603. X        m0 = hex_area(0.0,0.0,0.0,0.0,radius) * mtabscale;
  604. X        m1 = hex_area(0.0,1.0,0.0,0.0,radius) * mtabscale;
  605. X        m2 = hex_area(0.0,-1.0,0.0,0.0,radius) * mtabscale;
  606. X        h0 = hex_area(0.0,0.0,radius,0.0,radius) * tabscale;
  607. X        h1 = hex_area(1.0,1.0,radius,0.0,radius) * tabscale;
  608. X        h2 = hex_area(1.0,0.0,radius,0.0,radius) * tabscale;
  609. X        h3 = hex_area(1.0,-1.0,radius,0.0,radius) * tabscale;
  610. X
  611. X        for (i=0; i <= MXIVAL; i++)
  612. X            {
  613. X            double fi;
  614. X            fi = (double)i;
  615. X            V0[i] = ROUND(fi * v0);
  616. X            V1[i] = ROUND(fi * v1);
  617. X            V2[i] = ROUND(fi * v2);
  618. X            V3[i] = ROUND(fi * v3);
  619. X            M0[i] = ROUND(fi * m0);
  620. X            M1[i] = ROUND(fi * m1);
  621. X            M2[i] = ROUND(fi * m2);
  622. X            H0[i] = ROUND(fi * h0);
  623. X            H1[i] = ROUND(fi * h1);
  624. X            H2[i] = ROUND(fi * h2);
  625. X            H3[i] = ROUND(fi * h3);
  626. X            }
  627. X        /* set up alpha fraction lookup table used on big/small */
  628. X        for (i=0; i < (NOIVAL * 8); i++)
  629. X            {
  630. X            ALFRAC[i] = ROUND((double)i * alphafraction);
  631. X            }
  632. X        }
  633. X    return alpharange;
  634. X    }
  635. X    
  636. X
  637. X/* Core pixel processing function - hand it 3x3 pixels and return result. */
  638. X/* Mean filter */
  639. Xint
  640. Xatfilt0(p)
  641. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  642. X    {
  643. X    int retv;
  644. X    /* map to scaled hexagon values */
  645. X    retv = M0[p[0]] + M1[p[3]] + M2[p[7]];
  646. X    retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  647. X    retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  648. X    retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  649. X    retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  650. X    retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  651. X    retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  652. X    return UNSCALE(retv);
  653. X    }
  654. X
  655. X/* Mean of 5 - 7 middle values */
  656. Xint
  657. Xatfilt1(p)
  658. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  659. X    {
  660. X    int h0,h1,h2,h3,h4,h5,h6;    /* hexagon values    2 3   */
  661. X                                /*                  1 0 4  */
  662. X                                /*                   6 5   */
  663. X    int big,small;
  664. X    /* map to scaled hexagon values */
  665. X    h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  666. X    h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  667. X    h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  668. X    h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  669. X    h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  670. X    h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  671. X    h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  672. X    /* sum values and also discover the largest and smallest */
  673. X    big = small = h0;
  674. X#define CHECK(xx) \
  675. X    h0 += xx; \
  676. X    if (xx > big) \
  677. X        big = xx; \
  678. X    else if (xx < small) \
  679. X        small = xx;
  680. X    CHECK(h1)
  681. X    CHECK(h2)
  682. X    CHECK(h3)
  683. X    CHECK(h4)
  684. X    CHECK(h5)
  685. X    CHECK(h6)
  686. X#undef CHECK
  687. X    /* Compute mean of middle 5-7 values */
  688. X    return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]);
  689. X    }
  690. X
  691. X/* Mean of 3 - 5 middle values */
  692. Xint
  693. Xatfilt2(p)
  694. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  695. X    {
  696. X    int h0,h1,h2,h3,h4,h5,h6;    /* hexagon values    2 3   */
  697. X                                /*                  1 0 4  */
  698. X                                /*                   6 5   */
  699. X    int big0,big1,small0,small1;
  700. X    /* map to scaled hexagon values */
  701. X    h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  702. X    h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  703. X    h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  704. X    h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  705. X    h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  706. X    h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  707. X    h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  708. X    /* sum values and also discover the 2 largest and 2 smallest */
  709. X    big0 = small0 = h0;
  710. X    small1 = MAXINT;
  711. X    big1 = 0;
  712. X#define CHECK(xx) \
  713. X    h0 += xx; \
  714. X    if (xx > big1) \
  715. X        { \
  716. X        if (xx > big0) \
  717. X            { \
  718. X            big1 = big0; \
  719. X            big0 = xx; \
  720. X            } \
  721. X        else \
  722. X            big1 = xx; \
  723. X        } \
  724. X    if (xx < small1) \
  725. X        { \
  726. X        if (xx < small0) \
  727. X            { \
  728. X            small1 = small0; \
  729. X            small0 = xx; \
  730. X            } \
  731. X        else \
  732. X            small1 = xx; \
  733. X        }
  734. X    CHECK(h1)
  735. X    CHECK(h2)
  736. X    CHECK(h3)
  737. X    CHECK(h4)
  738. X    CHECK(h5)
  739. X    CHECK(h6)
  740. X#undef CHECK
  741. X    /* Compute mean of middle 3-5 values */
  742. X    return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]);
  743. X    }
  744. X
  745. X/* Mean of 1 - 3 middle values. If only 1 value, then this is a median filter. */
  746. Xint
  747. Xatfilt3(p)
  748. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  749. X    {
  750. X    int h0,h1,h2,h3,h4,h5,h6;    /* hexagon values    2 3   */
  751. X                                /*                  1 0 4  */
  752. X                                /*                   6 5   */
  753. X    int big0,big1,big2,small0,small1,small2;
  754. X    /* map to scaled hexagon values */
  755. X    h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  756. X    h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  757. X    h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  758. X    h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  759. X    h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  760. X    h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  761. X    h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  762. X    /* sum values and also discover the 3 largest and 3 smallest */
  763. X    big0 = small0 = h0;
  764. X    small1 = small2 = MAXINT;
  765. X    big1 = big2 = 0;
  766. X#define CHECK(xx) \
  767. X    h0 += xx; \
  768. X    if (xx > big2) \
  769. X        { \
  770. X        if (xx > big1) \
  771. X            { \
  772. X            if (xx > big0) \
  773. X                { \
  774. X                big2 = big1; \
  775. X                big1 = big0; \
  776. X                big0 = xx; \
  777. X                } \
  778. X            else \
  779. X                { \
  780. X                big2 = big1; \
  781. X                big1 = xx; \
  782. X                } \
  783. X            } \
  784. X        else \
  785. X            big2 = xx; \
  786. X        } \
  787. X    if (xx < small2) \
  788. X        { \
  789. X        if (xx < small1) \
  790. X            { \
  791. X            if (xx < small0) \
  792. X                { \
  793. X                small2 = small1; \
  794. X                small1 = small0; \
  795. X                small0 = xx; \
  796. X                } \
  797. X            else \
  798. X                { \
  799. X                small2 = small1; \
  800. X                small1 = xx; \
  801. X                } \
  802. X            } \
  803. X        else \
  804. X            small2 = xx; \
  805. X        }
  806. X    CHECK(h1)
  807. X    CHECK(h2)
  808. X    CHECK(h3)
  809. X    CHECK(h4)
  810. X    CHECK(h5)
  811. X    CHECK(h6)
  812. X#undef CHECK
  813. X    /* Compute mean of middle 1-3 values */
  814. X    return  UNSCALE(h0 -big0 -big1 -small0 -small1 -ALFRAC[(big2 + small2)>>SCALEB]);
  815. X    }
  816. X
  817. X/* Edge enhancement */
  818. X/* notice we use the global omaxval */
  819. Xint
  820. Xatfilt4(p)
  821. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  822. X    {
  823. X    int hav;
  824. X    /* map to scaled hexagon values and compute enhance value */
  825. X    hav = M0[p[0]] + M1[p[3]] + M2[p[7]];
  826. X    hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  827. X    hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  828. X    hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  829. X    hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  830. X    hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  831. X    hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  832. X    if (hav < 0)
  833. X        hav = 0;
  834. X    hav = UNSCALE(hav);
  835. X    if (hav > omaxval)
  836. X        hav = omaxval;
  837. X    return hav;
  838. X    }
  839. X
  840. X/* Optimal estimation - do smoothing in inverse proportion */
  841. X/* to the local variance. */
  842. X/* notice we use the globals noisevariance and omaxval*/
  843. Xint
  844. Xatfilt5(p)
  845. Xint *p;        /* 9 pixel values from 3x3 neighbors */
  846. X    {
  847. X    int mean,variance,temp;
  848. X    int h0,h1,h2,h3,h4,h5,h6;    /* hexagon values    2 3   */
  849. X                                /*                  1 0 4  */
  850. X                                /*                   6 5   */
  851. X    /* map to scaled hexagon values */
  852. X    h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
  853. X    h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
  854. X    h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
  855. X    h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
  856. X    h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
  857. X    h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
  858. X    h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
  859. X    mean = h0 + h1 + h2 + h3 + h4 + h5 + h6;
  860. X    mean = AVEDIV[SCTOCSC(mean)];    /* compute scaled mean by dividing by 7 */
  861. X    temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)];     /* compute scaled variance */
  862. X    temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* and rescale to keep */
  863. X    temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* within 32 bit limits */
  864. X    temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  865. X    temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  866. X    temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
  867. X    temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];    /* (temp = h0 - mean) */
  868. X    if (variance != 0)    /* avoid possible divide by 0 */
  869. X        temp = mean + (variance * temp) / (variance + noisevariance);    /* optimal estimate */
  870. X    else temp = h0;
  871. X    if (temp < 0)
  872. X        temp = 0;
  873. X    temp = RUNSCALE(temp);
  874. X    if (temp > omaxval)
  875. X        temp = omaxval;
  876. X    return temp;
  877. X    }
  878. X
  879. X/* ************************************************** */
  880. X/* Hexagon intersecting square area functions */
  881. X/* Compute the area of the intersection of a triangle */
  882. X/* and a rectangle */
  883. X
  884. Xdouble triang_area();
  885. Xdouble rectang_area();
  886. X
  887. X/* Triangle orientation is per geometric axes (not graphical axies) */
  888. X
  889. X#define NW 0    /* North west triangle /| */
  890. X#define NE 1    /* North east triangle |\ */
  891. X#define SW 2    /* South west triangle \| */
  892. X#define SE 3    /* South east triangle |/ */
  893. X#define STH 2
  894. X#define EST 1
  895. X
  896. X#define SWAPI(a,b) (t = a, a = -b, b = -t)
  897. X
  898. X/* compute the area of overlap of a hexagon diameter d, */
  899. X/* centered at hx,hy, with a unit square of center sx,sy. */
  900. Xdouble
  901. Xhex_area(sx,sy,hx,hy,d)
  902. Xdouble sx,sy;    /* square center */
  903. Xdouble hx,hy,d;    /* hexagon center and diameter */
  904. X    {
  905. X    double hx0,hx1,hx2,hy0,hy1,hy2,hy3;
  906. X    double sx0,sx1,sy0,sy1;
  907. X
  908. X    /* compute square co-ordinates */
  909. X    sx0 = sx - 0.5;
  910. X    sy0 = sy - 0.5;
  911. X    sx1 = sx + 0.5;
  912. X    sy1 = sy + 0.5;
  913. X    /* compute hexagon co-ordinates */
  914. X    hx0 = hx - d/2.0;
  915. X    hx1 = hx;
  916. X    hx2 = hx + d/2.0;
  917. X    hy0 = hy - 0.5773502692 * d;     /* d / sqrt(3) */
  918. X    hy1 = hy - 0.2886751346 * d;     /* d / sqrt(12) */
  919. X    hy2 = hy + 0.2886751346 * d;     /* d / sqrt(12) */
  920. X    hy3 = hy + 0.5773502692 * d;     /* d / sqrt(3) */
  921. X
  922. X    return 
  923. X        triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) + 
  924. X        triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) + 
  925. X        rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) + 
  926. X        triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) + 
  927. X        triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE);
  928. X    }
  929. X
  930. Xdouble
  931. Xtriang_area(rx0,ry0,rx1,ry1,tx0,ty0,tx1,ty1,tt)
  932. Xdouble rx0,ry0,rx1,ry1;        /* rectangle boundaries */
  933. Xdouble tx0,ty0,tx1,ty1;        /* triangle boundaries */
  934. Xint tt;                        /* triangle type */
  935. X    {
  936. X    double a,b,c,d;
  937. X    double lx0,ly0,lx1,ly1;
  938. X    /* Convert everything to a NW triangle */
  939. X    if (tt & STH)
  940. X        {
  941. X        double t;
  942. X        SWAPI(ry0,ry1);
  943. X        SWAPI(ty0,ty1);
  944. X        }
  945. X    if (tt & EST)
  946. X        {
  947. X        double t;
  948. X        SWAPI(rx0,rx1);
  949. X        SWAPI(tx0,tx1);
  950. X        }
  951. X    /* Compute overlapping box */
  952. X    if (tx0 > rx0)
  953. X        rx0 = tx0;
  954. X    if (ty0 > ry0)
  955. X        ry0 = ty0;
  956. X    if (tx1 < rx1)
  957. X        rx1 = tx1;
  958. X    if (ty1 < ry1)
  959. X        ry1 = ty1;
  960. X    if (rx1 <= rx0 || ry1 <= ry0)
  961. X        return 0.0;
  962. X    /* Need to compute diagonal line intersection with the box */
  963. X    /* First compute co-efficients to formulas x = a + by and y = c + dx */
  964. X    b = (tx1 - tx0)/(ty1 - ty0);
  965. X    a = tx0 - b * ty0;
  966. X    d = (ty1 - ty0)/(tx1 - tx0);
  967. X    c = ty0 - d * tx0;
  968. X
  969. X    /* compute top or right intersection */
  970. X    tt = 0;
  971. X    ly1 = ry1; 
  972. X    lx1 = a + b * ly1;
  973. X    if (lx1 <= rx0)
  974. X        return (rx1 - rx0) * (ry1 - ry0);
  975. X    else if (lx1 > rx1)    /* could be right hand side */
  976. X        {
  977. X        lx1 = rx1;
  978. X        ly1 = c + d * lx1;
  979. X        if (ly1 <= ry0)
  980. X            return (rx1 - rx0) * (ry1 - ry0);
  981. X        tt = 1;    /* right hand side intersection */
  982. X        }
  983. X    /* compute left or bottom intersection */
  984. X    lx0 = rx0;
  985. X    ly0 = c + d * lx0;
  986. X    if (ly0 >= ry1)
  987. X        return (rx1 - rx0) * (ry1 - ry0);
  988. X    else if (ly0 < ry0)    /* could be right hand side */
  989. X        {
  990. X        ly0 = ry0; 
  991. X        lx0 = a + b * ly0;
  992. X        if (lx0 >= rx1)
  993. X            return (rx1 - rx0) * (ry1 - ry0);
  994. X        tt |= 2;    /* bottom intersection */
  995. X        }
  996. X    
  997. X    if (tt == 0)    /* top and left intersection */
  998. X        {    /* rectangle minus triangle */
  999. X        return ((rx1 - rx0) * (ry1 - ry0))
  1000. X             - (0.5 * (lx1 - rx0) * (ry1 - ly0));
  1001. X        }
  1002. X    else if (tt == 1)    /* right and left intersection */
  1003. X        {
  1004. X        return ((rx1 - rx0) * (ly0 - ry0))
  1005. X             + (0.5 * (rx1 - rx0) * (ly1 - ly0));
  1006. X        }
  1007. X    else if (tt == 2)    /* top and bottom intersection */
  1008. X        {
  1009. X        return ((rx1 - lx1) * (ry1 - ry0))
  1010. X             + (0.5 * (lx1 - lx0) * (ry1 - ry0));
  1011. X        }
  1012. X    else /* tt == 3 */    /* right and bottom intersection */
  1013. X        {    /* triangle */
  1014. X        return (0.5 * (rx1 - lx0) * (ly1 - ry0));
  1015. X        }
  1016. X    }
  1017. X
  1018. X/* Compute rectangle area */
  1019. Xdouble
  1020. Xrectang_area(rx0,ry0,rx1,ry1,tx0,ty0,tx1,ty1)
  1021. Xdouble rx0,ry0,rx1,ry1;        /* rectangle boundaries */
  1022. Xdouble tx0,ty0,tx1,ty1;        /* rectangle boundaries */
  1023. X    {
  1024. X    /* Compute overlapping box */
  1025. X    if (tx0 > rx0)
  1026. X        rx0 = tx0;
  1027. X    if (ty0 > ry0)
  1028. X        ry0 = ty0;
  1029. X    if (tx1 < rx1)
  1030. X        rx1 = tx1;
  1031. X    if (ty1 < ry1)
  1032. X        ry1 = ty1;
  1033. X    if (rx1 <= rx0 || ry1 <= ry0)
  1034. X        return 0.0;
  1035. X    return (rx1 - rx0) * (ry1 - ry0);
  1036. X    }
  1037. X
  1038. END_OF_FILE
  1039.   if test 25912 -ne `wc -c <'pnmnlfilt.c'`; then
  1040.     echo shar: \"'pnmnlfilt.c'\" unpacked with wrong size!
  1041.   fi
  1042.   # end of 'pnmnlfilt.c'
  1043. fi
  1044. echo shar: End of archive 1 \(of 1\).
  1045. cp /dev/null ark1isdone
  1046. MISSING=""
  1047. for I in 1 ; do
  1048.     if test ! -f ark${I}isdone ; then
  1049.     MISSING="${MISSING} ${I}"
  1050.     fi
  1051. done
  1052. if test "${MISSING}" = "" ; then
  1053.     echo You have the archive.
  1054.     rm -f ark[1-9]isdone
  1055. else
  1056.     echo You still must unpack the following archives:
  1057.     echo "        " ${MISSING}
  1058. fi
  1059. exit 0
  1060. exit 0 # Just in case...
  1061.