src/px/pf3.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 *  pf3.c - routines for gaussian and box filtering
 *
 *     8/13/86
 */

#include  "standard.h"

#include  "color.h"

#define  RSCA           1.13    /* square-radius multiplier: sqrt(4/PI) */
#define  TEPS           0.2     /* threshold proximity goal */
#define  REPS           0.1     /* radius 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 int  wrapfilt;           /* wrap filter horizontally? */

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 */

extern double  (*ourbright)();  /* brightness computation function */

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 = 111*CHECKRAD*CHECKRAD/(rad*rad);
        gausstable = (float *)malloc(gtabsiz*sizeof(float));
        if (gausstable == NULL)
                goto memerr;
        d = x_c*y_r*0.25/(rad*rad);
        gausstable[0] = exp(-d);
        for (x = 1; x < gtabsiz; x++)
                if (x*0.1 <= d)
                        gausstable[x] = gausstable[0];
                else
                        gausstable[x] = exp(-x*0.1);
        if (obarsize == 0)
                return;
                                        /* compute integral of filter */
        gaussN = PI*d*exp(-d);                  /* plateau */
        for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
                gaussN += 0.1*2.0*PI*d*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);
        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, offs;
        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+.5) : (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++) {
                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
                        if (offs && !wrapfilt)
                                continue;
                        d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
                        if (d < -0.5) continue;
                        if (d >= 0.5) break;
                        wsum++;
                        addcolor(csum, scan[x+offs]);
                }
        }
        if (wsum > 1) {
                d = 1.0/wsum;
                scalecolor(csum, d);
        }
}


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, offs;
        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++) {
                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
                        if (offs && !wrapfilt)
                                continue;
                        dx = (x_c*(x+.5) - (c+.5))/rad;
                        weight = lookgauss(dx*dx + dy*dy);
                        wsum += weight;
                        copycolor(ctmp, scan[x+offs]);
                        scalecolor(ctmp, weight);
                        addcolor(csum, ctmp);
                }
        }
        weight = 1.0/wsum;
        scalecolor(csum, weight);
}


dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
int  xcent, ycent;
int  ccent, rcent;
{
        double  d;
        int  r, y, offs;
        register int  c, x;
        register float  *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 >= nrows) break;
                gscan = greybar[(rcent+r)%obarsize];
                for (c = -orad; c <= orad; c++) {
                        offs = ccent+c < 0 ? ncols :
                                        ccent+c >= ncols ? -ncols : 0;
                        if (offs && !wrapfilt)
                                continue;
                        x = ringndx[c*c + r*r];
                        if (x < 0) continue;
                        ringsum[x] += gscan[ccent+c+offs];
                        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+.5) : (double)(y - ycent);
                if (d < -0.5) continue;
                if (d >= 0.5) break;
                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
                        offs = x < 0 ? xres : x >= xres ? -xres : 0;
                        if (offs && !wrapfilt)
                                continue;
                        d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
                        if (d < -0.5) continue;
                        if (d >= 0.5) break;
                        sumans(x, y, rcent, ccent,
                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
                }
        }
}


double
pickfilt(p0)                    /* find filter multiplier for p0 */
double  p0;
{
        double  m = 1.0;
        double  t, num, denom, avg, wsum;
        double  mlimit[2];
        int  ilimit = 4.0/TEPS;
        register int  i;
                                /* iterative search for m */
        mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
        do {
                                        /* compute grey weighted average */
                i = RSCA*CHECKRAD*rad*m + .5;
                if (i > orad) i = orad;
                avg = wsum = 0.0;
                while (i--) {
                        t = (i+.5)/(m*rad);
                        t = lookgauss(t*t);
                        avg += t*ringsum[i];
                        wsum += t*ringwt[i];
                }
                if (avg < 1e-20)                /* zero inclusive average */
                        return(1.0);
                avg /= wsum;
                                        /* check threshold */
                denom = m*m/gausstable[0] - p0/avg;
                if (denom <= FTINY) {           /* zero exclusive average */
                        if (m >= mlimit[1]-REPS)
                                break;
                        m = mlimit[1];
                        continue;
                }
                num = p0/avg - 1.0;
                if (num < 0.0) num = -num;
                t = num/denom;
                if (t <= thresh) {
                        if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
                                break;
                } else {
                        if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
                                break;
                }
                t = m;                  /* remember current m */
                                        /* next guesstimate */
                m = sqrt(gausstable[0]*(num/thresh + p0/avg));
                if (m < t) {            /* bound it */
                        if (m <= mlimit[0]+FTINY)
                                m = 0.5*(mlimit[0] + t);
                        mlimit[1] = t;
                } else {
                        if (m >= mlimit[1]-FTINY)
                                m = 0.5*(mlimit[1] + t);
                        mlimit[0] = t;
                }
        } 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  dy2, dx;
        COLOR  pval, ctmp;
        int  ksiz, r, offs;
        double  pc, pr, norm;
        register int  i, c;
        register COLOR  *scan;
        /*
         * This normalization method fails at the picture borders because
         * a different number of input pixels contribute there.
         */
        scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
        copycolor(pval, scan[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;
                dy2 = (pr - (r+.5))/(m*rad);
                dy2 *= dy2;
                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
                        if (!wrapfilt) {
                                if (c < 0) continue;
                                if (c >= ncols) break;
                        }
                        dx = (pc - (c+.5))/(m*rad);
                        norm += warr[i++] = lookgauss(dx*dx + dy2);
                }
        }
        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++) {
                        offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
                        if (offs && !wrapfilt)
                                continue;
                        copycolor(ctmp, pval);
                        dx = norm*warr[i++];
                        scalecolor(ctmp, dx);
                        addcolor(scan[c+offs], ctmp);
                }
        }
}
Revision:	2.14
Committed:	Sat Feb 22 02:07:27 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.13:	+22 -18 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id$";
3	#endif
4	/*
5	* pf3.c - routines for gaussian and box filtering
6	*
7	* 8/13/86
8	*/
9
10	#include "standard.h"
11
12	#include "color.h"
13
14	#define RSCA 1.13 /* square-radius multiplier: sqrt(4/PI) */
15	#define TEPS 0.2 /* threshold proximity goal */
16	#define REPS 0.1 /* radius proximity goal */
17
18	extern double CHECKRAD; /* radius over which gaussian is summed */
19
20	extern double rad; /* output pixel radius for filtering */
21
22	extern double thresh; /* maximum contribution for subpixel */
23
24	extern int nrows; /* number of rows for output */
25	extern int ncols; /* number of columns for output */
26
27	extern int xres, yres; /* resolution of input */
28
29	extern double x_c, y_r; /* conversion factors */
30
31	extern int xrad; /* x search radius */
32	extern int yrad; /* y search radius */
33	extern int xbrad; /* x box size */
34	extern int ybrad; /* y box size */
35
36	extern int barsize; /* size of input scan bar */
37	extern COLOR *scanin; / input scan bar */
38	extern COLOR scanout; / output scan line */
39	extern COLOR *scoutbar; / output scan bar (if thresh > 0) */
40	extern float *greybar; / grey-averaged input values */
41	extern int obarsize; /* size of output scan bar */
42	extern int orad; /* output window radius */
43
44	extern int wrapfilt; /* wrap filter horizontally? */
45
46	extern char *progname;
47
48	float gausstable; / gauss lookup table */
49
50	float ringsum; / sum of ring values */
51	short ringwt; / weight (count) of ring values */
52	short ringndx; / ring index table */
53	float warr; / array of pixel weights */
54
55	extern double (ourbright)(); / brightness computation function */
56
57	double pickfilt();
58
59	#define lookgauss(x) gausstable[(int)(10.*(x)+.5)]
60
61
62	initmask() /* initialize gaussian lookup table */
63	{
64	int gtabsiz;
65	double gaussN;
66	double d;
67	register int x;
68
69	gtabsiz = 111CHECKRADCHECKRAD/(rad*rad);
70	gausstable = (float )malloc(gtabsizsizeof(float));
71	if (gausstable == NULL)
72	goto memerr;
73	d = x_cy_r0.25/(rad*rad);
74	gausstable[0] = exp(-d);
75	for (x = 1; x < gtabsiz; x++)
76	if (x*0.1 <= d)
77	gausstable[x] = gausstable[0];
78	else
79	gausstable[x] = exp(-x*0.1);
80	if (obarsize == 0)
81	return;
82	/* compute integral of filter */
83	gaussN = PIdexp(-d); /* plateau */
84	for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
85	gaussN += 0.12.0PIdexp(-d*d);
86	/* normalize filter */
87	gaussN = x_cy_r/(radrad*gaussN);
88	for (x = 0; x < gtabsiz; x++)
89	gausstable[x] *= gaussN;
90	/* create ring averages table */
91	ringndx = (short )malloc((2oradorad+1)sizeof(short));
92	if (ringndx == NULL)
93	goto memerr;
94	for (x = 2oradorad+1; --x > orad*orad; )
95	ringndx[x] = -1;
96	do
97	ringndx[x] = sqrt((double)x);
98	while (x--);
99	ringsum = (float )malloc((orad+1)sizeof(float));
100	ringwt = (short )malloc((orad+1)sizeof(short));
101	warr = (float )malloc(obarsizeobarsize*sizeof(float));
102	if (ringsum == NULL \| ringwt == 0 \| warr == NULL)
103	goto memerr;
104	return;
105	memerr:
106	fprintf(stderr, "%s: out of memory in initmask\n", progname);
107	quit(1);
108	}
109
110
111	dobox(csum, xcent, ycent, c, r) /* simple box filter */
112	COLOR csum;
113	int xcent, ycent;
114	int c, r;
115	{
116	int wsum;
117	double d;
118	int y;
119	register int x, offs;
120	register COLOR *scan;
121
122	wsum = 0;
123	setcolor(csum, 0.0, 0.0, 0.0);
124	for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
125	if (y < 0) continue;
126	if (y >= yres) break;
127	d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
128	if (d < -0.5) continue;
129	if (d >= 0.5) break;
130	scan = scanin[y%barsize];
131	for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
132	offs = x < 0 ? xres : x >= xres ? -xres : 0;
133	if (offs && !wrapfilt)
134	continue;
135	d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
136	if (d < -0.5) continue;
137	if (d >= 0.5) break;
138	wsum++;
139	addcolor(csum, scan[x+offs]);
140	}
141	}
142	if (wsum > 1) {
143	d = 1.0/wsum;
144	scalecolor(csum, d);
145	}
146	}
147
148
149	dogauss(csum, xcent, ycent, c, r) /* gaussian filter */
150	COLOR csum;
151	int xcent, ycent;
152	int c, r;
153	{
154	double dy, dx, weight, wsum;
155	COLOR ctmp;
156	int y;
157	register int x, offs;
158	register COLOR *scan;
159
160	wsum = FTINY;
161	setcolor(csum, 0.0, 0.0, 0.0);
162	for (y = ycent-yrad; y <= ycent+yrad; y++) {
163	if (y < 0) continue;
164	if (y >= yres) break;
165	dy = (y_r*(y+.5) - (r+.5))/rad;
166	scan = scanin[y%barsize];
167	for (x = xcent-xrad; x <= xcent+xrad; x++) {
168	offs = x < 0 ? xres : x >= xres ? -xres : 0;
169	if (offs && !wrapfilt)
170	continue;
171	dx = (x_c*(x+.5) - (c+.5))/rad;
172	weight = lookgauss(dxdx + dydy);
173	wsum += weight;
174	copycolor(ctmp, scan[x+offs]);
175	scalecolor(ctmp, weight);
176	addcolor(csum, ctmp);
177	}
178	}
179	weight = 1.0/wsum;
180	scalecolor(csum, weight);
181	}
182
183
184	dothresh(xcent, ycent, ccent, rcent) /* gaussian threshold filter */
185	int xcent, ycent;
186	int ccent, rcent;
187	{
188	double d;
189	int r, y, offs;
190	register int c, x;
191	register float *gscan;
192	/* compute ring sums */
193	bzero((char )ringsum, (orad+1)sizeof(float));
194	bzero((char )ringwt, (orad+1)sizeof(short));
195	for (r = -orad; r <= orad; r++) {
196	if (rcent+r < 0) continue;
197	if (rcent+r >= nrows) break;
198	gscan = greybar[(rcent+r)%obarsize];
199	for (c = -orad; c <= orad; c++) {
200	offs = ccent+c < 0 ? ncols :
201	ccent+c >= ncols ? -ncols : 0;
202	if (offs && !wrapfilt)
203	continue;
204	x = ringndx[cc + rr];
205	if (x < 0) continue;
206	ringsum[x] += gscan[ccent+c+offs];
207	ringwt[x]++;
208	}
209	}
210	/* filter each subpixel */
211	for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
212	if (y < 0) continue;
213	if (y >= yres) break;
214	d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
215	if (d < -0.5) continue;
216	if (d >= 0.5) break;
217	for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
218	offs = x < 0 ? xres : x >= xres ? -xres : 0;
219	if (offs && !wrapfilt)
220	continue;
221	d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
222	if (d < -0.5) continue;
223	if (d >= 0.5) break;
224	sumans(x, y, rcent, ccent,
225	pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
226	}
227	}
228	}
229
230
231	double
232	pickfilt(p0) /* find filter multiplier for p0 */
233	double p0;
234	{
235	double m = 1.0;
236	double t, num, denom, avg, wsum;
237	double mlimit[2];
238	int ilimit = 4.0/TEPS;
239	register int i;
240	/* iterative search for m */
241	mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
242	do {
243	/* compute grey weighted average */
244	i = RSCACHECKRADrad*m + .5;
245	if (i > orad) i = orad;
246	avg = wsum = 0.0;
247	while (i--) {
248	t = (i+.5)/(m*rad);
249	t = lookgauss(t*t);
250	avg += t*ringsum[i];
251	wsum += t*ringwt[i];
252	}
253	if (avg < 1e-20) /* zero inclusive average */
254	return(1.0);
255	avg /= wsum;
256	/* check threshold */
257	denom = m*m/gausstable[0] - p0/avg;
258	if (denom <= FTINY) { /* zero exclusive average */
259	if (m >= mlimit[1]-REPS)
260	break;
261	m = mlimit[1];
262	continue;
263	}
264	num = p0/avg - 1.0;
265	if (num < 0.0) num = -num;
266	t = num/denom;
267	if (t <= thresh) {
268	if (m <= mlimit[0]+REPS \|\| (thresh-t)/thresh <= TEPS)
269	break;
270	} else {
271	if (m >= mlimit[1]-REPS \|\| (t-thresh)/thresh <= TEPS)
272	break;
273	}
274	t = m; /* remember current m */
275	/* next guesstimate */
276	m = sqrt(gausstable[0]*(num/thresh + p0/avg));
277	if (m < t) { /* bound it */
278	if (m <= mlimit[0]+FTINY)
279	m = 0.5*(mlimit[0] + t);
280	mlimit[1] = t;
281	} else {
282	if (m >= mlimit[1]-FTINY)
283	m = 0.5*(mlimit[1] + t);
284	mlimit[0] = t;
285	}
286	} while (--ilimit > 0);
287	return(m);
288	}
289
290
291	sumans(px, py, rcent, ccent, m) /* sum input pixel to output */
292	int px, py;
293	int rcent, ccent;
294	double m;
295	{
296	double dy2, dx;
297	COLOR pval, ctmp;
298	int ksiz, r, offs;
299	double pc, pr, norm;
300	register int i, c;
301	register COLOR *scan;
302	/*
303	* This normalization method fails at the picture borders because
304	* a different number of input pixels contribute there.
305	*/
306	scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
307	copycolor(pval, scan[px]);
308	pc = x_c*(px+.5);
309	pr = y_r*(py+.5);
310	ksiz = CHECKRADmrad + 1;
311	if (ksiz > orad) ksiz = orad;
312	/* compute normalization */
313	norm = 0.0;
314	i = 0;
315	for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
316	if (r < 0) continue;
317	if (r >= nrows) break;
318	dy2 = (pr - (r+.5))/(m*rad);
319	dy2 *= dy2;
320	for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
321	if (!wrapfilt) {
322	if (c < 0) continue;
323	if (c >= ncols) break;
324	}
325	dx = (pc - (c+.5))/(m*rad);
326	norm += warr[i++] = lookgauss(dx*dx + dy2);
327	}
328	}
329	norm = 1.0/norm;
330	if (x_c < 1.0) norm *= x_c;
331	if (y_r < 1.0) norm *= y_r;
332	/* sum pixels */
333	i = 0;
334	for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
335	if (r < 0) continue;
336	if (r >= nrows) break;
337	scan = scoutbar[r%obarsize];
338	for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
339	offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
340	if (offs && !wrapfilt)
341	continue;
342	copycolor(ctmp, pval);
343	dx = norm*warr[i++];
344	scalecolor(ctmp, dx);
345	addcolor(scan[c+offs], ctmp);
346	}
347	}
348	}