ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhpict2.c
Revision: 3.9
Committed: Wed Mar 10 16:15:19 1999 UTC (25 years ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.8: +52 -68 lines
Log Message:
reduced NNEIGH default from 7 to 5
changed random sampling function to use sampled Voronoi regions

File Contents

# User Rev Content
1 gwlarson 3.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 gwlarson 3.6 #include "random.h"
14 gwlarson 3.1
15 gwlarson 3.2 #ifndef DEPS
16 gwlarson 3.3 #define DEPS 0.02 /* depth epsilon */
17 gwlarson 3.2 #endif
18 gwlarson 3.3 #ifndef PEPS
19     #define PEPS 0.04 /* pixel value epsilon */
20     #endif
21 gwlarson 3.2 #ifndef MAXRAD
22     #define MAXRAD 64 /* maximum kernel radius */
23     #endif
24     #ifndef NNEIGH
25 gwlarson 3.9 #define NNEIGH 5 /* find this many neighbors */
26 gwlarson 3.2 #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 gwlarson 3.3 static short isqrttab[MAXRAD2]; /* integer square root table */
47 gwlarson 3.2
48 gwlarson 3.9 #define isqrt(i2) (isqrttab[i2])
49 gwlarson 3.3
50 gwlarson 3.1 extern VIEW myview; /* current output view */
51     extern COLOR *mypixel; /* pixels being rendered */
52     extern float *myweight; /* weights (used to compute final pixels) */
53 gwlarson 3.2 extern float *mydepth; /* depth values (visibility culling) */
54 gwlarson 3.1 extern int hres, vres; /* current horizontal and vertical res. */
55    
56    
57 gwlarson 3.2 pixBeam(bp, hb) /* render a particular beam */
58 gwlarson 3.1 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 gwlarson 3.2 int n;
67     register int4 p;
68 gwlarson 3.1
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 gwlarson 3.2 prox = d*d/DOT(ip,ip); /* cos(diff_angle)^32 */
80     prox *= prox; prox *= prox; prox *= prox; prox *= prox;
81 gwlarson 3.1 } else {
82     if (myview.type == VT_PAR || myview.vaft > FTINY)
83     continue; /* inf. off view */
84     VSUM(wp, myview.vp, rdir, FHUGE);
85 gwlarson 3.2 prox = 1.;
86 gwlarson 3.1 }
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 gwlarson 3.2 continue; /* not exact for VT_PER */
96 gwlarson 3.1 p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
97 gwlarson 3.2 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 gwlarson 3.1 colr_color(col, rv->v);
106 gwlarson 3.2 scalecolor(col, prox);
107 gwlarson 3.1 addcolor(mypixel[p], col);
108 gwlarson 3.2 myweight[p] += prox;
109     mydepth[p] = ip[2];
110 gwlarson 3.1 }
111 gwlarson 3.2 }
112    
113    
114     int
115 gwlarson 3.9 kill_occl(h, v, nl, nd, n) /* check for occlusion errors */
116 gwlarson 3.2 int h, v;
117 gwlarson 3.3 register short nl[NNEIGH][2];
118 gwlarson 3.9 int nd[NNEIGH];
119 gwlarson 3.2 int n;
120     {
121     short forequad[2][2];
122     int d;
123 gwlarson 3.9 register int i;
124     register int4 p;
125 gwlarson 3.2
126 gwlarson 3.5 if (n <= 0) {
127     #ifdef DEBUG
128     error(WARNING, "neighborless sample in kill_occl");
129     #endif
130 gwlarson 3.2 return(1);
131 gwlarson 3.5 }
132 gwlarson 3.3 p = v*hres + h;
133 gwlarson 3.2 forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
134     for (i = n; i--; ) {
135 gwlarson 3.9 d = isqrt(nd[i]);
136 gwlarson 3.3 if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
137 gwlarson 3.2 forequad[nl[i][0]<h][nl[i][1]<v] = 1;
138     }
139 gwlarson 3.3 if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
140     setcolor(mypixel[p], 0., 0., 0.);
141     myweight[p] = 0.; /* occupancy reset afterwards */
142 gwlarson 3.2 }
143     return(1);
144     }
145    
146    
147     int
148 gwlarson 3.9 smooth_samp(h, v, nl, nd, n) /* grow sample point smoothly */
149 gwlarson 3.2 int h, v;
150     register short nl[NNEIGH][2];
151 gwlarson 3.9 int nd[NNEIGH];
152 gwlarson 3.2 int n;
153     {
154 gwlarson 3.3 int dis[NNEIGH], ndis;
155 gwlarson 3.2 COLOR mykern[MAXRAD2];
156 gwlarson 3.9 int maxr2;
157 gwlarson 3.5 double d;
158 gwlarson 3.9 register int4 p;
159     register int r2;
160 gwlarson 3.3 int i, r, maxr, h2, v2;
161 gwlarson 3.2
162     if (n <= 0)
163     return(1);
164     p = v*hres + h; /* build kernel values */
165 gwlarson 3.9 maxr2 = nd[n-1];
166 gwlarson 3.2 DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
167 gwlarson 3.3 maxr = isqrt(maxr2);
168     for (v2 = 1; v2 <= maxr; v2++)
169     for (h2 = 0; h2 <= v2; h2++) {
170     r2 = h2*h2 + v2*v2;
171     if (r2 > maxr2) break;
172     copycolor(mykern[r2], mypixel[p]);
173     scalecolor(mykern[r2], pixWeight[r2]);
174     }
175     ndis = 0; /* find discontinuities */
176     for (i = n; i--; ) {
177 gwlarson 3.9 r = isqrt(nd[i]);
178 gwlarson 3.3 d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
179     d = d>=1. ? d-1. : 1.-d;
180     if (d > r*DEPS || bigdiff(mypixel[p],
181     mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
182     dis[ndis++] = i;
183 gwlarson 3.2 }
184 gwlarson 3.3 /* stamp out that kernel */
185 gwlarson 3.2 for (v2 = v-maxr; v2 <= v+maxr; v2++) {
186 gwlarson 3.3 if (v2 < 0) v2 = 0;
187     else if (v2 >= vres) break;
188 gwlarson 3.2 for (h2 = h-maxr; h2 <= h+maxr; h2++) {
189 gwlarson 3.3 if (h2 < 0) h2 = 0;
190     else if (h2 >= hres) break;
191     r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
192 gwlarson 3.2 if (r2 > maxr2) continue;
193     if (CHK4(pixFlags, v2*hres+h2))
194     continue; /* occupied */
195 gwlarson 3.3 for (i = ndis; i--; ) {
196     r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
197     (v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
198     if (r < r2) break;
199     }
200     if (i >= 0) continue; /* outside edge */
201 gwlarson 3.2 addcolor(mypixel[v2*hres+h2], mykern[r2]);
202 gwlarson 3.5 myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
203 gwlarson 3.2 }
204     }
205     return(1);
206     }
207    
208    
209 gwlarson 3.6 int
210 gwlarson 3.9 random_samp(h, v, nl, nd, n, rf) /* gather samples randomly */
211 gwlarson 3.6 int h, v;
212     register short nl[NNEIGH][2];
213 gwlarson 3.9 int nd[NNEIGH];
214 gwlarson 3.6 int n;
215 gwlarson 3.7 double *rf;
216 gwlarson 3.2 {
217 gwlarson 3.9 float rnt[NNEIGH];
218     double rvar;
219     register int4 p, pn;
220     register int ni;
221 gwlarson 3.6
222     if (n <= 0)
223     return(1);
224 gwlarson 3.9 p = v*hres + h;
225     if (*rf <= FTINY) /* straight Voronoi regions */
226     ni = 0;
227     else { /* weighted choice */
228     DCHECK(nd[n-1]>=MAXRAD2, CONSISTENCY, "out of range neighbor");
229     rnt[0] = pixWeight[nd[0]];
230     for (ni = 1; ni < n; ni++)
231     rnt[ni] = rnt[ni-1] + pixWeight[nd[ni]];
232     rvar = rnt[n-1]*pow(frandom(), 1. / *rf);
233     for (ni = 0; rvar > rnt[ni]+FTINY; ni++)
234     ;
235 gwlarson 3.6 }
236 gwlarson 3.9 pn = nl[ni][1]*hres + nl[ni][0];
237     addcolor(mypixel[p], mypixel[pn]);
238     myweight[p] += myweight[pn];
239 gwlarson 3.6 return(1);
240     }
241    
242    
243     pixFinish(ransamp) /* done with beams -- compute pixel values */
244 gwlarson 3.7 double ransamp;
245 gwlarson 3.6 {
246 gwlarson 3.5 if (pixWeight[0] <= FTINY)
247     init_wfunc(); /* initialize weighting function */
248 gwlarson 3.2 reset_flags(); /* set occupancy flags */
249 gwlarson 3.9 meet_neighbors(1,kill_occl,NULL); /* identify occlusion errors */
250 gwlarson 3.2 reset_flags(); /* reset occupancy flags */
251 gwlarson 3.7 if (ransamp >= 0.) /* spread samples over image */
252 gwlarson 3.9 meet_neighbors(0,random_samp,&ransamp);
253 gwlarson 3.6 else
254 gwlarson 3.9 meet_neighbors(1,smooth_samp,NULL);
255 gwlarson 3.2 free((char *)pixFlags); /* free pixel flags */
256     pixFlags = NULL;
257     }
258    
259    
260     reset_flags() /* allocate/set/reset occupancy flags */
261     {
262 gwlarson 3.4 register int4 p;
263 gwlarson 3.2
264     if (pixFlags == NULL) {
265     pixFlags = (int4 *)calloc(FL4NELS(hres*vres), sizeof(int4));
266     CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
267     } else
268     CLR4ALL(pixFlags, hres*vres);
269     for (p = hres*vres; p--; )
270     if (myweight[p] > FTINY)
271     SET4(pixFlags, p);
272 gwlarson 3.4 }
273    
274    
275     init_wfunc() /* initialize weighting function */
276     {
277 gwlarson 3.9 register int r2;
278 gwlarson 3.4 register double d;
279    
280 gwlarson 3.6 for (r2 = MAXRAD2; --r2; ) {
281     d = sqrt((double)r2);
282     pixWeight[r2] = G0NORM/d;
283     isqrttab[r2] = d + 0.99;
284     }
285 gwlarson 3.4 pixWeight[0] = 1.;
286     isqrttab[0] = 0;
287 gwlarson 3.2 }
288    
289    
290     int
291 gwlarson 3.9 findneigh(nl, nd, h, v, rnl) /* find NNEIGH neighbors for pixel */
292 gwlarson 3.2 short nl[NNEIGH][2];
293 gwlarson 3.9 int nd[NNEIGH];
294 gwlarson 3.2 int h, v;
295     register short (*rnl)[NNEIGH];
296     {
297     int nn = 0;
298 gwlarson 3.9 int d, n, hoff;
299 gwlarson 3.2 register int h2, n2;
300    
301 gwlarson 3.3 nd[NNEIGH-1] = MAXRAD2;
302 gwlarson 3.9 for (hoff = 0; hoff < hres; hoff = (hoff<=0) - hoff) {
303 gwlarson 3.2 h2 = h + hoff;
304     if (h2 < 0 | h2 >= hres)
305     continue;
306 gwlarson 3.3 if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
307 gwlarson 3.2 break;
308     for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
309     d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
310 gwlarson 3.9 if (d == 0 | d >= nd[NNEIGH-1])
311 gwlarson 3.2 continue;
312 gwlarson 3.3 if (nn < NNEIGH) /* insert neighbor */
313     nn++;
314     for (n2 = nn; n2--; ) {
315 gwlarson 3.2 if (!n2 || d >= nd[n2-1]) {
316     nd[n2] = d;
317     nl[n2][0] = h2;
318     nl[n2][1] = rnl[h2][n];
319     break;
320     }
321     nd[n2] = nd[n2-1];
322     nl[n2][0] = nl[n2-1][0];
323     nl[n2][1] = nl[n2-1][1];
324     }
325     }
326     }
327     return(nn);
328     }
329    
330    
331 gwlarson 3.9 meet_neighbors(occ, nf, dp) /* run through samples and their neighbors */
332     int occ;
333 gwlarson 3.2 int (*nf)();
334 gwlarson 3.7 char *dp;
335 gwlarson 3.2 {
336     short ln[NNEIGH][2];
337 gwlarson 3.9 int nd[NNEIGH];
338 gwlarson 3.2 int h, v, n, v2;
339     register short (*rnl)[NNEIGH];
340     /* initialize bottom row list */
341     rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
342     CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
343     for (h = 0; h < hres; h++) {
344     for (n = v = 0; v < vres; v++)
345     if (CHK4(pixFlags, v*hres+h)) {
346     rnl[h][n++] = v;
347     if (n >= NNEIGH)
348     break;
349     }
350     while (n < NNEIGH)
351     rnl[h][n++] = NINF;
352     }
353     v = 0; /* do each row */
354     for ( ; ; ) {
355     for (h = 0; h < hres; h++) {
356 gwlarson 3.9 if (!CHK4(pixFlags, v*hres+h) != !occ)
357     continue; /* occupancy mismatch */
358     /* find neighbors */
359     n = findneigh(ln, nd, h, v, rnl);
360     /* call on neighbors */
361     (*nf)(h, v, ln, nd, n, dp);
362 gwlarson 3.2 }
363     if (++v >= vres) /* reinitialize row list */
364     break;
365     for (h = 0; h < hres; h++)
366     for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
367     if (v2 - v > v - rnl[h][0])
368     break; /* not close enough */
369     if (CHK4(pixFlags, v2*hres+h)) {
370     for (n = 0; n < NNEIGH-1; n++)
371     rnl[h][n] = rnl[h][n+1];
372     rnl[h][NNEIGH-1] = v2;
373     }
374     }
375     }
376     free((char *)rnl); /* free row list */
377 gwlarson 3.1 }