ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhpict2.c
(Generate patch)

Comparing ray/src/hd/rhpict2.c (file contents):
Revision 3.1 by gwlarson, Thu Mar 4 10:30:04 1999 UTC vs.
Revision 3.10 by greg, Sat Feb 22 02:07:25 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1999 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Rendering routines for rhpict.
6   */
7  
8   #include "holo.h"
9   #include "view.h"
10 + #include "random.h"
11  
12 + #ifndef DEPS
13 + #define DEPS            0.02    /* depth epsilon */
14 + #endif
15 + #ifndef PEPS
16 + #define PEPS            0.04    /* pixel value epsilon */
17 + #endif
18 + #ifndef MAXRAD
19 + #define MAXRAD          64      /* maximum kernel radius */
20 + #endif
21 + #ifndef NNEIGH
22 + #define NNEIGH          5       /* find this many neighbors */
23 + #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 + static short    isqrttab[MAXRAD2];      /* integer square root table */
44 +
45 + #define isqrt(i2)       (isqrttab[i2])
46 +
47   extern VIEW     myview;         /* current output view */
48   extern COLOR    *mypixel;       /* pixels being rendered */
49   extern float    *myweight;      /* weights (used to compute final pixels) */
50 + extern float    *mydepth;       /* depth values (visibility culling) */
51   extern int      hres, vres;     /* current horizontal and vertical res. */
52  
53  
54 < render_beam(bp, hb)             /* render a particular beam */
54 > pixBeam(bp, hb)                 /* render a particular beam */
55   BEAM    *bp;
56   register HDBEAMI        *hb;
57   {
# Line 26 | Line 60 | register HDBEAMI       *hb;
60          FVECT   rorg, rdir, wp, ip;
61          double  d, prox;
62          COLOR   col;
63 <        int     n, p;
63 >        int     n;
64 >        register int4   p;
65  
66          if (!hdbcoord(gc, hb->h, hb->b))
67                  error(CONSISTENCY, "bad beam in render_beam");
# Line 38 | Line 73 | register HDBEAMI       *hb;
73                          VSUM(wp, rorg, rdir, d);
74                          VSUB(ip, wp, myview.vp);
75                          d = DOT(ip,rdir);
76 <                        prox = 1. - d*d/DOT(ip,ip);     /* sin(diff_angle)^2 */
76 >                        prox = d*d/DOT(ip,ip);  /* cos(diff_angle)^32 */
77 >                        prox *= prox; prox *= prox; prox *= prox; prox *= prox;
78                  } else {
79                          if (myview.type == VT_PAR || myview.vaft > FTINY)
80                                  continue;               /* inf. off view */
81                          VSUM(wp, myview.vp, rdir, FHUGE);
82 <                        prox = 0.;
82 >                        prox = 1.;
83                  }
84                  viewloc(ip, &myview, wp);       /* frustum clipping */
85                  if (ip[2] < 0.)
# Line 53 | Line 89 | register HDBEAMI       *hb;
89                  if (ip[1] < 0. || ip[1] >= 1.)
90                          continue;
91                  if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
92 <                        continue;
57 <                                                /* we're in! */
92 >                        continue;               /* not exact for VT_PER */
93                  p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
94 +                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                  colr_color(col, rv->v);
103 +                scalecolor(col, prox);
104                  addcolor(mypixel[p], col);
105 <                myweight[p] += 1.0;
105 >                myweight[p] += prox;
106 >                mydepth[p] = ip[2];
107          }
108 + }
109 +
110 +
111 + int
112 + kill_occl(h, v, nl, nd, n)              /* check for occlusion errors */
113 + int     h, v;
114 + register short  nl[NNEIGH][2];
115 + int     nd[NNEIGH];
116 + int     n;
117 + {
118 +        short   forequad[2][2];
119 +        int     d;
120 +        register int    i;
121 +        register int4   p;
122 +
123 +        if (n <= 0) {
124 + #ifdef DEBUG
125 +                error(WARNING, "neighborless sample in kill_occl");
126 + #endif
127 +                return(1);
128 +        }
129 +        p = v*hres + h;
130 +        forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
131 +        for (i = n; i--; ) {
132 +                d = isqrt(nd[i]);
133 +                if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
134 +                        forequad[nl[i][0]<h][nl[i][1]<v] = 1;
135 +        }
136 +        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 +        }
140 +        return(1);
141 + }
142 +
143 +
144 + int
145 + smooth_samp(h, v, nl, nd, n)    /* grow sample point smoothly */
146 + int     h, v;
147 + register short  nl[NNEIGH][2];
148 + int     nd[NNEIGH];
149 + int     n;
150 + {
151 +        int     dis[NNEIGH], ndis;
152 +        COLOR   mykern[MAXRAD2];
153 +        int     maxr2;
154 +        double  d;
155 +        register int4   p;
156 +        register int    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 = nd[n-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 +                r = isqrt(nd[i]);
175 +                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 +        }
181 +                                                /* stamp out that kernel */
182 +        for (v2 = v-maxr; v2 <= v+maxr; v2++) {
183 +                if (v2 < 0) v2 = 0;
184 +                else if (v2 >= vres) break;
185 +                for (h2 = h-maxr; h2 <= h+maxr; h2++) {
186 +                        if (h2 < 0) h2 = 0;
187 +                        else if (h2 >= hres) break;
188 +                        r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
189 +                        if (r2 > maxr2) continue;
190 +                        if (CHK4(pixFlags, v2*hres+h2))
191 +                                continue;       /* occupied */
192 +                        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 +                        addcolor(mypixel[v2*hres+h2], mykern[r2]);
199 +                        myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
200 +                }
201 +        }
202 +        return(1);
203 + }
204 +
205 +
206 + int
207 + random_samp(h, v, nl, nd, n, rf)        /* gather samples randomly */
208 + int     h, v;
209 + register short  nl[NNEIGH][2];
210 + int     nd[NNEIGH];
211 + int     n;
212 + double  *rf;
213 + {
214 +        float   rnt[NNEIGH];
215 +        double  rvar;
216 +        register int4   p, pn;
217 +        register int    ni;
218 +
219 +        if (n <= 0)
220 +                return(1);
221 +        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 +        }
233 +        pn = nl[ni][1]*hres + nl[ni][0];
234 +        addcolor(mypixel[p], mypixel[pn]);
235 +        myweight[p] += myweight[pn];
236 +        return(1);
237 + }
238 +
239 +
240 + pixFinish(ransamp)              /* done with beams -- compute pixel values */
241 + double  ransamp;
242 + {
243 +        if (pixWeight[0] <= FTINY)
244 +                init_wfunc();           /* initialize weighting function */
245 +        reset_flags();                  /* set occupancy flags */
246 +        meet_neighbors(1,kill_occl,NULL); /* identify occlusion errors */
247 +        reset_flags();                  /* reset occupancy flags */
248 +        if (ransamp >= 0.)              /* spread samples over image */
249 +                meet_neighbors(0,random_samp,&ransamp);
250 +        else
251 +                meet_neighbors(1,smooth_samp,NULL);
252 +        free((void *)pixFlags);         /* free pixel flags */
253 +        pixFlags = NULL;
254 + }
255 +
256 +
257 + reset_flags()                   /* allocate/set/reset occupancy flags */
258 + {
259 +        register int4   p;
260 +
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 + }
270 +
271 +
272 + init_wfunc()                    /* initialize weighting function */
273 + {
274 +        register int    r2;
275 +        register double d;
276 +
277 +        for (r2 = MAXRAD2; --r2; ) {
278 +                d = sqrt((double)r2);
279 +                pixWeight[r2] = G0NORM/d;
280 +                isqrttab[r2] = d + 0.99;
281 +        }
282 +        pixWeight[0] = 1.;
283 +        isqrttab[0] = 0;
284 + }
285 +
286 +
287 + int
288 + findneigh(nl, nd, h, v, rnl)    /* find NNEIGH neighbors for pixel */
289 + short   nl[NNEIGH][2];
290 + int     nd[NNEIGH];
291 + int     h, v;
292 + register short  (*rnl)[NNEIGH];
293 + {
294 +        int     nn = 0;
295 +        int     d, n, hoff;
296 +        register int    h2, n2;
297 +
298 +        nd[NNEIGH-1] = MAXRAD2;
299 +        for (hoff = 0; hoff < hres; hoff = (hoff<=0) - hoff) {
300 +                h2 = h + hoff;
301 +                if (h2 < 0 | h2 >= hres)
302 +                        continue;
303 +                if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
304 +                        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 +                        if (d == 0 | d >= nd[NNEIGH-1])
308 +                                continue;
309 +                        if (nn < NNEIGH)        /* insert neighbor */
310 +                                nn++;
311 +                        for (n2 = nn; n2--; )   {
312 +                                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 + meet_neighbors(occ, nf, dp)     /* run through samples and their neighbors */
329 + int     occ;
330 + int     (*nf)();
331 + char    *dp;
332 + {
333 +        short   ln[NNEIGH][2];
334 +        int     nd[NNEIGH];
335 +        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 +                        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 +                }
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 +        free((void *)rnl);              /* free row list */
374   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines