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.6 by greg, Wed Jun 23 12:17:09 1993 UTC vs.
Revision 2.7 by greg, Fri Jun 25 15:01:20 1993 UTC

#	Line 10 | Line 10 | static char SCCSid[] = "$SunId$ LBL";
10		*     8/13/86
11		*/
12
13	<	#include  <stdio.h>
13	>	#include  "standard.h"
14
15	–	#include  <math.h>
16	–
15		#include  "color.h"
16
17	<	#define  FTINY          1e-7
17	>	#define  TEPS           0.25    /* threshold proximity goal */
18
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
26	–	extern int  boxfilt;            /* do box filtering? */
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 window size */
33	<	extern int  yrad;               /* y window size */
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
47		float  *gausstable;             /* gauss lookup table */
48
49	+	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		#define  lookgauss(x)           gausstable[(int)(10.*(x)+.5)]
57
58
59		initmask()                      /* initialize gaussian lookup table */
60		{
61	<	extern char  *malloc();
61	>	int  gtabsiz;
62	>	double  gaussN;
63		double  d;
64		register int  x;
65
66	<	gausstable = (float *)malloc(100*sizeof(float));
67	<	if (gausstable == NULL) {
68	<	fprintf(stderr, "%s: out of memory in initmask\n", progname);
69	<	quit(1);
70	<	}
57	<	d = x_c*y_r*0.25/rad/rad;
66	>	gtabsiz = 150*CHECKRAD;
67	>	gausstable = (float *)malloc(gtabsiz*sizeof(float));
68	>	if (gausstable == NULL)
69	>	goto memerr;
70	>	d = x_c*y_r*0.25/(rad*rad);
71		gausstable[0] = exp(-d)/sqrt(d);
72	<	for (x = 1; x < 100; x++)
72	>	for (x = 1; x < gtabsiz; x++)
73		if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
74		gausstable[x] = gausstable[0];
75	+	if (obarsize == 0)
76	+	return;
77	+	/* compute integral of filter */
78	+	gaussN = PI*sqrt(d)*exp(-d);            /* plateau */
79	+	for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1)
80	+	gaussN += 0.1*2.0*PI*exp(-d*d);
81	+	/* normalize filter */
82	+	gaussN = x_c*y_r/(rad*rad*gaussN);
83	+	for (x = 0; x < gtabsiz; x++)
84	+	gausstable[x] *= gaussN;
85	+	/* create ring averages table */
86	+	ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
87	+	if (ringndx == NULL)
88	+	goto memerr;
89	+	for (x = 2*orad*orad+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(obarsize*obarsize*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		}
104
105
#	Line 72 | Line 113 | int  c, r;
113		int  y;
114		register int  x;
115		register COLOR  *scan;
116	<
116	>
117		wsum = 0;
118		setcolor(csum, 0.0, 0.0, 0.0);
119	<	for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
119	>	for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
120		if (y < 0) continue;
121		if (y >= yres) break;
122		d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
123	<	if (d > 0.5+FTINY || d < -0.5-FTINY)
124	<	continue;
123	>	if (d < -0.5) continue;
124	>	if (d >= 0.5) break;
125		scan = scanin[y%barsize];
126	<	for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
126	>	for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
127		if (x < 0) continue;
128		if (x >= xres) break;
129		d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
130	<	if (d > 0.5+FTINY || d < -0.5-FTINY)
131	<	continue;
130	>	if (d < -0.5) continue;
131	>	if (d >= 0.5) break;
132		wsum++;
133		addcolor(csum, scan[x]);
134		}
#	Line 127 | Line 168 | int  c, r;
168		}
169		}
170		scalecolor(csum, 1.0/wsum);
171	+	}
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[c*c + r*r];
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.25*CHECKRAD*rad*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]/(m*m*avg);
245	+	if (t <= thresh && (m <= 1.0+FTINY ||
246	+	(thresh-t)/thresh <= TEPS))
247	+	break;
248	+	if (t > thresh && (m*CHECKRAD*rad >= 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 (m*CHECKRAD*rad > 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 = CHECKRAD*m*rad + 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(dx*dx + dy*dy);
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		}

Diff Legend

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