ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/findglare.c
Revision: 1.25
Committed: Tue Nov 12 13:21:01 1991 UTC (32 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.24: +1 -1 lines
Log Message:
fixed bug in floating point comparison

File Contents

# Content
1 /* Copyright (c) 1991 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Find glare sources in a scene or image.
9 *
10 * Greg Ward March 1991
11 */
12
13 #include "glare.h"
14
15 #define FEQ(a,b) ((a)-(b)<=FTINY&&(b)-(a)<=FTINY)
16 #define VEQ(v1,v2) (FEQ((v1)[0],(v2)[0])&&FEQ((v1)[1],(v2)[1]) \
17 &&FEQ((v1)[2],(v2)[2]))
18
19 char *rtargv[32] = {"rtrace", "-h-", "-ov", "-fff"};
20 int rtargc = 4;
21
22 VIEW ourview = STDVIEW; /* our view */
23 VIEW pictview = STDVIEW; /* picture view */
24 VIEW leftview, rightview; /* leftmost and rightmost views */
25
26 char *picture = NULL; /* picture file name */
27 char *octree = NULL; /* octree file name */
28
29 int verbose = 0; /* verbose reporting */
30 char *progname; /* global argv[0] */
31
32 double threshold = 0.; /* glare threshold */
33
34 int sampdens = SAMPDENS; /* sample density */
35 ANGLE glarang[180] = {AEND}; /* glare calculation angles */
36 int nglarangs = 0;
37 double maxtheta; /* maximum angle (in radians) */
38 int hsize; /* horizontal size */
39
40 struct illum *indirect; /* array of indirect illuminances */
41
42 long npixinvw; /* number of pixels in view */
43 long npixmiss; /* number of pixels missed */
44
45
46 main(argc, argv)
47 int argc;
48 char *argv[];
49 {
50 int combine = 1;
51 int gotview = 0;
52 int rval, i;
53 char *err;
54
55 progname = argv[0];
56 /* process options */
57 for (i = 1; i < argc && argv[i][0] == '-'; i++) {
58 rval = getviewopt(&ourview, argc-i, argv+i);
59 if (rval >= 0) {
60 i += rval;
61 gotview++;
62 continue;
63 }
64 switch (argv[i][1]) {
65 case 't':
66 threshold = atof(argv[++i]);
67 break;
68 case 'r':
69 sampdens = atoi(argv[++i])/2;
70 break;
71 case 'v':
72 if (argv[i][2] == '\0') {
73 verbose++;
74 break;
75 }
76 if (argv[i][2] != 'f')
77 goto userr;
78 rval = viewfile(argv[++i], &ourview, 0, 0);
79 if (rval < 0) {
80 fprintf(stderr,
81 "%s: cannot open view file \"%s\"\n",
82 progname, argv[i]);
83 exit(1);
84 } else if (rval == 0) {
85 fprintf(stderr,
86 "%s: bad view file \"%s\"\n",
87 progname, argv[i]);
88 exit(1);
89 } else
90 gotview++;
91 break;
92 case 'g':
93 if (argv[i][2] != 'a')
94 goto userr;
95 if (setscan(glarang, argv[++i]) < 0) {
96 fprintf(stderr, "%s: bad angle spec \"%s\"\n",
97 progname, argv[i]);
98 exit(1);
99 }
100 break;
101 case 'p':
102 picture = argv[++i];
103 break;
104 case 'c':
105 combine = !combine;
106 break;
107 case 'd':
108 case 'l':
109 rtargv[rtargc++] = argv[i];
110 rtargv[rtargc++] = argv[++i];
111 break;
112 case 'a':
113 rtargv[rtargc++] = argv[i];
114 if (argv[i][2] == 'v') {
115 rtargv[rtargc++] = argv[++i];
116 rtargv[rtargc++] = argv[++i];
117 }
118 rtargv[rtargc++] = argv[++i];
119 break;
120 default:
121 goto userr;
122 }
123 }
124 /* get octree */
125 if (i < argc-1)
126 goto userr;
127 if (i == argc-1) {
128 rtargv[rtargc++] = octree = argv[i];
129 rtargv[rtargc] = NULL;
130 }
131 /* get view */
132 if (picture != NULL) {
133 rval = viewfile(picture, &pictview, 0, 0);
134 if (rval < 0) {
135 fprintf(stderr, "%s: cannot open picture file \"%s\"\n",
136 progname, picture);
137 exit(1);
138 } else if (rval == 0) {
139 fprintf(stderr,
140 "%s: cannot get view from picture \"%s\"\n",
141 progname, picture);
142 exit(1);
143 }
144 if (pictview.type == VT_PAR) {
145 fprintf(stderr, "%s: %s: cannot use parallel view\n",
146 progname, picture);
147 exit(1);
148 }
149 if ((err = setview(&pictview)) != NULL) {
150 fprintf(stderr, "%s: %s\n", picture, err);
151 exit(1);
152 }
153 }
154 if (!gotview) {
155 if (picture == NULL) {
156 fprintf(stderr, "%s: must have view or picture\n",
157 progname);
158 exit(1);
159 }
160 copystruct(&ourview, &pictview);
161 } else if (picture != NULL && !VEQ(ourview.vp, pictview.vp)) {
162 fprintf(stderr, "%s: picture must have same viewpoint\n",
163 progname);
164 exit(1);
165 }
166 ourview.type = VT_HEM;
167 ourview.horiz = ourview.vert = 180.0;
168 ourview.hoff = ourview.voff = 0.0;
169 fvsum(ourview.vdir, ourview.vdir, ourview.vup,
170 -DOT(ourview.vdir,ourview.vup));
171 if ((err = setview(&ourview)) != NULL) {
172 fprintf(stderr, "%s: %s\n", progname, err);
173 exit(1);
174 }
175 if (octree == NULL && picture == NULL) {
176 fprintf(stderr,
177 "%s: must specify at least one of picture or octree\n",
178 progname);
179 exit(1);
180 }
181 init(); /* initialize program */
182 if (threshold <= FTINY)
183 comp_thresh(); /* compute glare threshold */
184 analyze(); /* analyze view */
185 if (combine)
186 absorb_specks(); /* eliminate tiny sources */
187 cleanup(); /* tidy up */
188 /* print header */
189 printargs(argc, argv, stdout);
190 fputs(VIEWSTR, stdout);
191 fprintview(&ourview, stdout);
192 printf("\n");
193 fputformat("ascii", stdout);
194 printf("\n");
195 printsources(); /* print glare sources */
196 printillum(); /* print illuminances */
197 exit(0);
198 userr:
199 fprintf(stderr,
200 "Usage: %s [view options][-ga angles][-p picture][[rtrace options] octree]\n",
201 progname);
202 exit(1);
203 }
204
205
206 int
207 angcmp(ap1, ap2) /* compare two angles */
208 ANGLE *ap1, *ap2;
209 {
210 register int a1, a2;
211
212 a1 = *ap1;
213 a2 = *ap2;
214 if (a1 == a2) {
215 fprintf(stderr, "%s: duplicate glare angle (%d)\n",
216 progname, a1);
217 exit(1);
218 }
219 return(a1-a2);
220 }
221
222
223 init() /* initialize global variables */
224 {
225 double d;
226 register int i;
227
228 if (verbose)
229 fprintf(stderr, "%s: initializing data structures...\n",
230 progname);
231 /* set direction vectors */
232 for (i = 0; glarang[i] != AEND; i++)
233 ;
234 qsort(glarang, i, sizeof(ANGLE), angcmp);
235 if (i > 0 && (glarang[0] <= 0 || glarang[i-1] > 180)) {
236 fprintf(stderr, "%s: glare angles must be between 1 and 180\n",
237 progname);
238 exit(1);
239 }
240 nglarangs = i;
241 /* nglardirs = 2*nglarangs + 1; */
242 /* vsize = sampdens - 1; */
243 if (nglarangs > 0)
244 maxtheta = (PI/180.)*glarang[nglarangs-1];
245 else
246 maxtheta = 0.0;
247 hsize = hlim(0) + sampdens - 1;
248 if (hsize > (int)(PI*sampdens))
249 hsize = PI*sampdens;
250 indirect = (struct illum *)calloc(nglardirs, sizeof(struct illum));
251 if (indirect == NULL)
252 memerr("indirect illuminances");
253 npixinvw = npixmiss = 0L;
254 copystruct(&leftview, &ourview);
255 copystruct(&rightview, &ourview);
256 spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta);
257 spinvector(rightview.vdir, ourview.vdir, ourview.vup, -maxtheta);
258 setview(&leftview);
259 setview(&rightview);
260 indirect[nglarangs].lcos =
261 indirect[nglarangs].rcos = cos(maxtheta);
262 indirect[nglarangs].rsin =
263 -(indirect[nglarangs].lsin = sin(maxtheta));
264 indirect[nglarangs].theta = 0.0;
265 for (i = 0; i < nglarangs; i++) {
266 d = (glarang[nglarangs-1] - glarang[i])*(PI/180.);
267 indirect[nglarangs-i-1].lcos =
268 indirect[nglarangs+i+1].rcos = cos(d);
269 indirect[nglarangs+i+1].rsin =
270 -(indirect[nglarangs-i-1].lsin = sin(d));
271 d = (glarang[nglarangs-1] + glarang[i])*(PI/180.);
272 indirect[nglarangs-i-1].rcos =
273 indirect[nglarangs+i+1].lcos = cos(d);
274 indirect[nglarangs-i-1].rsin =
275 -(indirect[nglarangs+i+1].lsin = sin(d));
276 indirect[nglarangs-i-1].theta = (PI/180.)*glarang[i];
277 indirect[nglarangs+i+1].theta = -(PI/180.)*glarang[i];
278 }
279 /* open picture */
280 if (picture != NULL) {
281 if (verbose)
282 fprintf(stderr, "%s: opening picture file \"%s\"...\n",
283 progname, picture);
284 open_pict(picture);
285 }
286 /* start rtrace */
287 if (octree != NULL) {
288 if (verbose) {
289 fprintf(stderr,
290 "%s: starting luminance calculation...\n\t",
291 progname);
292 printargs(rtargc, rtargv, stderr);
293 }
294 fork_rtrace(rtargv);
295 }
296 }
297
298
299 cleanup() /* close files, wait for children */
300 {
301 if (verbose)
302 fprintf(stderr, "%s: cleaning up... \n", progname);
303 if (picture != NULL)
304 close_pict();
305 if (octree != NULL)
306 done_rtrace();
307 if (npixinvw < 100*npixmiss)
308 fprintf(stderr, "%s: warning -- missing %ld%% of samples\n",
309 progname, 100L*npixmiss/npixinvw);
310 }
311
312
313 compdir(vd, x, y) /* compute direction for x,y */
314 FVECT vd;
315 int x, y;
316 {
317 int hl;
318 FVECT org; /* dummy variable */
319
320 hl = hlim(y);
321 if (x <= -hl) { /* left region */
322 if (x <= -hl-sampdens)
323 return(-1);
324 return(viewray(org, vd, &leftview,
325 (double)(x+hl)/(2*sampdens)+.5,
326 (double)y/(2*sampdens)+.5));
327 }
328 if (x >= hl) { /* right region */
329 if (x >= hl+sampdens)
330 return(-1);
331 return(viewray(org, vd, &rightview,
332 (double)(x-hl)/(2*sampdens)+.5,
333 (double)y/(2*sampdens)+.5));
334 }
335 /* central region */
336 if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0)
337 return(-1);
338 spinvector(vd, vd, ourview.vup, h_theta(x,y));
339 return(0);
340 }
341
342
343 double
344 pixsize(x, y) /* return the solid angle of pixel at (x,y) */
345 int x, y;
346 {
347 register int hl, xo;
348 double disc;
349
350 hl = hlim(y);
351 if (x < -hl)
352 xo = x+hl;
353 else if (x > hl)
354 xo = x-hl;
355 else
356 xo = 0;
357 disc = 1. - (double)(xo*xo + y*y)/(sampdens*sampdens);
358 if (disc <= FTINY)
359 return(0.);
360 return(1./(sampdens*sampdens*sqrt(disc)));
361 }
362
363
364 memerr(s) /* malloc failure */
365 char *s;
366 {
367 fprintf(stderr, "%s: out of memory for %s\n", progname, s);
368 exit(1);
369 }
370
371
372 printsources() /* print out glare sources */
373 {
374 register struct source *sp;
375
376 printf("BEGIN glare source\n");
377 for (sp = donelist; sp != NULL; sp = sp->next)
378 printf("\t%f %f %f\t%f\t%f\n",
379 sp->dir[0], sp->dir[1], sp->dir[2],
380 sp->dom, sp->brt);
381 printf("END glare source\n");
382 }
383
384
385 printillum() /* print out indirect illuminances */
386 {
387 register int i;
388
389 printf("BEGIN indirect illuminance\n");
390 for (i = 0; i < nglardirs; i++)
391 if (indirect[i].n > FTINY)
392 printf("\t%.0f\t%f\n", (180.0/PI)*indirect[i].theta,
393 PI * indirect[i].sum / indirect[i].n);
394 printf("END indirect illuminance\n");
395 }