ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhdisp3.c
Revision: 3.12
Committed: Sat Feb 22 02:07:24 2003 UTC (21 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 3.11: +2 -5 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

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