ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcdraw.c
Revision: 2.5
Committed: Sat Feb 22 02:07:29 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +62 -6 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 * Draw small sources into image in case we missed them.
6 *
7 * External symbols declared in ray.h
8 */
9
10 /* ====================================================================
11 * The Radiance Software License, Version 1.0
12 *
13 * Copyright (c) 1990 - 2002 The Regents of the University of California,
14 * through Lawrence Berkeley National Laboratory. All rights reserved.
15 *
16 * Redistribution and use in source and binary forms, with or without
17 * modification, are permitted provided that the following conditions
18 * are met:
19 *
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 *
23 * 2. Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions and the following disclaimer in
25 * the documentation and/or other materials provided with the
26 * distribution.
27 *
28 * 3. The end-user documentation included with the redistribution,
29 * if any, must include the following acknowledgment:
30 * "This product includes Radiance software
31 * (http://radsite.lbl.gov/)
32 * developed by the Lawrence Berkeley National Laboratory
33 * (http://www.lbl.gov/)."
34 * Alternately, this acknowledgment may appear in the software itself,
35 * if and wherever such third-party acknowledgments normally appear.
36 *
37 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
38 * and "The Regents of the University of California" must
39 * not be used to endorse or promote products derived from this
40 * software without prior written permission. For written
41 * permission, please contact [email protected].
42 *
43 * 5. Products derived from this software may not be called "Radiance",
44 * nor may "Radiance" appear in their name, without prior written
45 * permission of Lawrence Berkeley National Laboratory.
46 *
47 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
48 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
49 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50 * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
51 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
52 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
53 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
54 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
55 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
56 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
57 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
58 * SUCH DAMAGE.
59 * ====================================================================
60 *
61 * This software consists of voluntary contributions made by many
62 * individuals on behalf of Lawrence Berkeley National Laboratory. For more
63 * information on Lawrence Berkeley National Laboratory, please see
64 * <http://www.lbl.gov/>.
65 */
66
67 #include "ray.h"
68
69 #include "view.h"
70
71 #include "source.h"
72
73
74 #define CLIP_ABOVE 1
75 #define CLIP_BELOW 2
76 #define CLIP_RIGHT 3
77 #define CLIP_LEFT 4
78
79 #define MAXVERT 10
80
81 typedef struct splist {
82 struct splist *next; /* next source in list */
83 int sn; /* source number */
84 short nv; /* number of vertices */
85 FLOAT vl[3][2]; /* vertex array (last) */
86 } SPLIST; /* source polygon list */
87
88 extern VIEW ourview; /* our view parameters */
89 extern int hres, vres; /* our image resolution */
90 static SPLIST *sphead = NULL; /* our list of source polys */
91
92
93 static int
94 inregion(p, cv, crit) /* check if vertex is in region */
95 FLOAT p[2];
96 double cv;
97 int crit;
98 {
99 switch (crit) {
100 case CLIP_ABOVE:
101 return(p[1] < cv);
102 case CLIP_BELOW:
103 return(p[1] >= cv);
104 case CLIP_RIGHT:
105 return(p[0] < cv);
106 case CLIP_LEFT:
107 return(p[0] >= cv);
108 }
109 return(-1);
110 }
111
112
113 static
114 clipregion(a, b, cv, crit, r) /* find intersection with boundary */
115 register FLOAT a[2], b[2];
116 double cv;
117 int crit;
118 FLOAT r[2]; /* return value */
119 {
120 switch (crit) {
121 case CLIP_ABOVE:
122 case CLIP_BELOW:
123 r[1] = cv;
124 r[0] = a[0] + (cv-a[1])/(b[1]-a[1])*(b[0]-a[0]);
125 return;
126 case CLIP_RIGHT:
127 case CLIP_LEFT:
128 r[0] = cv;
129 r[1] = a[1] + (cv-a[0])/(b[0]-a[0])*(b[1]-a[1]);
130 return;
131 }
132 }
133
134
135 static int
136 hp_clip_poly(vl, nv, cv, crit, vlo) /* clip polygon to half-plane */
137 FLOAT vl[][2];
138 int nv;
139 double cv;
140 int crit;
141 FLOAT vlo[][2]; /* return value */
142 {
143 FLOAT *s, *p;
144 register int j, nvo;
145
146 s = vl[nv-1];
147 nvo = 0;
148 for (j = 0; j < nv; j++) {
149 p = vl[j];
150 if (inregion(p, cv, crit)) {
151 if (!inregion(s, cv, crit))
152 clipregion(s, p, cv, crit, vlo[nvo++]);
153 vlo[nvo][0] = p[0]; vlo[nvo++][1] = p[1];
154 } else if (inregion(s, cv, crit))
155 clipregion(s, p, cv, crit, vlo[nvo++]);
156 s = p;
157 }
158 return(nvo);
159 }
160
161
162 static int
163 box_clip_poly(vl, nv, xl, xr, yb, ya, vlo) /* clip polygon to box */
164 FLOAT vl[MAXVERT][2];
165 int nv;
166 double xl, xr, yb, ya;
167 FLOAT vlo[MAXVERT][2]; /* return value */
168 {
169 FLOAT vlt[MAXVERT][2];
170 int nvt, nvo;
171
172 nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
173 nvo = hp_clip_poly(vlt, nvt, ya, CLIP_ABOVE, vlo);
174 nvt = hp_clip_poly(vlo, nvo, xl, CLIP_LEFT, vlt);
175 nvo = hp_clip_poly(vlt, nvt, xr, CLIP_RIGHT, vlo);
176
177 return(nvo);
178 }
179
180
181 static double
182 minw2(vl, nv, ar2) /* compute square of minimum width */
183 FLOAT vl[][2];
184 int nv;
185 double ar2;
186 {
187 double d2, w2, w2min, w2max;
188 register FLOAT *p0, *p1, *p2;
189 int i, j;
190 /* find minimum for all widths */
191 w2min = FHUGE;
192 p0 = vl[nv-1];
193 for (i = 0; i < nv; i++) { /* for each edge */
194 p1 = vl[i];
195 d2 = (p1[0]-p0[0])*(p1[0]-p0[0]) +
196 (p1[1]-p0[1])*(p1[1]-p0[1])*ar2;
197 w2max = 0.; /* find maximum for this side */
198 for (j = 1; j < nv-1; j++) {
199 p2 = vl[(i+j)%nv];
200 w2 = (p1[0]-p0[0])*(p2[1]-p0[1]) -
201 (p1[1]-p0[1])*(p2[0]-p0[0]);
202 w2 = w2*w2*ar2/d2; /* triangle height squared */
203 if (w2 > w2max)
204 w2max = w2;
205 }
206 if (w2max < w2min) /* global min. based on local max.'s */
207 w2min = w2max;
208 p0 = p1;
209 }
210 return(w2min);
211 }
212
213
214 static
215 convex_center(vl, nv, cv) /* compute center of convex polygon */
216 register FLOAT vl[][2];
217 int nv;
218 FLOAT cv[2]; /* return value */
219 {
220 register int i;
221 /* simple average (suboptimal) */
222 cv[0] = cv[1] = 0.;
223 for (i = 0; i < nv; i++) {
224 cv[0] += vl[i][0];
225 cv[1] += vl[i][1];
226 }
227 cv[0] /= (double)nv;
228 cv[1] /= (double)nv;
229 }
230
231
232 static double
233 poly_area(vl, nv) /* compute area of polygon */
234 register FLOAT vl[][2];
235 int nv;
236 {
237 double a;
238 FLOAT v0[2], v1[2];
239 register int i;
240
241 a = 0.;
242 v0[0] = vl[1][0] - vl[0][0];
243 v0[1] = vl[1][1] - vl[0][1];
244 for (i = 2; i < nv; i++) {
245 v1[0] = vl[i][0] - vl[0][0];
246 v1[1] = vl[i][1] - vl[0][1];
247 a += v0[0]*v1[1] - v0[1]*v1[0];
248 v0[0] = v1[0]; v0[1] = v1[1];
249 }
250 return(a * (a >= 0. ? .5 : -.5));
251 }
252
253
254 static int
255 convex_hull(vl, nv, vlo) /* compute polygon's convex hull */
256 FLOAT vl[][2];
257 int nv;
258 FLOAT vlo[][2]; /* return value */
259 {
260 int nvo, nvt;
261 FLOAT vlt[MAXVERT][2];
262 double voa, vta;
263 register int i, j;
264 /* start with original polygon */
265 for (i = nvo = nv; i--; ) {
266 vlo[i][0] = vl[i][0]; vlo[i][1] = vl[i][1];
267 }
268 voa = poly_area(vlo, nvo); /* compute its area */
269 for (i = 0; i < nvo; i++) { /* for each output vertex */
270 for (j = 0; j < i; j++) {
271 vlt[j][0] = vlo[j][0]; vlt[j][1] = vlo[j][1];
272 }
273 nvt = nvo - 1; /* make poly w/o vertex */
274 for (j = i; j < nvt; j++) {
275 vlt[j][0] = vlo[j+1][0]; vlt[j][1] = vlo[j+1][1];
276 }
277 vta = poly_area(vlt, nvt);
278 if (vta >= voa) { /* is simpler poly bigger? */
279 voa = vta; /* then use it */
280 for (j = nvo = nvt; j--; ) {
281 vlo[j][0] = vlt[j][0]; vlo[j][1] = vlt[j][1];
282 }
283 i--; /* next adjust */
284 }
285 }
286 return(nvo);
287 }
288
289
290 static
291 spinsert(sn, vl, nv) /* insert new source polygon */
292 int sn;
293 FLOAT vl[][2];
294 int nv;
295 {
296 register SPLIST *spn;
297 register int i;
298
299 if (nv < 3)
300 return;
301 if (nv > 3)
302 spn = (SPLIST *)malloc(sizeof(SPLIST)+sizeof(FLOAT)*2*(nv-3));
303 else
304 spn = (SPLIST *)malloc(sizeof(SPLIST));
305 if (spn == NULL)
306 error(SYSTEM, "out of memory in spinsert");
307 spn->sn = sn;
308 for (i = spn->nv = nv; i--; ) {
309 spn->vl[i][0] = vl[i][0]; spn->vl[i][1] = vl[i][1];
310 }
311 spn->next = sphead; /* push onto global list */
312 sphead = spn;
313 }
314
315
316 int
317 sourcepoly(sn, sp) /* compute image polygon for source */
318 int sn;
319 FLOAT sp[MAXVERT][2];
320 {
321 static char cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
322 {0,1,3,7,6,4},{0,1,5,7,6,2},
323 {0,2,6,7,5,1},{0,4,6,7,3,1},
324 {0,2,3,7,5,4},{1,5,4,6,2,3}};
325 register SRCREC *s = source + sn;
326 FVECT ap, ip;
327 FLOAT pt[6][2];
328 int dir;
329 register int i, j;
330
331 if (s->sflags & (SDISTANT|SFLAT)) {
332 if (s->sflags & SDISTANT && ourview.type == VT_PAR)
333 return(0); /* all or nothing case */
334 if (s->sflags & SFLAT) {
335 for (i = 0; i < 3; i++)
336 ap[i] = s->sloc[i] - ourview.vp[i];
337 if (DOT(ap, s->snorm) >= 0.)
338 return(0); /* source faces away */
339 }
340 for (j = 0; j < 4; j++) { /* four corners */
341 for (i = 0; i < 3; i++) {
342 ap[i] = s->sloc[i];
343 if (j==1|j==2) ap[i] += s->ss[SU][i];
344 else ap[i] -= s->ss[SU][i];
345 if (j==2|j==3) ap[i] += s->ss[SV][i];
346 else ap[i] -= s->ss[SV][i];
347 if (s->sflags & SDISTANT) {
348 ap[i] *= 1. + ourview.vfore;
349 ap[i] += ourview.vp[i];
350 }
351 }
352 viewloc(ip, &ourview, ap); /* find image point */
353 if (ip[2] <= 0.)
354 return(0); /* in front of view */
355 sp[j][0] = ip[0]; sp[j][1] = ip[1];
356 }
357 return(4);
358 }
359 /* identify furthest corner */
360 for (i = 0; i < 3; i++)
361 ap[i] = s->sloc[i] - ourview.vp[i];
362 dir = (DOT(ap,s->ss[SU])>0.) |
363 (DOT(ap,s->ss[SV])>0.)<<1 |
364 (DOT(ap,s->ss[SW])>0.)<<2 ;
365 /* order vertices based on this */
366 for (j = 0; j < 6; j++) {
367 for (i = 0; i < 3; i++) {
368 ap[i] = s->sloc[i];
369 if (cubeord[dir][j] & 1) ap[i] += s->ss[SU][i];
370 else ap[i] -= s->ss[SU][i];
371 if (cubeord[dir][j] & 2) ap[i] += s->ss[SV][i];
372 else ap[i] -= s->ss[SV][i];
373 if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
374 else ap[i] -= s->ss[SW][i];
375 }
376 viewloc(ip, &ourview, ap); /* find image point */
377 if (ip[2] <= 0.)
378 return(0); /* in front of view */
379 pt[j][0] = ip[0]; pt[j][1] = ip[1];
380 }
381 return(convex_hull(pt, 6, sp)); /* make sure it's convex */
382 }
383
384
385 /* initialize by finding sources smaller than rad */
386 init_drawsources(rad)
387 int rad; /* source sample size */
388 {
389 FLOAT spoly[MAXVERT][2];
390 int nsv;
391 register SPLIST *sp;
392 register int i;
393 /* free old source list if one */
394 for (sp = sphead; sp != NULL; sp = sphead) {
395 sphead = sp->next;
396 free((void *)sp);
397 }
398 /* loop through all sources */
399 for (i = nsources; i--; ) {
400 /* compute image polygon for source */
401 if (!(nsv = sourcepoly(i, spoly)))
402 continue;
403 /* clip to image boundaries */
404 if (!(nsv = box_clip_poly(spoly, nsv, 0., 1., 0., 1., spoly)))
405 continue;
406 /* big enough for standard sampling? */
407 if (minw2(spoly, nsv, ourview.vn2/ourview.hn2) >
408 (double)rad*rad/hres/hres)
409 continue;
410 /* OK, add to our list */
411 spinsert(i, spoly, nsv);
412 }
413 }
414
415 void /* add sources smaller than rad to computed subimage */
416 drawsources(pic, zbf, x0, xsiz, y0, ysiz)
417 COLOR *pic[]; /* subimage pixel value array */
418 float *zbf[]; /* subimage distance array (opt.) */
419 int x0, xsiz, y0, ysiz; /* origin and size of subimage */
420 {
421 FLOAT spoly[MAXVERT][2], ppoly[MAXVERT][2];
422 int nsv, npv;
423 int xmin, xmax, ymin, ymax, x, y;
424 FLOAT cxy[2];
425 double w;
426 RAY sr;
427 register SPLIST *sp;
428 register int i;
429 /* check each source in our list */
430 for (sp = sphead; sp != NULL; sp = sp->next) {
431 /* clip source poly to subimage */
432 nsv = box_clip_poly(sp->vl, sp->nv,
433 (double)x0/hres, (double)(x0+xsiz)/hres,
434 (double)y0/vres, (double)(y0+ysiz)/vres, spoly);
435 if (!nsv)
436 continue;
437 /* find common subimage (BBox) */
438 xmin = x0 + xsiz; xmax = x0;
439 ymin = y0 + ysiz; ymax = y0;
440 for (i = 0; i < nsv; i++) {
441 if ((double)xmin/hres > spoly[i][0])
442 xmin = spoly[i][0]*hres + FTINY;
443 if ((double)xmax/hres < spoly[i][0])
444 xmax = spoly[i][0]*hres - FTINY;
445 if ((double)ymin/vres > spoly[i][1])
446 ymin = spoly[i][1]*vres + FTINY;
447 if ((double)ymax/vres < spoly[i][1])
448 ymax = spoly[i][1]*vres - FTINY;
449 }
450 /* evaluate each pixel in BBox */
451 for (y = ymin; y <= ymax; y++)
452 for (x = xmin; x <= xmax; x++) {
453 /* subarea for pixel */
454 npv = box_clip_poly(spoly, nsv,
455 (double)x/hres, (x+1.)/hres,
456 (double)y/vres, (y+1.)/vres,
457 ppoly);
458 if (!npv)
459 continue; /* no overlap */
460 convex_center(ppoly, npv, cxy);
461 if ((sr.rmax = viewray(sr.rorg,sr.rdir,&ourview,
462 cxy[0],cxy[1])) < -FTINY)
463 continue; /* not in view */
464 if (source[sp->sn].sflags & SSPOT &&
465 spotout(&sr, source[sp->sn].sl.s))
466 continue; /* outside spot */
467 rayorigin(&sr, NULL, SHADOW, 1.0);
468 sr.rsrc = sp->sn;
469 rayvalue(&sr); /* compute value */
470 if (bright(sr.rcol) <= FTINY)
471 continue; /* missed/blocked */
472 /* modify pixel */
473 if (zbf[y-y0] != NULL &&
474 sr.rt < 0.999*zbf[y-y0][x-x0])
475 zbf[y-y0][x-x0] = sr.rt;
476 else if (!bigdiff(sr.rcol, pic[y-y0][x-x0],
477 0.001)) /* source sample */
478 setcolor(pic[y-y0][x-x0], 0., 0., 0.);
479 w = poly_area(ppoly, npv) * hres * vres;
480 scalecolor(sr.rcol, w);
481 scalecolor(pic[y-y0][x-x0], 1.-w);
482 addcolor(pic[y-y0][x-x0], sr.rcol);
483 }
484 }
485 }