ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhpict2.c
Revision: 3.8
Committed: Tue Mar 9 15:08:22 1999 UTC (25 years ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.7: +12 -11 lines
Log Message:
final cut on random resampling algorithm -- fairer choice algorithm

File Contents

# Content
1 /* Copyright (c) 1999 Silicon Graphics, Inc. */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ SGI";
5 #endif
6
7 /*
8 * Rendering routines for rhpict.
9 */
10
11 #include "holo.h"
12 #include "view.h"
13 #include "random.h"
14
15 #ifndef DEPS
16 #define DEPS 0.02 /* depth epsilon */
17 #endif
18 #ifndef PEPS
19 #define PEPS 0.04 /* pixel value epsilon */
20 #endif
21 #ifndef MAXRAD
22 #define MAXRAD 64 /* maximum kernel radius */
23 #endif
24 #ifndef NNEIGH
25 #define NNEIGH 7 /* find this many neighbors */
26 #endif
27
28 #define NINF 16382
29
30 #define MAXRAD2 (MAXRAD*MAXRAD+1)
31
32 #define G0NORM 0.286 /* ground zero normalization (1/x integral) */
33
34 #ifndef FL4OP
35 #define FL4OP(f,i,op) ((f)[(i)>>5] op (1L<<((i)&0x1f)))
36 #define CHK4(f,i) FL4OP(f,i,&)
37 #define SET4(f,i) FL4OP(f,i,|=)
38 #define CLR4(f,i) FL4OP(f,i,&=~)
39 #define TGL4(f,i) FL4OP(f,i,^=)
40 #define FL4NELS(n) (((n)+0x1f)>>5)
41 #define CLR4ALL(f,n) bzero((char *)(f),FL4NELS(n)*sizeof(int4))
42 #endif
43
44 static int4 *pixFlags; /* pixel occupancy flags */
45 static float pixWeight[MAXRAD2]; /* pixel weighting function */
46 static short isqrttab[MAXRAD2]; /* integer square root table */
47
48 #define isqrt(i2) ((int)isqrttab[(int)(i2)])
49
50 extern VIEW myview; /* current output view */
51 extern COLOR *mypixel; /* pixels being rendered */
52 extern float *myweight; /* weights (used to compute final pixels) */
53 extern float *mydepth; /* depth values (visibility culling) */
54 extern int hres, vres; /* current horizontal and vertical res. */
55
56
57 pixBeam(bp, hb) /* render a particular beam */
58 BEAM *bp;
59 register HDBEAMI *hb;
60 {
61 GCOORD gc[2];
62 register RAYVAL *rv;
63 FVECT rorg, rdir, wp, ip;
64 double d, prox;
65 COLOR col;
66 int n;
67 register int4 p;
68
69 if (!hdbcoord(gc, hb->h, hb->b))
70 error(CONSISTENCY, "bad beam in render_beam");
71 for (n = bp->nrm, rv = hdbray(bp); n--; rv++) {
72 /* reproject each sample */
73 hdray(rorg, rdir, hb->h, gc, rv->r);
74 if (rv->d < DCINF) {
75 d = hddepth(hb->h, rv->d);
76 VSUM(wp, rorg, rdir, d);
77 VSUB(ip, wp, myview.vp);
78 d = DOT(ip,rdir);
79 prox = d*d/DOT(ip,ip); /* cos(diff_angle)^32 */
80 prox *= prox; prox *= prox; prox *= prox; prox *= prox;
81 } else {
82 if (myview.type == VT_PAR || myview.vaft > FTINY)
83 continue; /* inf. off view */
84 VSUM(wp, myview.vp, rdir, FHUGE);
85 prox = 1.;
86 }
87 viewloc(ip, &myview, wp); /* frustum clipping */
88 if (ip[2] < 0.)
89 continue;
90 if (ip[0] < 0. || ip[0] >= 1.)
91 continue;
92 if (ip[1] < 0. || ip[1] >= 1.)
93 continue;
94 if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
95 continue; /* not exact for VT_PER */
96 p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
97 if (mydepth[p] > FTINY) { /* check depth */
98 if (ip[2] > mydepth[p]*(1.+DEPS))
99 continue;
100 if (ip[2] < mydepth[p]*(1.-DEPS)) {
101 setcolor(mypixel[p], 0., 0., 0.);
102 myweight[p] = 0.;
103 }
104 }
105 colr_color(col, rv->v);
106 scalecolor(col, prox);
107 addcolor(mypixel[p], col);
108 myweight[p] += prox;
109 mydepth[p] = ip[2];
110 }
111 }
112
113
114 int
115 kill_occl(h, v, nl, n) /* check for occlusion errors */
116 int h, v;
117 register short nl[NNEIGH][2];
118 int n;
119 {
120 short forequad[2][2];
121 int d;
122 register int4 i, p;
123
124 if (n <= 0) {
125 #ifdef DEBUG
126 error(WARNING, "neighborless sample in kill_occl");
127 #endif
128 return(1);
129 }
130 p = v*hres + h;
131 forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
132 for (i = n; i--; ) {
133 d = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
134 d = isqrt(d);
135 if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
136 forequad[nl[i][0]<h][nl[i][1]<v] = 1;
137 }
138 if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
139 setcolor(mypixel[p], 0., 0., 0.);
140 myweight[p] = 0.; /* occupancy reset afterwards */
141 }
142 return(1);
143 }
144
145
146 int
147 smooth_samp(h, v, nl, n) /* grow sample point smoothly */
148 int h, v;
149 register short nl[NNEIGH][2];
150 int n;
151 {
152 int dis[NNEIGH], ndis;
153 COLOR mykern[MAXRAD2];
154 int4 maxr2;
155 double d;
156 register int4 p, r2;
157 int i, r, maxr, h2, v2;
158
159 if (n <= 0)
160 return(1);
161 p = v*hres + h; /* build kernel values */
162 maxr2 = (h-nl[n-1][0])*(h-nl[n-1][0]) + (v-nl[n-1][1])*(v-nl[n-1][1]);
163 DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
164 maxr = isqrt(maxr2);
165 for (v2 = 1; v2 <= maxr; v2++)
166 for (h2 = 0; h2 <= v2; h2++) {
167 r2 = h2*h2 + v2*v2;
168 if (r2 > maxr2) break;
169 copycolor(mykern[r2], mypixel[p]);
170 scalecolor(mykern[r2], pixWeight[r2]);
171 }
172 ndis = 0; /* find discontinuities */
173 for (i = n; i--; ) {
174 r2 = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
175 r = isqrt(r2);
176 d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
177 d = d>=1. ? d-1. : 1.-d;
178 if (d > r*DEPS || bigdiff(mypixel[p],
179 mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
180 dis[ndis++] = i;
181 }
182 /* stamp out that kernel */
183 for (v2 = v-maxr; v2 <= v+maxr; v2++) {
184 if (v2 < 0) v2 = 0;
185 else if (v2 >= vres) break;
186 for (h2 = h-maxr; h2 <= h+maxr; h2++) {
187 if (h2 < 0) h2 = 0;
188 else if (h2 >= hres) break;
189 r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
190 if (r2 > maxr2) continue;
191 if (CHK4(pixFlags, v2*hres+h2))
192 continue; /* occupied */
193 for (i = ndis; i--; ) {
194 r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
195 (v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
196 if (r < r2) break;
197 }
198 if (i >= 0) continue; /* outside edge */
199 addcolor(mypixel[v2*hres+h2], mykern[r2]);
200 myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
201 }
202 }
203 return(1);
204 }
205
206
207 int
208 random_samp(h, v, nl, n, rf) /* grow sample point noisily */
209 int h, v;
210 register short nl[NNEIGH][2];
211 int n;
212 double *rf;
213 {
214 double tv = *rf * (1.-1./NNEIGH);
215 int4 maxr2;
216 register int4 p, r2;
217 register int i;
218 int maxr, h2, v2;
219 COLOR ctmp;
220
221 if (n <= 0)
222 return(1);
223 p = v*hres + h; /* compute kernel radius */
224 maxr2 = (h-nl[n-1][0])*(h-nl[n-1][0]) + (v-nl[n-1][1])*(v-nl[n-1][1]);
225 DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
226 maxr = isqrt(maxr2);
227 /* sample kernel */
228 for (v2 = v-maxr; v2 <= v+maxr; v2++) {
229 if (v2 < 0) v2 = 0;
230 else if (v2 >= vres) break;
231 for (h2 = h-maxr; h2 <= h+maxr; h2++) {
232 if (h2 < 0) h2 = 0;
233 else if (h2 >= hres) break;
234 r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
235 if (r2 > maxr2) continue;
236 if (CHK4(pixFlags, v2*hres+h2))
237 continue; /* occupied */
238 if (frandom() < tv) { /* use neighbor instead? */
239 i = random() % n;
240 r2 = nl[i][1]*hres + nl[i][0];
241 copycolor(ctmp, mypixel[r2]);
242 r2 = (h2-nl[i][0])*(h2-nl[i][0]) +
243 (v2-nl[i][1])*(v2-nl[i][1]);
244 if (r2 < MAXRAD2) {
245 scalecolor(ctmp, pixWeight[r2]);
246 addcolor(mypixel[v2*hres+h2], ctmp);
247 myweight[v2*hres+h2] += pixWeight[r2] *
248 myweight[nl[i][1]*hres+nl[i][0]];
249 continue;
250 }
251 }
252 /* use central sample */
253 copycolor(ctmp, mypixel[p]);
254 scalecolor(ctmp, pixWeight[r2]);
255 addcolor(mypixel[v2*hres+h2], ctmp);
256 myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
257 }
258 }
259 return(1);
260 }
261
262
263 pixFinish(ransamp) /* done with beams -- compute pixel values */
264 double ransamp;
265 {
266 if (pixWeight[0] <= FTINY)
267 init_wfunc(); /* initialize weighting function */
268 reset_flags(); /* set occupancy flags */
269 meet_neighbors(kill_occl,NULL); /* eliminate occlusion errors */
270 reset_flags(); /* reset occupancy flags */
271 if (ransamp >= 0.) /* spread samples over image */
272 meet_neighbors(random_samp,&ransamp);
273 else
274 meet_neighbors(smooth_samp,NULL);
275 free((char *)pixFlags); /* free pixel flags */
276 pixFlags = NULL;
277 }
278
279
280 reset_flags() /* allocate/set/reset occupancy flags */
281 {
282 register int4 p;
283
284 if (pixFlags == NULL) {
285 pixFlags = (int4 *)calloc(FL4NELS(hres*vres), sizeof(int4));
286 CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
287 } else
288 CLR4ALL(pixFlags, hres*vres);
289 for (p = hres*vres; p--; )
290 if (myweight[p] > FTINY)
291 SET4(pixFlags, p);
292 }
293
294
295 init_wfunc() /* initialize weighting function */
296 {
297 register int4 r2;
298 register double d;
299
300 for (r2 = MAXRAD2; --r2; ) {
301 d = sqrt((double)r2);
302 pixWeight[r2] = G0NORM/d;
303 isqrttab[r2] = d + 0.99;
304 }
305 pixWeight[0] = 1.;
306 isqrttab[0] = 0;
307 }
308
309
310 int
311 findneigh(nl, h, v, rnl) /* find NNEIGH neighbors for pixel */
312 short nl[NNEIGH][2];
313 int h, v;
314 register short (*rnl)[NNEIGH];
315 {
316 int nn = 0;
317 int4 d, nd[NNEIGH];
318 int n, hoff;
319 register int h2, n2;
320
321 nd[NNEIGH-1] = MAXRAD2;
322 for (hoff = 1; hoff < hres; hoff = (hoff<0) - hoff) {
323 h2 = h + hoff;
324 if (h2 < 0 | h2 >= hres)
325 continue;
326 if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
327 break;
328 for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
329 d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
330 if (d >= nd[NNEIGH-1])
331 continue;
332 if (nn < NNEIGH) /* insert neighbor */
333 nn++;
334 for (n2 = nn; n2--; ) {
335 if (!n2 || d >= nd[n2-1]) {
336 nd[n2] = d;
337 nl[n2][0] = h2;
338 nl[n2][1] = rnl[h2][n];
339 break;
340 }
341 nd[n2] = nd[n2-1];
342 nl[n2][0] = nl[n2-1][0];
343 nl[n2][1] = nl[n2-1][1];
344 }
345 }
346 }
347 return(nn);
348 }
349
350
351 meet_neighbors(nf, dp) /* run through samples and their neighbors */
352 int (*nf)();
353 char *dp;
354 {
355 short ln[NNEIGH][2];
356 int h, v, n, v2;
357 register short (*rnl)[NNEIGH];
358 /* initialize bottom row list */
359 rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
360 CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
361 for (h = 0; h < hres; h++) {
362 for (n = v = 0; v < vres; v++)
363 if (CHK4(pixFlags, v*hres+h)) {
364 rnl[h][n++] = v;
365 if (n >= NNEIGH)
366 break;
367 }
368 while (n < NNEIGH)
369 rnl[h][n++] = NINF;
370 }
371 v = 0; /* do each row */
372 for ( ; ; ) {
373 for (h = 0; h < hres; h++) {
374 if (!CHK4(pixFlags, v*hres+h))
375 continue; /* no one home */
376 n = findneigh(ln, h, v, rnl);
377 (*nf)(h, v, ln, n, dp); /* call on neighbors */
378 }
379 if (++v >= vres) /* reinitialize row list */
380 break;
381 for (h = 0; h < hres; h++)
382 for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
383 if (v2 - v > v - rnl[h][0])
384 break; /* not close enough */
385 if (CHK4(pixFlags, v2*hres+h)) {
386 for (n = 0; n < NNEIGH-1; n++)
387 rnl[h][n] = rnl[h][n+1];
388 rnl[h][NNEIGH-1] = v2;
389 }
390 }
391 }
392 free((char *)rnl); /* free row list */
393 }