ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo2.c
Revision: 3.28
Committed: Thu Jan 1 11:21:55 2004 UTC (20 years, 3 months ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1
Changes since 3.27: +47 -28 lines
Log Message:
Ansification and prototypes.

File Contents

# User Rev Content
1 gregl 3.1 #ifndef lint
2 schorsch 3.28 static const char RCSid[] = "$Id: rholo2.c,v 3.27 2003/07/27 22:12:02 schorsch Exp $";
3 gregl 3.1 #endif
4     /*
5     * Rtrace support routines for holodeck rendering
6     */
7    
8 schorsch 3.28 #include <time.h>
9    
10 gregl 3.1 #include "rholo.h"
11 gregl 3.5 #include "paths.h"
12 gregl 3.1 #include "random.h"
13    
14    
15 gwlarson 3.16 VIEWPOINT myeye; /* target view position */
16    
17 gwlarson 3.19 struct gclim {
18     HOLO *hp; /* holodeck pointer */
19     GCOORD gc; /* grid cell */
20     FVECT egp; /* eye grid point */
21 gwlarson 3.20 double erg2; /* mean square eye grid range */
22 gwlarson 3.19 double gmin[2], gmax[2]; /* grid coordinate limits */
23     }; /* a grid coordinate range */
24 gwlarson 3.16
25 schorsch 3.28 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 gwlarson 3.19
31 schorsch 3.28 static void
32     initeyelim( /* initialize grid coordinate limits */
33     register struct gclim *gcl,
34     register HOLO *hp,
35     GCOORD *gc
36     )
37 gwlarson 3.19 {
38 schorsch 3.23 register RREAL *v;
39 gwlarson 3.19 register int i;
40    
41 gwlarson 3.20 if (hp != NULL) {
42     hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
43     gcl->erg2 = 0;
44     for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
45     gcl->erg2 += DOT(v,v);
46     gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
47 gwlarson 3.19 }
48 gwlarson 3.20 if (gc != NULL)
49 schorsch 3.26 gcl->gc = *gc;
50 gwlarson 3.20 gcl->gmin[0] = gcl->gmin[1] = FHUGE;
51     gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
52 gwlarson 3.19 }
53    
54    
55 schorsch 3.28 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 gwlarson 3.19 {
64 gwlarson 3.20 FVECT gp, ab;
65 gwlarson 3.21 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 gwlarson 3.20 i = gc->w>>1;
73 gwlarson 3.21 gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
74 gwlarson 3.20 gp[hdwg0[gc->w]] = gc->i[0] + r0;
75     gp[hdwg1[gc->w]] = gc->i[1] + r1;
76 gwlarson 3.19 VSUB(ab, gcl->egp, gp);
77 gwlarson 3.21 ab2 = DOT(ab, ab);
78     gw = gcl->gc.w>>1;
79     if ((i==gw ? ab[gw]*ab[gw] : ab2) <= gcl->erg2 + FTINY) {
80 gwlarson 3.19 gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
81     gcl->gmax[0] = gcl->gmax[1] = FHUGE;
82 gwlarson 3.21 return; /* too close (to wall) */
83 gwlarson 3.19 }
84 gwlarson 3.21 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 gwlarson 3.20 }
114 gwlarson 3.21 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 schorsch 3.27 (gc->w&1) ^ (root[n] < gp[gc->w>>1])) {
120 gwlarson 3.21 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 schorsch 3.27 (gc->w&1) ^ (root[n] < gp[gc->w>>1]))
146 gwlarson 3.21 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 gwlarson 3.19 }
154     }
155    
156    
157     static int
158 schorsch 3.28 clipeyelim( /* clip eye limits to grid cell */
159     register short rrng[2][2],
160     register struct gclim *gcl
161     )
162 gwlarson 3.19 {
163     int incell = 1;
164     register int i;
165    
166     for (i = 0; i < 2; i++) {
167     if (gcl->gmin[i] < gcl->gc.i[i])
168     gcl->gmin[i] = gcl->gc.i[i];
169     if (gcl->gmax[i] > gcl->gc.i[i]+1)
170     gcl->gmax[i] = gcl->gc.i[i]+1;
171 gwlarson 3.20 if (gcl->gmax[i] > gcl->gmin[i]) {
172 gwlarson 3.19 rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
173     (1.-FTINY);
174     rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
175     (1.-FTINY) - rrng[i][0];
176 gwlarson 3.20 } else
177     rrng[i][0] = rrng[i][1] = 0;
178     incell &= rrng[i][1] > 0;
179 gwlarson 3.19 }
180     return(incell);
181     }
182    
183    
184 schorsch 3.28 extern void
185     packrays( /* pack ray origins and directions */
186     register float *rod,
187     register PACKET *p
188     )
189 gregl 3.1 {
190 gwlarson 3.21 #if 0
191     double dist2sum = 0.;
192     FVECT vt;
193     #endif
194 gwlarson 3.19 int nretries = p->nr + 2;
195     struct gclim eyelim;
196     short rrng0[2][2], rrng1[2][2];
197     int useyelim;
198 gregl 3.3 GCOORD gc[2];
199 gwlarson 3.19 FVECT ro, rd;
200     double d;
201     register int i;
202 gregl 3.1
203     if (!hdbcoord(gc, hdlist[p->hd], p->bi))
204     error(CONSISTENCY, "bad beam index in packrays");
205 gwlarson 3.19 if ((useyelim = myeye.rng > FTINY)) {
206 gwlarson 3.20 initeyelim(&eyelim, hdlist[p->hd], gc);
207 gwlarson 3.21 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 gwlarson 3.19 }
215     for (i = 0; i < p->nr; i++) {
216     retry:
217     if (useyelim) {
218 gwlarson 3.20 initeyelim(&eyelim, NULL, gc+1);
219 gwlarson 3.19 p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
220     + rrng0[0][0];
221     p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
222     + rrng0[1][0];
223 gwlarson 3.20 groweyelim(&eyelim, gc,
224     (1./256.)*(p->ra[i].r[0][0]+.5),
225 gwlarson 3.21 (1./256.)*(p->ra[i].r[0][1]+.5), 1);
226 gwlarson 3.19 if (!clipeyelim(rrng1, &eyelim)) {
227 gwlarson 3.21 useyelim = nretries-- > 0;
228     #ifdef DEBUG
229     if (!useyelim)
230     error(WARNING,
231     "exceeded retry limit in packrays");
232     #endif
233 gwlarson 3.19 goto retry;
234     }
235     p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
236     + rrng1[0][0];
237     p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
238     + rrng1[1][0];
239     } else {
240     p->ra[i].r[0][0] = frandom() * 256.;
241     p->ra[i].r[0][1] = frandom() * 256.;
242     p->ra[i].r[1][0] = frandom() * 256.;
243     p->ra[i].r[1][1] = frandom() * 256.;
244     }
245     d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
246 gwlarson 3.21 #if 0
247     VSUM(vt, ro, rd, d);
248     dist2sum += dist2line(myeye.vpt, ro, vt);
249     #endif
250 gregl 3.1 if (p->offset != NULL) {
251 gwlarson 3.16 if (!vdef(OBSTRUCTIONS))
252     d *= frandom(); /* random offset */
253 gregl 3.8 VSUM(ro, ro, rd, d); /* advance ray */
254 gwlarson 3.19 p->offset[i] = d;
255 gregl 3.1 }
256 gwlarson 3.19 VCOPY(rod, ro);
257     rod += 3;
258     VCOPY(rod, rd);
259     rod += 3;
260 gregl 3.1 }
261 gwlarson 3.21 #if 0
262     fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
263     p->nr + 2 - nretries);
264     #endif
265 gregl 3.1 }
266    
267    
268 schorsch 3.28 extern void
269     donerays( /* encode finished ray computations */
270     register PACKET *p,
271     register float *rvl
272     )
273 gregl 3.1 {
274     double d;
275     register int i;
276    
277     for (i = 0; i < p->nr; i++) {
278     setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
279     d = rvl[3];
280     if (p->offset != NULL)
281     d += p->offset[i];
282     p->ra[i].d = hdcode(hdlist[p->hd], d);
283     rvl += 4;
284     }
285 gregl 3.4 p->nc += p->nr;
286 gregl 3.5 }
287    
288    
289 schorsch 3.28 extern int
290     done_rtrace(void) /* clean up and close rtrace calculation */
291 gregl 3.5 {
292     int status;
293     /* already closed? */
294     if (!nprocs)
295 greg 3.25 return(0);
296 gregl 3.5 /* flush beam queue */
297     done_packets(flush_queue());
298 gregl 3.7 /* sync holodeck */
299     hdsync(NULL, 1);
300 gregl 3.5 /* close rtrace */
301     if ((status = end_rtrace()))
302     error(WARNING, "bad exit status from rtrace");
303 gregl 3.9 if (vdef(REPORT)) { /* report time */
304 gregl 3.10 eputs("rtrace process closed\n");
305 gregl 3.5 report(0);
306 gregl 3.9 }
307 gregl 3.5 return(status); /* return status */
308     }
309    
310    
311 schorsch 3.28 extern void
312     new_rtrace(void) /* restart rtrace calculation */
313 gregl 3.5 {
314     char combuf[128];
315    
316     if (nprocs > 0) /* already running? */
317     return;
318 gregl 3.6 starttime = time(NULL); /* reset start time and counts */
319     npacksdone = nraysdone = 0L;
320     if (vdef(TIME)) /* reset end time */
321     endtime = starttime + vflt(TIME)*3600. + .5;
322 gregl 3.5 if (vdef(RIF)) { /* rerun rad to update octree */
323     sprintf(combuf, "rad -v 0 -s -w %s", vval(RIF));
324     if (system(combuf))
325     error(WARNING, "error running rad");
326     }
327     if (start_rtrace() < 1) /* start rtrace */
328     error(WARNING, "cannot restart rtrace");
329 gregl 3.9 else if (vdef(REPORT)) {
330     eputs("rtrace process restarted\n");
331 gregl 3.5 report(0);
332 gregl 3.9 }
333 gregl 3.5 }
334    
335    
336 schorsch 3.28 extern int
337     getradfile(void) /* run rad and get needed variables */
338 gregl 3.5 {
339 gwlarson 3.12 static short mvar[] = {OCTREE,EYESEP,-1};
340 gregl 3.5 static char tf1[] = TEMPLATE;
341     char tf2[64];
342     char combuf[256];
343 schorsch 3.28 char *pippt = NULL;
344 gregl 3.5 register int i;
345     register char *cp;
346     /* check if rad file specified */
347     if (!vdef(RIF))
348 gregl 3.11 return(0);
349 gregl 3.5 /* create rad command */
350     mktemp(tf1);
351     sprintf(tf2, "%s.rif", tf1);
352     sprintf(combuf,
353     "rad -v 0 -s -e -w %s OPTFILE=%s | egrep '^[ \t]*(NOMATCH",
354     vval(RIF), tf1);
355     cp = combuf;
356     while (*cp){
357     if (*cp == '|') pippt = cp;
358     cp++;
359     } /* match unset variables */
360     for (i = 0; mvar[i] >= 0; i++)
361     if (!vdef(mvar[i])) {
362     *cp++ = '|';
363     strcpy(cp, vnam(mvar[i]));
364     while (*cp) cp++;
365     pippt = NULL;
366     }
367     if (pippt != NULL)
368 schorsch 3.24 strcpy(pippt, "> " NULL_DEVICE); /* nothing to match */
369 gregl 3.5 else
370     sprintf(cp, ")[ \t]*=' > %s", tf2);
371 gwlarson 3.13 #ifdef DEBUG
372     wputs(combuf); wputs("\n");
373     #endif
374     system(combuf); /* ignore exit code */
375 gregl 3.5 if (pippt == NULL) {
376     loadvars(tf2); /* load variables */
377     unlink(tf2);
378     }
379     rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */
380     unlink(tf1); /* clean up */
381 gregl 3.11 return(1);
382 gregl 3.6 }
383    
384    
385 schorsch 3.28 extern void
386     report( /* report progress so far */
387     time_t t
388     )
389 gregl 3.6 {
390     static time_t seconds2go = 1000000;
391    
392     if (t == 0L)
393     t = time(NULL);
394     sprintf(errmsg, "%ld packets (%ld rays) done after %.2f hours\n",
395     npacksdone, nraysdone, (t-starttime)/3600.);
396     eputs(errmsg);
397     if (seconds2go == 1000000)
398     seconds2go = vdef(REPORT) ? (long)(vflt(REPORT)*60. + .5) : 0L;
399     if (seconds2go)
400     reporttime = t + seconds2go;
401 gregl 3.1 }