ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcdraw.c
Revision: 2.1
Committed: Fri Mar 15 21:26:42 1996 UTC (28 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

File Contents

# User Rev Content
1 greg 2.1 /* Copyright (c) 1996 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * Draw small sources into image in case we missed them.
9     */
10    
11     #include "ray.h"
12    
13     #include "view.h"
14    
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    
26     static int
27     inregion(p, cv, crit) /* check if vertex is in region */
28     FLOAT p[2];
29     double cv;
30     int crit;
31     {
32     switch (crit) {
33     case CLIP_ABOVE:
34     return(p[1] < cv);
35     case CLIP_BELOW:
36     return(p[1] >= cv);
37     case CLIP_RIGHT:
38     return(p[0] < cv);
39     case CLIP_LEFT:
40     return(p[0] >= cv);
41     }
42     return(-1);
43     }
44    
45    
46     static
47     clipregion(a, b, cv, crit, r) /* find intersection with boundary */
48     register FLOAT a[2], b[2];
49     double cv;
50     int crit;
51     FLOAT r[2]; /* return value */
52     {
53     switch (crit) {
54     case CLIP_ABOVE:
55     case CLIP_BELOW:
56     r[1] = cv;
57     r[0] = a[0] + (cv-a[1])/(b[1]-a[1])*(b[0]-a[0]);
58     return;
59     case CLIP_RIGHT:
60     case CLIP_LEFT:
61     r[0] = cv;
62     r[1] = a[1] + (cv-a[0])/(b[0]-a[0])*(b[1]-a[1]);
63     return;
64     }
65     }
66    
67    
68     static int
69     hp_clip_poly(vl, nv, cv, crit, vlo) /* clip polygon to half-plane */
70     FLOAT vl[][2];
71     int nv;
72     double cv;
73     int crit;
74     FLOAT vlo[][2]; /* return value */
75     {
76     FLOAT *s, *p;
77     register int j, nvo;
78    
79     s = vl[nv-1];
80     nvo = 0;
81     for (j = 0; j < nv; j++) {
82     p = vl[j];
83     if (inregion(p, cv, crit)) {
84     if (!inregion(s, cv, crit))
85     clipregion(s, p, cv, crit, vlo[nvo++]);
86     vlo[nvo][0] = p[0]; vlo[nvo++][1] = p[1];
87     } else if (inregion(s, cv, crit))
88     clipregion(s, p, cv, crit, vlo[nvo++]);
89     s = p;
90     }
91     return(nvo);
92     }
93    
94    
95     static int
96     box_clip_poly(vl, nv, xl, xr, yb, ya, vlo) /* clip polygon to box */
97     FLOAT vl[MAXVERT][2];
98     int nv;
99     double xl, xr, yb, ya;
100     FLOAT vlo[MAXVERT][2]; /* return value */
101     {
102     FLOAT vlt[MAXVERT][2];
103     int nvt, nvo;
104    
105     nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
106     nvo = hp_clip_poly(vlt, nvt, ya, CLIP_ABOVE, vlo);
107     nvt = hp_clip_poly(vlo, nvo, xl, CLIP_LEFT, vlt);
108     nvo = hp_clip_poly(vlt, nvt, xr, CLIP_RIGHT, vlo);
109    
110     return(nvo);
111     }
112    
113    
114     static double
115     minw2(vl, nv, ar2) /* compute square of minimum width */
116     FLOAT vl[][2];
117     int nv;
118     double ar2;
119     {
120     double d2, w2, w2min, w2max;
121     register FLOAT *p0, *p1, *p2;
122     int i, j;
123     /* find minimum for all widths */
124     w2min = FHUGE;
125     p0 = vl[nv-1];
126     for (i = 0; i < nv; i++) { /* for each edge */
127     p1 = vl[i];
128     d2 = (p1[0]-p0[0])*(p1[0]-p0[0]) +
129     (p1[1]-p0[1])*(p1[1]-p0[1])*ar2;
130     w2max = 0.; /* find maximum for this side */
131     for (j = 1; j < nv-1; j++) {
132     p2 = vl[(i+j)%nv];
133     w2 = (p1[0]-p0[0])*(p2[1]-p0[1]) -
134     (p1[1]-p0[1])*(p2[0]-p0[0]);
135     w2 = w2*w2*ar2/d2; /* triangle height squared */
136     if (w2 > w2max)
137     w2max = w2;
138     }
139     if (w2max < w2min) /* global min. based on local max.'s */
140     w2min = w2max;
141     p0 = p1;
142     }
143     return(w2min);
144     }
145    
146    
147     static
148     convex_center(vl, nv, cv) /* compute center of convex polygon */
149     register FLOAT vl[][2];
150     int nv;
151     FLOAT cv[2]; /* return value */
152     {
153     register int i;
154     /* simple average (suboptimal) */
155     cv[0] = cv[1] = 0.;
156     for (i = 0; i < nv; i++) {
157     cv[0] += vl[i][0];
158     cv[1] += vl[i][1];
159     }
160     cv[0] /= (double)nv;
161     cv[1] /= (double)nv;
162     }
163    
164    
165     static double
166     poly_area(vl, nv) /* compute area of polygon */
167     register FLOAT vl[][2];
168     int nv;
169     {
170     double a;
171     FLOAT v0[2], v1[2];
172     register int i;
173    
174     a = 0.;
175     v0[0] = vl[1][0] - vl[0][0];
176     v0[1] = vl[1][1] - vl[0][1];
177     for (i = 2; i < nv; i++) {
178     v1[0] = vl[i][0] - vl[0][0];
179     v1[1] = vl[i][1] - vl[0][1];
180     a += v0[0]*v1[1] - v0[1]*v1[0];
181     v0[0] = v1[0]; v0[1] = v1[1];
182     }
183     return(a * (a >= 0. ? .5 : -.5));
184     }
185    
186    
187     static int
188     convex_hull(vl, nv, vlo) /* compute polygon's convex hull */
189     FLOAT vl[][2];
190     int nv;
191     FLOAT vlo[][2]; /* return value */
192     {
193     int nvo, nvt;
194     FLOAT vlt[MAXVERT][2];
195     double voa, vta;
196     register int i, j;
197     /* start with original polygon */
198     for (i = nvo = nv; i--; ) {
199     vlo[i][0] = vl[i][0]; vlo[i][1] = vl[i][1];
200     }
201     voa = poly_area(vlo, nvo); /* compute its area */
202     for (i = 0; i < nvo; i++) { /* for each output vertex */
203     for (j = 0; j < i; j++) {
204     vlt[j][0] = vlo[j][0]; vlt[j][1] = vlo[j][1];
205     }
206     nvt = nvo - 1; /* make poly w/o vertex */
207     for (j = i; j < nvt; j++) {
208     vlt[j][0] = vlo[j+1][0]; vlt[j][1] = vlo[j+1][1];
209     }
210     vta = poly_area(vlt, nvt);
211     if (vta >= voa) { /* is simpler poly bigger? */
212     voa = vta; /* then use it */
213     for (j = nvo = nvt; j--; ) {
214     vlo[j][0] = vlt[j][0]; vlo[j][1] = vlt[j][1];
215     }
216     i--; /* next adjust */
217     }
218     }
219     return(nvo);
220     }
221    
222    
223     int
224     sourcepoly(vw, sn, sp) /* compute image polygon for source */
225     VIEW *vw;
226     int sn;
227     FLOAT sp[MAXVERT][2];
228     {
229     static char cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
230     {0,1,3,7,6,4},{0,1,5,7,6,2},
231     {0,2,6,7,5,1},{0,4,6,7,3,1},
232     {0,2,3,7,5,4},{1,5,4,6,2,3}};
233     register SRCREC *s = source + sn;
234     FVECT ap, ip;
235     FLOAT pt[6][2];
236     int dir;
237     register int i, j;
238    
239     if (s->sflags & (SDISTANT|SFLAT)) {
240     if (s->sflags & SDISTANT && vw->type == VT_PAR)
241     return(0); /* all or nothing case */
242     if (s->sflags & SFLAT) {
243     for (i = 0; i < 3; i++)
244     ap[i] = s->sloc[i] - vw->vp[i];
245     if (DOT(ap, s->snorm) >= 0.)
246     return(0); /* source faces away */
247     }
248     for (j = 0; j < 4; j++) { /* four corners */
249     for (i = 0; i < 3; i++) {
250     ap[i] = s->sloc[i];
251     if (j==1|j==2) ap[i] += s->ss[SU][i];
252     else ap[i] -= s->ss[SU][i];
253     if (j==2|j==3) ap[i] += s->ss[SV][i];
254     else ap[i] -= s->ss[SV][i];
255     if (s->sflags & SDISTANT) {
256     ap[i] *= 1. + vw->vfore;
257     ap[i] += vw->vp[i];
258     }
259     }
260     viewloc(ip, vw, ap); /* find image point */
261     if (ip[2] <= 0.)
262     return(0); /* in front of view */
263     sp[j][0] = ip[0]; sp[j][1] = ip[1];
264     }
265     return(4);
266     }
267     /* identify furthest corner */
268     for (i = 0; i < 3; i++)
269     ap[i] = s->sloc[i] - vw->vp[i];
270     dir = (DOT(ap,s->ss[SU])>0.) |
271     (DOT(ap,s->ss[SV])>0.)<<1 |
272     (DOT(ap,s->ss[SW])>0.)<<2 ;
273     /* order vertices based on this */
274     for (j = 0; j < 6; j++) {
275     for (i = 0; i < 3; i++) {
276     ap[i] = s->sloc[i];
277     if (cubeord[dir][j] & 1) ap[i] += s->ss[SU][i];
278     else ap[i] -= s->ss[SU][i];
279     if (cubeord[dir][j] & 2) ap[i] += s->ss[SV][i];
280     else ap[i] -= s->ss[SV][i];
281     if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
282     else ap[i] -= s->ss[SW][i];
283     }
284     viewloc(ip, vw, ap); /* find image point */
285     if (ip[2] <= 0.)
286     return(0); /* in front of view */
287     pt[j][0] = ip[0]; pt[j][1] = ip[1];
288     }
289     return(convex_hull(pt, 6, sp)); /* make sure it's convex */
290     }
291    
292    
293     /* add sources smaller than rad to computed subimage */
294     drawsources(vw, xr, yr, pic, zbf, x0, xsiz, y0, ysiz, rad)
295     VIEW *vw; /* full image view */
296     int xr, yr; /* full image dimensions */
297     COLOR *pic[]; /* subimage pixel value array */
298     float *zbf[]; /* subimage distance array (opt.) */
299     int x0, xsiz, y0, ysiz; /* origin and size of subimage */
300     int rad; /* source sample size */
301     {
302     int sn;
303     FLOAT spoly[MAXVERT][2], ppoly[MAXVERT][2];
304     int nsv, npv;
305     int xmin, xmax, ymin, ymax, x, y, i;
306     FLOAT cxy[2];
307     double pa;
308     RAY sr;
309     /* loop through all sources */
310     for (sn = 0; sn < nsources; sn++) {
311     /* compute image polygon for source */
312     if (!(nsv = sourcepoly(vw, sn, spoly)))
313     continue;
314     /* clip source poly to subimage */
315     nsv = box_clip_poly(spoly, nsv,
316     (double)x0/xr, (double)(x0+xsiz)/xr,
317     (double)y0/yr, (double)(y0+ysiz)/yr, spoly);
318     if (!nsv)
319     continue;
320     /* is it big so sampling found it? */
321     if (minw2(spoly, nsv, vw->vn2/vw->hn2) > (double)rad*rad/xr/xr)
322     continue;
323     /* find common subimage (BBox) */
324     xmin = x0 + xsiz; xmax = x0;
325     ymin = y0 + ysiz; ymax = y0;
326     for (i = 0; i < nsv; i++) {
327     if ((double)xmin/xr > spoly[i][0])
328     xmin = spoly[i][0]*xr + FTINY;
329     if ((double)xmax/xr < spoly[i][0])
330     xmax = spoly[i][0]*xr - FTINY;
331     if ((double)ymin/yr > spoly[i][1])
332     ymin = spoly[i][1]*yr + FTINY;
333     if ((double)ymax/yr < spoly[i][1])
334     ymax = spoly[i][1]*yr - FTINY;
335     }
336     /* evaluate each pixel in BBox */
337     for (y = ymin; y <= ymax; y++)
338     for (x = xmin; x <= xmax; x++) {
339     /* subarea for pixel */
340     npv = box_clip_poly(spoly, nsv,
341     (double)x/xr, (x+1.)/xr,
342     (double)y/yr, (y+1.)/yr, ppoly);
343     if (!npv)
344     continue; /* no overlap */
345     convex_center(ppoly, npv, cxy);
346     if ((sr.rmax = viewray(sr.rorg, sr.rdir, vw,
347     cxy[0], cxy[1])) < -FTINY)
348     continue; /* not in view */
349     if (source[sn].sflags & SSPOT &&
350     spotout(sr, source[sn].sl.s))
351     continue; /* outside spot */
352     rayorigin(&sr, NULL, SHADOW, 1.0);
353     sr.rsrc = sn;
354     rayvalue(&sr); /* compute value */
355     if (bright(sr.rcol) <= FTINY)
356     continue; /* missed/blocked */
357     /* modify pixel */
358     if (zbf[y-y0] != NULL &&
359     sr.rt < zbf[y-y0][x-x0])
360     zbf[y-y0][x-x0] = sr.rt;
361     pa = poly_area(ppoly, npv);
362     scalecolor(sr.rcol, pa*xr*yr);
363     scalecolor(pic[y-y0][x-x0], (1.-pa*xr*yr));
364     addcolor(pic[y-y0][x-x0], sr.rcol);
365     }
366     }
367     }