ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo2.c
Revision: 3.26
Committed: Mon Jul 21 22:30:18 2003 UTC (20 years, 8 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 3.25: +2 -2 lines
Log Message:
Eliminated copystruct() macro, which is unnecessary in ANSI.
Reduced ambiguity warnings for nested if/if/else clauses.

File Contents

# User Rev Content
1 gregl 3.1 #ifndef lint
2 schorsch 3.26 static const char RCSid[] = "$Id: rholo2.c,v 3.25 2003/07/07 17:21:51 greg Exp $";
3 gregl 3.1 #endif
4     /*
5     * Rtrace support routines for holodeck rendering
6     */
7    
8     #include "rholo.h"
9 gregl 3.5 #include "paths.h"
10 gregl 3.1 #include "random.h"
11    
12    
13 gwlarson 3.16 VIEWPOINT myeye; /* target view position */
14    
15 gwlarson 3.19 struct gclim {
16     HOLO *hp; /* holodeck pointer */
17     GCOORD gc; /* grid cell */
18     FVECT egp; /* eye grid point */
19 gwlarson 3.20 double erg2; /* mean square eye grid range */
20 gwlarson 3.19 double gmin[2], gmax[2]; /* grid coordinate limits */
21     }; /* a grid coordinate range */
22 gwlarson 3.16
23 gwlarson 3.19
24     static
25 gwlarson 3.20 initeyelim(gcl, hp, gc) /* initialize grid coordinate limits */
26 gwlarson 3.19 register struct gclim *gcl;
27 gwlarson 3.20 register HOLO *hp;
28 gwlarson 3.19 GCOORD *gc;
29     {
30 schorsch 3.23 register RREAL *v;
31 gwlarson 3.19 register int i;
32    
33 gwlarson 3.20 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 gwlarson 3.19 }
40 gwlarson 3.20 if (gc != NULL)
41 schorsch 3.26 gcl->gc = *gc;
42 gwlarson 3.20 gcl->gmin[0] = gcl->gmin[1] = FHUGE;
43     gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
44 gwlarson 3.19 }
45    
46    
47     static
48 gwlarson 3.21 groweyelim(gcl, gc, r0, r1, tight) /* grow grid limits about eye point */
49 gwlarson 3.19 register struct gclim *gcl;
50 gwlarson 3.20 GCOORD *gc;
51     double r0, r1;
52 gwlarson 3.21 int tight;
53 gwlarson 3.19 {
54 gwlarson 3.20 FVECT gp, ab;
55 gwlarson 3.21 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 gwlarson 3.20 i = gc->w>>1;
63 gwlarson 3.21 gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
64 gwlarson 3.20 gp[hdwg0[gc->w]] = gc->i[0] + r0;
65     gp[hdwg1[gc->w]] = gc->i[1] + r1;
66 gwlarson 3.19 VSUB(ab, gcl->egp, gp);
67 gwlarson 3.21 ab2 = DOT(ab, ab);
68     gw = gcl->gc.w>>1;
69     if ((i==gw ? ab[gw]*ab[gw] : ab2) <= gcl->erg2 + FTINY) {
70 gwlarson 3.19 gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
71     gcl->gmax[0] = gcl->gmax[1] = FHUGE;
72 gwlarson 3.21 return; /* too close (to wall) */
73 gwlarson 3.19 }
74 gwlarson 3.21 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 gwlarson 3.20 }
104 gwlarson 3.21 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 gwlarson 3.19 }
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 gwlarson 3.20 if (gcl->gmax[i] > gcl->gmin[i]) {
161 gwlarson 3.19 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 gwlarson 3.20 } else
166     rrng[i][0] = rrng[i][1] = 0;
167     incell &= rrng[i][1] > 0;
168 gwlarson 3.19 }
169     return(incell);
170     }
171    
172    
173 gregl 3.1 packrays(rod, p) /* pack ray origins and directions */
174 gwlarson 3.15 register float *rod;
175 gregl 3.1 register PACKET *p;
176     {
177 gwlarson 3.21 #if 0
178     double dist2sum = 0.;
179     FVECT vt;
180     #endif
181 gwlarson 3.19 int nretries = p->nr + 2;
182     struct gclim eyelim;
183     short rrng0[2][2], rrng1[2][2];
184     int useyelim;
185 gregl 3.3 GCOORD gc[2];
186 gwlarson 3.19 FVECT ro, rd;
187     double d;
188     register int i;
189 gregl 3.1
190     if (!hdbcoord(gc, hdlist[p->hd], p->bi))
191     error(CONSISTENCY, "bad beam index in packrays");
192 gwlarson 3.19 if ((useyelim = myeye.rng > FTINY)) {
193 gwlarson 3.20 initeyelim(&eyelim, hdlist[p->hd], gc);
194 gwlarson 3.21 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 gwlarson 3.19 }
202     for (i = 0; i < p->nr; i++) {
203     retry:
204     if (useyelim) {
205 gwlarson 3.20 initeyelim(&eyelim, NULL, gc+1);
206 gwlarson 3.19 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 gwlarson 3.20 groweyelim(&eyelim, gc,
211     (1./256.)*(p->ra[i].r[0][0]+.5),
212 gwlarson 3.21 (1./256.)*(p->ra[i].r[0][1]+.5), 1);
213 gwlarson 3.19 if (!clipeyelim(rrng1, &eyelim)) {
214 gwlarson 3.21 useyelim = nretries-- > 0;
215     #ifdef DEBUG
216     if (!useyelim)
217     error(WARNING,
218     "exceeded retry limit in packrays");
219     #endif
220 gwlarson 3.19 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 gwlarson 3.21 #if 0
234     VSUM(vt, ro, rd, d);
235     dist2sum += dist2line(myeye.vpt, ro, vt);
236     #endif
237 gregl 3.1 if (p->offset != NULL) {
238 gwlarson 3.16 if (!vdef(OBSTRUCTIONS))
239     d *= frandom(); /* random offset */
240 gregl 3.8 VSUM(ro, ro, rd, d); /* advance ray */
241 gwlarson 3.19 p->offset[i] = d;
242 gregl 3.1 }
243 gwlarson 3.19 VCOPY(rod, ro);
244     rod += 3;
245     VCOPY(rod, rd);
246     rod += 3;
247 gregl 3.1 }
248 gwlarson 3.21 #if 0
249     fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
250     p->nr + 2 - nretries);
251     #endif
252 gregl 3.1 }
253    
254    
255     donerays(p, rvl) /* encode finished ray computations */
256     register PACKET *p;
257     register float *rvl;
258     {
259     double d;
260     register int i;
261    
262     for (i = 0; i < p->nr; i++) {
263     setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
264     d = rvl[3];
265     if (p->offset != NULL)
266     d += p->offset[i];
267     p->ra[i].d = hdcode(hdlist[p->hd], d);
268     rvl += 4;
269     }
270 gregl 3.4 p->nc += p->nr;
271 gregl 3.5 }
272    
273    
274     int
275     done_rtrace() /* clean up and close rtrace calculation */
276     {
277     int status;
278     /* already closed? */
279     if (!nprocs)
280 greg 3.25 return(0);
281 gregl 3.5 /* flush beam queue */
282     done_packets(flush_queue());
283 gregl 3.7 /* sync holodeck */
284     hdsync(NULL, 1);
285 gregl 3.5 /* close rtrace */
286     if ((status = end_rtrace()))
287     error(WARNING, "bad exit status from rtrace");
288 gregl 3.9 if (vdef(REPORT)) { /* report time */
289 gregl 3.10 eputs("rtrace process closed\n");
290 gregl 3.5 report(0);
291 gregl 3.9 }
292 gregl 3.5 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 gregl 3.6 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 gregl 3.5 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 gregl 3.9 else if (vdef(REPORT)) {
314     eputs("rtrace process restarted\n");
315 gregl 3.5 report(0);
316 gregl 3.9 }
317 gregl 3.5 }
318    
319    
320     getradfile() /* run rad and get needed variables */
321     {
322 gwlarson 3.12 static short mvar[] = {OCTREE,EYESEP,-1};
323 gregl 3.5 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 gregl 3.11 return(0);
332 gregl 3.5 /* 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 schorsch 3.24 strcpy(pippt, "> " NULL_DEVICE); /* nothing to match */
352 gregl 3.5 else
353     sprintf(cp, ")[ \t]*=' > %s", tf2);
354 gwlarson 3.13 #ifdef DEBUG
355     wputs(combuf); wputs("\n");
356     #endif
357     system(combuf); /* ignore exit code */
358 gregl 3.5 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 gregl 3.11 return(1);
365 gregl 3.6 }
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 gregl 3.1 }