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 sythesized 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
52 #include "pm_c_util.h"
53 #include "pnm.h"
54
55 struct cmdlineInfo {
56 const char * inputFileName;
57 double alpha;
58 double radius;
59 };
60
61
62 static void
parseCommandLine(int argc,char ** argv,struct cmdlineInfo * const cmdlineP)63 parseCommandLine(int argc,
64 char ** argv,
65 struct cmdlineInfo * const cmdlineP) {
66
67 if (argc-1 < 2)
68 pm_error("You must specify at least two arguments: alpha and radius. "
69 "You specified %u", argc-1);
70
71 if (sscanf(argv[1], "%lf", &cmdlineP->alpha) != 1)
72 pm_error("Invalid alpha (1st) argument '%s'. "
73 "Must be a decimal number",
74 argv[1]);
75
76 if (sscanf( argv[2], "%lf", &cmdlineP->radius ) != 1)
77 pm_error("Invalid radius (2nd) argument '%s'. "
78 "Must be a decimal number",
79 argv[2]);
80
81 if ((cmdlineP->alpha > -0.1 && cmdlineP->alpha < 0.0) ||
82 (cmdlineP->alpha > 0.5 && cmdlineP->alpha < 1.0))
83 pm_error( "Alpha must be in range 0.0 <= alpha <= 0.5 "
84 "for alpha trimmed mean. "
85 "You specified %f", cmdlineP->alpha);
86 if (cmdlineP->alpha > 2.0)
87 pm_error("Alpha must be in range 1.0 <= cmdlineP->alpha <= 2.0 "
88 "for optimal estimation. You specified %f",
89 cmdlineP->alpha);
90
91 if (cmdlineP->alpha < -0.9 ||
92 (cmdlineP->alpha > -0.1 && cmdlineP->alpha < 0.0))
93 pm_error( "Alpha must be in range -0.9 <= alpha <= -0.1 "
94 "for edge enhancement. You specified %f",
95 cmdlineP->alpha);
96
97 if (cmdlineP->radius < 1.0/3 || cmdlineP->radius > 1.0)
98 pm_error("Radius must be in range 1/3 <= radius <= 1. "
99 "You specified %f", cmdlineP->radius);
100
101 if (argc-1 < 3)
102 cmdlineP->inputFileName = "-";
103 else
104 cmdlineP->inputFileName = argv[3];
105
106 if (argc-1 > 3)
107 pm_error("Too many arguments: %u. The most allowed are 3: alpha, "
108 "radius, and file name", argc-1);
109 }
110
111
112 /* MXIVAL is the maximum input sample value we can handle.
113 It is limited by our willingness to allocate storage in various arrays
114 that are indexed by sample values.
115
116 We use PPM_MAXMAXVAL because that used to be the maximum possible
117 sample value in the format, and most images still limit themselves to
118 this value.
119 */
120
121 #define MXIVAL PPM_MAXMAXVAL
122
123 xelval omaxval;
124 /* global so that pixel processing code can get at it quickly */
125 int noisevariance;
126 /* global so that pixel processing code can get at it quickly */
127
128 /*
129 * Declared static here rather than passing a jillion options in the call to
130 * do_one_frame(). Also it makes a huge amount of sense to only malloc the
131 * row buffers once instead of for each frame (with the corresponding free'ing
132 * of course).
133 */
134 static xel *irows[3];
135 static xel *irow0, *irow1, *irow2, *orow;
136 static int rows, cols, format, oformat, row, col;
137 static int (*atfunc)(int *);
138 static xelval maxval;
139
140 #define NOIVAL (MXIVAL + 1) /* number of possible input values */
141
142 #define SCALEB 8 /* scale bits */
143 #define SCALE (1 << SCALEB) /* scale factor */
144 #define MXSVAL (MXIVAL * SCALE) /* maximum scaled values */
145
146 #define CSCALEB 2 /* coarse scale bits */
147 #define CSCALE (1 << CSCALEB) /* coarse scale factor */
148 #define MXCSVAL (MXIVAL * CSCALE) /* maximum coarse scaled values */
149 #define NOCSVAL (MXCSVAL + 1) /* number of coarse scaled values */
150 #define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB)) /* scaled to coarse scaled */
151 #define CSCTOSC(x) ((x) << (SCALEB - CSCALEB)) /* course scaled to scaled */
152
153 #ifndef MAXINT
154 # define MAXINT 0x7fffffff /* assume this is a 32 bit machine */
155 #endif
156
157 /* round and scale floating point to scaled integer */
158 #define ROUNDSCALE(x) ((int)(((x) * (double)SCALE) + 0.5))
159 /* round and un-scale scaled integer value */
160 #define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB)
161 /* rounded un-scale */
162 #define UNSCALE(x) ((x) >> SCALEB)
163
164 static double
sqr(double const arg)165 sqr(double const arg) {
166 return arg * arg;
167 }
168
169
170
171 /* We restrict radius to the values: 0.333333 <= radius <= 1.0 */
172 /* so that no fewer and no more than a 3x3 grid of pixels around */
173 /* the pixel in question needs to be read. Given this, we only */
174 /* need 3 or 4 weightings per hexagon, as follows: */
175 /* _ _ */
176 /* Vertical hex: |_|_| 1 2 */
177 /* |X|_| 0 3 */
178 /* _ */
179 /* _ _|_| 1 */
180 /* Middle hex: |_| 1 Horizontal hex: |X|_| 0 2 */
181 /* |X| 0 |_| 3 */
182 /* |_| 2 */
183
184 /* all filters */
185 int V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL]; /* vertical hex */
186 int M0[NOIVAL],M1[NOIVAL],M2[NOIVAL]; /* middle hex */
187 int H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL]; /* horizontal hex */
188
189 /* alpha trimmed and edge enhancement only */
190 int ALFRAC[NOIVAL * 8]; /* fractional alpha divider table */
191
192 /* optimal estimation only */
193 int AVEDIV[7 * NOCSVAL]; /* divide by 7 to give average value */
194 int SQUARE[2 * NOCSVAL]; /* scaled square lookup table */
195
196 /* ************************************************** *
197 Hexagon intersecting square area functions
198 Compute the area of the intersection of a triangle
199 and a rectangle
200 ************************************************** */
201
202 /* Triangle orientation is per geometric axes (not graphical axies) */
203
204 #define NW 0 /* North west triangle /| */
205 #define NE 1 /* North east triangle |\ */
206 #define SW 2 /* South west triangle \| */
207 #define SE 3 /* South east triangle |/ */
208 #define STH 2
209 #define EST 1
210
211 #define SWAPI(a,b) (t = a, a = -b, b = -t)
212
213 static double
triang_area(double rx0,double ry0,double rx1,double ry1,double tx0,double ty0,double tx1,double ty1,int tt)214 triang_area(double rx0, double ry0, double rx1, double ry1,
215 double tx0, double ty0, double tx1, double ty1,
216 int tt) {
217 /* rx0,ry0,rx1,ry1: rectangle boundaries */
218 /* tx0,ty0,tx1,ty1: triangle boundaries */
219 /* tt: triangle type */
220
221 double a,b,c,d;
222 double lx0,ly0,lx1,ly1;
223
224 /* Convert everything to a NW triangle */
225 if (tt & STH) {
226 double t;
227 SWAPI(ry0,ry1);
228 SWAPI(ty0,ty1);
229 }
230 if (tt & EST) {
231 double t;
232 SWAPI(rx0,rx1);
233 SWAPI(tx0,tx1);
234 }
235 /* Compute overlapping box */
236 if (tx0 > rx0)
237 rx0 = tx0;
238 if (ty0 > ry0)
239 ry0 = ty0;
240 if (tx1 < rx1)
241 rx1 = tx1;
242 if (ty1 < ry1)
243 ry1 = ty1;
244 if (rx1 <= rx0 || ry1 <= ry0)
245 return 0.0;
246
247 /* Need to compute diagonal line intersection with the box */
248 /* First compute co-efficients to formulas x = a + by and y = c + dx */
249 b = (tx1 - tx0)/(ty1 - ty0);
250 a = tx0 - b * ty0;
251 d = (ty1 - ty0)/(tx1 - tx0);
252 c = ty0 - d * tx0;
253
254 /* compute top or right intersection */
255 tt = 0;
256 ly1 = ry1;
257 lx1 = a + b * ly1;
258 if (lx1 <= rx0)
259 return (rx1 - rx0) * (ry1 - ry0);
260 else if (lx1 > rx1) {
261 /* could be right hand side */
262 lx1 = rx1;
263 ly1 = c + d * lx1;
264 if (ly1 <= ry0)
265 return (rx1 - rx0) * (ry1 - ry0);
266 tt = 1; /* right hand side intersection */
267 }
268 /* compute left or bottom intersection */
269 lx0 = rx0;
270 ly0 = c + d * lx0;
271 if (ly0 >= ry1)
272 return (rx1 - rx0) * (ry1 - ry0);
273 else if (ly0 < ry0) {
274 /* could be right hand side */
275 ly0 = ry0;
276 lx0 = a + b * ly0;
277 if (lx0 >= rx1)
278 return (rx1 - rx0) * (ry1 - ry0);
279 tt |= 2; /* bottom intersection */
280 }
281
282 if (tt == 0) {
283 /* top and left intersection */
284 /* rectangle minus triangle */
285 return ((rx1 - rx0) * (ry1 - ry0))
286 - (0.5 * (lx1 - rx0) * (ry1 - ly0));
287 } else if (tt == 1) {
288 /* right and left intersection */
289 return ((rx1 - rx0) * (ly0 - ry0))
290 + (0.5 * (rx1 - rx0) * (ly1 - ly0));
291 } else if (tt == 2) {
292 /* top and bottom intersection */
293 return ((rx1 - lx1) * (ry1 - ry0))
294 + (0.5 * (lx1 - lx0) * (ry1 - ry0));
295 } else {
296 /* tt == 3 */
297 /* right and bottom intersection */
298 /* triangle */
299 return (0.5 * (rx1 - lx0) * (ly1 - ry0));
300 }
301 }
302
303
304
305 static double
rectang_area(double rx0,double ry0,double rx1,double ry1,double tx0,double ty0,double tx1,double ty1)306 rectang_area(double rx0, double ry0, double rx1, double ry1,
307 double tx0, double ty0, double tx1, double ty1) {
308 /* Compute rectangle area */
309 /* rx0,ry0,rx1,ry1: rectangle boundaries */
310 /* tx0,ty0,tx1,ty1: rectangle boundaries */
311
312 /* Compute overlapping box */
313 if (tx0 > rx0)
314 rx0 = tx0;
315 if (ty0 > ry0)
316 ry0 = ty0;
317 if (tx1 < rx1)
318 rx1 = tx1;
319 if (ty1 < ry1)
320 ry1 = ty1;
321 if (rx1 <= rx0 || ry1 <= ry0)
322 return 0.0;
323 return (rx1 - rx0) * (ry1 - ry0);
324 }
325
326
327
328
329 static double
hex_area(double sx,double sy,double hx,double hy,double d)330 hex_area(double sx, double sy, double hx, double hy, double d) {
331 /* compute the area of overlap of a hexagon diameter d, */
332 /* centered at hx,hy, with a unit square of center sx,sy. */
333 /* sx,sy: square center */
334 /* hx,hy,d: hexagon center and diameter */
335
336 double hx0,hx1,hx2,hy0,hy1,hy2,hy3;
337 double sx0,sx1,sy0,sy1;
338
339 /* compute square co-ordinates */
340 sx0 = sx - 0.5;
341 sy0 = sy - 0.5;
342 sx1 = sx + 0.5;
343 sy1 = sy + 0.5;
344
345 /* compute hexagon co-ordinates */
346 hx0 = hx - d/2.0;
347 hx1 = hx;
348 hx2 = hx + d/2.0;
349 hy0 = hy - 0.5773502692 * d; /* d / sqrt(3) */
350 hy1 = hy - 0.2886751346 * d; /* d / sqrt(12) */
351 hy2 = hy + 0.2886751346 * d; /* d / sqrt(12) */
352 hy3 = hy + 0.5773502692 * d; /* d / sqrt(3) */
353
354 return
355 triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) +
356 triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) +
357 rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) +
358 triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) +
359 triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE);
360 }
361
362
363
364
365 static void
setupAvediv(void)366 setupAvediv(void) {
367
368 unsigned int i;
369
370 for (i=0; i < (7 * NOCSVAL); ++i) {
371 /* divide scaled value by 7 lookup */
372 AVEDIV[i] = CSCTOSC(i)/7; /* scaled divide by 7 */
373 }
374
375 }
376
377
378
379
380 static void
setupSquare(void)381 setupSquare(void) {
382
383 unsigned int i;
384
385 for (i=0; i < (2 * NOCSVAL); ++i) {
386 /* compute square and rescale by (val >> (2 * SCALEB + 2)) table */
387 int const val = CSCTOSC(i - NOCSVAL);
388 /* NOCSVAL offset to cope with -ve input values */
389 SQUARE[i] = (val * val) >> (2 * SCALEB + 2);
390 }
391 }
392
393
394
395
396 static void
setup1(double const alpha,double const radius,double const maxscale,int * const alpharangeP,double * const meanscaleP,double * const mmeanscaleP,double * const alphafractionP,int * const noisevarianceP)397 setup1(double const alpha,
398 double const radius,
399 double const maxscale,
400 int * const alpharangeP,
401 double * const meanscaleP,
402 double * const mmeanscaleP,
403 double * const alphafractionP,
404 int * const noisevarianceP) {
405
406
407 setupAvediv();
408 setupSquare();
409
410 if (alpha >= 0.0 && alpha <= 0.5) {
411 /* alpha trimmed mean */
412 double const noinmean = ((0.5 - alpha) * 12.0) + 1.0;
413 /* number of elements (out of a possible 7) used in the mean */
414
415 *mmeanscaleP = *meanscaleP = maxscale/noinmean;
416 if (alpha == 0.0) {
417 /* mean filter */
418 *alpharangeP = 0;
419 *alphafractionP = 0.0; /* not used */
420 } else if (alpha < (1.0/6.0)) {
421 /* mean of 5 to 7 middle values */
422 *alpharangeP = 1;
423 *alphafractionP = (7.0 - noinmean)/2.0;
424 } else if (alpha < (1.0/3.0)) {
425 /* mean of 3 to 5 middle values */
426 *alpharangeP = 2;
427 *alphafractionP = (5.0 - noinmean)/2.0;
428 } else {
429 /* mean of 1 to 3 middle values */
430 /* alpha == 0.5 == median filter */
431 *alpharangeP = 3;
432 *alphafractionP = (3.0 - noinmean)/2.0;
433 }
434 } else if (alpha >= 1.0 && alpha <= 2.0) {
435 /* optimal estimation - alpha controls noise variance threshold. */
436 double const alphaNormalized = alpha - 1.0;
437 /* normalize it to 0.0 -> 1.0 */
438 double const noinmean = 7.0;
439 *alpharangeP = 5; /* edge enhancement function */
440 *mmeanscaleP = *meanscaleP = maxscale; /* compute scaled hex values */
441 *alphafractionP = 1.0/noinmean;
442 /* Set up 1:1 division lookup - not used */
443 *noisevarianceP = sqr(alphaNormalized * omaxval) / 8.0;
444 /* estimate of noise variance */
445 } else if (alpha >= -0.9 && alpha <= -0.1) {
446 /* edge enhancement function */
447 double const posAlpha = -alpha;
448 /* positive alpha value */
449 *alpharangeP = 4; /* edge enhancement function */
450 *meanscaleP = maxscale * (-posAlpha/((1.0 - posAlpha) * 7.0));
451 /* mean of 7 and scaled by -posAlpha/(1-posAlpha) */
452 *mmeanscaleP = maxscale * (1.0/(1.0 - posAlpha) + *meanscaleP);
453 /* middle pixel has 1/(1-posAlpha) as well */
454 *alphafractionP = 0.0; /* not used */
455 } else {
456 /* An entry condition on 'alpha' makes this impossible */
457 pm_error("INTERNAL ERROR: impossible alpha value: %f", alpha);
458 }
459 }
460
461
462
463
464 static void
setupAlfrac(double const alphafraction)465 setupAlfrac(double const alphafraction) {
466 /* set up alpha fraction lookup table used on big/small */
467
468 unsigned int i;
469
470 for (i=0; i < (NOIVAL * 8); ++i) {
471 ALFRAC[i] = ROUNDSCALE(i * alphafraction);
472 }
473 }
474
475
476
477
478 static void
setupPixelWeightingTables(double const radius,double const meanscale,double const mmeanscale)479 setupPixelWeightingTables(double const radius,
480 double const meanscale,
481 double const mmeanscale) {
482
483 /* Setup pixel weighting tables - note we pre-compute mean
484 division here too.
485 */
486 double const hexhoff = radius/2;
487 /* horizontal offset of vertical hex centers */
488 double const hexvoff = 3.0 * radius/sqrt(12.0);
489 /* vertical offset of vertical hex centers */
490
491 double const tabscale = meanscale / (radius * hexvoff);
492 double const mtabscale = mmeanscale / (radius * hexvoff);
493
494 /* scale tables to normalize by hexagon area, and number of
495 hexes used in mean
496 */
497 double const v0 =
498 hex_area(0.0, 0.0, hexhoff, hexvoff, radius) * tabscale;
499 double const v1 =
500 hex_area(0.0, 1.0, hexhoff, hexvoff, radius) * tabscale;
501 double const v2 =
502 hex_area(1.0, 1.0, hexhoff, hexvoff, radius) * tabscale;
503 double const v3 =
504 hex_area(1.0, 0.0, hexhoff, hexvoff, radius) * tabscale;
505 double const m0 =
506 hex_area(0.0, 0.0, 0.0, 0.0, radius) * mtabscale;
507 double const m1 =
508 hex_area(0.0, 1.0, 0.0, 0.0, radius) * mtabscale;
509 double const m2 =
510 hex_area(0.0, -1.0, 0.0, 0.0, radius) * mtabscale;
511 double const h0 =
512 hex_area(0.0, 0.0, radius, 0.0, radius) * tabscale;
513 double const h1 =
514 hex_area(1.0, 1.0, radius, 0.0, radius) * tabscale;
515 double const h2 =
516 hex_area(1.0, 0.0, radius, 0.0, radius) * tabscale;
517 double const h3 =
518 hex_area(1.0, -1.0, radius, 0.0, radius) * tabscale;
519
520 unsigned int i;
521
522 for (i=0; i <= MXIVAL; ++i) {
523 V0[i] = ROUNDSCALE(i * v0);
524 V1[i] = ROUNDSCALE(i * v1);
525 V2[i] = ROUNDSCALE(i * v2);
526 V3[i] = ROUNDSCALE(i * v3);
527 M0[i] = ROUNDSCALE(i * m0);
528 M1[i] = ROUNDSCALE(i * m1);
529 M2[i] = ROUNDSCALE(i * m2);
530 H0[i] = ROUNDSCALE(i * h0);
531 H1[i] = ROUNDSCALE(i * h1);
532 H2[i] = ROUNDSCALE(i * h2);
533 H3[i] = ROUNDSCALE(i * h3);
534 }
535 }
536
537
538
539
540 /* Table initialization function - return alpha range */
541 static int
atfilt_setup(double const alpha,double const radius,double const maxscale)542 atfilt_setup(double const alpha,
543 double const radius,
544 double const maxscale) {
545
546 int alpharange; /* alpha range value 0 - 5 */
547 double meanscale; /* scale for finding mean */
548 double mmeanscale; /* scale for finding mean - midle hex */
549 double alphafraction;
550 /* fraction of next largest/smallest to subtract from sum */
551
552 setup1(alpha, radius, maxscale,
553 &alpharange, &meanscale, &mmeanscale, &alphafraction,
554 &noisevariance);
555
556 setupAlfrac(alphafraction);
557
558 setupPixelWeightingTables(radius, meanscale, mmeanscale);
559
560 return alpharange;
561 }
562
563
564
565 static int
atfilt0(int * p)566 atfilt0(int * p) {
567 /* Core pixel processing function - hand it 3x3 pixels and return result. */
568 /* Mean filter */
569 /* 'p' is 9 pixel values from 3x3 neighbors */
570 int retv;
571 /* map to scaled hexagon values */
572 retv = M0[p[0]] + M1[p[3]] + M2[p[7]];
573 retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
574 retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
575 retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
576 retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
577 retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
578 retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
579 return UNSCALE(retv);
580 }
581
582 #define CHECK(xx) {\
583 h0 += xx; \
584 if (xx > big) \
585 big = xx; \
586 else if (xx < small) \
587 small = xx; }
588
589 static int
atfilt1(int * p)590 atfilt1(int * p) {
591 /* Mean of 5 - 7 middle values */
592 /* 'p' is 9 pixel values from 3x3 neighbors */
593
594 int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
595 /* 1 0 4 */
596 /* 6 5 */
597 int big,small;
598 /* map to scaled hexagon values */
599 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
600 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
601 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
602 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
603 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
604 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
605 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
606 /* sum values and also discover the largest and smallest */
607 big = small = h0;
608 CHECK(h1);
609 CHECK(h2);
610 CHECK(h3);
611 CHECK(h4);
612 CHECK(h5);
613 CHECK(h6);
614 /* Compute mean of middle 5-7 values */
615 return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]);
616 }
617 #undef CHECK
618
619 #define CHECK(xx) {\
620 h0 += xx; \
621 if (xx > big1) {\
622 if (xx > big0) {\
623 big1 = big0; \
624 big0 = xx; \
625 } else \
626 big1 = xx; \
627 } \
628 if (xx < small1) {\
629 if (xx < small0) {\
630 small1 = small0; \
631 small0 = xx; \
632 } else \
633 small1 = xx; \
634 }\
635 }
636
637
638 static int
atfilt2(int * p)639 atfilt2(int *p) {
640 /* Mean of 3 - 5 middle values */
641 /* 'p' is 9 pixel values from 3x3 neighbors */
642 int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
643 /* 1 0 4 */
644 /* 6 5 */
645 int big0,big1,small0,small1;
646 /* map to scaled hexagon values */
647 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
648 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
649 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
650 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
651 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
652 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
653 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
654 /* sum values and also discover the 2 largest and 2 smallest */
655 big0 = small0 = h0;
656 small1 = MAXINT;
657 big1 = 0;
658 CHECK(h1);
659 CHECK(h2);
660 CHECK(h3);
661 CHECK(h4);
662 CHECK(h5);
663 CHECK(h6);
664 /* Compute mean of middle 3-5 values */
665 return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]);
666 }
667
668 #undef CHECK
669
670 #define CHECK(xx) {\
671 h0 += xx; \
672 if (xx > big2) \
673 { \
674 if (xx > big1) \
675 { \
676 if (xx > big0) \
677 { \
678 big2 = big1; \
679 big1 = big0; \
680 big0 = xx; \
681 } \
682 else \
683 { \
684 big2 = big1; \
685 big1 = xx; \
686 } \
687 } \
688 else \
689 big2 = xx; \
690 } \
691 if (xx < small2) \
692 { \
693 if (xx < small1) \
694 { \
695 if (xx < small0) \
696 { \
697 small2 = small1; \
698 small1 = small0; \
699 small0 = xx; \
700 } \
701 else \
702 { \
703 small2 = small1; \
704 small1 = xx; \
705 } \
706 } \
707 else \
708 small2 = xx; \
709 }}
710
711 static int
atfilt3(int * p)712 atfilt3(int *p) {
713 /* Mean of 1 - 3 middle values. If only 1 value, then this is a median
714 filter.
715 */
716 /* 'p' is pixel values from 3x3 neighbors */
717 int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
718 /* 1 0 4 */
719 /* 6 5 */
720 int big0,big1,big2,small0,small1,small2;
721 /* map to scaled hexagon values */
722 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
723 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
724 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
725 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
726 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
727 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
728 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
729 /* sum values and also discover the 3 largest and 3 smallest */
730 big0 = small0 = h0;
731 small1 = small2 = MAXINT;
732 big1 = big2 = 0;
733 CHECK(h1);
734 CHECK(h2);
735 CHECK(h3);
736 CHECK(h4);
737 CHECK(h5);
738 CHECK(h6);
739 /* Compute mean of middle 1-3 values */
740 return UNSCALE(h0 -big0 -big1 -small0 -small1
741 -ALFRAC[(big2 + small2)>>SCALEB]);
742 }
743 #undef CHECK
744
745 static int
atfilt4(int * p)746 atfilt4(int *p) {
747 /* Edge enhancement */
748 /* notice we use the global omaxval */
749 /* 'p' is 9 pixel values from 3x3 neighbors */
750
751 int hav;
752 /* map to scaled hexagon values and compute enhance value */
753 hav = M0[p[0]] + M1[p[3]] + M2[p[7]];
754 hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
755 hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
756 hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
757 hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
758 hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
759 hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
760 if (hav < 0)
761 hav = 0;
762 hav = UNSCALE(hav);
763 if (hav > omaxval)
764 hav = omaxval;
765 return hav;
766 }
767
768 static int
atfilt5(int * p)769 atfilt5(int *p) {
770 /* Optimal estimation - do smoothing in inverse proportion */
771 /* to the local variance. */
772 /* notice we use the globals noisevariance and omaxval*/
773 /* 'p' is 9 pixel values from 3x3 neighbors */
774
775 int mean,variance,temp;
776 int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
777 /* 1 0 4 */
778 /* 6 5 */
779 /* map to scaled hexagon values */
780 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
781 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
782 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
783 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
784 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
785 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
786 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
787 mean = h0 + h1 + h2 + h3 + h4 + h5 + h6;
788 mean = AVEDIV[SCTOCSC(mean)]; /* compute scaled mean by dividing by 7 */
789 temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)];
790 /* compute scaled variance */
791 temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
792 /* and rescale to keep */
793 temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
794 /* within 32 bit limits */
795 temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
796 temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
797 temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
798 temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
799 /* (temp = h0 - mean) */
800 if (variance != 0) /* avoid possible divide by 0 */
801 temp = mean + (variance * temp) / (variance + noisevariance);
802 /* optimal estimate */
803 else temp = h0;
804 if (temp < 0)
805 temp = 0;
806 temp = RUNSCALE(temp);
807 if (temp > omaxval)
808 temp = omaxval;
809 return temp;
810 }
811
812
813
814 static void
do_one_frame(FILE * const ifP)815 do_one_frame(FILE * const ifP) {
816
817 pnm_writepnminit( stdout, cols, rows, omaxval, oformat, 0 );
818
819 if ( PNM_FORMAT_TYPE(oformat) == PPM_TYPE ) {
820 int pr[9],pg[9],pb[9]; /* 3x3 neighbor pixel values */
821 int r,g,b;
822
823 for ( row = 0; row < rows; row++ ) {
824 int po,no; /* offsets for left and right columns in 3x3 */
825 xel *ip0, *ip1, *ip2, *op;
826
827 if (row == 0) {
828 irow0 = irow1;
829 pnm_readpnmrow( ifP, irow1, cols, maxval, format );
830 }
831 if (row == (rows-1))
832 irow2 = irow1;
833 else
834 pnm_readpnmrow( ifP, irow2, cols, maxval, format );
835
836 for (col = cols-1,po= col>0?1:0,no=0,
837 ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
838 col >= 0;
839 col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0) {
840 /* grab 3x3 pixel values */
841 pr[0] = PPM_GETR( *ip1 );
842 pg[0] = PPM_GETG( *ip1 );
843 pb[0] = PPM_GETB( *ip1 );
844 pr[1] = PPM_GETR( *(ip1-no) );
845 pg[1] = PPM_GETG( *(ip1-no) );
846 pb[1] = PPM_GETB( *(ip1-no) );
847 pr[5] = PPM_GETR( *(ip1+po) );
848 pg[5] = PPM_GETG( *(ip1+po) );
849 pb[5] = PPM_GETB( *(ip1+po) );
850 pr[3] = PPM_GETR( *(ip2) );
851 pg[3] = PPM_GETG( *(ip2) );
852 pb[3] = PPM_GETB( *(ip2) );
853 pr[2] = PPM_GETR( *(ip2-no) );
854 pg[2] = PPM_GETG( *(ip2-no) );
855 pb[2] = PPM_GETB( *(ip2-no) );
856 pr[4] = PPM_GETR( *(ip2+po) );
857 pg[4] = PPM_GETG( *(ip2+po) );
858 pb[4] = PPM_GETB( *(ip2+po) );
859 pr[6] = PPM_GETR( *(ip0+po) );
860 pg[6] = PPM_GETG( *(ip0+po) );
861 pb[6] = PPM_GETB( *(ip0+po) );
862 pr[8] = PPM_GETR( *(ip0-no) );
863 pg[8] = PPM_GETG( *(ip0-no) );
864 pb[8] = PPM_GETB( *(ip0-no) );
865 pr[7] = PPM_GETR( *(ip0) );
866 pg[7] = PPM_GETG( *(ip0) );
867 pb[7] = PPM_GETB( *(ip0) );
868 r = (*atfunc)(pr);
869 g = (*atfunc)(pg);
870 b = (*atfunc)(pb);
871 PPM_ASSIGN( *op, r, g, b );
872 }
873 pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
874 if (irow1 == irows[2]) {
875 irow1 = irows[0];
876 irow2 = irows[1];
877 irow0 = irows[2];
878 } else if (irow1 == irows[1]) {
879 irow2 = irows[0];
880 irow0 = irows[1];
881 irow1 = irows[2];
882 }
883 else /* must be at irows[0] */
884 {
885 irow0 = irows[0];
886 irow1 = irows[1];
887 irow2 = irows[2];
888 }
889 }
890 } else {
891 /* Else must be PGM */
892 int p[9]; /* 3x3 neighbor pixel values */
893 int pv;
894 int promote;
895
896 /* we scale maxval to omaxval */
897 promote = ( PNM_FORMAT_TYPE(format) != PNM_FORMAT_TYPE(oformat) );
898
899 for ( row = 0; row < rows; row++ ) {
900 int po,no; /* offsets for left and right columns in 3x3 */
901 xel *ip0, *ip1, *ip2, *op;
902
903 if (row == 0) {
904 irow0 = irow1;
905 pnm_readpnmrow( ifP, irow1, cols, maxval, format );
906 if ( promote )
907 pnm_promoteformatrow( irow1, cols, maxval,
908 format, maxval, oformat );
909 }
910 if (row == (rows-1))
911 irow2 = irow1;
912 else {
913 pnm_readpnmrow( ifP, irow2, cols, maxval, format );
914 if ( promote )
915 pnm_promoteformatrow( irow2, cols, maxval,
916 format, maxval, oformat );
917 }
918
919 for (col = cols-1,po= col>0?1:0,no=0,
920 ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
921 col >= 0;
922 col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0) {
923 /* grab 3x3 pixel values */
924 p[0] = PNM_GET1( *ip1 );
925 p[1] = PNM_GET1( *(ip1-no) );
926 p[5] = PNM_GET1( *(ip1+po) );
927 p[3] = PNM_GET1( *(ip2) );
928 p[2] = PNM_GET1( *(ip2-no) );
929 p[4] = PNM_GET1( *(ip2+po) );
930 p[6] = PNM_GET1( *(ip0+po) );
931 p[8] = PNM_GET1( *(ip0-no) );
932 p[7] = PNM_GET1( *(ip0) );
933 pv = (*atfunc)(p);
934 PNM_ASSIGN1( *op, pv );
935 }
936 pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
937 if (irow1 == irows[2]) {
938 irow1 = irows[0];
939 irow2 = irows[1];
940 irow0 = irows[2];
941 } else if (irow1 == irows[1]) {
942 irow2 = irows[0];
943 irow0 = irows[1];
944 irow1 = irows[2];
945 } else {
946 /* must be at irows[0] */
947 irow0 = irows[0];
948 irow1 = irows[1];
949 irow2 = irows[2];
950 }
951 }
952 }
953 }
954
955
956
957 static void
verifySame(unsigned int const imageSeq,int const imageCols,int const imageRows,xelval const imageMaxval,int const imageFormat,int const cols,int const rows,xelval const maxval,int const format)958 verifySame(unsigned int const imageSeq,
959 int const imageCols, int const imageRows,
960 xelval const imageMaxval, int const imageFormat,
961 int const cols, int const rows,
962 xelval const maxval, int const format) {
963 /*----------------------------------------------------------------------------
964 Issue error message and exit the program if the imageXXX arguments don't
965 match the XXX arguments.
966 -----------------------------------------------------------------------------*/
967 if (imageCols != cols)
968 pm_error("Width of Image %u (%d) is not the same as Image 0 (%d)",
969 imageSeq, imageCols, cols);
970 if (imageRows != rows)
971 pm_error("Height of Image %u (%d) is not the same as Image 0 (%d)",
972 imageSeq, imageRows, rows);
973 if (imageMaxval != maxval)
974 pm_error("Maxval of Image %u (%u) is not the same as Image 0 (%u)",
975 imageSeq, imageMaxval, maxval);
976 if (imageFormat != format)
977 pm_error("Format of Image %u is not the same as Image 0",
978 imageSeq);
979 }
980
981
982
983 int (*atfuncs[6]) (int *) = {atfilt0,atfilt1,atfilt2,atfilt3,atfilt4,atfilt5};
984
985
986
987
988 int
main(int argc,char * argv[])989 main(int argc, char *argv[]) {
990
991 FILE * ifP;
992 struct cmdlineInfo cmdline;
993 int eof; /* We've hit the end of the input stream */
994 unsigned int imageSeq; /* Sequence number of image, starting from 0 */
995
996 pnm_init(&argc, argv);
997
998 parseCommandLine(argc, argv, &cmdline);
999
1000 ifP = pm_openr(cmdline.inputFileName);
1001
1002 pnm_readpnminit(ifP, &cols, &rows, &maxval, &format);
1003
1004 if (maxval > MXIVAL)
1005 pm_error("The maxval of the input image (%d) is too large.\n"
1006 "This program's limit is %d.",
1007 maxval, MXIVAL);
1008
1009 oformat = PNM_FORMAT_TYPE(format);
1010 /* force output to max precision without forcing new 2-byte format */
1011 omaxval = MIN(maxval, PPM_MAXMAXVAL);
1012
1013 atfunc = atfuncs[atfilt_setup(cmdline.alpha, cmdline.radius,
1014 (double)omaxval/(double)maxval)];
1015
1016 if (oformat < PGM_TYPE) {
1017 oformat = RPGM_FORMAT;
1018 pm_message("promoting file to PGM");
1019 }
1020
1021 orow = pnm_allocrow(cols);
1022 irows[0] = pnm_allocrow(cols);
1023 irows[1] = pnm_allocrow(cols);
1024 irows[2] = pnm_allocrow(cols);
1025 irow0 = irows[0];
1026 irow1 = irows[1];
1027 irow2 = irows[2];
1028
1029 eof = FALSE; /* We're already in the middle of the first image */
1030 imageSeq = 0;
1031 while (!eof) {
1032 do_one_frame(ifP);
1033 pm_nextimage(ifP, &eof);
1034 if (!eof) {
1035 /* Read and validate header of next image */
1036 int imageCols, imageRows;
1037 xelval imageMaxval;
1038 int imageFormat;
1039
1040 ++imageSeq;
1041 pnm_readpnminit(ifP, &imageCols, &imageRows,
1042 &imageMaxval, &imageFormat);
1043 verifySame(imageSeq,
1044 imageCols, imageRows, imageMaxval, imageFormat,
1045 cols, rows, maxval, format);
1046 }
1047 }
1048
1049 pnm_freerow(irow0);
1050 pnm_freerow(irow1);
1051 pnm_freerow(irow2);
1052 pnm_freerow(orow);
1053 pm_close(ifP);
1054
1055 return 0;
1056 }
1057
1058
1059
1060