ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/px/pf3.c

Comparing ray/src/px/pf3.c (file contents):
Revision 2.9 by greg, Fri Jul 9 10:27:52 1993 UTC vs.
Revision 2.21 by greg, Fri Dec 8 17:56:26 2023 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1992 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		*  pf3.c - routines for gaussian and box filtering
6		*
7		*     8/13/86
8		*/
9
10	<	#include  "standard.h"
10	>	#include  "pfilt.h"
11
12	<	#include  "color.h"
16	<
12	>	#define  RSCA           1.13    /* square-radius multiplier: sqrt(4/PI) */
13		#define  TEPS           0.2     /* threshold proximity goal */
14	+	#define  REPS           0.1     /* radius proximity goal */
15
19	–	extern double  CHECKRAD;        /* radius over which gaussian is summed */
20	–
21	–	extern double  rad;             /* output pixel radius for filtering */
22	–
23	–	extern double  thresh;          /* maximum contribution for subpixel */
24	–
25	–	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	–	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	–
37	–	extern int  barsize;            /* size of input scan bar */
38	–	extern COLOR  **scanin;         /* input scan bar */
39	–	extern COLOR  *scanout;         /* output scan line */
40	–	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	–
45	–	extern char  *progname;
46	–
16		float  *gausstable;             /* gauss lookup table */
17
18		float  *ringsum;                /* sum of ring values */
#	Line 51 | Line 20 | short  *ringwt;                        /* weight (count) of ring values */
20		short  *ringndx;                /* ring index table */
21		float  *warr;                   /* array of pixel weights */
22
23	<	double  pickfilt();
23	>	#define  lookgauss(x)           gausstable[(int)(20.*(x)+.5)]
24
25	<	#define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
25	>	static double pickfilt(double  p0);
26	>	static void sumans(int  px, int  py, int  rcent, int  ccent, double  m);
27
28
29	<	initmask()                      /* initialize gaussian lookup table */
29	>	void
30	>	initmask(void)                  /* initialize gaussian lookup table */
31		{
32		int  gtabsiz;
33		double  gaussN;
34		double  d;
35	<	register int  x;
35	>	int  x;
36
37	<	gtabsiz = 150*CHECKRAD;
37	>	gtabsiz = 444*CHECKRAD*CHECKRAD;
38		gausstable = (float *)malloc(gtabsiz*sizeof(float));
39		if (gausstable == NULL)
40		goto memerr;
41		d = x_c*y_r*0.25/(rad*rad);
42	<	gausstable[0] = exp(-d)/sqrt(d);
42	>	gausstable[0] = exp(-d);
43		for (x = 1; x < gtabsiz; x++)
44	<	if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
44	>	if (x*0.05 <= d)
45		gausstable[x] = gausstable[0];
46	+	else
47	+	gausstable[x] = exp(-x*0.05);
48		if (obarsize == 0)
49		return;
50		/* compute integral of filter */
51	<	gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
52	<	for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
53	<	gaussN += 0.1*2.0*PI*exp(-d*d);
51	>	gaussN = PI*d*exp(-d);                  /* plateau */
52	>	for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
53	>	gaussN += 0.1*2.0*PI*d*exp(-d*d);
54		/* normalize filter */
55		gaussN = x_c*y_r/(rad*rad*gaussN);
56		for (x = 0; x < gtabsiz; x++)
#	Line 94 | Line 67 | initmask()                     /* initialize gaussian lookup table */
67		ringsum = (float *)malloc((orad+1)*sizeof(float));
68		ringwt = (short *)malloc((orad+1)*sizeof(short));
69		warr = (float *)malloc(obarsize*obarsize*sizeof(float));
70	<	if (ringsum == NULL | ringwt == 0 | warr == NULL)
70	>	if ((ringsum == NULL) | (ringwt == 0) | (warr == NULL))
71		goto memerr;
72		return;
73		memerr:
#	Line 103 | Line 76 | memerr:
76		}
77
78
79	<	dobox(csum, xcent, ycent, c, r)                 /* simple box filter */
80	<	COLOR  csum;
81	<	int  xcent, ycent;
82	<	int  c, r;
79	>	void
80	>	dobox(                  /* simple box filter */
81	>	SCOLOR  csum,
82	>	int  xcent,
83	>	int  ycent,
84	>	int  c,
85	>	int  r
86	>	)
87		{
88		int  wsum;
89		double  d;
90		int  y;
91	<	register int  x;
92	<	register COLOR  *scan;
91	>	int  x, offs;
92	>	COLORV  *scan, *scp;
93
94		wsum = 0;
95	<	setcolor(csum, 0.0, 0.0, 0.0);
95	>	scolorblack(csum);
96		for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
97		if (y < 0) continue;
98		if (y >= yres) break;
99	<	d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
99	>	d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
100		if (d < -0.5) continue;
101		if (d >= 0.5) break;
102		scan = scanin[y%barsize];
103		for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
104	<	if (x < 0) continue;
105	<	if (x >= xres) break;
106	<	d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
104	>	offs = x < 0 ? xres : x >= xres ? -xres : 0;
105	>	if (offs && !wrapfilt)
106	>	continue;
107	>	d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
108		if (d < -0.5) continue;
109		if (d >= 0.5) break;
110		wsum++;
111	<	addcolor(csum, scan[x]);
111	>	scp = scan + (x+offs)*NCSAMP;
112	>	saddscolor(csum, scp);
113		}
114		}
115	<	if (wsum > 1)
116	<	scalecolor(csum, 1.0/wsum);
115	>	if (wsum > 1) {
116	>	d = 1.0/wsum;
117	>	scalescolor(csum, d);
118	>	}
119		}
120
121
122	<	dogauss(csum, xcent, ycent, c, r)               /* gaussian filter */
123	<	COLOR  csum;
124	<	int  xcent, ycent;
125	<	int  c, r;
122	>	void
123	>	dogauss(                /* gaussian filter */
124	>	SCOLOR  csum,
125	>	int  xcent,
126	>	int  ycent,
127	>	int  c,
128	>	int  r
129	>	)
130		{
131		double  dy, dx, weight, wsum;
132	<	COLOR  ctmp;
132	>	SCOLOR  ctmp;
133		int  y;
134	<	register int  x;
135	<	register COLOR  *scan;
134	>	int  x, offs;
135	>	COLORV  *scan;
136
137		wsum = FTINY;
138	<	setcolor(csum, 0.0, 0.0, 0.0);
138	>	scolorblack(csum);
139		for (y = ycent-yrad; y <= ycent+yrad; y++) {
140		if (y < 0) continue;
141		if (y >= yres) break;
142		dy = (y_r*(y+.5) - (r+.5))/rad;
143		scan = scanin[y%barsize];
144		for (x = xcent-xrad; x <= xcent+xrad; x++) {
145	<	if (x < 0) continue;
146	<	if (x >= xres) break;
145	>	offs = x < 0 ? xres : x >= xres ? -xres : 0;
146	>	if (offs && !wrapfilt)
147	>	continue;
148		dx = (x_c*(x+.5) - (c+.5))/rad;
149		weight = lookgauss(dx*dx + dy*dy);
150		wsum += weight;
151	<	copycolor(ctmp, scan[x]);
152	<	scalecolor(ctmp, weight);
153	<	addcolor(csum, ctmp);
151	>	copyscolor(ctmp, scan+(x+offs)*NCSAMP);
152	>	scalescolor(ctmp, weight);
153	>	saddscolor(csum, ctmp);
154		}
155		}
156	<	scalecolor(csum, 1.0/wsum);
156	>	weight = 1.0/wsum;
157	>	scalescolor(csum, weight);
158		}
159
160
161	<	dothresh(xcent, ycent, ccent, rcent)    /* gaussian threshold filter */
162	<	int  xcent, ycent;
163	<	int  ccent, rcent;
161	>	void
162	>	dothresh(       /* gaussian threshold filter */
163	>	int  xcent,
164	>	int  ycent,
165	>	int  ccent,
166	>	int  rcent
167	>	)
168		{
169		double  d;
170	<	int  r, y;
171	<	register int  c, x;
172	<	register float  *gscan;
182	<	#define pval gscan
170	>	int  r, y, offs;
171	>	int  c, x;
172	>	float  *gscan;
173		/* compute ring sums */
174	<	bzero((char *)ringsum, (orad+1)*sizeof(float));
175	<	bzero((char *)ringwt, (orad+1)*sizeof(short));
174	>	memset((char *)ringsum, '\0', (orad+1)*sizeof(float));
175	>	memset((char *)ringwt, '\0', (orad+1)*sizeof(short));
176		for (r = -orad; r <= orad; r++) {
177		if (rcent+r < 0) continue;
178		if (rcent+r >= nrows) break;
179		gscan = greybar[(rcent+r)%obarsize];
180		for (c = -orad; c <= orad; c++) {
181	<	if (ccent+c < 0) continue;
182	<	if (ccent+c >= ncols) break;
181	>	offs = ccent+c < 0 ? ncols :
182	>	ccent+c >= ncols ? -ncols : 0;
183	>	if (offs && !wrapfilt)
184	>	continue;
185		x = ringndx[c*c + r*r];
186		if (x < 0) continue;
187	<	ringsum[x] += gscan[ccent+c];
187	>	ringsum[x] += gscan[ccent+c+offs];
188		ringwt[x]++;
189		}
190		}
#	Line 200 | Line 192 | int  ccent, rcent;
192		for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
193		if (y < 0) continue;
194		if (y >= yres) break;
195	<	d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent);
195	>	d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
196		if (d < -0.5) continue;
197		if (d >= 0.5) break;
198		for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
199	<	if (x < 0) continue;
200	<	if (x >= xres) break;
201	<	d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
199	>	offs = x < 0 ? xres : x >= xres ? -xres : 0;
200	>	if (offs && !wrapfilt)
201	>	continue;
202	>	d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
203		if (d < -0.5) continue;
204		if (d >= 0.5) break;
205	<	pval = scanin[y%barsize][x];
206	<	sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
205	>	sumans(x, y, rcent, ccent,
206	>	pickfilt((*ourbright)(scanin[y%barsize]+(x+offs)*NCSAMP)));
207		}
208		}
216	–	#undef pval
209		}
210
211
212	<	double
213	<	pickfilt(p0)                    /* find filter multiplier for p0 */
214	<	double  p0;
212	>	static double
213	>	pickfilt(                       /* find filter multiplier for p0 */
214	>	double  p0
215	>	)
216		{
217		double  m = 1.0;
218	<	double  t, avg, wsum;
219	<	int  ilimit = 4/TEPS;
220	<	register int  i;
218	>	double  t, num, denom, avg, wsum;
219	>	double  mlimit[2];
220	>	int  ilimit = 4.0/TEPS;
221	>	int  i;
222		/* iterative search for m */
223	+	mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
224		do {
225		/* compute grey weighted average */
226	<	i = 1.25*CHECKRAD*rad*m + .5;
226	>	i = RSCA*CHECKRAD*rad*m + .5;
227		if (i > orad) i = orad;
228		avg = wsum = 0.0;
229		while (i--) {
#	Line 237 | Line 232 | double  p0;
232		avg += t*ringsum[i];
233		wsum += t*ringwt[i];
234		}
235	+	if (avg < 1e-20)                /* zero inclusive average */
236	+	return(1.0);
237		avg /= wsum;
238		/* check threshold */
239	<	t = p0 - avg;
240	<	if (t < 0.0) t = -t;
241	<	t *= gausstable[0]/(m*m*avg);
242	<	if (t <= thresh && (m <= 1.0+FTINY ||
243	<	(thresh-t)/thresh <= TEPS))
244	<	break;
245	<	if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
246	<	(t-thresh)/thresh <= TEPS))
247	<	break;
239	>	denom = m*m/gausstable[0] - p0/avg;
240	>	if (denom <= FTINY) {           /* zero exclusive average */
241	>	if (m >= mlimit[1]-REPS)
242	>	break;
243	>	m = mlimit[1];
244	>	continue;
245	>	}
246	>	num = p0/avg - 1.0;
247	>	if (num < 0.0) num = -num;
248	>	t = num/denom;
249	>	if (t <= thresh) {
250	>	if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
251	>	break;
252	>	} else {
253	>	if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
254	>	break;
255	>	}
256	>	t = m;                  /* remember current m */
257		/* next guesstimate */
258	<	m *= sqrt(t/thresh);
259	<	if (m < 1.0) m = 1.0;
260	<	else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
258	>	m = sqrt(gausstable[0]*(num/thresh + p0/avg));
259	>	if (m < t) {            /* bound it */
260	>	if (m <= mlimit[0]+FTINY)
261	>	m = 0.5*(mlimit[0] + t);
262	>	mlimit[1] = t;
263	>	} else {
264	>	if (m >= mlimit[1]-FTINY)
265	>	m = 0.5*(mlimit[1] + t);
266	>	mlimit[0] = t;
267	>	}
268		} while (--ilimit > 0);
269		return(m);
270		}
271
272
273	<	sumans(px, py, rcent, ccent, m)         /* sum input pixel to output */
274	<	int  px, py;
275	<	int  rcent, ccent;
276	<	double  m;
273	>	static void
274	>	sumans(         /* sum input pixel to output */
275	>	int  px,
276	>	int  py,
277	>	int  rcent,
278	>	int  ccent,
279	>	double  m
280	>	)
281		{
282	<	double  dy, dx;
283	<	COLOR  pval, ctmp;
284	<	int  ksiz, r;
282	>	double  dy2, dx;
283	>	SCOLOR  pval, ctmp;
284	>	int  ksiz, r, offs;
285		double  pc, pr, norm;
286	<	register int  i, c;
287	<	register COLOR  *scan;
288	<
289	<	copycolor(pval, scanin[py%barsize][px]);
286	>	int  i, c;
287	>	COLORV  *scan, *scp;
288	>	/*
289	>	* This normalization method fails at the picture borders because
290	>	* a different number of input pixels contribute there.
291	>	*/
292	>	scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0)*NCSAMP;
293	>	copyscolor(pval, scan+px*NCSAMP);
294		pc = x_c*(px+.5);
295		pr = y_r*(py+.5);
296		ksiz = CHECKRAD*m*rad + 1;
#	Line 280 | Line 301 | double  m;
301		for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
302		if (r < 0) continue;
303		if (r >= nrows) break;
304	<	dy = (pr - (r+.5))/(m*rad);
304	>	dy2 = (pr - (r+.5))/(m*rad);
305	>	dy2 *= dy2;
306		for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
307	<	if (c < 0) continue;
308	<	if (c >= ncols) break;
307	>	if (!wrapfilt) {
308	>	if (c < 0) continue;
309	>	if (c >= ncols) break;
310	>	}
311		dx = (pc - (c+.5))/(m*rad);
312	<	norm += warr[i++] = lookgauss(dx*dx + dy*dy);
312	>	norm += warr[i++] = lookgauss(dx*dx + dy2);
313		}
314		}
315		norm = 1.0/norm;
#	Line 298 | Line 322 | double  m;
322		if (r >= nrows) break;
323		scan = scoutbar[r%obarsize];
324		for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
325	<	if (c < 0) continue;
326	<	if (c >= ncols) break;
327	<	copycolor(ctmp, pval);
325	>	offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
326	>	if ((offs != 0) & !wrapfilt)
327	>	continue;
328	>	copyscolor(ctmp, pval);
329		dx = norm*warr[i++];
330	<	scalecolor(ctmp, dx);
331	<	addcolor(scan[c], ctmp);
330	>	scalescolor(ctmp, dx);
331	>	scp = scan + (c+offs)*NCSAMP;
332	>	saddscolor(scp, ctmp);
333		}
334		}
335		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)