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

Comparing ray/src/hd/rholo2.c (file contents):
Revision 3.20 by gwlarson, Thu Dec 3 15:21:05 1998 UTC vs.
Revision 3.28 by schorsch, Thu Jan 1 11:21:55 2004 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Rtrace support routines for holodeck rendering
6   */
7  
8 + #include <time.h>
9 +
10   #include "rholo.h"
11   #include "paths.h"
12   #include "random.h"
# Line 23 | Line 22 | struct gclim {
22          double  gmin[2], gmax[2];       /* grid coordinate limits */
23   };                              /* a grid coordinate range */
24  
25 + static void initeyelim(struct gclim     *gcl, HOLO      *hp, GCOORD     *gc);
26 + static void groweyelim(struct gclim *gcl, GCOORD *gc,
27 +                double r0, double r1, int tight);
28 + static int clipeyelim(short     rrng[2][2], struct gclim        *gcl);
29  
30 < static
31 < initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
32 < register struct gclim   *gcl;
33 < register HOLO   *hp;
34 < GCOORD  *gc;
30 >
31 > static void
32 > initeyelim(             /* initialize grid coordinate limits */
33 >        register struct gclim   *gcl,
34 >        register HOLO   *hp,
35 >        GCOORD  *gc
36 > )
37   {
38 <        register FLOAT  *v;
38 >        register RREAL  *v;
39          register int    i;
40  
41          if (hp != NULL) {
# Line 41 | Line 46 | GCOORD *gc;
46                  gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
47          }
48          if (gc != NULL)
49 <                copystruct(&gcl->gc, gc);
49 >                gcl->gc = *gc;
50          gcl->gmin[0] = gcl->gmin[1] = FHUGE;
51          gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
52   }
53  
54  
55 < static
56 < groweyelim(gcl, gc, r0, r1)     /* grow grid limits about eye point */
57 < register struct gclim   *gcl;
58 < GCOORD  *gc;
59 < double  r0, r1;
55 > static void
56 > groweyelim(     /* grow grid limits about eye point */
57 >        register struct gclim   *gcl,
58 >        GCOORD  *gc,
59 >        double  r0,
60 >        double  r1,
61 >        int     tight
62 > )
63   {
64          FVECT   gp, ab;
65 <        double  vlen, plen, dv0, dv1;
66 <        double  rd2, dwall, gpos;
67 <        int     eyeout;
68 <        register int    i, g0, g1;
69 <
65 >        double  ab2, od, cfact;
66 >        double  sqcoef[3], ctcoef[3], licoef[3], cnst;
67 >        int     gw, gi[2];
68 >        double  wallpos, a, b, c, d, e, f;
69 >        double  root[2], yex;
70 >        int     n, i, j, nex;
71 >                                                /* point/view cone */
72          i = gc->w>>1;
73 <        if (gc->w&1)
64 <                eyeout = (gp[i] = gcl->hp->grid[i]) < gcl->egp[i];
65 <        else
66 <                eyeout = (gp[i] = 0) > gcl->egp[i];
73 >        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
74          gp[hdwg0[gc->w]] = gc->i[0] + r0;
75          gp[hdwg1[gc->w]] = gc->i[1] + r1;
76          VSUB(ab, gcl->egp, gp);
77 <        rd2 = DOT(ab,ab);
78 <        if (rd2 <= gcl->erg2) {
77 >        ab2 = DOT(ab, ab);
78 >        gw = gcl->gc.w>>1;
79 >        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
80                  gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
81                  gcl->gmax[0] = gcl->gmax[1] = FHUGE;
82 <                return;
82 >                return;                 /* too close (to wall) */
83          }
84 <        rd2 = gcl->erg2 / rd2;
85 <        vlen = 1. - rd2;
86 <        plen = sqrt(rd2 * vlen);
87 <        g0 = gcl->gc.w>>1;
88 <        dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0];
89 <        for (i = 0; i < 4; i++) {
90 <                if (i == 2)
91 <                        plen = -plen;
92 <                g1 = (g0+(i&1)+1)%3;
93 <                dv0 = vlen*ab[g0] + plen*ab[g1];
94 <                dv1 = vlen*ab[g1] - plen*ab[g0];
95 <                if ((dv0 < 0 ^ dwall < 0 ^ eyeout) ||
96 <                                (dv0 <= FTINY && dv0 >= -FTINY)) {
97 <                        if (eyeout)
98 <                                dv1 = -dv1;
99 <                        if (dv1 > FTINY)
100 <                                gcl->gmax[i&1] = FHUGE;
101 <                        else if (dv1 < -FTINY)
102 <                                gcl->gmin[i&1] = -FHUGE;
103 <                } else {
104 <                        gpos = gp[g1] + dv1*dwall/dv0;
105 <                        if (gpos < gcl->gmin[i&1])
106 <                                gcl->gmin[i&1] = gpos;
107 <                        if (gpos > gcl->gmax[i&1])
108 <                                gcl->gmax[i&1] = gpos;
84 >        ab2 = 1./ab2;                           /* 1/norm2(ab) */
85 >        od = DOT(gp, ab);                       /* origin dot direction */
86 >        cfact = 1./(1. - ab2*gcl->erg2);        /* tan^2 + 1 of cone angle */
87 >        for (i = 0; i < 3; i++) {               /* compute cone equation */
88 >                sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
89 >                ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
90 >                licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
91 >        }
92 >        cnst = cfact*od*od*ab2 - DOT(gp,gp);
93 >        /*
94 >         * CONE:        sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
95 >         *              + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
96 >         *              + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
97 >         */
98 >                                /* equation for conic section in plane */
99 >        gi[0] = hdwg0[gcl->gc.w];
100 >        gi[1] = hdwg1[gcl->gc.w];
101 >        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
102 >        a = sqcoef[gi[0]];                                      /* x2 */
103 >        b = ctcoef[gi[0]];                                      /* xy */
104 >        c = sqcoef[gi[1]];                                      /* y2 */
105 >        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
106 >        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
107 >        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
108 >        for (i = 0; i < 2; i++) {
109 >                if (i) {                /* swap x and y coefficients */
110 >                        register double t;
111 >                        t = a; a = c; c = t;
112 >                        t = d; d = e; e = t;
113                  }
114 +                nex = 0;                /* check global extrema */
115 +                n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
116 +                                d*(c*d-b*e) + f*b*b);
117 +                while (n-- > 0) {
118 +                        if (gc->w>>1 == gi[i] &&
119 +                                        (gc->w&1) ^ (root[n] < gp[gc->w>>1])) {
120 +                                if (gc->w&1)
121 +                                        gcl->gmin[i] = -FHUGE;
122 +                                else
123 +                                        gcl->gmax[i] = FHUGE;
124 +                                nex++;
125 +                                continue;               /* hyperbolic */
126 +                        }
127 +                        if (tight) {
128 +                                yex = (-2.*a*root[n] - d)/b;
129 +                                if (yex < gcl->gc.i[1-i] ||
130 +                                                yex > gcl->gc.i[1-i]+1)
131 +                                        continue;       /* outside cell */
132 +                        }
133 +                        if (root[n] < gcl->gmin[i])
134 +                                gcl->gmin[i] = root[n];
135 +                        if (root[n] > gcl->gmax[i])
136 +                                gcl->gmax[i] = root[n];
137 +                        nex++;
138 +                }
139 +                                        /* check local extrema */
140 +                for (j = nex < 2 ? 2 : 0; j--; ) {
141 +                        yex = gcl->gc.i[1-i] + j;
142 +                        n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
143 +                        while (n-- > 0) {
144 +                                if (gc->w>>1 == gi[i] &&
145 +                                        (gc->w&1) ^ (root[n] < gp[gc->w>>1]))
146 +                                        continue;
147 +                                if (root[n] < gcl->gmin[i])
148 +                                        gcl->gmin[i] = root[n];
149 +                                if (root[n] > gcl->gmax[i])
150 +                                        gcl->gmax[i] = root[n];
151 +                        }
152 +                }
153          }
154   }
155  
156  
157   static int
158 < clipeyelim(rrng, gcl)           /* clip eye limits to grid cell */
159 < register short  rrng[2][2];
160 < register struct gclim   *gcl;
158 > clipeyelim(             /* clip eye limits to grid cell */
159 >        register short  rrng[2][2],
160 >        register struct gclim   *gcl
161 > )
162   {
163          int     incell = 1;
164          register int    i;
# Line 129 | Line 181 | register struct gclim  *gcl;
181   }
182  
183  
184 < packrays(rod, p)                /* pack ray origins and directions */
185 < register float  *rod;
186 < register PACKET *p;
184 > extern void
185 > packrays(               /* pack ray origins and directions */
186 >        register float  *rod,
187 >        register PACKET *p
188 > )
189   {
190 + #if 0
191 +        double  dist2sum = 0.;
192 +        FVECT   vt;
193 + #endif
194          int     nretries = p->nr + 2;
195          struct gclim    eyelim;
196          short   rrng0[2][2], rrng1[2][2];
# Line 146 | Line 204 | register PACKET        *p;
204                  error(CONSISTENCY, "bad beam index in packrays");
205          if ((useyelim = myeye.rng > FTINY)) {
206                  initeyelim(&eyelim, hdlist[p->hd], gc);
207 <                groweyelim(&eyelim, gc+1, 0., 0.);
208 <                groweyelim(&eyelim, gc+1, 1., 1.);
209 <                useyelim &= clipeyelim(rrng0, &eyelim);
207 >                groweyelim(&eyelim, gc+1, 0., 0., 0);
208 >                groweyelim(&eyelim, gc+1, 1., 1., 0);
209 >                useyelim = clipeyelim(rrng0, &eyelim);
210 > #ifdef DEBUG
211 >                if (!useyelim)
212 >                        error(WARNING, "no eye overlap in packrays");
213 > #endif
214          }
215          for (i = 0; i < p->nr; i++) {
216          retry:
# Line 160 | Line 222 | register PACKET        *p;
222                                                  + rrng0[1][0];
223                          groweyelim(&eyelim, gc,
224                                          (1./256.)*(p->ra[i].r[0][0]+.5),
225 <                                        (1./256.)*(p->ra[i].r[0][1]+.5));
225 >                                        (1./256.)*(p->ra[i].r[0][1]+.5), 1);
226                          if (!clipeyelim(rrng1, &eyelim)) {
227 <                                useyelim &= nretries-- > 0;
227 >                                useyelim = nretries-- > 0;
228 > #ifdef DEBUG
229 >                                if (!useyelim)
230 >                                        error(WARNING,
231 >                                        "exceeded retry limit in packrays");
232 > #endif
233                                  goto retry;
234                          }
235                          p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
# Line 176 | Line 243 | register PACKET        *p;
243                          p->ra[i].r[1][1] = frandom() * 256.;
244                  }
245                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
246 + #if 0
247 +                VSUM(vt, ro, rd, d);
248 +                dist2sum += dist2line(myeye.vpt, ro, vt);
249 + #endif
250                  if (p->offset != NULL) {
251                          if (!vdef(OBSTRUCTIONS))
252                                  d *= frandom();         /* random offset */
# Line 187 | Line 258 | register PACKET        *p;
258                  VCOPY(rod, rd);
259                  rod += 3;
260          }
261 + #if 0
262 +        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
263 +                        p->nr + 2 - nretries);
264 + #endif
265   }
266  
267  
268 < donerays(p, rvl)                /* encode finished ray computations */
269 < register PACKET *p;
270 < register float  *rvl;
268 > extern void
269 > donerays(               /* encode finished ray computations */
270 >        register PACKET *p,
271 >        register float  *rvl
272 > )
273   {
274          double  d;
275          register int    i;
# Line 209 | Line 286 | register float *rvl;
286   }
287  
288  
289 < int
290 < done_rtrace()                   /* clean up and close rtrace calculation */
289 > extern int
290 > done_rtrace(void)                       /* clean up and close rtrace calculation */
291   {
292          int     status;
293                                          /* already closed? */
294          if (!nprocs)
295 <                return;
295 >                return(0);
296                                          /* flush beam queue */
297          done_packets(flush_queue());
298                                          /* sync holodeck */
# Line 231 | Line 308 | done_rtrace()                  /* clean up and close rtrace calculati
308   }
309  
310  
311 < new_rtrace()                    /* restart rtrace calculation */
311 > extern void
312 > new_rtrace(void)                        /* restart rtrace calculation */
313   {
314          char    combuf[128];
315  
# Line 255 | Line 333 | new_rtrace()                   /* restart rtrace calculation */
333   }
334  
335  
336 < getradfile()                    /* run rad and get needed variables */
336 > extern int
337 > getradfile(void)                        /* run rad and get needed variables */
338   {
339          static short    mvar[] = {OCTREE,EYESEP,-1};
340          static char     tf1[] = TEMPLATE;
341          char    tf2[64];
342          char    combuf[256];
343 <        char    *pippt;
343 >        char    *pippt = NULL;
344          register int    i;
345          register char   *cp;
346                                          /* check if rad file specified */
# Line 286 | Line 365 | getradfile()                   /* run rad and get needed variables */
365                          pippt = NULL;
366                  }
367          if (pippt != NULL)
368 <                strcpy(pippt, "> /dev/null");   /* nothing to match */
368 >                strcpy(pippt, "> " NULL_DEVICE);        /* nothing to match */
369          else
370                  sprintf(cp, ")[ \t]*=' > %s", tf2);
371   #ifdef DEBUG
# Line 303 | Line 382 | getradfile()                   /* run rad and get needed variables */
382   }
383  
384  
385 < report(t)                       /* report progress so far */
386 < time_t  t;
385 > extern void
386 > report(                 /* report progress so far */
387 >        time_t  t
388 > )
389   {
390          static time_t   seconds2go = 1000000;
391  

Diff Legend

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