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.3 by gwlarson, Mon Mar 8 14:09:11 1999 UTC

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

Diff Legend

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