ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo2.c
Revision: 3.19
Committed: Tue Dec 1 15:47:05 1998 UTC (25 years, 4 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.18: +162 -44 lines
Log Message:
eliminated rejection sampling, but problems remain

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     double erg[2]; /* eye range in wall grid coords */
23     double gmin[2], gmax[2]; /* grid coordinate limits */
24     }; /* a grid coordinate range */
25 gwlarson 3.16
26 gwlarson 3.19
27     static
28     initeyelim(gcl, hd, gc) /* initialize grid coordinate limits */
29     register struct gclim *gcl;
30     int hd;
31     GCOORD *gc;
32     {
33     register FLOAT *v;
34     register int i;
35    
36     gcl->hp = hdlist[hd];
37     copystruct(&gcl->gc, gc);
38     hdgrid(gcl->egp, gcl->hp, myeye.vpt);
39     for (i = 0; i < 2; i++) {
40     v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
41     gcl->erg[i] = myeye.rng * VLEN(v);
42     gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
43     }
44     }
45    
46    
47     static
48     groweyelim(gcl, gp) /* grow grid limits about eye point */
49     register struct gclim *gcl;
50     FVECT gp;
51     {
52     FVECT ab;
53     double l2, d, mult, wg;
54     register int i, g;
55    
56     VSUB(ab, gcl->egp, gp);
57     l2 = DOT(ab,ab);
58     if (l2 <= gcl->erg[0]*gcl->erg[1]) {
59     gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
60     gcl->gmax[0] = gcl->gmax[1] = FHUGE;
61     return;
62     }
63     mult = gp[g = gcl->gc.w>>1];
64     if (gcl->gc.w&1)
65     mult -= gcl->hp->grid[g];
66     if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
67     mult /= -ab[g];
68     else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
69     mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
70     ab[hdwg0[gcl->gc.w]];
71     else
72     mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
73     ab[hdwg1[gcl->gc.w]];
74     for (i = 0; i < 2; i++) {
75     g = ((gcl->gc.w>>1)+i+1)%3;
76     wg = gp[g] + mult*ab[g];
77     d = mult*gcl->erg[i];
78     if (d < 0.) d = -d;
79     if (wg - d < gcl->gmin[i])
80     gcl->gmin[i] = wg - d;
81     if (wg + d > gcl->gmax[i])
82     gcl->gmax[i] = wg + d;
83     }
84     }
85    
86    
87     static int
88     clipeyelim(rrng, gcl) /* clip eye limits to grid cell */
89     register short rrng[2][2];
90     register struct gclim *gcl;
91     {
92     int incell = 1;
93     register int i;
94    
95     for (i = 0; i < 2; i++) {
96     if (gcl->gmin[i] < gcl->gc.i[i])
97     gcl->gmin[i] = gcl->gc.i[i];
98     if (gcl->gmax[i] > gcl->gc.i[i]+1)
99     gcl->gmax[i] = gcl->gc.i[i]+1;
100     if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
101     rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
102     (1.-FTINY);
103     rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
104     (1.-FTINY) - rrng[i][0];
105     incell &= rrng[i][1] > 0;
106     }
107     }
108     return(incell);
109     }
110    
111    
112 gregl 3.1 packrays(rod, p) /* pack ray origins and directions */
113 gwlarson 3.15 register float *rod;
114 gregl 3.1 register PACKET *p;
115     {
116 gwlarson 3.19 #define gp ro
117     #ifdef DEBUG
118     double dist2sum = 0.;
119     FVECT vt;
120     #endif
121     int nretries = p->nr + 2;
122     struct gclim eyelim;
123     short rrng0[2][2], rrng1[2][2];
124     int useyelim;
125 gregl 3.3 GCOORD gc[2];
126 gwlarson 3.19 FVECT ro, rd;
127     double d;
128     register int i;
129 gregl 3.1
130     if (!hdbcoord(gc, hdlist[p->hd], p->bi))
131     error(CONSISTENCY, "bad beam index in packrays");
132 gwlarson 3.19 if ((useyelim = myeye.rng > FTINY)) {
133     initeyelim(&eyelim, p->hd, gc);
134     gp[gc[1].w>>1] = gc[1].w&1 ?
135     hdlist[p->hd]->grid[gc[1].w>>1] : 0;
136     gp[hdwg0[gc[1].w]] = gc[1].i[0];
137     gp[hdwg1[gc[1].w]] = gc[1].i[1];
138     groweyelim(&eyelim, gp);
139     gp[hdwg0[gc[1].w]]++;
140     gp[hdwg1[gc[1].w]]++;
141     groweyelim(&eyelim, gp);
142     useyelim &= clipeyelim(rrng0, &eyelim);
143     }
144     for (i = 0; i < p->nr; i++) {
145     retry:
146     if (useyelim) {
147     p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
148     + rrng0[0][0];
149     p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
150     + rrng0[1][0];
151     initeyelim(&eyelim, p->hd, gc+1);
152     gp[gc[0].w>>1] = gc[0].w&1 ?
153     hdlist[p->hd]->grid[gc[0].w>>1] : 0;
154     gp[hdwg0[gc[0].w]] = gc[0].i[0] +
155     (1./256.)*(p->ra[i].r[0][0]+.5);
156     gp[hdwg1[gc[0].w]] = gc[0].i[1] +
157     (1./256.)*(p->ra[i].r[0][1]+.5);
158     groweyelim(&eyelim, gp);
159     if (!clipeyelim(rrng1, &eyelim)) {
160     useyelim &= nretries-- > 0;
161     #ifdef DEBUG
162     if (!useyelim)
163     error(WARNING, "exceeded retry limit in packrays");
164     #endif
165     goto retry;
166     }
167     p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
168     + rrng1[0][0];
169     p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
170     + rrng1[1][0];
171     } else {
172     p->ra[i].r[0][0] = frandom() * 256.;
173     p->ra[i].r[0][1] = frandom() * 256.;
174     p->ra[i].r[1][0] = frandom() * 256.;
175     p->ra[i].r[1][1] = frandom() * 256.;
176     }
177     d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
178     #ifdef DEBUG
179     VSUM(vt, ro, rd, d);
180     dist2sum += dist2line(myeye.vpt, ro, vt);
181     #endif
182 gregl 3.1 if (p->offset != NULL) {
183 gwlarson 3.16 if (!vdef(OBSTRUCTIONS))
184     d *= frandom(); /* random offset */
185 gregl 3.8 VSUM(ro, ro, rd, d); /* advance ray */
186 gwlarson 3.19 p->offset[i] = d;
187 gregl 3.1 }
188 gwlarson 3.19 VCOPY(rod, ro);
189     rod += 3;
190     VCOPY(rod, rd);
191     rod += 3;
192 gregl 3.1 }
193 gwlarson 3.16 #ifdef DEBUG
194 gwlarson 3.19 fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
195 gwlarson 3.16 #endif
196 gwlarson 3.19 #undef gp
197 gregl 3.1 }
198    
199    
200     donerays(p, rvl) /* encode finished ray computations */
201     register PACKET *p;
202     register float *rvl;
203     {
204     double d;
205     register int i;
206    
207     for (i = 0; i < p->nr; i++) {
208     setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
209     d = rvl[3];
210     if (p->offset != NULL)
211     d += p->offset[i];
212     p->ra[i].d = hdcode(hdlist[p->hd], d);
213     rvl += 4;
214     }
215 gregl 3.4 p->nc += p->nr;
216 gregl 3.5 }
217    
218    
219     int
220     done_rtrace() /* clean up and close rtrace calculation */
221     {
222     int status;
223     /* already closed? */
224     if (!nprocs)
225     return;
226     /* flush beam queue */
227     done_packets(flush_queue());
228 gregl 3.7 /* sync holodeck */
229     hdsync(NULL, 1);
230 gregl 3.5 /* close rtrace */
231     if ((status = end_rtrace()))
232     error(WARNING, "bad exit status from rtrace");
233 gregl 3.9 if (vdef(REPORT)) { /* report time */
234 gregl 3.10 eputs("rtrace process closed\n");
235 gregl 3.5 report(0);
236 gregl 3.9 }
237 gregl 3.5 return(status); /* return status */
238     }
239    
240    
241     new_rtrace() /* restart rtrace calculation */
242     {
243     char combuf[128];
244    
245     if (nprocs > 0) /* already running? */
246     return;
247 gregl 3.6 starttime = time(NULL); /* reset start time and counts */
248     npacksdone = nraysdone = 0L;
249     if (vdef(TIME)) /* reset end time */
250     endtime = starttime + vflt(TIME)*3600. + .5;
251 gregl 3.5 if (vdef(RIF)) { /* rerun rad to update octree */
252     sprintf(combuf, "rad -v 0 -s -w %s", vval(RIF));
253     if (system(combuf))
254     error(WARNING, "error running rad");
255     }
256     if (start_rtrace() < 1) /* start rtrace */
257     error(WARNING, "cannot restart rtrace");
258 gregl 3.9 else if (vdef(REPORT)) {
259     eputs("rtrace process restarted\n");
260 gregl 3.5 report(0);
261 gregl 3.9 }
262 gregl 3.5 }
263    
264    
265     getradfile() /* run rad and get needed variables */
266     {
267 gwlarson 3.12 static short mvar[] = {OCTREE,EYESEP,-1};
268 gregl 3.5 static char tf1[] = TEMPLATE;
269     char tf2[64];
270     char combuf[256];
271     char *pippt;
272     register int i;
273     register char *cp;
274     /* check if rad file specified */
275     if (!vdef(RIF))
276 gregl 3.11 return(0);
277 gregl 3.5 /* create rad command */
278     mktemp(tf1);
279     sprintf(tf2, "%s.rif", tf1);
280     sprintf(combuf,
281     "rad -v 0 -s -e -w %s OPTFILE=%s | egrep '^[ \t]*(NOMATCH",
282     vval(RIF), tf1);
283     cp = combuf;
284     while (*cp){
285     if (*cp == '|') pippt = cp;
286     cp++;
287     } /* match unset variables */
288     for (i = 0; mvar[i] >= 0; i++)
289     if (!vdef(mvar[i])) {
290     *cp++ = '|';
291     strcpy(cp, vnam(mvar[i]));
292     while (*cp) cp++;
293     pippt = NULL;
294     }
295     if (pippt != NULL)
296     strcpy(pippt, "> /dev/null"); /* nothing to match */
297     else
298     sprintf(cp, ")[ \t]*=' > %s", tf2);
299 gwlarson 3.13 #ifdef DEBUG
300     wputs(combuf); wputs("\n");
301     #endif
302     system(combuf); /* ignore exit code */
303 gregl 3.5 if (pippt == NULL) {
304     loadvars(tf2); /* load variables */
305     unlink(tf2);
306     }
307     rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */
308     unlink(tf1); /* clean up */
309 gregl 3.11 return(1);
310 gregl 3.6 }
311    
312    
313     report(t) /* report progress so far */
314     time_t t;
315     {
316     static time_t seconds2go = 1000000;
317    
318     if (t == 0L)
319     t = time(NULL);
320     sprintf(errmsg, "%ld packets (%ld rays) done after %.2f hours\n",
321     npacksdone, nraysdone, (t-starttime)/3600.);
322     eputs(errmsg);
323     if (seconds2go == 1000000)
324     seconds2go = vdef(REPORT) ? (long)(vflt(REPORT)*60. + .5) : 0L;
325     if (seconds2go)
326     reporttime = t + seconds2go;
327 gregl 3.1 }