ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/hd/rholo2.c

Comparing ray/src/hd/rholo2.c (file contents):
Revision 3.19 by gwlarson, Tue Dec 1 15:47:05 1998 UTC vs.
Revision 3.23 by schorsch, Thu Jun 26 00:58:10 2003 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1998 Silicon Graphics, Inc. */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ SGI";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Rtrace support routines for holodeck rendering
6		*/
#	Line 19 | Line 16 | struct gclim {
16		HOLO    *hp;                    /* holodeck pointer */
17		GCOORD  gc;                     /* grid cell */
18		FVECT   egp;                    /* eye grid point */
19	<	double  erg[2];                 /* eye range in wall grid coords */
19	>	double  erg2;                   /* mean square eye grid range */
20		double  gmin[2], gmax[2];       /* grid coordinate limits */
21		};                              /* a grid coordinate range */
22
23
24		static
25	<	initeyelim(gcl, hd, gc)         /* initialize grid coordinate limits */
25	>	initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
26		register struct gclim   *gcl;
27	<	int     hd;
27	>	register HOLO   *hp;
28		GCOORD  *gc;
29		{
30	<	register FLOAT  *v;
30	>	register RREAL  *v;
31		register int    i;
32
33	<	gcl->hp = hdlist[hd];
34	<	copystruct(&gcl->gc, gc);
35	<	hdgrid(gcl->egp, gcl->hp, myeye.vpt);
36	<	for (i = 0; i < 2; i++) {
37	<	v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
38	<	gcl->erg[i] = myeye.rng * VLEN(v);
42	<	gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
33	>	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		}
40	+	if (gc != NULL)
41	+	copystruct(&gcl->gc, gc);
42	+	gcl->gmin[0] = gcl->gmin[1] = FHUGE;
43	+	gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
44		}
45
46
47		static
48	<	groweyelim(gcl, gp)             /* grow grid limits about eye point */
48	>	groweyelim(gcl, gc, r0, r1, tight)      /* grow grid limits about eye point */
49		register struct gclim   *gcl;
50	<	FVECT   gp;
50	>	GCOORD  *gc;
51	>	double  r0, r1;
52	>	int     tight;
53		{
54	<	FVECT   ab;
55	<	double  l2, d, mult, wg;
56	<	register int    i, g;
57	<
54	>	FVECT   gp, ab;
55	>	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	>	i = gc->w>>1;
63	>	gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
64	>	gp[hdwg0[gc->w]] = gc->i[0] + r0;
65	>	gp[hdwg1[gc->w]] = gc->i[1] + r1;
66		VSUB(ab, gcl->egp, gp);
67	<	l2 = DOT(ab,ab);
68	<	if (l2 <= gcl->erg[0]*gcl->erg[1]) {
67	>	ab2 = DOT(ab, ab);
68	>	gw = gcl->gc.w>>1;
69	>	if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
70		gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
71		gcl->gmax[0] = gcl->gmax[1] = FHUGE;
72	<	return;
72	>	return;                 /* too close (to wall) */
73		}
74	<	mult = gp[g = gcl->gc.w>>1];
75	<	if (gcl->gc.w&1)
76	<	mult -= gcl->hp->grid[g];
77	<	if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
78	<	mult /= -ab[g];
79	<	else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
80	<	mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
81	<	ab[hdwg0[gcl->gc.w]];
82	<	else
83	<	mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
84	<	ab[hdwg1[gcl->gc.w]];
74	>	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	<	g = ((gcl->gc.w>>1)+i+1)%3;
100	<	wg = gp[g] + mult*ab[g];
101	<	d = mult*gcl->erg[i];
102	<	if (d < 0.) d = -d;
103	<	if (wg - d < gcl->gmin[i])
104	<	gcl->gmin[i] = wg - d;
105	<	if (wg + d > gcl->gmax[i])
106	<	gcl->gmax[i] = wg + d;
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	>	}
104	>	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		}
144		}
145
#	Line 97 | Line 157 | register struct gclim  *gcl;
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	<	if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
160	>	if (gcl->gmax[i] > gcl->gmin[i]) {
161		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	<	incell &= rrng[i][1] > 0;
166	<	}
165	>	} else
166	>	rrng[i][0] = rrng[i][1] = 0;
167	>	incell &= rrng[i][1] > 0;
168		}
169		return(incell);
170		}
#	Line 113 | Line 174 | packrays(rod, p)               /* pack ray origins and directions *
174		register float  *rod;
175		register PACKET *p;
176		{
177	<	#define gp      ro
178	<	#ifdef DEBUG
118	<	double dist2sum = 0.;
177	>	#if 0
178	>	double  dist2sum = 0.;
179		FVECT   vt;
180		#endif
181		int     nretries = p->nr + 2;
#	Line 130 | Line 190 | register PACKET        *p;
190		if (!hdbcoord(gc, hdlist[p->hd], p->bi))
191		error(CONSISTENCY, "bad beam index in packrays");
192		if ((useyelim = myeye.rng > FTINY)) {
193	<	initeyelim(&eyelim, p->hd, gc);
194	<	gp[gc[1].w>>1] = gc[1].w&1 ?
195	<	hdlist[p->hd]->grid[gc[1].w>>1] : 0;
196	<	gp[hdwg0[gc[1].w]] = gc[1].i[0];
197	<	gp[hdwg1[gc[1].w]] = gc[1].i[1];
198	<	groweyelim(&eyelim, gp);
199	<	gp[hdwg0[gc[1].w]]++;
200	<	gp[hdwg1[gc[1].w]]++;
141	<	groweyelim(&eyelim, gp);
142	<	useyelim &= clipeyelim(rrng0, &eyelim);
193	>	initeyelim(&eyelim, hdlist[p->hd], gc);
194	>	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		}
202		for (i = 0; i < p->nr; i++) {
203		retry:
204		if (useyelim) {
205	+	initeyelim(&eyelim, NULL, gc+1);
206		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	<	initeyelim(&eyelim, p->hd, gc+1);
211	<	gp[gc[0].w>>1] = gc[0].w&1 ?
212	<	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);
210	>	groweyelim(&eyelim, gc,
211	>	(1./256.)*(p->ra[i].r[0][0]+.5),
212	>	(1./256.)*(p->ra[i].r[0][1]+.5), 1);
213		if (!clipeyelim(rrng1, &eyelim)) {
214	<	useyelim &= nretries-- > 0;
214	>	useyelim = nretries-- > 0;
215		#ifdef DEBUG
216		if (!useyelim)
217	<	error(WARNING, "exceeded retry limit in packrays");
217	>	error(WARNING,
218	>	"exceeded retry limit in packrays");
219		#endif
220		goto retry;
221		}
#	Line 175 | Line 230 | register PACKET        *p;
230		p->ra[i].r[1][1] = frandom() * 256.;
231		}
232		d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
233	<	#ifdef DEBUG
233	>	#if 0
234		VSUM(vt, ro, rd, d);
235		dist2sum += dist2line(myeye.vpt, ro, vt);
236		#endif
#	Line 190 | Line 245 | register PACKET        *p;
245		VCOPY(rod, rd);
246		rod += 3;
247		}
248	<	#ifdef DEBUG
249	<	fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
248	>	#if 0
249	>	fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
250	>	p->nr + 2 - nretries);
251		#endif
196	–	#undef gp
252		}
253
254

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)