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 |
} |