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.4 by gregl, Mon Dec 1 16:34:36 1997 UTC vs.
Revision 3.22 by greg, Sat Feb 22 02:07:25 2003 UTC

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

Diff Legend

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