ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo2.c
Revision: 3.21
Committed: Mon Dec 7 16:56:08 1998 UTC (25 years, 4 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.20: +108 -43 lines
Log Message:
went to more sophisticated cone extrema for view vicinity ray samples

File Contents

# User Rev Content
1 gwlarson 3.12 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 gregl 3.1
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7     /*
8     * Rtrace support routines for holodeck rendering
9     */
10    
11     #include "rholo.h"
12 gregl 3.5 #include "paths.h"
13 gregl 3.1 #include "random.h"
14    
15    
16 gwlarson 3.16 VIEWPOINT myeye; /* target view position */
17    
18 gwlarson 3.19 struct gclim {
19     HOLO *hp; /* holodeck pointer */
20     GCOORD gc; /* grid cell */
21     FVECT egp; /* eye grid point */
22 gwlarson 3.20 double erg2; /* mean square eye grid range */
23 gwlarson 3.19 double gmin[2], gmax[2]; /* grid coordinate limits */
24     }; /* a grid coordinate range */
25 gwlarson 3.16
26 gwlarson 3.19
27     static
28 gwlarson 3.20 initeyelim(gcl, hp, gc) /* initialize grid coordinate limits */
29 gwlarson 3.19 register struct gclim *gcl;
30 gwlarson 3.20 register HOLO *hp;
31 gwlarson 3.19 GCOORD *gc;
32     {
33     register FLOAT *v;
34     register int i;
35    
36 gwlarson 3.20 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 gwlarson 3.19 }
43 gwlarson 3.20 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 gwlarson 3.19 }
48    
49    
50     static
51 gwlarson 3.21 groweyelim(gcl, gc, r0, r1, tight) /* grow grid limits about eye point */
52 gwlarson 3.19 register struct gclim *gcl;
53 gwlarson 3.20 GCOORD *gc;
54     double r0, r1;
55 gwlarson 3.21 int tight;
56 gwlarson 3.19 {
57 gwlarson 3.20 FVECT gp, ab;
58 gwlarson 3.21 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 gwlarson 3.20 i = gc->w>>1;
66 gwlarson 3.21 gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
67 gwlarson 3.20 gp[hdwg0[gc->w]] = gc->i[0] + r0;
68     gp[hdwg1[gc->w]] = gc->i[1] + r1;
69 gwlarson 3.19 VSUB(ab, gcl->egp, gp);
70 gwlarson 3.21 ab2 = DOT(ab, ab);
71     gw = gcl->gc.w>>1;
72     if ((i==gw ? ab[gw]*ab[gw] : ab2) <= gcl->erg2 + FTINY) {
73 gwlarson 3.19 gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
74     gcl->gmax[0] = gcl->gmax[1] = FHUGE;
75 gwlarson 3.21 return; /* too close (to wall) */
76 gwlarson 3.19 }
77 gwlarson 3.21 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 gwlarson 3.20 }
107 gwlarson 3.21 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 gwlarson 3.19 }
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 gwlarson 3.20 if (gcl->gmax[i] > gcl->gmin[i]) {
164 gwlarson 3.19 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 gwlarson 3.20 } else
169     rrng[i][0] = rrng[i][1] = 0;
170     incell &= rrng[i][1] > 0;
171 gwlarson 3.19 }
172     return(incell);
173     }
174    
175    
176 gregl 3.1 packrays(rod, p) /* pack ray origins and directions */
177 gwlarson 3.15 register float *rod;
178 gregl 3.1 register PACKET *p;
179     {
180 gwlarson 3.21 #if 0
181     double dist2sum = 0.;
182     FVECT vt;
183     #endif
184 gwlarson 3.19 int nretries = p->nr + 2;
185     struct gclim eyelim;
186     short rrng0[2][2], rrng1[2][2];
187     int useyelim;
188 gregl 3.3 GCOORD gc[2];
189 gwlarson 3.19 FVECT ro, rd;
190     double d;
191     register int i;
192 gregl 3.1
193     if (!hdbcoord(gc, hdlist[p->hd], p->bi))
194     error(CONSISTENCY, "bad beam index in packrays");
195 gwlarson 3.19 if ((useyelim = myeye.rng > FTINY)) {
196 gwlarson 3.20 initeyelim(&eyelim, hdlist[p->hd], gc);
197 gwlarson 3.21 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 gwlarson 3.19 }
205     for (i = 0; i < p->nr; i++) {
206     retry:
207     if (useyelim) {
208 gwlarson 3.20 initeyelim(&eyelim, NULL, gc+1);
209 gwlarson 3.19 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 gwlarson 3.20 groweyelim(&eyelim, gc,
214     (1./256.)*(p->ra[i].r[0][0]+.5),
215 gwlarson 3.21 (1./256.)*(p->ra[i].r[0][1]+.5), 1);
216 gwlarson 3.19 if (!clipeyelim(rrng1, &eyelim)) {
217 gwlarson 3.21 useyelim = nretries-- > 0;
218     #ifdef DEBUG
219     if (!useyelim)
220     error(WARNING,
221     "exceeded retry limit in packrays");
222     #endif
223 gwlarson 3.19 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 gwlarson 3.21 #if 0
237     VSUM(vt, ro, rd, d);
238     dist2sum += dist2line(myeye.vpt, ro, vt);
239     #endif
240 gregl 3.1 if (p->offset != NULL) {
241 gwlarson 3.16 if (!vdef(OBSTRUCTIONS))
242     d *= frandom(); /* random offset */
243 gregl 3.8 VSUM(ro, ro, rd, d); /* advance ray */
244 gwlarson 3.19 p->offset[i] = d;
245 gregl 3.1 }
246 gwlarson 3.19 VCOPY(rod, ro);
247     rod += 3;
248     VCOPY(rod, rd);
249     rod += 3;
250 gregl 3.1 }
251 gwlarson 3.21 #if 0
252     fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
253     p->nr + 2 - nretries);
254     #endif
255 gregl 3.1 }
256    
257    
258     donerays(p, rvl) /* encode finished ray computations */
259     register PACKET *p;
260     register float *rvl;
261     {
262     double d;
263     register int i;
264    
265     for (i = 0; i < p->nr; i++) {
266     setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
267     d = rvl[3];
268     if (p->offset != NULL)
269     d += p->offset[i];
270     p->ra[i].d = hdcode(hdlist[p->hd], d);
271     rvl += 4;
272     }
273 gregl 3.4 p->nc += p->nr;
274 gregl 3.5 }
275    
276    
277     int
278     done_rtrace() /* clean up and close rtrace calculation */
279     {
280     int status;
281     /* already closed? */
282     if (!nprocs)
283     return;
284     /* flush beam queue */
285     done_packets(flush_queue());
286 gregl 3.7 /* sync holodeck */
287     hdsync(NULL, 1);
288 gregl 3.5 /* close rtrace */
289     if ((status = end_rtrace()))
290     error(WARNING, "bad exit status from rtrace");
291 gregl 3.9 if (vdef(REPORT)) { /* report time */
292 gregl 3.10 eputs("rtrace process closed\n");
293 gregl 3.5 report(0);
294 gregl 3.9 }
295 gregl 3.5 return(status); /* return status */
296     }
297    
298    
299     new_rtrace() /* restart rtrace calculation */
300     {
301     char combuf[128];
302    
303     if (nprocs > 0) /* already running? */
304     return;
305 gregl 3.6 starttime = time(NULL); /* reset start time and counts */
306     npacksdone = nraysdone = 0L;
307     if (vdef(TIME)) /* reset end time */
308     endtime = starttime + vflt(TIME)*3600. + .5;
309 gregl 3.5 if (vdef(RIF)) { /* rerun rad to update octree */
310     sprintf(combuf, "rad -v 0 -s -w %s", vval(RIF));
311     if (system(combuf))
312     error(WARNING, "error running rad");
313     }
314     if (start_rtrace() < 1) /* start rtrace */
315     error(WARNING, "cannot restart rtrace");
316 gregl 3.9 else if (vdef(REPORT)) {
317     eputs("rtrace process restarted\n");
318 gregl 3.5 report(0);
319 gregl 3.9 }
320 gregl 3.5 }
321    
322    
323     getradfile() /* run rad and get needed variables */
324     {
325 gwlarson 3.12 static short mvar[] = {OCTREE,EYESEP,-1};
326 gregl 3.5 static char tf1[] = TEMPLATE;
327     char tf2[64];
328     char combuf[256];
329     char *pippt;
330     register int i;
331     register char *cp;
332     /* check if rad file specified */
333     if (!vdef(RIF))
334 gregl 3.11 return(0);
335 gregl 3.5 /* create rad command */
336     mktemp(tf1);
337     sprintf(tf2, "%s.rif", tf1);
338     sprintf(combuf,
339     "rad -v 0 -s -e -w %s OPTFILE=%s | egrep '^[ \t]*(NOMATCH",
340     vval(RIF), tf1);
341     cp = combuf;
342     while (*cp){
343     if (*cp == '|') pippt = cp;
344     cp++;
345     } /* match unset variables */
346     for (i = 0; mvar[i] >= 0; i++)
347     if (!vdef(mvar[i])) {
348     *cp++ = '|';
349     strcpy(cp, vnam(mvar[i]));
350     while (*cp) cp++;
351     pippt = NULL;
352     }
353     if (pippt != NULL)
354     strcpy(pippt, "> /dev/null"); /* nothing to match */
355     else
356     sprintf(cp, ")[ \t]*=' > %s", tf2);
357 gwlarson 3.13 #ifdef DEBUG
358     wputs(combuf); wputs("\n");
359     #endif
360     system(combuf); /* ignore exit code */
361 gregl 3.5 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 gregl 3.11 return(1);
368 gregl 3.6 }
369    
370    
371     report(t) /* report progress so far */
372     time_t t;
373     {
374     static time_t seconds2go = 1000000;
375    
376     if (t == 0L)
377     t = time(NULL);
378     sprintf(errmsg, "%ld packets (%ld rays) done after %.2f hours\n",
379     npacksdone, nraysdone, (t-starttime)/3600.);
380     eputs(errmsg);
381     if (seconds2go == 1000000)
382     seconds2go = vdef(REPORT) ? (long)(vflt(REPORT)*60. + .5) : 0L;
383     if (seconds2go)
384     reporttime = t + seconds2go;
385 gregl 3.1 }