src/px/pf3.c

/* Copyright (c) 1992 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  pf3.c - routines for gaussian and box filtering
 *
 *     8/13/86
 */

#include  "standard.h"

#include  "color.h"

#define  TEPS           0.25    /* threshold proximity goal */

extern double  CHECKRAD;        /* radius over which gaussian is summed */

extern double  rad;             /* output pixel radius for filtering */

extern double  thresh;          /* maximum contribution for subpixel */

extern int  nrows;              /* number of rows for output */
extern int  ncols;              /* number of columns for output */

extern int  xres, yres;         /* resolution of input */

extern double  x_c, y_r;        /* conversion factors */

extern int  xrad;               /* x search radius */
extern int  yrad;               /* y search radius */
extern int  xbrad;              /* x box size */
extern int  ybrad;              /* y box size */

extern int  barsize;            /* size of input scan bar */
extern COLOR  **scanin;         /* input scan bar */
extern COLOR  *scanout;         /* output scan line */
extern COLOR  **scoutbar;       /* output scan bar (if thresh > 0) */
extern float  **greybar;        /* grey-averaged input values */
extern int  obarsize;           /* size of output scan bar */
extern int  orad;               /* output window radius */

extern char  *progname;

float  *gausstable;             /* gauss lookup table */

float  *ringsum;                /* sum of ring values */
short  *ringwt;                 /* weight (count) of ring values */
short  *ringndx;                /* ring index table */
float  *warr;                   /* array of pixel weights */

double  pickfilt();

#define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]


initmask()                      /* initialize gaussian lookup table */
{
        int  gtabsiz;
        double  gaussN;
        double  d;
        register int  x;

        gtabsiz = 150*CHECKRAD;
        gausstable = (float *)malloc(gtabsiz*sizeof(float));
        if (gausstable == NULL)
                goto memerr;
        d = x_c*y_r*0.25/(rad*rad);
        gausstable[0] = exp(-d)/sqrt(d);
        for (x = 1; x < gtabsiz; x++)
                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
                        gausstable[x] = gausstable[0];
        if (obarsize == 0)
                return;
                                        /* compute integral of filter */
        gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
        for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
                gaussN += 0.1*2.0*PI*exp(-d*d);
                                        /* normalize filter */
        gaussN = x_c*y_r/(rad*rad*gaussN);
        for (x = 0; x < gtabsiz; x++)
                gausstable[x] *= gaussN;
                                        /* create ring averages table */
        ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
        if (ringndx == NULL)
                goto memerr;
        for (x = 2*orad*orad+1; --x > orad*orad; )
                ringndx[x] = -1;
        do
                ringndx[x] = sqrt((double)x) + 0.5;
        while (x--);
        ringsum = (float *)malloc((orad+1)*sizeof(float));
        ringwt = (short *)malloc((orad+1)*sizeof(short));
        warr = (float *)malloc(obarsize*obarsize*sizeof(float));
        if (ringsum == NULL | ringwt == 0 | warr == NULL)
                goto memerr;
        return;
memerr:
        fprintf(stderr, "%s: out of memory in initmask\n", progname);
        quit(1);
}


dobox(csum, xcent, ycent, c, r)                 /* simple box filter */
COLOR  csum;
int  xcent, ycent;
int  c, r;
{
        int  wsum;
        double  d;
        int  y;
        register int  x;
        register COLOR  *scan;
        
        wsum = 0;
        setcolor(csum, 0.0, 0.0, 0.0);
        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
                if (y < 0) continue;
                if (y >= yres) break;
                d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
                if (d < -0.5) continue;
                if (d >= 0.5) break;
                scan = scanin[y%barsize];
                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
                        if (x < 0) continue;
                        if (x >= xres) break;
                        d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
                        if (d < -0.5) continue;
                        if (d >= 0.5) break;
                        wsum++;
                        addcolor(csum, scan[x]);
                }
        }
        if (wsum > 1)
                scalecolor(csum, 1.0/wsum);
}


dogauss(csum, xcent, ycent, c, r)               /* gaussian filter */
COLOR  csum;
int  xcent, ycent;
int  c, r;
{
        double  dy, dx, weight, wsum;
        COLOR  ctmp;
        int  y;
        register int  x;
        register COLOR  *scan;

        wsum = FTINY;
        setcolor(csum, 0.0, 0.0, 0.0);
        for (y = ycent-yrad; y <= ycent+yrad; y++) {
                if (y < 0) continue;
                if (y >= yres) break;
                dy = (y_r*(y+.5) - (r+.5))/rad;
                scan = scanin[y%barsize];
                for (x = xcent-xrad; x <= xcent+xrad; x++) {
                        if (x < 0) continue;
                        if (x >= xres) break;
                        dx = (x_c*(x+.5) - (c+.5))/rad;
                        weight = lookgauss(dx*dx + dy*dy);
                        wsum += weight;
                        copycolor(ctmp, scan[x]);
                        scalecolor(ctmp, weight);
                        addcolor(csum, ctmp);
                }
        }
        scalecolor(csum, 1.0/wsum);
}


dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
int  xcent, ycent;
int  ccent, rcent;
{
        double  d;
        int  r, y;
        register int  c, x;
        register float  *gscan;
#define pval gscan
                                        /* compute ring sums */
        bzero((char *)ringsum, (orad+1)*sizeof(float));
        bzero((char *)ringwt, (orad+1)*sizeof(short));
        for (r = -orad; r <= orad; r++) {
                if (rcent+r < 0) continue;
                if (rcent+r >= ncols) break;
                gscan = greybar[(rcent+r)%obarsize];
                for (c = -orad; c <= orad; c++) {
                        if (ccent+c < 0) continue;
                        if (ccent+c >= ncols) break;
                        x = ringndx[c*c + r*r];
                        if (x < 0) continue;
                        ringsum[x] += gscan[ccent+c];
                        ringwt[x]++;
                }
        }
                                        /* filter each subpixel */
        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
                if (y < 0) continue;
                if (y >= yres) break;
                d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
                if (d < -0.5) continue;
                if (d >= 0.5) break;
                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
                        if (x < 0) continue;
                        if (x >= xres) break;
                        d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
                        if (d < -0.5) continue;
                        if (d >= 0.5) break;
                        pval = scanin[y%barsize][x];
                        sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
                }
        }
#undef pval
}


double
pickfilt(p0)                    /* find filter multiplier for p0 */
double  p0;
{
        double  m = 1.0;
        double  t, avg, wsum;
        int  ilimit = 5/TEPS;
        register int  i;
                                /* iterative search for m */
        do {
                                        /* compute grey weighted average */
                i = 1.25*CHECKRAD*rad*m + .5;
                if (i > orad) i = orad+1;
                avg = wsum = 0.0;
                while (i--) {
                        t = i/(m*rad);
                        t = lookgauss(t*t);
                        avg += t*ringsum[i];
                        wsum += t*ringwt[i];
                }
                avg /= wsum;
                                        /* check threshold */
                t = p0 - avg;
                if (t < 0.0) t = -t;
                t *= gausstable[0]/(m*m*avg);
                if (t <= thresh && (m <= 1.0+FTINY ||
                                (thresh-t)/thresh <= TEPS))
                        break;
                if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
                                (t-thresh)/thresh <= TEPS))
                        break;
                                        /* next guesstimate */
                m *= sqrt(t/thresh);
                if (m < 1.0) m = 1.0;
                else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
        } while (--ilimit > 0);
        return(m);
}


sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
int  px, py;
int  rcent, ccent;
double  m;
{
        double  dy, dx;
        COLOR  pval, ctmp;
        int  ksiz, r;
        double  pc, pr, norm;
        register int  i, c;
        register COLOR  *scan;

        copycolor(pval, scanin[py%barsize][px]);
        pc = x_c*(px+.5);
        pr = y_r*(py+.5);
        ksiz = CHECKRAD*m*rad + 1;
        if (ksiz > orad) ksiz = orad;
                                                /* compute normalization */
        norm = 0.0;
        i = 0;
        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
                if (r < 0) continue;
                if (r >= nrows) break;
                dy = (pr - (r+.5))/(m*rad);
                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
                        if (c < 0) continue;
                        if (c >= ncols) break;
                        dx = (pc - (c+.5))/(m*rad);
                        norm += warr[i++] = lookgauss(dx*dx + dy*dy);
                }
        }
        norm = 1.0/norm;
        if (x_c < 1.0) norm *= x_c;
        if (y_r < 1.0) norm *= y_r;
                                                /* sum pixels */
        i = 0;
        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
                if (r < 0) continue;
                if (r >= nrows) break;
                scan = scoutbar[r%obarsize];
                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
                        if (c < 0) continue;
                        if (c >= ncols) break;
                        copycolor(ctmp, pval);
                        dx = norm*warr[i++];
                        scalecolor(ctmp, dx);
                        addcolor(scan[c], ctmp);
                }
        }
}
Revision:	2.7
Committed:	Fri Jun 25 15:01:20 1993 UTC (30 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.6:	+202 -23 lines
Log Message:	added -m option for sample smoothing
#	User	Rev	Content
1	greg	2.2	/* Copyright (c) 1992 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* pf3.c - routines for gaussian and box filtering
9			*
10			* 8/13/86
11			*/
12
13	greg	2.7	#include "standard.h"
14	greg	1.1
15			#include "color.h"
16
17	greg	2.7	#define TEPS 0.25 /* threshold proximity goal */
18	greg	1.1
19	greg	2.7	extern double CHECKRAD; /* radius over which gaussian is summed */
20
21	greg	1.1	extern double rad; /* output pixel radius for filtering */
22
23	greg	2.7	extern double thresh; /* maximum contribution for subpixel */
24
25	greg	1.1	extern int nrows; /* number of rows for output */
26			extern int ncols; /* number of columns for output */
27
28			extern int xres, yres; /* resolution of input */
29
30			extern double x_c, y_r; /* conversion factors */
31
32	greg	2.7	extern int xrad; /* x search radius */
33			extern int yrad; /* y search radius */
34			extern int xbrad; /* x box size */
35			extern int ybrad; /* y box size */
36	greg	1.1
37			extern int barsize; /* size of input scan bar */
38			extern COLOR *scanin; / input scan bar */
39			extern COLOR scanout; / output scan line */
40	greg	2.7	extern COLOR *scoutbar; / output scan bar (if thresh > 0) */
41			extern float *greybar; / grey-averaged input values */
42			extern int obarsize; /* size of output scan bar */
43			extern int orad; /* output window radius */
44	greg	1.1
45			extern char *progname;
46
47	greg	2.5	float gausstable; / gauss lookup table */
48	greg	1.1
49	greg	2.7	float ringsum; / sum of ring values */
50			short ringwt; / weight (count) of ring values */
51			short ringndx; / ring index table */
52			float warr; / array of pixel weights */
53
54			double pickfilt();
55
56	greg	2.6	#define lookgauss(x) gausstable[(int)(10.*(x)+.5)]
57	greg	1.1
58
59			initmask() /* initialize gaussian lookup table */
60			{
61	greg	2.7	int gtabsiz;
62			double gaussN;
63	greg	2.5	double d;
64	greg	1.1	register int x;
65
66	greg	2.7	gtabsiz = 150*CHECKRAD;
67			gausstable = (float )malloc(gtabsizsizeof(float));
68			if (gausstable == NULL)
69			goto memerr;
70			d = x_cy_r0.25/(rad*rad);
71	greg	2.5	gausstable[0] = exp(-d)/sqrt(d);
72	greg	2.7	for (x = 1; x < gtabsiz; x++)
73	greg	2.5	if ((gausstable[x] = exp(-x0.1)/sqrt(x0.1)) > gausstable[0])
74			gausstable[x] = gausstable[0];
75	greg	2.7	if (obarsize == 0)
76			return;
77			/* compute integral of filter */
78			gaussN = PIsqrt(d)exp(-d); /* plateau */
79			for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
80			gaussN += 0.12.0PIexp(-dd);
81			/* normalize filter */
82			gaussN = x_cy_r/(radrad*gaussN);
83			for (x = 0; x < gtabsiz; x++)
84			gausstable[x] *= gaussN;
85			/* create ring averages table */
86			ringndx = (short )malloc((2oradorad+1)sizeof(short));
87			if (ringndx == NULL)
88			goto memerr;
89			for (x = 2oradorad+1; --x > orad*orad; )
90			ringndx[x] = -1;
91			do
92			ringndx[x] = sqrt((double)x) + 0.5;
93			while (x--);
94			ringsum = (float )malloc((orad+1)sizeof(float));
95			ringwt = (short )malloc((orad+1)sizeof(short));
96			warr = (float )malloc(obarsizeobarsize*sizeof(float));
97			if (ringsum == NULL \| ringwt == 0 \| warr == NULL)
98			goto memerr;
99			return;
100			memerr:
101			fprintf(stderr, "%s: out of memory in initmask\n", progname);
102			quit(1);
103	greg	1.1	}
104
105
106			dobox(csum, xcent, ycent, c, r) /* simple box filter */
107			COLOR csum;
108			int xcent, ycent;
109			int c, r;
110			{
111	greg	2.6	int wsum;
112			double d;
113			int y;
114	greg	1.1	register int x;
115	greg	2.2	register COLOR *scan;
116	greg	2.7
117	greg	1.1	wsum = 0;
118			setcolor(csum, 0.0, 0.0, 0.0);
119	greg	2.7	for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
120	greg	2.4	if (y < 0) continue;
121			if (y >= yres) break;
122			d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
123	greg	2.7	if (d < -0.5) continue;
124			if (d >= 0.5) break;
125	greg	1.1	scan = scanin[y%barsize];
126	greg	2.7	for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
127	greg	2.4	if (x < 0) continue;
128			if (x >= xres) break;
129			d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
130	greg	2.7	if (d < -0.5) continue;
131			if (d >= 0.5) break;
132	greg	1.1	wsum++;
133			addcolor(csum, scan[x]);
134			}
135			}
136			if (wsum > 1)
137			scalecolor(csum, 1.0/wsum);
138			}
139
140
141			dogauss(csum, xcent, ycent, c, r) /* gaussian filter */
142			COLOR csum;
143			int xcent, ycent;
144			int c, r;
145			{
146	greg	2.6	double dy, dx, weight, wsum;
147			COLOR ctmp;
148			int y;
149	greg	1.1	register int x;
150	greg	2.2	register COLOR *scan;
151	greg	1.1
152			wsum = FTINY;
153			setcolor(csum, 0.0, 0.0, 0.0);
154	greg	1.2	for (y = ycent-yrad; y <= ycent+yrad; y++) {
155	greg	2.4	if (y < 0) continue;
156			if (y >= yres) break;
157	greg	1.2	dy = (y_r*(y+.5) - (r+.5))/rad;
158	greg	1.1	scan = scanin[y%barsize];
159	greg	1.2	for (x = xcent-xrad; x <= xcent+xrad; x++) {
160	greg	2.4	if (x < 0) continue;
161			if (x >= xres) break;
162	greg	1.2	dx = (x_c*(x+.5) - (c+.5))/rad;
163	greg	2.6	weight = lookgauss(dxdx + dydy);
164	greg	1.1	wsum += weight;
165			copycolor(ctmp, scan[x]);
166			scalecolor(ctmp, weight);
167			addcolor(csum, ctmp);
168			}
169			}
170			scalecolor(csum, 1.0/wsum);
171	greg	2.7	}
172
173
174			dothresh(xcent, ycent, ccent, rcent) /* gaussian threshold filter */
175			int xcent, ycent;
176			int ccent, rcent;
177			{
178			double d;
179			int r, y;
180			register int c, x;
181			register float *gscan;
182			#define pval gscan
183			/* compute ring sums */
184			bzero((char )ringsum, (orad+1)sizeof(float));
185			bzero((char )ringwt, (orad+1)sizeof(short));
186			for (r = -orad; r <= orad; r++) {
187			if (rcent+r < 0) continue;
188			if (rcent+r >= ncols) break;
189			gscan = greybar[(rcent+r)%obarsize];
190			for (c = -orad; c <= orad; c++) {
191			if (ccent+c < 0) continue;
192			if (ccent+c >= ncols) break;
193			x = ringndx[cc + rr];
194			if (x < 0) continue;
195			ringsum[x] += gscan[ccent+c];
196			ringwt[x]++;
197			}
198			}
199			/* filter each subpixel */
200			for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
201			if (y < 0) continue;
202			if (y >= yres) break;
203			d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
204			if (d < -0.5) continue;
205			if (d >= 0.5) break;
206			for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
207			if (x < 0) continue;
208			if (x >= xres) break;
209			d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
210			if (d < -0.5) continue;
211			if (d >= 0.5) break;
212			pval = scanin[y%barsize][x];
213			sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
214			}
215			}
216			#undef pval
217			}
218
219
220			double
221			pickfilt(p0) /* find filter multiplier for p0 */
222			double p0;
223			{
224			double m = 1.0;
225			double t, avg, wsum;
226			int ilimit = 5/TEPS;
227			register int i;
228			/* iterative search for m */
229			do {
230			/* compute grey weighted average */
231			i = 1.25CHECKRADrad*m + .5;
232			if (i > orad) i = orad+1;
233			avg = wsum = 0.0;
234			while (i--) {
235			t = i/(m*rad);
236			t = lookgauss(t*t);
237			avg += t*ringsum[i];
238			wsum += t*ringwt[i];
239			}
240			avg /= wsum;
241			/* check threshold */
242			t = p0 - avg;
243			if (t < 0.0) t = -t;
244			t = gausstable[0]/(mm*avg);
245			if (t <= thresh && (m <= 1.0+FTINY \|\|
246			(thresh-t)/thresh <= TEPS))
247			break;
248			if (t > thresh && (mCHECKRADrad >= orad-FTINY \|\|
249			(t-thresh)/thresh <= TEPS))
250			break;
251			/* next guesstimate */
252			m *= sqrt(t/thresh);
253			if (m < 1.0) m = 1.0;
254			else if (mCHECKRADrad > orad) m = orad/rad/CHECKRAD;
255			} while (--ilimit > 0);
256			return(m);
257			}
258
259
260			sumans(px, py, rcent, ccent, m) /* sum input pixel to output */
261			int px, py;
262			int rcent, ccent;
263			double m;
264			{
265			double dy, dx;
266			COLOR pval, ctmp;
267			int ksiz, r;
268			double pc, pr, norm;
269			register int i, c;
270			register COLOR *scan;
271
272			copycolor(pval, scanin[py%barsize][px]);
273			pc = x_c*(px+.5);
274			pr = y_r*(py+.5);
275			ksiz = CHECKRADmrad + 1;
276			if (ksiz > orad) ksiz = orad;
277			/* compute normalization */
278			norm = 0.0;
279			i = 0;
280			for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
281			if (r < 0) continue;
282			if (r >= nrows) break;
283			dy = (pr - (r+.5))/(m*rad);
284			for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
285			if (c < 0) continue;
286			if (c >= ncols) break;
287			dx = (pc - (c+.5))/(m*rad);
288			norm += warr[i++] = lookgauss(dxdx + dydy);
289			}
290			}
291			norm = 1.0/norm;
292			if (x_c < 1.0) norm *= x_c;
293			if (y_r < 1.0) norm *= y_r;
294			/* sum pixels */
295			i = 0;
296			for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
297			if (r < 0) continue;
298			if (r >= nrows) break;
299			scan = scoutbar[r%obarsize];
300			for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
301			if (c < 0) continue;
302			if (c >= ncols) break;
303			copycolor(ctmp, pval);
304			dx = norm*warr[i++];
305			scalecolor(ctmp, dx);
306			addcolor(scan[c], ctmp);
307			}
308			}
309	greg	1.1	}