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.6 by gregl, Fri Dec 12 18:33:50 1997 UTC vs.
Revision 3.21 by gwlarson, Mon Dec 7 16:56:08 1998 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1997 Silicon Graphics, Inc. */
1 > /* Copyright (c) 1998 Silicon Graphics, Inc. */
2  
3   #ifndef lint
4   static char SCCSid[] = "$SunId$ SGI";
# Line 13 | Line 13 | static char SCCSid[] = "$SunId$ SGI";
13   #include "random.h"
14  
15  
16 + VIEWPOINT       myeye;          /* target view position */
17 +
18 + struct gclim {
19 +        HOLO    *hp;                    /* holodeck pointer */
20 +        GCOORD  gc;                     /* grid cell */
21 +        FVECT   egp;                    /* eye grid point */
22 +        double  erg2;                   /* mean square eye grid range */
23 +        double  gmin[2], gmax[2];       /* grid coordinate limits */
24 + };                              /* a grid coordinate range */
25 +
26 +
27 + static
28 + initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
29 + register struct gclim   *gcl;
30 + register HOLO   *hp;
31 + GCOORD  *gc;
32 + {
33 +        register FLOAT  *v;
34 +        register int    i;
35 +
36 +        if (hp != NULL) {
37 +                hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
38 +                gcl->erg2 = 0;
39 +                for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
40 +                        gcl->erg2 += DOT(v,v);
41 +                gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
42 +        }
43 +        if (gc != NULL)
44 +                copystruct(&gcl->gc, gc);
45 +        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
46 +        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
47 + }
48 +
49 +
50 + static
51 + groweyelim(gcl, gc, r0, r1, tight)      /* grow grid limits about eye point */
52 + register struct gclim   *gcl;
53 + GCOORD  *gc;
54 + double  r0, r1;
55 + int     tight;
56 + {
57 +        FVECT   gp, ab;
58 +        double  ab2, od, cfact;
59 +        double  sqcoef[3], ctcoef[3], licoef[3], cnst;
60 +        int     gw, gi[2];
61 +        double  wallpos, a, b, c, d, e, f;
62 +        double  root[2], yex;
63 +        int     n, i, j, nex;
64 +                                                /* point/view cone */
65 +        i = gc->w>>1;
66 +        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
67 +        gp[hdwg0[gc->w]] = gc->i[0] + r0;
68 +        gp[hdwg1[gc->w]] = gc->i[1] + r1;
69 +        VSUB(ab, gcl->egp, gp);
70 +        ab2 = DOT(ab, ab);
71 +        gw = gcl->gc.w>>1;
72 +        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
73 +                gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
74 +                gcl->gmax[0] = gcl->gmax[1] = FHUGE;
75 +                return;                 /* too close (to wall) */
76 +        }
77 +        ab2 = 1./ab2;                           /* 1/norm2(ab) */
78 +        od = DOT(gp, ab);                       /* origin dot direction */
79 +        cfact = 1./(1. - ab2*gcl->erg2);        /* tan^2 + 1 of cone angle */
80 +        for (i = 0; i < 3; i++) {               /* compute cone equation */
81 +                sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
82 +                ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
83 +                licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
84 +        }
85 +        cnst = cfact*od*od*ab2 - DOT(gp,gp);
86 +        /*
87 +         * CONE:        sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
88 +         *              + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
89 +         *              + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
90 +         */
91 +                                /* equation for conic section in plane */
92 +        gi[0] = hdwg0[gcl->gc.w];
93 +        gi[1] = hdwg1[gcl->gc.w];
94 +        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
95 +        a = sqcoef[gi[0]];                                      /* x2 */
96 +        b = ctcoef[gi[0]];                                      /* xy */
97 +        c = sqcoef[gi[1]];                                      /* y2 */
98 +        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
99 +        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
100 +        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
101 +        for (i = 0; i < 2; i++) {
102 +                if (i) {                /* swap x and y coefficients */
103 +                        register double t;
104 +                        t = a; a = c; c = t;
105 +                        t = d; d = e; e = t;
106 +                }
107 +                nex = 0;                /* check global extrema */
108 +                n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
109 +                                d*(c*d-b*e) + f*b*b);
110 +                while (n-- > 0) {
111 +                        if (gc->w>>1 == gi[i] &&
112 +                                        (gc->w&1) ^ root[n] < gp[gc->w>>1]) {
113 +                                if (gc->w&1)
114 +                                        gcl->gmin[i] = -FHUGE;
115 +                                else
116 +                                        gcl->gmax[i] = FHUGE;
117 +                                nex++;
118 +                                continue;               /* hyperbolic */
119 +                        }
120 +                        if (tight) {
121 +                                yex = (-2.*a*root[n] - d)/b;
122 +                                if (yex < gcl->gc.i[1-i] ||
123 +                                                yex > gcl->gc.i[1-i]+1)
124 +                                        continue;       /* outside cell */
125 +                        }
126 +                        if (root[n] < gcl->gmin[i])
127 +                                gcl->gmin[i] = root[n];
128 +                        if (root[n] > gcl->gmax[i])
129 +                                gcl->gmax[i] = root[n];
130 +                        nex++;
131 +                }
132 +                                        /* check local extrema */
133 +                for (j = nex < 2 ? 2 : 0; j--; ) {
134 +                        yex = gcl->gc.i[1-i] + j;
135 +                        n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
136 +                        while (n-- > 0) {
137 +                                if (gc->w>>1 == gi[i] &&
138 +                                        (gc->w&1) ^ root[n] < gp[gc->w>>1])
139 +                                        continue;
140 +                                if (root[n] < gcl->gmin[i])
141 +                                        gcl->gmin[i] = root[n];
142 +                                if (root[n] > gcl->gmax[i])
143 +                                        gcl->gmax[i] = root[n];
144 +                        }
145 +                }
146 +        }
147 + }
148 +
149 +
150 + static int
151 + clipeyelim(rrng, gcl)           /* clip eye limits to grid cell */
152 + register short  rrng[2][2];
153 + register struct gclim   *gcl;
154 + {
155 +        int     incell = 1;
156 +        register int    i;
157 +
158 +        for (i = 0; i < 2; i++) {
159 +                if (gcl->gmin[i] < gcl->gc.i[i])
160 +                        gcl->gmin[i] = gcl->gc.i[i];
161 +                if (gcl->gmax[i] > gcl->gc.i[i]+1)
162 +                        gcl->gmax[i] = gcl->gc.i[i]+1;
163 +                if (gcl->gmax[i] > gcl->gmin[i]) {
164 +                        rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
165 +                                        (1.-FTINY);
166 +                        rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
167 +                                        (1.-FTINY) - rrng[i][0];
168 +                } else
169 +                        rrng[i][0] = rrng[i][1] = 0;
170 +                incell &= rrng[i][1] > 0;
171 +        }
172 +        return(incell);
173 + }
174 +
175 +
176   packrays(rod, p)                /* pack ray origins and directions */
177   register float  *rod;
178   register PACKET *p;
179   {
180 <        static FVECT    ro, rd;
180 > #if 0
181 >        double  dist2sum = 0.;
182 >        FVECT   vt;
183 > #endif
184 >        int     nretries = p->nr + 2;
185 >        struct gclim    eyelim;
186 >        short   rrng0[2][2], rrng1[2][2];
187 >        int     useyelim;
188          GCOORD  gc[2];
189 <        int     ila[2], hsh;
190 <        double  d, sl[4];
189 >        FVECT   ro, rd;
190 >        double  d;
191          register int    i;
192  
193          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
194                  error(CONSISTENCY, "bad beam index in packrays");
195 <        ila[0] = p->hd; ila[1] = p->bi;
196 <        hsh = ilhash(ila,2) + p->nc;
195 >        if ((useyelim = myeye.rng > FTINY)) {
196 >                initeyelim(&eyelim, hdlist[p->hd], gc);
197 >                groweyelim(&eyelim, gc+1, 0., 0., 0);
198 >                groweyelim(&eyelim, gc+1, 1., 1., 0);
199 >                useyelim = clipeyelim(rrng0, &eyelim);
200 > #ifdef DEBUG
201 >                if (!useyelim)
202 >                        error(WARNING, "no eye overlap in packrays");
203 > #endif
204 >        }
205          for (i = 0; i < p->nr; i++) {
206 <                multisamp(sl, 4, urand(hsh+i));
207 <                p->ra[i].r[0][0] = sl[0] * 256.;
208 <                p->ra[i].r[0][1] = sl[1] * 256.;
209 <                p->ra[i].r[1][0] = sl[2] * 256.;
210 <                p->ra[i].r[1][1] = sl[3] * 256.;
206 >        retry:
207 >                if (useyelim) {
208 >                        initeyelim(&eyelim, NULL, gc+1);
209 >                        p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
210 >                                                + rrng0[0][0];
211 >                        p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
212 >                                                + rrng0[1][0];
213 >                        groweyelim(&eyelim, gc,
214 >                                        (1./256.)*(p->ra[i].r[0][0]+.5),
215 >                                        (1./256.)*(p->ra[i].r[0][1]+.5), 1);
216 >                        if (!clipeyelim(rrng1, &eyelim)) {
217 >                                useyelim = nretries-- > 0;
218 > #ifdef DEBUG
219 >                                if (!useyelim)
220 >                                        error(WARNING,
221 >                                        "exceeded retry limit in packrays");
222 > #endif
223 >                                goto retry;
224 >                        }
225 >                        p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
226 >                                                + rrng1[0][0];
227 >                        p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
228 >                                                + rrng1[1][0];
229 >                } else {
230 >                        p->ra[i].r[0][0] = frandom() * 256.;
231 >                        p->ra[i].r[0][1] = frandom() * 256.;
232 >                        p->ra[i].r[1][0] = frandom() * 256.;
233 >                        p->ra[i].r[1][1] = frandom() * 256.;
234 >                }
235                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
236 + #if 0
237 +                VSUM(vt, ro, rd, d);
238 +                dist2sum += dist2line(myeye.vpt, ro, vt);
239 + #endif
240                  if (p->offset != NULL) {
241 <                        VSUM(ro, ro, rd, d);            /* exterior only */
241 >                        if (!vdef(OBSTRUCTIONS))
242 >                                d *= frandom();         /* random offset */
243 >                        VSUM(ro, ro, rd, d);            /* advance ray */
244                          p->offset[i] = d;
245                  }
246                  VCOPY(rod, ro);
# Line 43 | Line 248 | register PACKET        *p;
248                  VCOPY(rod, rd);
249                  rod += 3;
250          }
251 + #if 0
252 +        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
253 +                        p->nr + 2 - nretries);
254 + #endif
255   }
256  
257  
# Line 72 | Line 281 | done_rtrace()                  /* clean up and close rtrace calculati
281                                          /* already closed? */
282          if (!nprocs)
283                  return;
75 –        wputs("closing rtrace process...\n");
284                                          /* flush beam queue */
285          done_packets(flush_queue());
286 +                                        /* sync holodeck */
287 +        hdsync(NULL, 1);
288                                          /* close rtrace */
289          if ((status = end_rtrace()))
290                  error(WARNING, "bad exit status from rtrace");
291 <        if (vdef(REPORT))               /* report time */
291 >        if (vdef(REPORT)) {             /* report time */
292 >                eputs("rtrace process closed\n");
293                  report(0);
294 +        }
295          return(status);                 /* return status */
296   }
297  
# Line 90 | Line 302 | new_rtrace()                   /* restart rtrace calculation */
302  
303          if (nprocs > 0)                 /* already running? */
304                  return;
93 –        wputs("restarting rtrace process...\n");
305          starttime = time(NULL);         /* reset start time and counts */
306          npacksdone = nraysdone = 0L;
307          if (vdef(TIME))                 /* reset end time */
# Line 102 | Line 313 | new_rtrace()                   /* restart rtrace calculation */
313          }
314          if (start_rtrace() < 1)         /* start rtrace */
315                  error(WARNING, "cannot restart rtrace");
316 <        else if (vdef(REPORT))
316 >        else if (vdef(REPORT)) {
317 >                eputs("rtrace process restarted\n");
318                  report(0);
319 +        }
320   }
321  
322  
323   getradfile()                    /* run rad and get needed variables */
324   {
325 <        static short    mvar[] = {OCTREE,-1};
325 >        static short    mvar[] = {OCTREE,EYESEP,-1};
326          static char     tf1[] = TEMPLATE;
327          char    tf2[64];
328          char    combuf[256];
# Line 118 | Line 331 | getradfile()                   /* run rad and get needed variables */
331          register char   *cp;
332                                          /* check if rad file specified */
333          if (!vdef(RIF))
334 <                return;
334 >                return(0);
335                                          /* create rad command */
336          mktemp(tf1);
337          sprintf(tf2, "%s.rif", tf1);
# Line 141 | Line 354 | getradfile()                   /* run rad and get needed variables */
354                  strcpy(pippt, "> /dev/null");   /* nothing to match */
355          else
356                  sprintf(cp, ")[ \t]*=' > %s", tf2);
357 <        if (system(combuf)) {
358 <                unlink(tf2);                    /* clean up */
359 <                unlink(tf1);
360 <                error(SYSTEM, "cannot execute rad command");
148 <        }
357 > #ifdef DEBUG
358 >        wputs(combuf); wputs("\n");
359 > #endif
360 >        system(combuf);                         /* ignore exit code */
361          if (pippt == NULL) {
362                  loadvars(tf2);                  /* load variables */
363                  unlink(tf2);
364          }
365          rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */
366          unlink(tf1);                    /* clean up */
367 +        return(1);
368   }
369  
370  

Diff Legend

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