ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhdisp3.c
Revision: 3.16
Committed: Tue Jun 8 19:48:30 2004 UTC (19 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1
Changes since 3.15: +1 -2 lines
Log Message:
Removed redundant #include's and fixed ordering on some headers

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rhdisp3.c,v 3.15 2004/01/01 11:21:55 schorsch Exp $";
3 #endif
4 /*
5 * Holodeck beam support for display process
6 */
7
8 #include "rholo.h"
9 #include "rhdisp.h"
10
11 struct cellist {
12 GCOORD *cl;
13 int n;
14 };
15
16
17 int
18 npixels(vp, hr, vr, hp, bi) /* compute appropriate nrays to evaluate */
19 register VIEW *vp;
20 int hr, vr;
21 HOLO *hp;
22 int bi;
23 {
24 VIEW vrev;
25 GCOORD gc[2];
26 FVECT cp[4], ip[4], pf, pb;
27 double af, ab, sf2, sb2, dfb2, df2, db2, penalty;
28 register int i;
29 /* special case */
30 if (hr <= 0 | vr <= 0)
31 return(0);
32 /* compute cell corners in image */
33 if (!hdbcoord(gc, hp, bi))
34 error(CONSISTENCY, "bad beam index in npixels");
35 hdcell(cp, hp, gc+1); /* find cell on front image */
36 for (i = 3; i--; ) /* compute front center */
37 pf[i] = 0.5*(cp[0][i] + cp[2][i]);
38 sf2 = 0.25*dist2(cp[0], cp[2]); /* compute half diagonal length */
39 for (i = 0; i < 4; i++) { /* compute visible quad */
40 viewloc(ip[i], vp, cp[i]);
41 if (ip[i][2] < 0.) {
42 af = 0;
43 goto getback;
44 }
45 ip[i][0] *= (double)hr; /* scale by resolution */
46 ip[i][1] *= (double)vr;
47 }
48 /* compute front area */
49 af = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
50 (ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
51 af += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
52 (ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
53 af *= af >= 0 ? 0.5 : -0.5;
54 getback:
55 vrev = *vp; /* compute reverse view */
56 for (i = 0; i < 3; i++) {
57 vrev.vdir[i] = -vp->vdir[i];
58 vrev.vup[i] = -vp->vup[i];
59 vrev.hvec[i] = -vp->hvec[i];
60 vrev.vvec[i] = -vp->vvec[i];
61 }
62 hdcell(cp, hp, gc); /* find cell on back image */
63 for (i = 3; i--; ) /* compute rear center */
64 pb[i] = 0.5*(cp[0][i] + cp[2][i]);
65 sb2 = 0.25*dist2(cp[0], cp[2]); /* compute half diagonal length */
66 for (i = 0; i < 4; i++) { /* compute visible quad */
67 viewloc(ip[i], &vrev, cp[i]);
68 if (ip[i][2] < 0.) {
69 ab = 0;
70 goto finish;
71 }
72 ip[i][0] *= (double)hr; /* scale by resolution */
73 ip[i][1] *= (double)vr;
74 }
75 /* compute back area */
76 ab = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
77 (ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
78 ab += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
79 (ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
80 ab *= ab >= 0 ? 0.5 : -0.5;
81 finish: /* compute penalty based on dist. sightline - viewpoint */
82 df2 = dist2(vp->vp, pf);
83 db2 = dist2(vp->vp, pb);
84 dfb2 = dist2(pf, pb);
85 penalty = dfb2 + df2 - db2;
86 penalty = df2 - 0.25*penalty*penalty/dfb2;
87 if (df2 > db2) penalty /= df2 <= dfb2 ? sb2 : sb2*df2/dfb2;
88 else penalty /= db2 <= dfb2 ? sf2 : sf2*db2/dfb2;
89 if (penalty < 1.) penalty = 1.;
90 /* round off smaller non-zero area */
91 if (ab <= FTINY || (af > FTINY && af <= ab))
92 return((int)(af/penalty + 0.5));
93 return((int)(ab/penalty + 0.5));
94 }
95
96
97 /*
98 * The ray directions that define the pyramid in visit_cells() needn't
99 * be normalized, but they must be given in clockwise order as seen
100 * from the pyramid's apex (origin).
101 * If no cell centers fall within the domain, the closest cell is visited.
102 */
103 int
104 visit_cells(orig, pyrd, hp, vf, dp) /* visit cells within a pyramid */
105 FVECT orig, pyrd[4]; /* pyramid ray directions in clockwise order */
106 register HOLO *hp;
107 int (*vf)();
108 char *dp;
109 {
110 int ncalls = 0, n = 0;
111 int inflags = 0;
112 FVECT gp, pn[4], lo, ld;
113 double po[4], lbeg, lend, d, t;
114 GCOORD gc, gc2[2];
115 register int i;
116 /* figure out whose side we're on */
117 hdgrid(gp, hp, orig);
118 for (i = 0; i < 3; i++) {
119 inflags |= (gp[i] > FTINY) << (i<<1);
120 inflags |= (gp[i] < hp->grid[i]-FTINY) << (i<<1 | 1);
121 }
122 /* compute pyramid planes */
123 for (i = 0; i < 4; i++) {
124 fcross(pn[i], pyrd[i], pyrd[(i+1)&03]);
125 po[i] = DOT(pn[i], orig);
126 }
127 /* traverse each wall */
128 for (gc.w = 0; gc.w < 6; gc.w++) {
129 if (!(inflags & 1<<gc.w)) /* origin on wrong side */
130 continue;
131 /* scanline algorithm */
132 for (gc.i[1] = hp->grid[hdwg1[gc.w]]; gc.i[1]--; ) {
133 /* compute scanline */
134 gp[gc.w>>1] = gc.w&1 ? hp->grid[gc.w>>1] : 0;
135 gp[hdwg0[gc.w]] = 0;
136 gp[hdwg1[gc.w]] = gc.i[1] + 0.5;
137 hdworld(lo, hp, gp);
138 gp[hdwg0[gc.w]] = 1;
139 hdworld(ld, hp, gp);
140 ld[0] -= lo[0]; ld[1] -= lo[1]; ld[2] -= lo[2];
141 /* find scanline limits */
142 lbeg = 0; lend = hp->grid[hdwg0[gc.w]];
143 for (i = 0; i < 4; i++) {
144 t = DOT(pn[i], lo) - po[i];
145 d = -DOT(pn[i], ld);
146 if (d > FTINY) { /* <- plane */
147 if ((t /= d) < lend)
148 lend = t;
149 } else if (d < -FTINY) { /* plane -> */
150 if ((t /= d) > lbeg)
151 lbeg = t;
152 } else if (t < 0) { /* outside */
153 lend = -1;
154 break;
155 }
156 }
157 if (lbeg >= lend)
158 continue;
159 i = lend + .5; /* visit cells on this scan */
160 for (gc.i[0] = lbeg + .5; gc.i[0] < i; gc.i[0]++) {
161 n += (*vf)(&gc, dp);
162 ncalls++;
163 }
164 }
165 }
166 if (ncalls) /* got one at least */
167 return(n);
168 /* else find closest cell */
169 VSUM(ld, pyrd[0], pyrd[1], 1.);
170 VSUM(ld, ld, pyrd[2], 1.);
171 VSUM(ld, ld, pyrd[3], 1.);
172 #if 0
173 if (normalize(ld) == 0.0) /* technically not necessary */
174 return(0);
175 #endif
176 d = hdinter(gc2, NULL, &t, hp, orig, ld);
177 if (d >= FHUGE || t <= 0.)
178 return(0);
179 return((*vf)(gc2+1, dp)); /* visit it */
180 }
181
182
183 sect_behind(hp, vp) /* check if section is "behind" viewpoint */
184 register HOLO *hp;
185 register VIEW *vp;
186 {
187 FVECT hcent;
188 /* compute holodeck section center */
189 VSUM(hcent, hp->orig, hp->xv[0], 0.5);
190 VSUM(hcent, hcent, hp->xv[1], 0.5);
191 VSUM(hcent, hcent, hp->xv[2], 0.5);
192 /* behind if center is behind */
193 return(DOT(vp->vdir,hcent) < DOT(vp->vdir,vp->vp));
194 }
195
196
197 viewpyramid(org, dir, hp, vp) /* compute view pyramid */
198 FVECT org, dir[4];
199 HOLO *hp;
200 VIEW *vp;
201 {
202 register int i;
203 /* check view type */
204 if (vp->type == VT_PAR)
205 return(0);
206 /* in front or behind? */
207 if (!sect_behind(hp, vp)) {
208 if (viewray(org, dir[0], vp, 0., 0.) < -FTINY)
209 return(0);
210 if (viewray(org, dir[1], vp, 0., 1.) < -FTINY)
211 return(0);
212 if (viewray(org, dir[2], vp, 1., 1.) < -FTINY)
213 return(0);
214 if (viewray(org, dir[3], vp, 1., 0.) < -FTINY)
215 return(0);
216 return(1);
217 } /* reverse pyramid */
218 if (viewray(org, dir[3], vp, 0., 0.) < -FTINY)
219 return(0);
220 if (viewray(org, dir[2], vp, 0., 1.) < -FTINY)
221 return(0);
222 if (viewray(org, dir[1], vp, 1., 1.) < -FTINY)
223 return(0);
224 if (viewray(org, dir[0], vp, 1., 0.) < -FTINY)
225 return(0);
226 for (i = 0; i < 3; i++) {
227 dir[0][i] = -dir[0][i];
228 dir[1][i] = -dir[1][i];
229 dir[2][i] = -dir[2][i];
230 dir[3][i] = -dir[3][i];
231 }
232 return(-1);
233 }
234
235
236 int
237 addcell(gcp, cl) /* add a cell to a list */
238 GCOORD *gcp;
239 register struct cellist *cl;
240 {
241 *(cl->cl+cl->n) = *gcp;
242 cl->n++;
243 return(1);
244 }
245
246
247 int
248 cellcmp(gcp1, gcp2) /* visit_cells() cell ordering */
249 register GCOORD *gcp1, *gcp2;
250 {
251 register int c;
252
253 if ((c = gcp1->w - gcp2->w))
254 return(c);
255 if ((c = gcp2->i[1] - gcp1->i[1])) /* wg1 is reverse-ordered */
256 return(c);
257 return(gcp1->i[0] - gcp2->i[0]);
258 }
259
260
261 GCOORD *
262 getviewcells(np, hp, vp) /* get ordered cell list for section view */
263 int *np; /* returned number of cells (negative if reversed) */
264 register HOLO *hp;
265 VIEW *vp;
266 {
267 FVECT org, dir[4];
268 int orient;
269 struct cellist cl;
270 /* compute view pyramid */
271 *np = 0;
272 orient = viewpyramid(org, dir, hp, vp);
273 if (!orient)
274 return(NULL);
275 /* allocate enough list space */
276 cl.n = 2*( hp->grid[0]*hp->grid[1] +
277 hp->grid[0]*hp->grid[2] +
278 hp->grid[1]*hp->grid[2] );
279 cl.cl = (GCOORD *)malloc(cl.n*sizeof(GCOORD));
280 if (cl.cl == NULL)
281 goto memerr;
282 cl.n = 0; /* add cells within pyramid */
283 visit_cells(org, dir, hp, addcell, (char *)&cl);
284 if (!cl.n) {
285 free((void *)cl.cl);
286 return(NULL);
287 }
288 *np = cl.n * orient;
289 #if 0
290 /* We're just going to free this memory in a moment, and list is
291 * sorted automatically by visit_cells(), so we don't need this.
292 */
293 /* optimize memory use */
294 cl.cl = (GCOORD *)realloc((void *)cl.cl, cl.n*sizeof(GCOORD));
295 if (cl.cl == NULL)
296 goto memerr;
297 /* sort the list */
298 qsort((char *)cl.cl, cl.n, sizeof(GCOORD), cellcmp);
299 #endif
300 return(cl.cl);
301 memerr:
302 error(SYSTEM, "out of memory in getviewcells");
303 }
304
305
306 extern void
307 gridlines( /* run through holodeck section grid lines */
308 void (*f)(FVECT wp[2])
309 )
310 {
311 register int hd, w, i;
312 int g0, g1;
313 FVECT wp[2], mov;
314 double d;
315 /* do each wall on each section */
316 for (hd = 0; hdlist[hd] != NULL; hd++)
317 for (w = 0; w < 6; w++) {
318 g0 = hdwg0[w];
319 g1 = hdwg1[w];
320 d = 1.0/hdlist[hd]->grid[g0];
321 mov[0] = d * hdlist[hd]->xv[g0][0];
322 mov[1] = d * hdlist[hd]->xv[g0][1];
323 mov[2] = d * hdlist[hd]->xv[g0][2];
324 if (w & 1) {
325 VSUM(wp[0], hdlist[hd]->orig,
326 hdlist[hd]->xv[w>>1], 1.);
327 VSUM(wp[0], wp[0], mov, 1.);
328 } else
329 VCOPY(wp[0], hdlist[hd]->orig);
330 VSUM(wp[1], wp[0], hdlist[hd]->xv[g1], 1.);
331 for (i = hdlist[hd]->grid[g0]; ; ) { /* g0 lines */
332 (*f)(wp);
333 if (!--i) break;
334 wp[0][0] += mov[0]; wp[0][1] += mov[1];
335 wp[0][2] += mov[2]; wp[1][0] += mov[0];
336 wp[1][1] += mov[1]; wp[1][2] += mov[2];
337 }
338 d = 1.0/hdlist[hd]->grid[g1];
339 mov[0] = d * hdlist[hd]->xv[g1][0];
340 mov[1] = d * hdlist[hd]->xv[g1][1];
341 mov[2] = d * hdlist[hd]->xv[g1][2];
342 if (w & 1)
343 VSUM(wp[0], hdlist[hd]->orig,
344 hdlist[hd]->xv[w>>1], 1.);
345 else
346 VSUM(wp[0], hdlist[hd]->orig, mov, 1.);
347 VSUM(wp[1], wp[0], hdlist[hd]->xv[g0], 1.);
348 for (i = hdlist[hd]->grid[g1]; ; ) { /* g1 lines */
349 (*f)(wp);
350 if (!--i) break;
351 wp[0][0] += mov[0]; wp[0][1] += mov[1];
352 wp[0][2] += mov[2]; wp[1][0] += mov[0];
353 wp[1][1] += mov[1]; wp[1][2] += mov[2];
354 }
355 }
356 }