ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhpict2.c
Revision: 3.10
Committed: Sat Feb 22 02:07:25 2003 UTC (21 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 3.9: +3 -6 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

File Contents

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