src/hd/rholo2.c

/* Copyright (c) 1998 Silicon Graphics, Inc. */

#ifndef lint
static char SCCSid[] = "$SunId$ SGI";
#endif

/*
 * Rtrace support routines for holodeck rendering
 */

#include "rholo.h"
#include "paths.h"
#include "random.h"


VIEWPOINT       myeye;          /* target view position */

struct gclim {
        HOLO    *hp;                    /* holodeck pointer */
        GCOORD  gc;                     /* grid cell */
        FVECT   egp;                    /* eye grid point */
        double  erg2;                   /* mean square eye grid range */
        double  gmin[2], gmax[2];       /* grid coordinate limits */
};                              /* a grid coordinate range */


static
initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
register struct gclim   *gcl;
register HOLO   *hp;
GCOORD  *gc;
{
        register FLOAT  *v;
        register int    i;

        if (hp != NULL) {
                hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
                gcl->erg2 = 0;
                for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
                        gcl->erg2 += DOT(v,v);
                gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
        }
        if (gc != NULL)
                copystruct(&gcl->gc, gc);
        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
}


static
groweyelim(gcl, gc, r0, r1)     /* grow grid limits about eye point */
register struct gclim   *gcl;
GCOORD  *gc;
double  r0, r1;
{
        FVECT   gp, ab;
        double  vlen, plen, dv0, dv1;
        double  rd2, dwall, gpos;
        int     eyeout;
        register int    i, g0, g1;

        i = gc->w>>1;
        if (gc->w&1)
                eyeout = (gp[i] = gcl->hp->grid[i]) < gcl->egp[i];
        else
                eyeout = (gp[i] = 0) > gcl->egp[i];
        gp[hdwg0[gc->w]] = gc->i[0] + r0;
        gp[hdwg1[gc->w]] = gc->i[1] + r1;
        VSUB(ab, gcl->egp, gp);
        rd2 = DOT(ab,ab);
        if (rd2 <= gcl->erg2) {
                gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
                gcl->gmax[0] = gcl->gmax[1] = FHUGE;
                return;
        }
        rd2 = gcl->erg2 / rd2;
        vlen = 1. - rd2;
        plen = sqrt(rd2 * vlen);
        g0 = gcl->gc.w>>1;
        dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0];
        for (i = 0; i < 4; i++) {
                if (i == 2)
                        plen = -plen;
                g1 = (g0+(i&1)+1)%3;
                dv0 = vlen*ab[g0] + plen*ab[g1];
                dv1 = vlen*ab[g1] - plen*ab[g0];
                if ((dv0 < 0 ^ dwall < 0 ^ eyeout) ||
                                (dv0 <= FTINY && dv0 >= -FTINY)) {
                        if (eyeout)
                                dv1 = -dv1;
                        if (dv1 > FTINY)
                                gcl->gmax[i&1] = FHUGE;
                        else if (dv1 < -FTINY)
                                gcl->gmin[i&1] = -FHUGE;
                } else {
                        gpos = gp[g1] + dv1*dwall/dv0;
                        if (gpos < gcl->gmin[i&1])
                                gcl->gmin[i&1] = gpos;
                        if (gpos > gcl->gmax[i&1])
                                gcl->gmax[i&1] = gpos;
                }
        }
}


static int
clipeyelim(rrng, gcl)           /* clip eye limits to grid cell */
register short  rrng[2][2];
register struct gclim   *gcl;
{
        int     incell = 1;
        register int    i;

        for (i = 0; i < 2; i++) {
                if (gcl->gmin[i] < gcl->gc.i[i])
                        gcl->gmin[i] = gcl->gc.i[i];
                if (gcl->gmax[i] > gcl->gc.i[i]+1)
                        gcl->gmax[i] = gcl->gc.i[i]+1;
                if (gcl->gmax[i] > gcl->gmin[i]) {
                        rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
                                        (1.-FTINY);
                        rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
                                        (1.-FTINY) - rrng[i][0];
                } else
                        rrng[i][0] = rrng[i][1] = 0;
                incell &= rrng[i][1] > 0;
        }
        return(incell);
}


packrays(rod, p)                /* pack ray origins and directions */
register float  *rod;
register PACKET *p;
{
        int     nretries = p->nr + 2;
        struct gclim    eyelim;
        short   rrng0[2][2], rrng1[2][2];
        int     useyelim;
        GCOORD  gc[2];
        FVECT   ro, rd;
        double  d;
        register int    i;

        if (!hdbcoord(gc, hdlist[p->hd], p->bi))
                error(CONSISTENCY, "bad beam index in packrays");
        if ((useyelim = myeye.rng > FTINY)) {
                initeyelim(&eyelim, hdlist[p->hd], gc);
                groweyelim(&eyelim, gc+1, 0., 0.);
                groweyelim(&eyelim, gc+1, 1., 1.);
                useyelim &= clipeyelim(rrng0, &eyelim);
        }
        for (i = 0; i < p->nr; i++) {
        retry:
                if (useyelim) {
                        initeyelim(&eyelim, NULL, gc+1);
                        p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
                                                + rrng0[0][0];
                        p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
                                                + rrng0[1][0];
                        groweyelim(&eyelim, gc,
                                        (1./256.)*(p->ra[i].r[0][0]+.5),
                                        (1./256.)*(p->ra[i].r[0][1]+.5));
                        if (!clipeyelim(rrng1, &eyelim)) {
                                useyelim &= nretries-- > 0;
                                goto retry;
                        }
                        p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
                                                + rrng1[0][0];
                        p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
                                                + rrng1[1][0];
                } else {
                        p->ra[i].r[0][0] = frandom() * 256.;
                        p->ra[i].r[0][1] = frandom() * 256.;
                        p->ra[i].r[1][0] = frandom() * 256.;
                        p->ra[i].r[1][1] = frandom() * 256.;
                }
                d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
                if (p->offset != NULL) {
                        if (!vdef(OBSTRUCTIONS))
                                d *= frandom();         /* random offset */
                        VSUM(ro, ro, rd, d);            /* advance ray */
                        p->offset[i] = d;
                }
                VCOPY(rod, ro);
                rod += 3;
                VCOPY(rod, rd);
                rod += 3;
        }
}


donerays(p, rvl)                /* encode finished ray computations */
register PACKET *p;
register float  *rvl;
{
        double  d;
        register int    i;

        for (i = 0; i < p->nr; i++) {
                setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
                d = rvl[3];
                if (p->offset != NULL)
                        d += p->offset[i];
                p->ra[i].d = hdcode(hdlist[p->hd], d);
                rvl += 4;
        }
        p->nc += p->nr;
}


int
done_rtrace()                   /* clean up and close rtrace calculation */
{
        int     status;
                                        /* already closed? */
        if (!nprocs)
                return;
                                        /* flush beam queue */
        done_packets(flush_queue());
                                        /* sync holodeck */
        hdsync(NULL, 1);
                                        /* close rtrace */
        if ((status = end_rtrace()))
                error(WARNING, "bad exit status from rtrace");
        if (vdef(REPORT)) {             /* report time */
                eputs("rtrace process closed\n");
                report(0);
        }
        return(status);                 /* return status */
}


new_rtrace()                    /* restart rtrace calculation */
{
        char    combuf[128];

        if (nprocs > 0)                 /* already running? */
                return;
        starttime = time(NULL);         /* reset start time and counts */
        npacksdone = nraysdone = 0L;
        if (vdef(TIME))                 /* reset end time */
                endtime = starttime + vflt(TIME)*3600. + .5;
        if (vdef(RIF)) {                /* rerun rad to update octree */
                sprintf(combuf, "rad -v 0 -s -w %s", vval(RIF));
                if (system(combuf))
                        error(WARNING, "error running rad");
        }
        if (start_rtrace() < 1)         /* start rtrace */
                error(WARNING, "cannot restart rtrace");
        else if (vdef(REPORT)) {
                eputs("rtrace process restarted\n");
                report(0);
        }
}


getradfile()                    /* run rad and get needed variables */
{
        static short    mvar[] = {OCTREE,EYESEP,-1};
        static char     tf1[] = TEMPLATE;
        char    tf2[64];
        char    combuf[256];
        char    *pippt;
        register int    i;
        register char   *cp;
                                        /* check if rad file specified */
        if (!vdef(RIF))
                return(0);
                                        /* create rad command */
        mktemp(tf1);
        sprintf(tf2, "%s.rif", tf1);
        sprintf(combuf,
                "rad -v 0 -s -e -w %s OPTFILE=%s | egrep '^[ \t]*(NOMATCH",
                        vval(RIF), tf1);
        cp = combuf;
        while (*cp){
                if (*cp == '|') pippt = cp;
                cp++;
        }                               /* match unset variables */
        for (i = 0; mvar[i] >= 0; i++)
                if (!vdef(mvar[i])) {
                        *cp++ = '|';
                        strcpy(cp, vnam(mvar[i]));
                        while (*cp) cp++;
                        pippt = NULL;
                }
        if (pippt != NULL)
                strcpy(pippt, "> /dev/null");   /* nothing to match */
        else
                sprintf(cp, ")[ \t]*=' > %s", tf2);
#ifdef DEBUG
        wputs(combuf); wputs("\n");
#endif
        system(combuf);                         /* ignore exit code */
        if (pippt == NULL) {
                loadvars(tf2);                  /* load variables */
                unlink(tf2);
        }
        rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */
        unlink(tf1);                    /* clean up */
        return(1);
}


report(t)                       /* report progress so far */
time_t  t;
{
        static time_t   seconds2go = 1000000;

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