ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcdraw.c
Revision: 2.16
Committed: Tue Apr 19 01:15:07 2005 UTC (18 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R1, rad4R0, rad3R8, rad3R9
Changes since 2.15: +2 -2 lines
Log Message:
Extensive changes to enable rtrace -oTW option for tracking ray contributions

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: srcdraw.c,v 2.15 2004/10/26 19:04:32 greg Exp $";
3 #endif
4 /*
5 * Draw small sources into image in case we missed them.
6 *
7 * External symbols declared in ray.h
8 */
9
10 #include "copyright.h"
11
12 #include "ray.h"
13 #include "view.h"
14 #include "otypes.h"
15 #include "source.h"
16
17
18 #define CLIP_ABOVE 1
19 #define CLIP_BELOW 2
20 #define CLIP_RIGHT 3
21 #define CLIP_LEFT 4
22
23 #define MAXVERT 10
24
25 typedef struct splist {
26 struct splist *next; /* next source in list */
27 int sn; /* source number */
28 short nv; /* number of vertices */
29 RREAL vl[3][2]; /* vertex array (last) */
30 } SPLIST; /* source polygon list */
31
32 extern VIEW ourview; /* our view parameters */
33 extern int hres, vres; /* our image resolution */
34 static SPLIST *sphead = NULL; /* our list of source polys */
35
36 static int inregion(RREAL p[2], double cv, int crit);
37 static void clipregion(RREAL a[2], RREAL b[2], double cv, int crit, RREAL r[2]);
38 static int hp_clip_poly(RREAL vl[][2], int nv, double cv, int crit,
39 RREAL vlo[][2]);
40 static int box_clip_poly(RREAL vl[MAXVERT][2], int nv,
41 double xl, double xr, double yb, double ya, RREAL vlo[MAXVERT][2]);
42 static double minw2(RREAL vl[][2], int nv, double ar2);
43 static void convex_center(RREAL vl[][2], int nv, RREAL cv[2]);
44 static double poly_area(RREAL vl[][2], int nv);
45 static int convex_hull(RREAL vl[][2], int nv, RREAL vlo[][2]);
46 static void spinsert(int sn, RREAL vl[][2], int nv);
47 static int sourcepoly(int sn, RREAL sp[MAXVERT][2]);
48
49
50 static int
51 inregion( /* check if vertex is in region */
52 RREAL p[2],
53 double cv,
54 int crit
55 )
56 {
57 switch (crit) {
58 case CLIP_ABOVE:
59 return(p[1] < cv);
60 case CLIP_BELOW:
61 return(p[1] >= cv);
62 case CLIP_RIGHT:
63 return(p[0] < cv);
64 case CLIP_LEFT:
65 return(p[0] >= cv);
66 }
67 return(-1);
68 }
69
70
71 static void
72 clipregion( /* find intersection with boundary */
73 register RREAL a[2],
74 register RREAL b[2],
75 double cv,
76 int crit,
77 RREAL r[2] /* return value */
78 )
79 {
80 switch (crit) {
81 case CLIP_ABOVE:
82 case CLIP_BELOW:
83 r[1] = cv;
84 r[0] = a[0] + (cv-a[1])/(b[1]-a[1])*(b[0]-a[0]);
85 return;
86 case CLIP_RIGHT:
87 case CLIP_LEFT:
88 r[0] = cv;
89 r[1] = a[1] + (cv-a[0])/(b[0]-a[0])*(b[1]-a[1]);
90 return;
91 }
92 }
93
94
95 static int
96 hp_clip_poly( /* clip polygon to half-plane */
97 RREAL vl[][2],
98 int nv,
99 double cv,
100 int crit,
101 RREAL vlo[][2] /* return value */
102 )
103 {
104 RREAL *s, *p;
105 register int j, nvo;
106
107 s = vl[nv-1];
108 nvo = 0;
109 for (j = 0; j < nv; j++) {
110 p = vl[j];
111 if (inregion(p, cv, crit)) {
112 if (!inregion(s, cv, crit))
113 clipregion(s, p, cv, crit, vlo[nvo++]);
114 vlo[nvo][0] = p[0]; vlo[nvo++][1] = p[1];
115 } else if (inregion(s, cv, crit))
116 clipregion(s, p, cv, crit, vlo[nvo++]);
117 s = p;
118 }
119 return(nvo);
120 }
121
122
123 static int
124 box_clip_poly( /* clip polygon to box */
125 RREAL vl[MAXVERT][2],
126 int nv,
127 double xl,
128 double xr,
129 double yb,
130 double ya,
131 RREAL vlo[MAXVERT][2] /* return value */
132 )
133 {
134 RREAL vlt[MAXVERT][2];
135 int nvt, nvo;
136
137 nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
138 nvo = hp_clip_poly(vlt, nvt, ya, CLIP_ABOVE, vlo);
139 nvt = hp_clip_poly(vlo, nvo, xl, CLIP_LEFT, vlt);
140 nvo = hp_clip_poly(vlt, nvt, xr, CLIP_RIGHT, vlo);
141
142 return(nvo);
143 }
144
145
146 static double
147 minw2( /* compute square of minimum width */
148 RREAL vl[][2],
149 int nv,
150 double ar2
151 )
152 {
153 double d2, w2, w2min, w2max;
154 register RREAL *p0, *p1, *p2;
155 int i, j;
156 /* find minimum for all widths */
157 w2min = FHUGE;
158 p0 = vl[nv-1];
159 for (i = 0; i < nv; i++) { /* for each edge */
160 p1 = vl[i];
161 d2 = (p1[0]-p0[0])*(p1[0]-p0[0]) +
162 (p1[1]-p0[1])*(p1[1]-p0[1])*ar2;
163 w2max = 0.; /* find maximum for this side */
164 for (j = 1; j < nv-1; j++) {
165 p2 = vl[(i+j)%nv];
166 w2 = (p1[0]-p0[0])*(p2[1]-p0[1]) -
167 (p1[1]-p0[1])*(p2[0]-p0[0]);
168 w2 = w2*w2*ar2/d2; /* triangle height squared */
169 if (w2 > w2max)
170 w2max = w2;
171 }
172 if (w2max < w2min) /* global min. based on local max.'s */
173 w2min = w2max;
174 p0 = p1;
175 }
176 return(w2min);
177 }
178
179
180 static void
181 convex_center( /* compute center of convex polygon */
182 register RREAL vl[][2],
183 int nv,
184 RREAL cv[2] /* return value */
185 )
186 {
187 register int i;
188 /* simple average (suboptimal) */
189 cv[0] = cv[1] = 0.;
190 for (i = 0; i < nv; i++) {
191 cv[0] += vl[i][0];
192 cv[1] += vl[i][1];
193 }
194 cv[0] /= (double)nv;
195 cv[1] /= (double)nv;
196 }
197
198
199 static double
200 poly_area( /* compute area of polygon */
201 register RREAL vl[][2],
202 int nv
203 )
204 {
205 double a;
206 RREAL v0[2], v1[2];
207 register int i;
208
209 a = 0.;
210 v0[0] = vl[1][0] - vl[0][0];
211 v0[1] = vl[1][1] - vl[0][1];
212 for (i = 2; i < nv; i++) {
213 v1[0] = vl[i][0] - vl[0][0];
214 v1[1] = vl[i][1] - vl[0][1];
215 a += v0[0]*v1[1] - v0[1]*v1[0];
216 v0[0] = v1[0]; v0[1] = v1[1];
217 }
218 return(a * (a >= 0. ? .5 : -.5));
219 }
220
221
222 static int
223 convex_hull( /* compute polygon's convex hull */
224 RREAL vl[][2],
225 int nv,
226 RREAL vlo[][2] /* return value */
227 )
228 {
229 int nvo, nvt;
230 RREAL vlt[MAXVERT][2];
231 double voa, vta;
232 register int i, j;
233 /* start with original polygon */
234 for (i = nvo = nv; i--; ) {
235 vlo[i][0] = vl[i][0]; vlo[i][1] = vl[i][1];
236 }
237 voa = poly_area(vlo, nvo); /* compute its area */
238 for (i = 0; i < nvo; i++) { /* for each output vertex */
239 for (j = 0; j < i; j++) {
240 vlt[j][0] = vlo[j][0]; vlt[j][1] = vlo[j][1];
241 }
242 nvt = nvo - 1; /* make poly w/o vertex */
243 for (j = i; j < nvt; j++) {
244 vlt[j][0] = vlo[j+1][0]; vlt[j][1] = vlo[j+1][1];
245 }
246 vta = poly_area(vlt, nvt);
247 if (vta >= voa) { /* is simpler poly bigger? */
248 voa = vta; /* then use it */
249 for (j = nvo = nvt; j--; ) {
250 vlo[j][0] = vlt[j][0]; vlo[j][1] = vlt[j][1];
251 }
252 i--; /* next adjust */
253 }
254 }
255 return(nvo);
256 }
257
258
259 static void
260 spinsert( /* insert new source polygon */
261 int sn,
262 RREAL vl[][2],
263 int nv
264 )
265 {
266 register SPLIST *spn;
267 register int i;
268
269 if (nv < 3)
270 return;
271 if (nv > 3)
272 spn = (SPLIST *)malloc(sizeof(SPLIST)+sizeof(RREAL)*2*(nv-3));
273 else
274 spn = (SPLIST *)malloc(sizeof(SPLIST));
275 if (spn == NULL)
276 error(SYSTEM, "out of memory in spinsert");
277 spn->sn = sn;
278 for (i = spn->nv = nv; i--; ) {
279 spn->vl[i][0] = vl[i][0]; spn->vl[i][1] = vl[i][1];
280 }
281 spn->next = sphead; /* push onto global list */
282 sphead = spn;
283 }
284
285
286 static int
287 sourcepoly( /* compute image polygon for source */
288 int sn,
289 RREAL sp[MAXVERT][2]
290 )
291 {
292 static short cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
293 {0,1,3,7,6,4},{0,1,5,7,6,2},
294 {0,2,6,7,5,1},{0,4,6,7,3,1},
295 {0,2,3,7,5,4},{1,5,4,6,2,3}};
296 register SRCREC *s = source + sn;
297 FVECT ap, ip;
298 RREAL pt[6][2];
299 int dir;
300 register int i, j;
301
302 if (s->sflags & (SDISTANT|SFLAT)) {
303 if (s->sflags & SDISTANT) {
304 if (ourview.type == VT_PAR)
305 return(0); /* all or nothing case */
306 if (s->srad >= 0.05)
307 return(0); /* should never be a problem */
308 }
309 if (s->sflags & SFLAT) {
310 for (i = 0; i < 3; i++)
311 ap[i] = s->sloc[i] - ourview.vp[i];
312 if (DOT(ap, s->snorm) >= 0.)
313 return(0); /* source faces away */
314 }
315 for (j = 0; j < 4; j++) { /* four corners */
316 for (i = 0; i < 3; i++) {
317 ap[i] = s->sloc[i];
318 if ((j==1)|(j==2)) ap[i] += s->ss[SU][i];
319 else ap[i] -= s->ss[SU][i];
320 if ((j==2)|(j==3)) ap[i] += s->ss[SV][i];
321 else ap[i] -= s->ss[SV][i];
322 if (s->sflags & SDISTANT) {
323 ap[i] *= 1. + ourview.vfore;
324 ap[i] += ourview.vp[i];
325 }
326 }
327 viewloc(ip, &ourview, ap); /* find image point */
328 if (ip[2] <= 0.)
329 return(0); /* in front of view */
330 sp[j][0] = ip[0]; sp[j][1] = ip[1];
331 }
332 return(4);
333 }
334 /* identify furthest corner */
335 for (i = 0; i < 3; i++)
336 ap[i] = s->sloc[i] - ourview.vp[i];
337 dir = (DOT(ap,s->ss[SU])>0.) |
338 (DOT(ap,s->ss[SV])>0.)<<1 |
339 (DOT(ap,s->ss[SW])>0.)<<2 ;
340 /* order vertices based on this */
341 for (j = 0; j < 6; j++) {
342 for (i = 0; i < 3; i++) {
343 ap[i] = s->sloc[i];
344 if (cubeord[dir][j] & 1) ap[i] += s->ss[SU][i];
345 else ap[i] -= s->ss[SU][i];
346 if (cubeord[dir][j] & 2) ap[i] += s->ss[SV][i];
347 else ap[i] -= s->ss[SV][i];
348 if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
349 else ap[i] -= s->ss[SW][i];
350 }
351 viewloc(ip, &ourview, ap); /* find image point */
352 if (ip[2] <= 0.)
353 return(0); /* in front of view */
354 pt[j][0] = ip[0]; pt[j][1] = ip[1];
355 }
356 return(convex_hull(pt, 6, sp)); /* make sure it's convex */
357 }
358
359
360 /* initialize by finding sources smaller than rad */
361 extern void
362 init_drawsources(
363 int rad /* source sample size */
364 )
365 {
366 RREAL spoly[MAXVERT][2];
367 int nsv;
368 register SPLIST *sp;
369 register int i;
370 /* free old source list if one */
371 for (sp = sphead; sp != NULL; sp = sphead) {
372 sphead = sp->next;
373 free((void *)sp);
374 }
375 /* loop through all sources */
376 for (i = nsources; i--; ) {
377 /* skip illum's */
378 if (findmaterial(source[i].so)->otype == MAT_ILLUM)
379 continue;
380 /* compute image polygon for source */
381 if (!(nsv = sourcepoly(i, spoly)))
382 continue;
383 /* clip to image boundaries */
384 if (!(nsv = box_clip_poly(spoly, nsv, 0., 1., 0., 1., spoly)))
385 continue;
386 /* big enough for standard sampling? */
387 if (minw2(spoly, nsv, ourview.vn2/ourview.hn2) >
388 (double)rad*rad/hres/hres)
389 continue;
390 /* OK, add to our list */
391 spinsert(i, spoly, nsv);
392 }
393 }
394
395 extern void /* add sources smaller than rad to computed subimage */
396 drawsources(
397 COLOR *pic[], /* subimage pixel value array */
398 float *zbf[], /* subimage distance array (opt.) */
399 int x0, /* origin and size of subimage */
400 int xsiz,
401 int y0,
402 int ysiz
403 )
404 {
405 RREAL spoly[MAXVERT][2], ppoly[MAXVERT][2];
406 int nsv, npv;
407 int xmin, xmax, ymin, ymax, x, y;
408 RREAL cxy[2];
409 double w;
410 RAY sr;
411 register SPLIST *sp;
412 register int i;
413 /* check each source in our list */
414 for (sp = sphead; sp != NULL; sp = sp->next) {
415 /* clip source poly to subimage */
416 nsv = box_clip_poly(sp->vl, sp->nv,
417 (double)x0/hres, (double)(x0+xsiz)/hres,
418 (double)y0/vres, (double)(y0+ysiz)/vres, spoly);
419 if (!nsv)
420 continue;
421 /* find common subimage (BBox) */
422 xmin = x0 + xsiz; xmax = x0;
423 ymin = y0 + ysiz; ymax = y0;
424 for (i = 0; i < nsv; i++) {
425 if ((double)xmin/hres > spoly[i][0])
426 xmin = spoly[i][0]*hres + FTINY;
427 if ((double)xmax/hres < spoly[i][0])
428 xmax = spoly[i][0]*hres - FTINY;
429 if ((double)ymin/vres > spoly[i][1])
430 ymin = spoly[i][1]*vres + FTINY;
431 if ((double)ymax/vres < spoly[i][1])
432 ymax = spoly[i][1]*vres - FTINY;
433 }
434 /* evaluate each pixel in BBox */
435 for (y = ymin; y <= ymax; y++)
436 for (x = xmin; x <= xmax; x++) {
437 /* subarea for pixel */
438 npv = box_clip_poly(spoly, nsv,
439 (double)x/hres, (x+1.)/hres,
440 (double)y/vres, (y+1.)/vres,
441 ppoly);
442 if (!npv)
443 continue; /* no overlap */
444 convex_center(ppoly, npv, cxy);
445 if ((sr.rmax = viewray(sr.rorg,sr.rdir,&ourview,
446 cxy[0],cxy[1])) < -FTINY)
447 continue; /* not in view */
448 if (source[sp->sn].sflags & SSPOT &&
449 spotout(&sr, source[sp->sn].sl.s))
450 continue; /* outside spot */
451 rayorigin(&sr, SHADOW, NULL, NULL);
452 sr.rsrc = sp->sn;
453 rayvalue(&sr); /* compute value */
454 if (bright(sr.rcol) <= FTINY)
455 continue; /* missed/blocked */
456 /* modify pixel */
457 w = poly_area(ppoly, npv) * hres * vres;
458 if (zbf[y-y0] != NULL &&
459 sr.rt < 0.99*zbf[y-y0][x-x0]) {
460 zbf[y-y0][x-x0] = sr.rt;
461 } else if (!bigdiff(sr.rcol, pic[y-y0][x-x0],
462 0.01)) { /* source sample */
463 scalecolor(pic[y-y0][x-x0], w);
464 continue;
465 }
466 scalecolor(sr.rcol, w);
467 scalecolor(pic[y-y0][x-x0], 1.-w);
468 addcolor(pic[y-y0][x-x0], sr.rcol);
469 }
470 }
471 }