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.18 by greg, Thu May 14 20:58:03 2020 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 <string.h>
9 +
10   #include "holo.h"
11   #include "view.h"
12 + #include "random.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          5       /* 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)    memset((char *)(f),'\0',FL4NELS(n)*sizeof(int32))
41 + #endif
42 +
43 + static int32    *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)       (isqrttab[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 + extern void pixFinish(double ransamp);
56 + extern void pixBeam(BEAM *bp, HDBEAMI *hb);
57  
58 < render_beam(bp, hb)             /* render a particular beam */
59 < BEAM    *bp;
60 < register HDBEAMI        *hb;
58 > typedef int (sampfunc)(int h, int v, short nl[NNEIGH][2], int nd[NNEIGH],
59 >                int n, double *rf);
60 > static sampfunc kill_occl;
61 > static sampfunc smooth_samp;
62 > static sampfunc random_samp;
63 > static void meet_neighbors(int occ, sampfunc *nf, double *dp);
64 > static void reset_flags(void);
65 > static void init_wfunc(void);
66 > static int findneigh(short nl[NNEIGH][2], int nd[NNEIGH], int h, int v,
67 >                short (*rnl)[NNEIGH]);
68 >
69 >
70 > void
71 > pixBeam(                        /* render a particular beam */
72 >        BEAM    *bp,
73 >        HDBEAMI *hb
74 > )
75   {
76          GCOORD  gc[2];
77 <        register RAYVAL *rv;
77 >        RAYVAL  *rv;
78          FVECT   rorg, rdir, wp, ip;
79          double  d, prox;
80          COLOR   col;
81 <        int     n, p;
81 >        int     n;
82 >        int32   p;
83  
84          if (!hdbcoord(gc, hb->h, hb->b))
85                  error(CONSISTENCY, "bad beam in render_beam");
# Line 38 | Line 91 | register HDBEAMI       *hb;
91                          VSUM(wp, rorg, rdir, d);
92                          VSUB(ip, wp, myview.vp);
93                          d = DOT(ip,rdir);
94 <                        prox = 1. - d*d/DOT(ip,ip);     /* sin(diff_angle)^2 */
94 >                        prox = d*d/DOT(ip,ip);  /* cos(diff_angle)^32 */
95 >                        prox *= prox; prox *= prox; prox *= prox; prox *= prox;
96                  } else {
97                          if (myview.type == VT_PAR || myview.vaft > FTINY)
98                                  continue;               /* inf. off view */
99                          VSUM(wp, myview.vp, rdir, FHUGE);
100 <                        prox = 0.;
100 >                        prox = 1.;
101                  }
102 <                viewloc(ip, &myview, wp);       /* frustum clipping */
103 <                if (ip[2] < 0.)
50 <                        continue;
51 <                if (ip[0] < 0. || ip[0] >= 1.)
52 <                        continue;
53 <                if (ip[1] < 0. || ip[1] >= 1.)
54 <                        continue;
102 >                if (viewloc(ip, &myview, wp) != VL_GOOD)
103 >                        continue;               /* frustum clipping */
104                  if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
105 <                        continue;
57 <                                                /* we're in! */
105 >                        continue;               /* not exact for VT_PER */
106                  p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
107 +                if (mydepth[p] > FTINY) {       /* check depth */
108 +                        if (ip[2] > mydepth[p]*(1.+DEPS))
109 +                                continue;
110 +                        if (ip[2] < mydepth[p]*(1.-DEPS)) {
111 +                                setcolor(mypixel[p], 0., 0., 0.);
112 +                                myweight[p] = 0.;
113 +                        }
114 +                }
115                  colr_color(col, rv->v);
116 +                scalecolor(col, prox);
117                  addcolor(mypixel[p], col);
118 <                myweight[p] += 1.0;
118 >                myweight[p] += prox;
119 >                mydepth[p] = ip[2];
120          }
121 + }
122 +
123 +
124 + static int
125 + kill_occl(              /* check for occlusion errors */
126 +        int     h,
127 +        int     v,
128 +        short   nl[NNEIGH][2],
129 +        int     nd[NNEIGH],
130 +        int     n,
131 +        double *rf
132 + )
133 + {
134 +        short   forequad[2][2];
135 +        int     d;
136 +        int     i;
137 +        int32   p;
138 +
139 +        if (n <= 0) {
140 + #ifdef DEBUG
141 +                error(WARNING, "neighborless sample in kill_occl");
142 + #endif
143 +                return(1);
144 +        }
145 +        p = v*hres + h;
146 +        forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
147 +        for (i = n; i--; ) {
148 +                d = isqrt(nd[i]);
149 +                if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
150 +                        forequad[nl[i][0]<h][nl[i][1]<v] = 1;
151 +        }
152 +        if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
153 +                setcolor(mypixel[p], 0., 0., 0.);
154 +                myweight[p] = 0.;       /* occupancy reset afterwards */
155 +        }
156 +        return(1);
157 + }
158 +
159 +
160 + static int
161 + smooth_samp(    /* grow sample point smoothly */
162 +        int     h,
163 +        int     v,
164 +        short   nl[NNEIGH][2],
165 +        int     nd[NNEIGH],
166 +        int     n,
167 +        double *rf
168 + )
169 + {
170 +        int     dis[NNEIGH], ndis;
171 +        COLOR   mykern[MAXRAD2];
172 +        int     maxr2;
173 +        double  d;
174 +        int32   p;
175 +        int     r2;
176 +        int     i, r, maxr, h2, v2;
177 +
178 +        if (n <= 0)
179 +                return(1);
180 +        p = v*hres + h;                         /* build kernel values */
181 +        maxr2 = nd[n-1];
182 +        DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
183 +        maxr = isqrt(maxr2);
184 +        for (v2 = 1; v2 <= maxr; v2++)
185 +                for (h2 = 0; h2 <= v2; h2++) {
186 +                        r2 = h2*h2 + v2*v2;
187 +                        if (r2 > maxr2) break;
188 +                        copycolor(mykern[r2], mypixel[p]);
189 +                        scalecolor(mykern[r2], pixWeight[r2]);
190 +                }
191 +        ndis = 0;                               /* find discontinuities */
192 +        for (i = n; i--; ) {
193 +                r = isqrt(nd[i]);
194 +                d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
195 +                d = d>=1. ? d-1. : 1.-d;
196 +                if (d > r*DEPS || bigdiff(mypixel[p],
197 +                                mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
198 +                        dis[ndis++] = i;
199 +        }
200 +                                                /* stamp out that kernel */
201 +        for (v2 = v-maxr; v2 <= v+maxr; v2++) {
202 +                if (v2 < 0) v2 = 0;
203 +                else if (v2 >= vres) break;
204 +                for (h2 = h-maxr; h2 <= h+maxr; h2++) {
205 +                        if (h2 < 0) h2 = 0;
206 +                        else if (h2 >= hres) break;
207 +                        r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
208 +                        if (r2 > maxr2) continue;
209 +                        if (CHK4(pixFlags, v2*hres+h2))
210 +                                continue;       /* occupied */
211 +                        for (i = ndis; i--; ) {
212 +                                r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
213 +                                        (v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
214 +                                if (r < r2) break;
215 +                        }
216 +                        if (i >= 0) continue;   /* outside edge */
217 +                        addcolor(mypixel[v2*hres+h2], mykern[r2]);
218 +                        myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
219 +                }
220 +        }
221 +        return(1);
222 + }
223 +
224 +
225 + static int
226 + random_samp(    /* gather samples randomly */
227 +        int     h,
228 +        int     v,
229 +        short   nl[NNEIGH][2],
230 +        int     nd[NNEIGH],
231 +        int     n,
232 +        double  *rf
233 + )
234 + {
235 +        float   rnt[NNEIGH];
236 +        double  rvar;
237 +        int32   p, pn;
238 +        int     ni;
239 +
240 +        if (n <= 0)
241 +                return(1);
242 +        p = v*hres + h;
243 +        if (*rf <= FTINY)               /* straight Voronoi regions */
244 +                ni = 0;
245 +        else {                          /* weighted choice */
246 +                DCHECK(nd[n-1]>=MAXRAD2, CONSISTENCY, "out of range neighbor");
247 +                rnt[0] = pixWeight[nd[0]];
248 +                for (ni = 1; ni < n; ni++)
249 +                        rnt[ni] = rnt[ni-1] + pixWeight[nd[ni]];
250 +                rvar = rnt[n-1]*pow(frandom(), 1. / *rf);
251 +                for (ni = 0; rvar > rnt[ni]+FTINY; ni++)
252 +                        ;
253 +        }
254 +        pn = nl[ni][1]*hres + nl[ni][0];
255 +        addcolor(mypixel[p], mypixel[pn]);
256 +        myweight[p] += myweight[pn];
257 +        return(1);
258 + }
259 +
260 +
261 + void
262 + pixFinish(              /* done with beams -- compute pixel values */
263 +        double  ransamp
264 + )
265 + {
266 +        if (pixWeight[0] <= FTINY)
267 +                init_wfunc();           /* initialize weighting function */
268 +        reset_flags();                  /* set occupancy flags */
269 +        meet_neighbors(1,kill_occl,NULL); /* identify occlusion errors */
270 +        reset_flags();                  /* reset occupancy flags */
271 +        if (ransamp >= 0.)              /* spread samples over image */
272 +                meet_neighbors(0,random_samp,&ransamp);
273 +        else
274 +                meet_neighbors(1,smooth_samp,NULL);
275 +        free((void *)pixFlags);         /* free pixel flags */
276 +        pixFlags = NULL;
277 + }
278 +
279 +
280 + static void
281 + reset_flags(void)                       /* allocate/set/reset occupancy flags */
282 + {
283 +        int32   p;
284 +
285 +        if (pixFlags == NULL) {
286 +                pixFlags = (int32 *)calloc(FL4NELS(hres*vres), sizeof(int32));
287 +                CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
288 +        } else
289 +                CLR4ALL(pixFlags, hres*vres);
290 +        for (p = hres*vres; p--; )
291 +                if (myweight[p] > FTINY)
292 +                        SET4(pixFlags, p);
293 + }
294 +
295 +
296 + static void
297 + init_wfunc(void)                        /* initialize weighting function */
298 + {
299 +        int     r2;
300 +        double  d;
301 +
302 +        for (r2 = MAXRAD2; --r2; ) {
303 +                d = sqrt((double)r2);
304 +                pixWeight[r2] = G0NORM/d;
305 +                isqrttab[r2] = d + 0.99;
306 +        }
307 +        pixWeight[0] = 1.;
308 +        isqrttab[0] = 0;
309 + }
310 +
311 +
312 + static int
313 + findneigh(      /* find NNEIGH neighbors for pixel */
314 +        short   nl[NNEIGH][2],
315 +        int     nd[NNEIGH],
316 +        int     h,
317 +        int     v,
318 +        short   (*rnl)[NNEIGH]
319 + )
320 + {
321 +        int     nn = 0;
322 +        int     d, n, hoff;
323 +        int     h2, n2;
324 +
325 +        nd[NNEIGH-1] = MAXRAD2;
326 +        for (hoff = 0; hoff < hres; hoff = (hoff<=0) - hoff) {
327 +                h2 = h + hoff;
328 +                if ((h2 < 0) | (h2 >= hres))
329 +                        continue;
330 +                if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
331 +                        break;
332 +                for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
333 +                        d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
334 +                        if ((d == 0) | (d >= nd[NNEIGH-1]))
335 +                                continue;
336 +                        if (nn < NNEIGH)        /* insert neighbor */
337 +                                nn++;
338 +                        for (n2 = nn; n2--; )   {
339 +                                if (!n2 || d >= nd[n2-1]) {
340 +                                        nd[n2] = d;
341 +                                        nl[n2][0] = h2;
342 +                                        nl[n2][1] = rnl[h2][n];
343 +                                        break;
344 +                                }
345 +                                nd[n2] = nd[n2-1];
346 +                                nl[n2][0] = nl[n2-1][0];
347 +                                nl[n2][1] = nl[n2-1][1];
348 +                        }
349 +                }
350 +        }
351 +        return(nn);
352 + }
353 +
354 +
355 + static void
356 + meet_neighbors( /* run through samples and their neighbors */
357 +        int     occ,
358 +        sampfunc *nf,
359 +        double  *dp
360 + )
361 + {
362 +        short   ln[NNEIGH][2];
363 +        int     nd[NNEIGH];
364 +        int     h, v, n, v2;
365 +        short   (*rnl)[NNEIGH];
366 +                                        /* initialize bottom row list */
367 +        rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
368 +        CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
369 +        for (h = 0; h < hres; h++) {
370 +                for (n = v = 0; v < vres; v++)
371 +                        if (CHK4(pixFlags, v*hres+h)) {
372 +                                rnl[h][n++] = v;
373 +                                if (n >= NNEIGH)
374 +                                        break;
375 +                        }
376 +                while (n < NNEIGH)
377 +                        rnl[h][n++] = NINF;
378 +        }
379 +        v = 0;                          /* do each row */
380 +        for ( ; ; ) {
381 +                for (h = 0; h < hres; h++) {
382 +                        if (!CHK4(pixFlags, v*hres+h) != !occ)
383 +                                continue;       /* occupancy mismatch */
384 +                                                /* find neighbors */
385 +                        n = findneigh(ln, nd, h, v, rnl);
386 +                                                /* call on neighbors */
387 +                        (*nf)(h, v, ln, nd, n, dp);
388 +                }
389 +                if (++v >= vres)                /* reinitialize row list */
390 +                        break;
391 +                for (h = 0; h < hres; h++)
392 +                        for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
393 +                                if (v2 - v > v - rnl[h][0])
394 +                                        break;          /* not close enough */
395 +                                if (CHK4(pixFlags, v2*hres+h)) {
396 +                                        for (n = 0; n < NNEIGH-1; n++)
397 +                                                rnl[h][n] = rnl[h][n+1];
398 +                                        rnl[h][NNEIGH-1] = v2;
399 +                                }
400 +                        }
401 +        }
402 +        free((void *)rnl);              /* free row list */
403   }

Diff Legend

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