|
Packit |
78deda |
/* pnmnlfilt.c - 4 in 1 (2 non-linear) filter
|
|
Packit |
78deda |
** - smooth an anyimage
|
|
Packit |
78deda |
** - do alpha trimmed mean filtering on an anyimage
|
|
Packit |
78deda |
** - do optimal estimation smoothing on an anyimage
|
|
Packit |
78deda |
** - do edge enhancement on an anyimage
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** Version 1.0
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** The implementation of an alpha-trimmed mean filter
|
|
Packit |
78deda |
** is based on the description in IEEE CG&A May 1990
|
|
Packit |
78deda |
** Page 23 by Mark E. Lee and Richard A. Redner.
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** The paper recommends using a hexagon sampling region around each
|
|
Packit |
78deda |
** pixel being processed, allowing an effective sub pixel radius to be
|
|
Packit |
78deda |
** specified. The hexagon values are sythesized by area sampling the
|
|
Packit |
78deda |
** rectangular pixels with a hexagon grid. The seven hexagon values
|
|
Packit |
78deda |
** obtained from the 3x3 pixel grid are used to compute the alpha
|
|
Packit |
78deda |
** trimmed mean. Note that an alpha value of 0.0 gives a conventional
|
|
Packit |
78deda |
** mean filter (where the radius controls the contribution of
|
|
Packit |
78deda |
** surrounding pixels), while a value of 0.5 gives a median filter.
|
|
Packit |
78deda |
** Although there are only seven values to trim from before finding
|
|
Packit |
78deda |
** the mean, the algorithm has been extended from that described in
|
|
Packit |
78deda |
** CG&A by using interpolation, to allow a continuous selection of
|
|
Packit |
78deda |
** alpha value between and including 0.0 to 0.5 The useful values
|
|
Packit |
78deda |
** for radius are between 0.3333333 (where the filter will have no
|
|
Packit |
78deda |
** effect because only one pixel is sampled), to 1.0, where all
|
|
Packit |
78deda |
** pixels in the 3x3 grid are sampled.
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** The optimal estimation filter is taken from an article "Converting Dithered
|
|
Packit |
78deda |
** Images Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November
|
|
Packit |
78deda |
** 1992, and this article references "Digital Image Enhancement andNoise Filtering by
|
|
Packit |
78deda |
** Use of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and
|
|
Packit |
78deda |
** Machine Intelligence, March 1980.
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** Also borrow the technique used in pgmenhance(1) to allow edge
|
|
Packit |
78deda |
** enhancement if the alpha value is negative.
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** Author:
|
|
Packit |
78deda |
** Graeme W. Gill, 30th Jan 1993
|
|
Packit |
78deda |
** graeme@labtam.oz.au
|
|
Packit |
78deda |
**
|
|
Packit |
78deda |
** Permission to use, copy, modify, and distribute this software and its
|
|
Packit |
78deda |
** documentation for any purpose and without fee is hereby granted, provided
|
|
Packit |
78deda |
** that the above copyright notice appear in all copies and that both that
|
|
Packit |
78deda |
** copyright notice and this permission notice appear in supporting
|
|
Packit |
78deda |
** documentation. This software is provided "as is" without express or
|
|
Packit |
78deda |
** implied warranty.
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#include <math.h>
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#include "pm_c_util.h"
|
|
Packit |
78deda |
#include "pnm.h"
|
|
Packit |
78deda |
|
|
Packit |
78deda |
struct cmdlineInfo {
|
|
Packit |
78deda |
const char * inputFileName;
|
|
Packit |
78deda |
double alpha;
|
|
Packit |
78deda |
double radius;
|
|
Packit |
78deda |
};
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
parseCommandLine(int argc,
|
|
Packit |
78deda |
char ** argv,
|
|
Packit |
78deda |
struct cmdlineInfo * const cmdlineP) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (argc-1 < 2)
|
|
Packit |
78deda |
pm_error("You must specify at least two arguments: alpha and radius. "
|
|
Packit |
78deda |
"You specified %u", argc-1);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (sscanf(argv[1], "%lf", &cmdlineP->alpha) != 1)
|
|
Packit |
78deda |
pm_error("Invalid alpha (1st) argument '%s'. "
|
|
Packit |
78deda |
"Must be a decimal number",
|
|
Packit |
78deda |
argv[1]);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (sscanf( argv[2], "%lf", &cmdlineP->radius ) != 1)
|
|
Packit |
78deda |
pm_error("Invalid radius (2nd) argument '%s'. "
|
|
Packit |
78deda |
"Must be a decimal number",
|
|
Packit |
78deda |
argv[2]);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if ((cmdlineP->alpha > -0.1 && cmdlineP->alpha < 0.0) ||
|
|
Packit |
78deda |
(cmdlineP->alpha > 0.5 && cmdlineP->alpha < 1.0))
|
|
Packit |
78deda |
pm_error( "Alpha must be in range 0.0 <= alpha <= 0.5 "
|
|
Packit |
78deda |
"for alpha trimmed mean. "
|
|
Packit |
78deda |
"You specified %f", cmdlineP->alpha);
|
|
Packit |
78deda |
if (cmdlineP->alpha > 2.0)
|
|
Packit |
78deda |
pm_error("Alpha must be in range 1.0 <= cmdlineP->alpha <= 2.0 "
|
|
Packit |
78deda |
"for optimal estimation. You specified %f",
|
|
Packit |
78deda |
cmdlineP->alpha);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (cmdlineP->alpha < -0.9 ||
|
|
Packit |
78deda |
(cmdlineP->alpha > -0.1 && cmdlineP->alpha < 0.0))
|
|
Packit |
78deda |
pm_error( "Alpha must be in range -0.9 <= alpha <= -0.1 "
|
|
Packit |
78deda |
"for edge enhancement. You specified %f",
|
|
Packit |
78deda |
cmdlineP->alpha);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (cmdlineP->radius < 1.0/3 || cmdlineP->radius > 1.0)
|
|
Packit |
78deda |
pm_error("Radius must be in range 1/3 <= radius <= 1. "
|
|
Packit |
78deda |
"You specified %f", cmdlineP->radius);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (argc-1 < 3)
|
|
Packit |
78deda |
cmdlineP->inputFileName = "-";
|
|
Packit |
78deda |
else
|
|
Packit |
78deda |
cmdlineP->inputFileName = argv[3];
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (argc-1 > 3)
|
|
Packit |
78deda |
pm_error("Too many arguments: %u. The most allowed are 3: alpha, "
|
|
Packit |
78deda |
"radius, and file name", argc-1);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* MXIVAL is the maximum input sample value we can handle.
|
|
Packit |
78deda |
It is limited by our willingness to allocate storage in various arrays
|
|
Packit |
78deda |
that are indexed by sample values.
|
|
Packit |
78deda |
|
|
Packit |
78deda |
We use PPM_MAXMAXVAL because that used to be the maximum possible
|
|
Packit |
78deda |
sample value in the format, and most images still limit themselves to
|
|
Packit |
78deda |
this value.
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define MXIVAL PPM_MAXMAXVAL
|
|
Packit |
78deda |
|
|
Packit |
78deda |
xelval omaxval;
|
|
Packit |
78deda |
/* global so that pixel processing code can get at it quickly */
|
|
Packit |
78deda |
int noisevariance;
|
|
Packit |
78deda |
/* global so that pixel processing code can get at it quickly */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/*
|
|
Packit |
78deda |
* Declared static here rather than passing a jillion options in the call to
|
|
Packit |
78deda |
* do_one_frame(). Also it makes a huge amount of sense to only malloc the
|
|
Packit |
78deda |
* row buffers once instead of for each frame (with the corresponding free'ing
|
|
Packit |
78deda |
* of course).
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
static xel *irows[3];
|
|
Packit |
78deda |
static xel *irow0, *irow1, *irow2, *orow;
|
|
Packit |
78deda |
static int rows, cols, format, oformat, row, col;
|
|
Packit |
78deda |
static int (*atfunc)(int *);
|
|
Packit |
78deda |
static xelval maxval;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define NOIVAL (MXIVAL + 1) /* number of possible input values */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define SCALEB 8 /* scale bits */
|
|
Packit |
78deda |
#define SCALE (1 << SCALEB) /* scale factor */
|
|
Packit |
78deda |
#define MXSVAL (MXIVAL * SCALE) /* maximum scaled values */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define CSCALEB 2 /* coarse scale bits */
|
|
Packit |
78deda |
#define CSCALE (1 << CSCALEB) /* coarse scale factor */
|
|
Packit |
78deda |
#define MXCSVAL (MXIVAL * CSCALE) /* maximum coarse scaled values */
|
|
Packit |
78deda |
#define NOCSVAL (MXCSVAL + 1) /* number of coarse scaled values */
|
|
Packit |
78deda |
#define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB)) /* scaled to coarse scaled */
|
|
Packit |
78deda |
#define CSCTOSC(x) ((x) << (SCALEB - CSCALEB)) /* course scaled to scaled */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#ifndef MAXINT
|
|
Packit |
78deda |
# define MAXINT 0x7fffffff /* assume this is a 32 bit machine */
|
|
Packit |
78deda |
#endif
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* round and scale floating point to scaled integer */
|
|
Packit |
78deda |
#define ROUNDSCALE(x) ((int)(((x) * (double)SCALE) + 0.5))
|
|
Packit |
78deda |
/* round and un-scale scaled integer value */
|
|
Packit |
78deda |
#define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB)
|
|
Packit |
78deda |
/* rounded un-scale */
|
|
Packit |
78deda |
#define UNSCALE(x) ((x) >> SCALEB)
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static double
|
|
Packit |
78deda |
sqr(double const arg) {
|
|
Packit |
78deda |
return arg * arg;
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* We restrict radius to the values: 0.333333 <= radius <= 1.0 */
|
|
Packit |
78deda |
/* so that no fewer and no more than a 3x3 grid of pixels around */
|
|
Packit |
78deda |
/* the pixel in question needs to be read. Given this, we only */
|
|
Packit |
78deda |
/* need 3 or 4 weightings per hexagon, as follows: */
|
|
Packit |
78deda |
/* _ _ */
|
|
Packit |
78deda |
/* Vertical hex: |_|_| 1 2 */
|
|
Packit |
78deda |
/* |X|_| 0 3 */
|
|
Packit |
78deda |
/* _ */
|
|
Packit |
78deda |
/* _ _|_| 1 */
|
|
Packit |
78deda |
/* Middle hex: |_| 1 Horizontal hex: |X|_| 0 2 */
|
|
Packit |
78deda |
/* |X| 0 |_| 3 */
|
|
Packit |
78deda |
/* |_| 2 */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* all filters */
|
|
Packit |
78deda |
int V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL]; /* vertical hex */
|
|
Packit |
78deda |
int M0[NOIVAL],M1[NOIVAL],M2[NOIVAL]; /* middle hex */
|
|
Packit |
78deda |
int H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL]; /* horizontal hex */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* alpha trimmed and edge enhancement only */
|
|
Packit |
78deda |
int ALFRAC[NOIVAL * 8]; /* fractional alpha divider table */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* optimal estimation only */
|
|
Packit |
78deda |
int AVEDIV[7 * NOCSVAL]; /* divide by 7 to give average value */
|
|
Packit |
78deda |
int SQUARE[2 * NOCSVAL]; /* scaled square lookup table */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* ************************************************** *
|
|
Packit |
78deda |
Hexagon intersecting square area functions
|
|
Packit |
78deda |
Compute the area of the intersection of a triangle
|
|
Packit |
78deda |
and a rectangle
|
|
Packit |
78deda |
************************************************** */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Triangle orientation is per geometric axes (not graphical axies) */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define NW 0 /* North west triangle /| */
|
|
Packit |
78deda |
#define NE 1 /* North east triangle |\ */
|
|
Packit |
78deda |
#define SW 2 /* South west triangle \| */
|
|
Packit |
78deda |
#define SE 3 /* South east triangle |/ */
|
|
Packit |
78deda |
#define STH 2
|
|
Packit |
78deda |
#define EST 1
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define SWAPI(a,b) (t = a, a = -b, b = -t)
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static double
|
|
Packit |
78deda |
triang_area(double rx0, double ry0, double rx1, double ry1,
|
|
Packit |
78deda |
double tx0, double ty0, double tx1, double ty1,
|
|
Packit |
78deda |
int tt) {
|
|
Packit |
78deda |
/* rx0,ry0,rx1,ry1: rectangle boundaries */
|
|
Packit |
78deda |
/* tx0,ty0,tx1,ty1: triangle boundaries */
|
|
Packit |
78deda |
/* tt: triangle type */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
double a,b,c,d;
|
|
Packit |
78deda |
double lx0,ly0,lx1,ly1;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Convert everything to a NW triangle */
|
|
Packit |
78deda |
if (tt & STH) {
|
|
Packit |
78deda |
double t;
|
|
Packit |
78deda |
SWAPI(ry0,ry1);
|
|
Packit |
78deda |
SWAPI(ty0,ty1);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
if (tt & EST) {
|
|
Packit |
78deda |
double t;
|
|
Packit |
78deda |
SWAPI(rx0,rx1);
|
|
Packit |
78deda |
SWAPI(tx0,tx1);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
/* Compute overlapping box */
|
|
Packit |
78deda |
if (tx0 > rx0)
|
|
Packit |
78deda |
rx0 = tx0;
|
|
Packit |
78deda |
if (ty0 > ry0)
|
|
Packit |
78deda |
ry0 = ty0;
|
|
Packit |
78deda |
if (tx1 < rx1)
|
|
Packit |
78deda |
rx1 = tx1;
|
|
Packit |
78deda |
if (ty1 < ry1)
|
|
Packit |
78deda |
ry1 = ty1;
|
|
Packit |
78deda |
if (rx1 <= rx0 || ry1 <= ry0)
|
|
Packit |
78deda |
return 0.0;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Need to compute diagonal line intersection with the box */
|
|
Packit |
78deda |
/* First compute co-efficients to formulas x = a + by and y = c + dx */
|
|
Packit |
78deda |
b = (tx1 - tx0)/(ty1 - ty0);
|
|
Packit |
78deda |
a = tx0 - b * ty0;
|
|
Packit |
78deda |
d = (ty1 - ty0)/(tx1 - tx0);
|
|
Packit |
78deda |
c = ty0 - d * tx0;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* compute top or right intersection */
|
|
Packit |
78deda |
tt = 0;
|
|
Packit |
78deda |
ly1 = ry1;
|
|
Packit |
78deda |
lx1 = a + b * ly1;
|
|
Packit |
78deda |
if (lx1 <= rx0)
|
|
Packit |
78deda |
return (rx1 - rx0) * (ry1 - ry0);
|
|
Packit |
78deda |
else if (lx1 > rx1) {
|
|
Packit |
78deda |
/* could be right hand side */
|
|
Packit |
78deda |
lx1 = rx1;
|
|
Packit |
78deda |
ly1 = c + d * lx1;
|
|
Packit |
78deda |
if (ly1 <= ry0)
|
|
Packit |
78deda |
return (rx1 - rx0) * (ry1 - ry0);
|
|
Packit |
78deda |
tt = 1; /* right hand side intersection */
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
/* compute left or bottom intersection */
|
|
Packit |
78deda |
lx0 = rx0;
|
|
Packit |
78deda |
ly0 = c + d * lx0;
|
|
Packit |
78deda |
if (ly0 >= ry1)
|
|
Packit |
78deda |
return (rx1 - rx0) * (ry1 - ry0);
|
|
Packit |
78deda |
else if (ly0 < ry0) {
|
|
Packit |
78deda |
/* could be right hand side */
|
|
Packit |
78deda |
ly0 = ry0;
|
|
Packit |
78deda |
lx0 = a + b * ly0;
|
|
Packit |
78deda |
if (lx0 >= rx1)
|
|
Packit |
78deda |
return (rx1 - rx0) * (ry1 - ry0);
|
|
Packit |
78deda |
tt |= 2; /* bottom intersection */
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (tt == 0) {
|
|
Packit |
78deda |
/* top and left intersection */
|
|
Packit |
78deda |
/* rectangle minus triangle */
|
|
Packit |
78deda |
return ((rx1 - rx0) * (ry1 - ry0))
|
|
Packit |
78deda |
- (0.5 * (lx1 - rx0) * (ry1 - ly0));
|
|
Packit |
78deda |
} else if (tt == 1) {
|
|
Packit |
78deda |
/* right and left intersection */
|
|
Packit |
78deda |
return ((rx1 - rx0) * (ly0 - ry0))
|
|
Packit |
78deda |
+ (0.5 * (rx1 - rx0) * (ly1 - ly0));
|
|
Packit |
78deda |
} else if (tt == 2) {
|
|
Packit |
78deda |
/* top and bottom intersection */
|
|
Packit |
78deda |
return ((rx1 - lx1) * (ry1 - ry0))
|
|
Packit |
78deda |
+ (0.5 * (lx1 - lx0) * (ry1 - ry0));
|
|
Packit |
78deda |
} else {
|
|
Packit |
78deda |
/* tt == 3 */
|
|
Packit |
78deda |
/* right and bottom intersection */
|
|
Packit |
78deda |
/* triangle */
|
|
Packit |
78deda |
return (0.5 * (rx1 - lx0) * (ly1 - ry0));
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static double
|
|
Packit |
78deda |
rectang_area(double rx0, double ry0, double rx1, double ry1,
|
|
Packit |
78deda |
double tx0, double ty0, double tx1, double ty1) {
|
|
Packit |
78deda |
/* Compute rectangle area */
|
|
Packit |
78deda |
/* rx0,ry0,rx1,ry1: rectangle boundaries */
|
|
Packit |
78deda |
/* tx0,ty0,tx1,ty1: rectangle boundaries */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Compute overlapping box */
|
|
Packit |
78deda |
if (tx0 > rx0)
|
|
Packit |
78deda |
rx0 = tx0;
|
|
Packit |
78deda |
if (ty0 > ry0)
|
|
Packit |
78deda |
ry0 = ty0;
|
|
Packit |
78deda |
if (tx1 < rx1)
|
|
Packit |
78deda |
rx1 = tx1;
|
|
Packit |
78deda |
if (ty1 < ry1)
|
|
Packit |
78deda |
ry1 = ty1;
|
|
Packit |
78deda |
if (rx1 <= rx0 || ry1 <= ry0)
|
|
Packit |
78deda |
return 0.0;
|
|
Packit |
78deda |
return (rx1 - rx0) * (ry1 - ry0);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static double
|
|
Packit |
78deda |
hex_area(double sx, double sy, double hx, double hy, double d) {
|
|
Packit |
78deda |
/* compute the area of overlap of a hexagon diameter d, */
|
|
Packit |
78deda |
/* centered at hx,hy, with a unit square of center sx,sy. */
|
|
Packit |
78deda |
/* sx,sy: square center */
|
|
Packit |
78deda |
/* hx,hy,d: hexagon center and diameter */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
double hx0,hx1,hx2,hy0,hy1,hy2,hy3;
|
|
Packit |
78deda |
double sx0,sx1,sy0,sy1;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* compute square co-ordinates */
|
|
Packit |
78deda |
sx0 = sx - 0.5;
|
|
Packit |
78deda |
sy0 = sy - 0.5;
|
|
Packit |
78deda |
sx1 = sx + 0.5;
|
|
Packit |
78deda |
sy1 = sy + 0.5;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* compute hexagon co-ordinates */
|
|
Packit |
78deda |
hx0 = hx - d/2.0;
|
|
Packit |
78deda |
hx1 = hx;
|
|
Packit |
78deda |
hx2 = hx + d/2.0;
|
|
Packit |
78deda |
hy0 = hy - 0.5773502692 * d; /* d / sqrt(3) */
|
|
Packit |
78deda |
hy1 = hy - 0.2886751346 * d; /* d / sqrt(12) */
|
|
Packit |
78deda |
hy2 = hy + 0.2886751346 * d; /* d / sqrt(12) */
|
|
Packit |
78deda |
hy3 = hy + 0.5773502692 * d; /* d / sqrt(3) */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
return
|
|
Packit |
78deda |
triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) +
|
|
Packit |
78deda |
triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) +
|
|
Packit |
78deda |
rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) +
|
|
Packit |
78deda |
triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) +
|
|
Packit |
78deda |
triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
setupAvediv(void) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
unsigned int i;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (i=0; i < (7 * NOCSVAL); ++i) {
|
|
Packit |
78deda |
/* divide scaled value by 7 lookup */
|
|
Packit |
78deda |
AVEDIV[i] = CSCTOSC(i)/7; /* scaled divide by 7 */
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
setupSquare(void) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
unsigned int i;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (i=0; i < (2 * NOCSVAL); ++i) {
|
|
Packit |
78deda |
/* compute square and rescale by (val >> (2 * SCALEB + 2)) table */
|
|
Packit |
78deda |
int const val = CSCTOSC(i - NOCSVAL);
|
|
Packit |
78deda |
/* NOCSVAL offset to cope with -ve input values */
|
|
Packit |
78deda |
SQUARE[i] = (val * val) >> (2 * SCALEB + 2);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
setup1(double const alpha,
|
|
Packit |
78deda |
double const radius,
|
|
Packit |
78deda |
double const maxscale,
|
|
Packit |
78deda |
int * const alpharangeP,
|
|
Packit |
78deda |
double * const meanscaleP,
|
|
Packit |
78deda |
double * const mmeanscaleP,
|
|
Packit |
78deda |
double * const alphafractionP,
|
|
Packit |
78deda |
int * const noisevarianceP) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
setupAvediv();
|
|
Packit |
78deda |
setupSquare();
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (alpha >= 0.0 && alpha <= 0.5) {
|
|
Packit |
78deda |
/* alpha trimmed mean */
|
|
Packit |
78deda |
double const noinmean = ((0.5 - alpha) * 12.0) + 1.0;
|
|
Packit |
78deda |
/* number of elements (out of a possible 7) used in the mean */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
*mmeanscaleP = *meanscaleP = maxscale/noinmean;
|
|
Packit |
78deda |
if (alpha == 0.0) {
|
|
Packit |
78deda |
/* mean filter */
|
|
Packit |
78deda |
*alpharangeP = 0;
|
|
Packit |
78deda |
*alphafractionP = 0.0; /* not used */
|
|
Packit |
78deda |
} else if (alpha < (1.0/6.0)) {
|
|
Packit |
78deda |
/* mean of 5 to 7 middle values */
|
|
Packit |
78deda |
*alpharangeP = 1;
|
|
Packit |
78deda |
*alphafractionP = (7.0 - noinmean)/2.0;
|
|
Packit |
78deda |
} else if (alpha < (1.0/3.0)) {
|
|
Packit |
78deda |
/* mean of 3 to 5 middle values */
|
|
Packit |
78deda |
*alpharangeP = 2;
|
|
Packit |
78deda |
*alphafractionP = (5.0 - noinmean)/2.0;
|
|
Packit |
78deda |
} else {
|
|
Packit |
78deda |
/* mean of 1 to 3 middle values */
|
|
Packit |
78deda |
/* alpha == 0.5 == median filter */
|
|
Packit |
78deda |
*alpharangeP = 3;
|
|
Packit |
78deda |
*alphafractionP = (3.0 - noinmean)/2.0;
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
} else if (alpha >= 1.0 && alpha <= 2.0) {
|
|
Packit |
78deda |
/* optimal estimation - alpha controls noise variance threshold. */
|
|
Packit |
78deda |
double const alphaNormalized = alpha - 1.0;
|
|
Packit |
78deda |
/* normalize it to 0.0 -> 1.0 */
|
|
Packit |
78deda |
double const noinmean = 7.0;
|
|
Packit |
78deda |
*alpharangeP = 5; /* edge enhancement function */
|
|
Packit |
78deda |
*mmeanscaleP = *meanscaleP = maxscale; /* compute scaled hex values */
|
|
Packit |
78deda |
*alphafractionP = 1.0/noinmean;
|
|
Packit |
78deda |
/* Set up 1:1 division lookup - not used */
|
|
Packit |
78deda |
*noisevarianceP = sqr(alphaNormalized * omaxval) / 8.0;
|
|
Packit |
78deda |
/* estimate of noise variance */
|
|
Packit |
78deda |
} else if (alpha >= -0.9 && alpha <= -0.1) {
|
|
Packit |
78deda |
/* edge enhancement function */
|
|
Packit |
78deda |
double const posAlpha = -alpha;
|
|
Packit |
78deda |
/* positive alpha value */
|
|
Packit |
78deda |
*alpharangeP = 4; /* edge enhancement function */
|
|
Packit |
78deda |
*meanscaleP = maxscale * (-posAlpha/((1.0 - posAlpha) * 7.0));
|
|
Packit |
78deda |
/* mean of 7 and scaled by -posAlpha/(1-posAlpha) */
|
|
Packit |
78deda |
*mmeanscaleP = maxscale * (1.0/(1.0 - posAlpha) + *meanscaleP);
|
|
Packit |
78deda |
/* middle pixel has 1/(1-posAlpha) as well */
|
|
Packit |
78deda |
*alphafractionP = 0.0; /* not used */
|
|
Packit |
78deda |
} else {
|
|
Packit |
78deda |
/* An entry condition on 'alpha' makes this impossible */
|
|
Packit |
78deda |
pm_error("INTERNAL ERROR: impossible alpha value: %f", alpha);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
setupAlfrac(double const alphafraction) {
|
|
Packit |
78deda |
/* set up alpha fraction lookup table used on big/small */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
unsigned int i;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (i=0; i < (NOIVAL * 8); ++i) {
|
|
Packit |
78deda |
ALFRAC[i] = ROUNDSCALE(i * alphafraction);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
setupPixelWeightingTables(double const radius,
|
|
Packit |
78deda |
double const meanscale,
|
|
Packit |
78deda |
double const mmeanscale) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Setup pixel weighting tables - note we pre-compute mean
|
|
Packit |
78deda |
division here too.
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
double const hexhoff = radius/2;
|
|
Packit |
78deda |
/* horizontal offset of vertical hex centers */
|
|
Packit |
78deda |
double const hexvoff = 3.0 * radius/sqrt(12.0);
|
|
Packit |
78deda |
/* vertical offset of vertical hex centers */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
double const tabscale = meanscale / (radius * hexvoff);
|
|
Packit |
78deda |
double const mtabscale = mmeanscale / (radius * hexvoff);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* scale tables to normalize by hexagon area, and number of
|
|
Packit |
78deda |
hexes used in mean
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
double const v0 =
|
|
Packit |
78deda |
hex_area(0.0, 0.0, hexhoff, hexvoff, radius) * tabscale;
|
|
Packit |
78deda |
double const v1 =
|
|
Packit |
78deda |
hex_area(0.0, 1.0, hexhoff, hexvoff, radius) * tabscale;
|
|
Packit |
78deda |
double const v2 =
|
|
Packit |
78deda |
hex_area(1.0, 1.0, hexhoff, hexvoff, radius) * tabscale;
|
|
Packit |
78deda |
double const v3 =
|
|
Packit |
78deda |
hex_area(1.0, 0.0, hexhoff, hexvoff, radius) * tabscale;
|
|
Packit |
78deda |
double const m0 =
|
|
Packit |
78deda |
hex_area(0.0, 0.0, 0.0, 0.0, radius) * mtabscale;
|
|
Packit |
78deda |
double const m1 =
|
|
Packit |
78deda |
hex_area(0.0, 1.0, 0.0, 0.0, radius) * mtabscale;
|
|
Packit |
78deda |
double const m2 =
|
|
Packit |
78deda |
hex_area(0.0, -1.0, 0.0, 0.0, radius) * mtabscale;
|
|
Packit |
78deda |
double const h0 =
|
|
Packit |
78deda |
hex_area(0.0, 0.0, radius, 0.0, radius) * tabscale;
|
|
Packit |
78deda |
double const h1 =
|
|
Packit |
78deda |
hex_area(1.0, 1.0, radius, 0.0, radius) * tabscale;
|
|
Packit |
78deda |
double const h2 =
|
|
Packit |
78deda |
hex_area(1.0, 0.0, radius, 0.0, radius) * tabscale;
|
|
Packit |
78deda |
double const h3 =
|
|
Packit |
78deda |
hex_area(1.0, -1.0, radius, 0.0, radius) * tabscale;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
unsigned int i;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (i=0; i <= MXIVAL; ++i) {
|
|
Packit |
78deda |
V0[i] = ROUNDSCALE(i * v0);
|
|
Packit |
78deda |
V1[i] = ROUNDSCALE(i * v1);
|
|
Packit |
78deda |
V2[i] = ROUNDSCALE(i * v2);
|
|
Packit |
78deda |
V3[i] = ROUNDSCALE(i * v3);
|
|
Packit |
78deda |
M0[i] = ROUNDSCALE(i * m0);
|
|
Packit |
78deda |
M1[i] = ROUNDSCALE(i * m1);
|
|
Packit |
78deda |
M2[i] = ROUNDSCALE(i * m2);
|
|
Packit |
78deda |
H0[i] = ROUNDSCALE(i * h0);
|
|
Packit |
78deda |
H1[i] = ROUNDSCALE(i * h1);
|
|
Packit |
78deda |
H2[i] = ROUNDSCALE(i * h2);
|
|
Packit |
78deda |
H3[i] = ROUNDSCALE(i * h3);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* Table initialization function - return alpha range */
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt_setup(double const alpha,
|
|
Packit |
78deda |
double const radius,
|
|
Packit |
78deda |
double const maxscale) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int alpharange; /* alpha range value 0 - 5 */
|
|
Packit |
78deda |
double meanscale; /* scale for finding mean */
|
|
Packit |
78deda |
double mmeanscale; /* scale for finding mean - midle hex */
|
|
Packit |
78deda |
double alphafraction;
|
|
Packit |
78deda |
/* fraction of next largest/smallest to subtract from sum */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
setup1(alpha, radius, maxscale,
|
|
Packit |
78deda |
&alpharange, &meanscale, &mmeanscale, &alphafraction,
|
|
Packit |
78deda |
&noisevariance);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
setupAlfrac(alphafraction);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
setupPixelWeightingTables(radius, meanscale, mmeanscale);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
return alpharange;
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt0(int * p) {
|
|
Packit |
78deda |
/* Core pixel processing function - hand it 3x3 pixels and return result. */
|
|
Packit |
78deda |
/* Mean filter */
|
|
Packit |
78deda |
/* 'p' is 9 pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
int retv;
|
|
Packit |
78deda |
/* map to scaled hexagon values */
|
|
Packit |
78deda |
retv = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
return UNSCALE(retv);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define CHECK(xx) {\
|
|
Packit |
78deda |
h0 += xx; \
|
|
Packit |
78deda |
if (xx > big) \
|
|
Packit |
78deda |
big = xx; \
|
|
Packit |
78deda |
else if (xx < small) \
|
|
Packit |
78deda |
small = xx; }
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt1(int * p) {
|
|
Packit |
78deda |
/* Mean of 5 - 7 middle values */
|
|
Packit |
78deda |
/* 'p' is 9 pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
|
|
Packit |
78deda |
/* 1 0 4 */
|
|
Packit |
78deda |
/* 6 5 */
|
|
Packit |
78deda |
int big,small;
|
|
Packit |
78deda |
/* map to scaled hexagon values */
|
|
Packit |
78deda |
h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
/* sum values and also discover the largest and smallest */
|
|
Packit |
78deda |
big = small = h0;
|
|
Packit |
78deda |
CHECK(h1);
|
|
Packit |
78deda |
CHECK(h2);
|
|
Packit |
78deda |
CHECK(h3);
|
|
Packit |
78deda |
CHECK(h4);
|
|
Packit |
78deda |
CHECK(h5);
|
|
Packit |
78deda |
CHECK(h6);
|
|
Packit |
78deda |
/* Compute mean of middle 5-7 values */
|
|
Packit |
78deda |
return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
#undef CHECK
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define CHECK(xx) {\
|
|
Packit |
78deda |
h0 += xx; \
|
|
Packit |
78deda |
if (xx > big1) {\
|
|
Packit |
78deda |
if (xx > big0) {\
|
|
Packit |
78deda |
big1 = big0; \
|
|
Packit |
78deda |
big0 = xx; \
|
|
Packit |
78deda |
} else \
|
|
Packit |
78deda |
big1 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
if (xx < small1) {\
|
|
Packit |
78deda |
if (xx < small0) {\
|
|
Packit |
78deda |
small1 = small0; \
|
|
Packit |
78deda |
small0 = xx; \
|
|
Packit |
78deda |
} else \
|
|
Packit |
78deda |
small1 = xx; \
|
|
Packit |
78deda |
}\
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt2(int *p) {
|
|
Packit |
78deda |
/* Mean of 3 - 5 middle values */
|
|
Packit |
78deda |
/* 'p' is 9 pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
|
|
Packit |
78deda |
/* 1 0 4 */
|
|
Packit |
78deda |
/* 6 5 */
|
|
Packit |
78deda |
int big0,big1,small0,small1;
|
|
Packit |
78deda |
/* map to scaled hexagon values */
|
|
Packit |
78deda |
h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
/* sum values and also discover the 2 largest and 2 smallest */
|
|
Packit |
78deda |
big0 = small0 = h0;
|
|
Packit |
78deda |
small1 = MAXINT;
|
|
Packit |
78deda |
big1 = 0;
|
|
Packit |
78deda |
CHECK(h1);
|
|
Packit |
78deda |
CHECK(h2);
|
|
Packit |
78deda |
CHECK(h3);
|
|
Packit |
78deda |
CHECK(h4);
|
|
Packit |
78deda |
CHECK(h5);
|
|
Packit |
78deda |
CHECK(h6);
|
|
Packit |
78deda |
/* Compute mean of middle 3-5 values */
|
|
Packit |
78deda |
return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#undef CHECK
|
|
Packit |
78deda |
|
|
Packit |
78deda |
#define CHECK(xx) {\
|
|
Packit |
78deda |
h0 += xx; \
|
|
Packit |
78deda |
if (xx > big2) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
if (xx > big1) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
if (xx > big0) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
big2 = big1; \
|
|
Packit |
78deda |
big1 = big0; \
|
|
Packit |
78deda |
big0 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
else \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
big2 = big1; \
|
|
Packit |
78deda |
big1 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
else \
|
|
Packit |
78deda |
big2 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
if (xx < small2) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
if (xx < small1) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
if (xx < small0) \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
small2 = small1; \
|
|
Packit |
78deda |
small1 = small0; \
|
|
Packit |
78deda |
small0 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
else \
|
|
Packit |
78deda |
{ \
|
|
Packit |
78deda |
small2 = small1; \
|
|
Packit |
78deda |
small1 = xx; \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
} \
|
|
Packit |
78deda |
else \
|
|
Packit |
78deda |
small2 = xx; \
|
|
Packit |
78deda |
}}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt3(int *p) {
|
|
Packit |
78deda |
/* Mean of 1 - 3 middle values. If only 1 value, then this is a median
|
|
Packit |
78deda |
filter.
|
|
Packit |
78deda |
*/
|
|
Packit |
78deda |
/* 'p' is pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
|
|
Packit |
78deda |
/* 1 0 4 */
|
|
Packit |
78deda |
/* 6 5 */
|
|
Packit |
78deda |
int big0,big1,big2,small0,small1,small2;
|
|
Packit |
78deda |
/* map to scaled hexagon values */
|
|
Packit |
78deda |
h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
/* sum values and also discover the 3 largest and 3 smallest */
|
|
Packit |
78deda |
big0 = small0 = h0;
|
|
Packit |
78deda |
small1 = small2 = MAXINT;
|
|
Packit |
78deda |
big1 = big2 = 0;
|
|
Packit |
78deda |
CHECK(h1);
|
|
Packit |
78deda |
CHECK(h2);
|
|
Packit |
78deda |
CHECK(h3);
|
|
Packit |
78deda |
CHECK(h4);
|
|
Packit |
78deda |
CHECK(h5);
|
|
Packit |
78deda |
CHECK(h6);
|
|
Packit |
78deda |
/* Compute mean of middle 1-3 values */
|
|
Packit |
78deda |
return UNSCALE(h0 -big0 -big1 -small0 -small1
|
|
Packit |
78deda |
-ALFRAC[(big2 + small2)>>SCALEB]);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
#undef CHECK
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt4(int *p) {
|
|
Packit |
78deda |
/* Edge enhancement */
|
|
Packit |
78deda |
/* notice we use the global omaxval */
|
|
Packit |
78deda |
/* 'p' is 9 pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int hav;
|
|
Packit |
78deda |
/* map to scaled hexagon values and compute enhance value */
|
|
Packit |
78deda |
hav = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
if (hav < 0)
|
|
Packit |
78deda |
hav = 0;
|
|
Packit |
78deda |
hav = UNSCALE(hav);
|
|
Packit |
78deda |
if (hav > omaxval)
|
|
Packit |
78deda |
hav = omaxval;
|
|
Packit |
78deda |
return hav;
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static int
|
|
Packit |
78deda |
atfilt5(int *p) {
|
|
Packit |
78deda |
/* Optimal estimation - do smoothing in inverse proportion */
|
|
Packit |
78deda |
/* to the local variance. */
|
|
Packit |
78deda |
/* notice we use the globals noisevariance and omaxval*/
|
|
Packit |
78deda |
/* 'p' is 9 pixel values from 3x3 neighbors */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int mean,variance,temp;
|
|
Packit |
78deda |
int h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
|
|
Packit |
78deda |
/* 1 0 4 */
|
|
Packit |
78deda |
/* 6 5 */
|
|
Packit |
78deda |
/* map to scaled hexagon values */
|
|
Packit |
78deda |
h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
|
|
Packit |
78deda |
h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
|
|
Packit |
78deda |
h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
|
|
Packit |
78deda |
h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
|
|
Packit |
78deda |
h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
|
|
Packit |
78deda |
h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
|
|
Packit |
78deda |
h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
|
|
Packit |
78deda |
mean = h0 + h1 + h2 + h3 + h4 + h5 + h6;
|
|
Packit |
78deda |
mean = AVEDIV[SCTOCSC(mean)]; /* compute scaled mean by dividing by 7 */
|
|
Packit |
78deda |
temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
/* compute scaled variance */
|
|
Packit |
78deda |
temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
/* and rescale to keep */
|
|
Packit |
78deda |
temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
/* within 32 bit limits */
|
|
Packit |
78deda |
temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
|
|
Packit |
78deda |
/* (temp = h0 - mean) */
|
|
Packit |
78deda |
if (variance != 0) /* avoid possible divide by 0 */
|
|
Packit |
78deda |
temp = mean + (variance * temp) / (variance + noisevariance);
|
|
Packit |
78deda |
/* optimal estimate */
|
|
Packit |
78deda |
else temp = h0;
|
|
Packit |
78deda |
if (temp < 0)
|
|
Packit |
78deda |
temp = 0;
|
|
Packit |
78deda |
temp = RUNSCALE(temp);
|
|
Packit |
78deda |
if (temp > omaxval)
|
|
Packit |
78deda |
temp = omaxval;
|
|
Packit |
78deda |
return temp;
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
do_one_frame(FILE * const ifP) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
pnm_writepnminit( stdout, cols, rows, omaxval, oformat, 0 );
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if ( PNM_FORMAT_TYPE(oformat) == PPM_TYPE ) {
|
|
Packit |
78deda |
int pr[9],pg[9],pb[9]; /* 3x3 neighbor pixel values */
|
|
Packit |
78deda |
int r,g,b;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for ( row = 0; row < rows; row++ ) {
|
|
Packit |
78deda |
int po,no; /* offsets for left and right colums in 3x3 */
|
|
Packit |
78deda |
xel *ip0, *ip1, *ip2, *op;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (row == 0) {
|
|
Packit |
78deda |
irow0 = irow1;
|
|
Packit |
78deda |
pnm_readpnmrow( ifP, irow1, cols, maxval, format );
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
if (row == (rows-1))
|
|
Packit |
78deda |
irow2 = irow1;
|
|
Packit |
78deda |
else
|
|
Packit |
78deda |
pnm_readpnmrow( ifP, irow2, cols, maxval, format );
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (col = cols-1,po= col>0?1:0,no=0,
|
|
Packit |
78deda |
ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
|
|
Packit |
78deda |
col >= 0;
|
|
Packit |
78deda |
col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0) {
|
|
Packit |
78deda |
/* grab 3x3 pixel values */
|
|
Packit |
78deda |
pr[0] = PPM_GETR( *ip1 );
|
|
Packit |
78deda |
pg[0] = PPM_GETG( *ip1 );
|
|
Packit |
78deda |
pb[0] = PPM_GETB( *ip1 );
|
|
Packit |
78deda |
pr[1] = PPM_GETR( *(ip1-no) );
|
|
Packit |
78deda |
pg[1] = PPM_GETG( *(ip1-no) );
|
|
Packit |
78deda |
pb[1] = PPM_GETB( *(ip1-no) );
|
|
Packit |
78deda |
pr[5] = PPM_GETR( *(ip1+po) );
|
|
Packit |
78deda |
pg[5] = PPM_GETG( *(ip1+po) );
|
|
Packit |
78deda |
pb[5] = PPM_GETB( *(ip1+po) );
|
|
Packit |
78deda |
pr[3] = PPM_GETR( *(ip2) );
|
|
Packit |
78deda |
pg[3] = PPM_GETG( *(ip2) );
|
|
Packit |
78deda |
pb[3] = PPM_GETB( *(ip2) );
|
|
Packit |
78deda |
pr[2] = PPM_GETR( *(ip2-no) );
|
|
Packit |
78deda |
pg[2] = PPM_GETG( *(ip2-no) );
|
|
Packit |
78deda |
pb[2] = PPM_GETB( *(ip2-no) );
|
|
Packit |
78deda |
pr[4] = PPM_GETR( *(ip2+po) );
|
|
Packit |
78deda |
pg[4] = PPM_GETG( *(ip2+po) );
|
|
Packit |
78deda |
pb[4] = PPM_GETB( *(ip2+po) );
|
|
Packit |
78deda |
pr[6] = PPM_GETR( *(ip0+po) );
|
|
Packit |
78deda |
pg[6] = PPM_GETG( *(ip0+po) );
|
|
Packit |
78deda |
pb[6] = PPM_GETB( *(ip0+po) );
|
|
Packit |
78deda |
pr[8] = PPM_GETR( *(ip0-no) );
|
|
Packit |
78deda |
pg[8] = PPM_GETG( *(ip0-no) );
|
|
Packit |
78deda |
pb[8] = PPM_GETB( *(ip0-no) );
|
|
Packit |
78deda |
pr[7] = PPM_GETR( *(ip0) );
|
|
Packit |
78deda |
pg[7] = PPM_GETG( *(ip0) );
|
|
Packit |
78deda |
pb[7] = PPM_GETB( *(ip0) );
|
|
Packit |
78deda |
r = (*atfunc)(pr);
|
|
Packit |
78deda |
g = (*atfunc)(pg);
|
|
Packit |
78deda |
b = (*atfunc)(pb);
|
|
Packit |
78deda |
PPM_ASSIGN( *op, r, g, b );
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
|
|
Packit |
78deda |
if (irow1 == irows[2]) {
|
|
Packit |
78deda |
irow1 = irows[0];
|
|
Packit |
78deda |
irow2 = irows[1];
|
|
Packit |
78deda |
irow0 = irows[2];
|
|
Packit |
78deda |
} else if (irow1 == irows[1]) {
|
|
Packit |
78deda |
irow2 = irows[0];
|
|
Packit |
78deda |
irow0 = irows[1];
|
|
Packit |
78deda |
irow1 = irows[2];
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
else /* must be at irows[0] */
|
|
Packit |
78deda |
{
|
|
Packit |
78deda |
irow0 = irows[0];
|
|
Packit |
78deda |
irow1 = irows[1];
|
|
Packit |
78deda |
irow2 = irows[2];
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
} else {
|
|
Packit |
78deda |
/* Else must be PGM */
|
|
Packit |
78deda |
int p[9]; /* 3x3 neighbor pixel values */
|
|
Packit |
78deda |
int pv;
|
|
Packit |
78deda |
int promote;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
/* we scale maxval to omaxval */
|
|
Packit |
78deda |
promote = ( PNM_FORMAT_TYPE(format) != PNM_FORMAT_TYPE(oformat) );
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for ( row = 0; row < rows; row++ ) {
|
|
Packit |
78deda |
int po,no; /* offsets for left and right colums in 3x3 */
|
|
Packit |
78deda |
xel *ip0, *ip1, *ip2, *op;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (row == 0) {
|
|
Packit |
78deda |
irow0 = irow1;
|
|
Packit |
78deda |
pnm_readpnmrow( ifP, irow1, cols, maxval, format );
|
|
Packit |
78deda |
if ( promote )
|
|
Packit |
78deda |
pnm_promoteformatrow( irow1, cols, maxval,
|
|
Packit |
78deda |
format, maxval, oformat );
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
if (row == (rows-1))
|
|
Packit |
78deda |
irow2 = irow1;
|
|
Packit |
78deda |
else {
|
|
Packit |
78deda |
pnm_readpnmrow( ifP, irow2, cols, maxval, format );
|
|
Packit |
78deda |
if ( promote )
|
|
Packit |
78deda |
pnm_promoteformatrow( irow2, cols, maxval,
|
|
Packit |
78deda |
format, maxval, oformat );
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
for (col = cols-1,po= col>0?1:0,no=0,
|
|
Packit |
78deda |
ip0=irow0,ip1=irow1,ip2=irow2,op=orow;
|
|
Packit |
78deda |
col >= 0;
|
|
Packit |
78deda |
col--,ip0++,ip1++,ip2++,op++, no |= 1,po = col!= 0 ? po : 0) {
|
|
Packit |
78deda |
/* grab 3x3 pixel values */
|
|
Packit |
78deda |
p[0] = PNM_GET1( *ip1 );
|
|
Packit |
78deda |
p[1] = PNM_GET1( *(ip1-no) );
|
|
Packit |
78deda |
p[5] = PNM_GET1( *(ip1+po) );
|
|
Packit |
78deda |
p[3] = PNM_GET1( *(ip2) );
|
|
Packit |
78deda |
p[2] = PNM_GET1( *(ip2-no) );
|
|
Packit |
78deda |
p[4] = PNM_GET1( *(ip2+po) );
|
|
Packit |
78deda |
p[6] = PNM_GET1( *(ip0+po) );
|
|
Packit |
78deda |
p[8] = PNM_GET1( *(ip0-no) );
|
|
Packit |
78deda |
p[7] = PNM_GET1( *(ip0) );
|
|
Packit |
78deda |
pv = (*atfunc)(p);
|
|
Packit |
78deda |
PNM_ASSIGN1( *op, pv );
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
pnm_writepnmrow( stdout, orow, cols, omaxval, oformat, 0 );
|
|
Packit |
78deda |
if (irow1 == irows[2]) {
|
|
Packit |
78deda |
irow1 = irows[0];
|
|
Packit |
78deda |
irow2 = irows[1];
|
|
Packit |
78deda |
irow0 = irows[2];
|
|
Packit |
78deda |
} else if (irow1 == irows[1]) {
|
|
Packit |
78deda |
irow2 = irows[0];
|
|
Packit |
78deda |
irow0 = irows[1];
|
|
Packit |
78deda |
irow1 = irows[2];
|
|
Packit |
78deda |
} else {
|
|
Packit |
78deda |
/* must be at irows[0] */
|
|
Packit |
78deda |
irow0 = irows[0];
|
|
Packit |
78deda |
irow1 = irows[1];
|
|
Packit |
78deda |
irow2 = irows[2];
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
static void
|
|
Packit |
78deda |
verifySame(unsigned int const imageSeq,
|
|
Packit |
78deda |
int const imageCols, int const imageRows,
|
|
Packit |
78deda |
xelval const imageMaxval, int const imageFormat,
|
|
Packit |
78deda |
int const cols, int const rows,
|
|
Packit |
78deda |
xelval const maxval, int const format) {
|
|
Packit |
78deda |
/*----------------------------------------------------------------------------
|
|
Packit |
78deda |
Issue error message and exit the program if the imageXXX arguments don't
|
|
Packit |
78deda |
match the XXX arguments.
|
|
Packit |
78deda |
-----------------------------------------------------------------------------*/
|
|
Packit |
78deda |
if (imageCols != cols)
|
|
Packit |
78deda |
pm_error("Width of Image %u (%d) is not the same as Image 0 (%d)",
|
|
Packit |
78deda |
imageSeq, imageCols, cols);
|
|
Packit |
78deda |
if (imageRows != rows)
|
|
Packit |
78deda |
pm_error("Height of Image %u (%d) is not the same as Image 0 (%d)",
|
|
Packit |
78deda |
imageSeq, imageRows, rows);
|
|
Packit |
78deda |
if (imageMaxval != maxval)
|
|
Packit |
78deda |
pm_error("Maxval of Image %u (%u) is not the same as Image 0 (%u)",
|
|
Packit |
78deda |
imageSeq, imageMaxval, maxval);
|
|
Packit |
78deda |
if (imageFormat != format)
|
|
Packit |
78deda |
pm_error("Format of Image %u is not the same as Image 0",
|
|
Packit |
78deda |
imageSeq);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int (*atfuncs[6]) (int *) = {atfilt0,atfilt1,atfilt2,atfilt3,atfilt4,atfilt5};
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
|
|
Packit |
78deda |
int
|
|
Packit |
78deda |
main(int argc, char *argv[]) {
|
|
Packit |
78deda |
|
|
Packit |
78deda |
FILE * ifP;
|
|
Packit |
78deda |
struct cmdlineInfo cmdline;
|
|
Packit |
78deda |
int eof; /* We've hit the end of the input stream */
|
|
Packit |
78deda |
unsigned int imageSeq; /* Sequence number of image, starting from 0 */
|
|
Packit |
78deda |
|
|
Packit |
78deda |
pnm_init(&argc, argv);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
parseCommandLine(argc, argv, &cmdline);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
ifP = pm_openr(cmdline.inputFileName);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
pnm_readpnminit(ifP, &cols, &rows, &maxval, &format);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (maxval > MXIVAL)
|
|
Packit |
78deda |
pm_error("The maxval of the input image (%d) is too large.\n"
|
|
Packit |
78deda |
"This program's limit is %d.",
|
|
Packit |
78deda |
maxval, MXIVAL);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
oformat = PNM_FORMAT_TYPE(format);
|
|
Packit |
78deda |
/* force output to max precision without forcing new 2-byte format */
|
|
Packit |
78deda |
omaxval = MIN(maxval, PPM_MAXMAXVAL);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
atfunc = atfuncs[atfilt_setup(cmdline.alpha, cmdline.radius,
|
|
Packit |
78deda |
(double)omaxval/(double)maxval)];
|
|
Packit |
78deda |
|
|
Packit |
78deda |
if (oformat < PGM_TYPE) {
|
|
Packit |
78deda |
oformat = RPGM_FORMAT;
|
|
Packit |
78deda |
pm_message("promoting file to PGM");
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
orow = pnm_allocrow(cols);
|
|
Packit |
78deda |
irows[0] = pnm_allocrow(cols);
|
|
Packit |
78deda |
irows[1] = pnm_allocrow(cols);
|
|
Packit |
78deda |
irows[2] = pnm_allocrow(cols);
|
|
Packit |
78deda |
irow0 = irows[0];
|
|
Packit |
78deda |
irow1 = irows[1];
|
|
Packit |
78deda |
irow2 = irows[2];
|
|
Packit |
78deda |
|
|
Packit |
78deda |
eof = FALSE; /* We're already in the middle of the first image */
|
|
Packit |
78deda |
imageSeq = 0;
|
|
Packit |
78deda |
while (!eof) {
|
|
Packit |
78deda |
do_one_frame(ifP);
|
|
Packit |
78deda |
pm_nextimage(ifP, &eof;;
|
|
Packit |
78deda |
if (!eof) {
|
|
Packit |
78deda |
/* Read and validate header of next image */
|
|
Packit |
78deda |
int imageCols, imageRows;
|
|
Packit |
78deda |
xelval imageMaxval;
|
|
Packit |
78deda |
int imageFormat;
|
|
Packit |
78deda |
|
|
Packit |
78deda |
++imageSeq;
|
|
Packit |
78deda |
pnm_readpnminit(ifP, &imageCols, &imageRows,
|
|
Packit |
78deda |
&imageMaxval, &imageFormat);
|
|
Packit |
78deda |
verifySame(imageSeq,
|
|
Packit |
78deda |
imageCols, imageRows, imageMaxval, imageFormat,
|
|
Packit |
78deda |
cols, rows, maxval, format);
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
}
|
|
Packit |
78deda |
|
|
Packit |
78deda |
pnm_freerow(irow0);
|
|
Packit |
78deda |
pnm_freerow(irow1);
|
|
Packit |
78deda |
pnm_freerow(irow2);
|
|
Packit |
78deda |
pnm_freerow(orow);
|
|
Packit |
78deda |
pm_close(ifP);
|
|
Packit |
78deda |
|
|
Packit |
78deda |
return 0;
|
|
Packit |
78deda |
}
|