ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcdraw.c
Revision: 2.4
Committed: Mon Oct 20 17:18:42 1997 UTC (26 years, 6 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 2.3: +8 -5 lines
Log Message:
checks to see if sample pixel hit source and changes to black
to avoid artifact of overbright pixels -- assumes relatively dark surround

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