src/util/findglare.c

/* Copyright (c) 1991 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Find glare sources in a scene or image.
 *
 *      Greg Ward       March 1991
 */

#include "glare.h"

#include "random.h"

#define FEQ(a,b)        ((a)-(b)<=FTINY&&(a)-(b)<=FTINY)
#define VEQ(v1,v2)      (FEQ((v1)[0],(v2)[0])&&FEQ((v1)[1],(v2)[1]) \
                                &&FEQ((v1)[2],(v2)[2]))

char    *rtargv[32] = {"rtrace", "-h", "-ov", "-fff"};
int     rtargc = 4;

VIEW    ourview = STDVIEW;              /* our view */
VIEW    pictview = STDVIEW;             /* picture view */
VIEW    leftview, rightview;            /* leftmost and rightmost views */

char    *picture = NULL;                /* picture file name */
char    *octree = NULL;                 /* octree file name */

int     verbose = 0;                    /* verbose reporting */
char    *progname;                      /* global argv[0] */

double  threshold = 0.;                 /* glare threshold */

int     sampdens = SAMPDENS;            /* sample density */
ANGLE   glarang[180] = {AEND};          /* glare calculation angles */
int     nglarangs = 0;
double  maxtheta;                       /* maximum angle (in radians) */
int     hsize;                          /* horizontal size */
int     hlim;                           /* central limit of horizontal */

struct illum    *indirect;              /* array of indirect illuminances */

long    npixinvw;                       /* number of pixels in view */
long    npixmiss;                       /* number of pixels missed */


main(argc, argv)
int     argc;
char    *argv[];
{
        int     gotview = 0;
        int     rval, i;
        char    *err;

        progname = argv[0];
                                        /* process options */
        for (i = 1; i < argc && argv[i][0] == '-'; i++) {
                rval = getviewopt(&ourview, argc-i, argv+i);
                if (rval >= 0) {
                        i += rval;
                        gotview++;
                        continue;
                }
                switch (argv[i][1]) {
                case 't':
                        threshold = atof(argv[++i]);
                        break;
                case 'r':
                        sampdens = atoi(argv[++i])/2;
                        break;
                case 'v':
                        if (argv[i][2] == '\0') {
                                verbose++;
                                break;
                        }
                        if (argv[i][2] != 'f')
                                goto userr;
                        rval = viewfile(argv[++i], &ourview, 0, 0);
                        if (rval < 0) {
                                fprintf(stderr,
                                "%s: cannot open view file \"%s\"\n",
                                                progname, argv[i]);
                                exit(1);
                        } else if (rval == 0) {
                                fprintf(stderr,
                                        "%s: bad view file \"%s\"\n",
                                                progname, argv[i]);
                                exit(1);
                        } else
                                gotview++;
                        break;
                case 'g':
                        if (argv[i][2] != 'a')
                                goto userr;
                        if (setscan(glarang, argv[++i]) < 0) {
                                fprintf(stderr, "%s: bad angle spec \"%s\"\n",
                                        progname, argv[i]);
                                exit(1);
                        }
                        break;
                case 'p':
                        picture = argv[++i];
                        break;
                case 'd':
                case 'l':
                        rtargv[rtargc++] = argv[i];
                        rtargv[rtargc++] = argv[++i];
                        break;
                case 'a':
                        rtargv[rtargc++] = argv[i];
                        if (argv[i][2] == 'v') {
                                rtargv[rtargc++] = argv[++i];
                                rtargv[rtargc++] = argv[++i];
                        }
                        rtargv[rtargc++] = argv[++i];
                        break;
                default:
                        goto userr;
                }
        }
                                                /* get octree */
        if (i < argc-1)
                goto userr;
        if (i == argc-1) {
                rtargv[rtargc++] = octree = argv[i];
                rtargv[rtargc] = NULL;
        }
                                                /* get view */
        if (picture != NULL) {
                rval = viewfile(picture, &pictview, 0, 0);
                if (rval < 0) {
                        fprintf(stderr, "%s: cannot open picture file \"%s\"\n",
                                        progname, picture);
                        exit(1);
                } else if (rval == 0) {
                        fprintf(stderr,
                                "%s: cannot get view from picture \"%s\"\n",
                                        progname, picture);
                        exit(1);
                }
                if (pictview.type == VT_PAR) {
                        fprintf(stderr, "%s: %s: cannot use parallel view\n",
                                        progname, picture);
                        exit(1);
                }
                if ((err = setview(&pictview)) != NULL) {
                        fprintf(stderr, "%s: %s\n", picture, err);
                        exit(1);
                }
        }
        if (!gotview) {
                if (picture == NULL) {
                        fprintf(stderr, "%s: must have view or picture\n",
                                        progname);
                        exit(1);
                }
                copystruct(&ourview, &pictview);
        } else if (picture != NULL && !VEQ(ourview.vp, pictview.vp)) {
                fprintf(stderr, "%s: picture must have same viewpoint\n",
                                progname);
                exit(1);
        }
        ourview.type = VT_HEM;
        ourview.horiz = ourview.vert = 180.0;
        ourview.hoff = ourview.voff = 0.0;
        fvsum(ourview.vdir, ourview.vdir, ourview.vup,
                        -DOT(ourview.vdir,ourview.vup));
        if ((err = setview(&ourview)) != NULL) {
                fprintf(stderr, "%s: %s\n", progname, err);
                exit(1);
        }
        if (octree == NULL && picture == NULL) {
                fprintf(stderr,
                "%s: must specify at least one of picture or octree\n",
                                progname);
                exit(1);
        }
        init();                                 /* initialize program */
        if (threshold <= FTINY)
                comp_thresh();                  /* compute glare threshold */
        analyze();                              /* analyze view */
        cleanup();                              /* tidy up */
                                                /* print header */
        printargs(argc, argv, stdout);
        fputs(VIEWSTR, stdout);
        fprintview(&ourview, stdout);
        printf("\n\n");
        printsources();                         /* print glare sources */
        printillum();                           /* print illuminances */
        exit(0);
userr:
        fprintf(stderr,
"Usage: %s [view options][-ga angles][-p picture][[rtrace options] octree]\n",
                        progname);
        exit(1);
}


init()                          /* initialize global variables */
{
        double  d;
        register int    i;

        if (verbose)
                fprintf(stderr, "%s: initializing data structures...\n",
                                progname);
                                                /* set direction vectors */
        for (i = 0; glarang[i] != AEND; i++)
                ;
        if (i > 0 && (glarang[0] <= 0 || glarang[i-1] >= 180)) {
                fprintf(stderr, "%s: glare angles must be between 1 and 179\n",
                                progname);
                exit(1);
        }
        nglarangs = i;
        /* nglardirs = 2*nglarangs + 1; */
        /* vsize = sampdens - 1; */
        if (nglarangs > 0)
                maxtheta = (PI/180.)*glarang[nglarangs-1];
        else
                maxtheta = 0.0;
        hlim = sampdens*maxtheta;
        hsize = hlim + sampdens - 1;
        if (hsize > (int)(PI*sampdens))
                hsize = PI*sampdens;
        indirect = (struct illum *)calloc(nglardirs, sizeof(struct illum));
        if (indirect == NULL)
                memerr("indirect illuminances");
        npixinvw = npixmiss = 0L;
        copystruct(&leftview, &ourview);
        copystruct(&rightview, &ourview);
        spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta);
        spinvector(rightview.vdir, ourview.vdir, ourview.vup, -maxtheta);
        setview(&leftview);
        setview(&rightview);
        indirect[nglarangs].lcos =
        indirect[nglarangs].rcos = cos(maxtheta);
        indirect[nglarangs].rsin =
        -(indirect[nglarangs].lsin = sin(maxtheta));
        indirect[nglarangs].theta = 0.0;
        for (i = 0; i < nglarangs; i++) {
                d = (glarang[nglarangs-1] - glarang[i])*(PI/180.);
                indirect[nglarangs-i-1].lcos =
                indirect[nglarangs+i+1].rcos = cos(d);
                indirect[nglarangs+i+1].rsin =
                -(indirect[nglarangs-i-1].lsin = sin(d));
                d = (glarang[nglarangs-1] + glarang[i])*(PI/180.);
                indirect[nglarangs-i-1].rcos =
                indirect[nglarangs+i+1].lcos = cos(d);
                indirect[nglarangs-i-1].rsin =
                -(indirect[nglarangs+i+1].lsin = sin(d));
                indirect[nglarangs-i-1].theta = (PI/180.)*glarang[i];
                indirect[nglarangs+i+1].theta = -(PI/180.)*glarang[i];
        }
                                                /* open picture */
        if (picture != NULL) {
                if (verbose)
                        fprintf(stderr, "%s: opening picture file \"%s\"...\n",
                                        progname, picture);
                open_pict(picture);
        }
                                                /* start rtrace */
        if (octree != NULL) {
                if (verbose) {
                        fprintf(stderr,
                                "%s: starting luminance calculation...\n\t",
                                        progname);
                        printargs(rtargc, rtargv, stderr);
                }
                fork_rtrace(rtargv);
        }
}


cleanup()                               /* close files, wait for children */
{
        if (verbose)
                fprintf(stderr, "%s: cleaning up...\n", progname);
        if (picture != NULL)
                close_pict();
        if (octree != NULL)
                done_rtrace();
        if (npixinvw < 100*npixmiss)
                fprintf(stderr, "%s: warning -- missing %ld%% of samples\n",
                                progname, 100L*npixmiss/npixinvw);
}


compdir(vd, x, y)                       /* compute direction for x,y */
FVECT   vd;
int     x, y;
{
        long    t;
        FVECT   org;                    /* dummy variable */

        if (x <= -hlim)                 /* left region */
                return(viewray(org, vd, &leftview,
                                (double)(x+hlim)/(2*sampdens)+.5,
                                (double)y/(2*sampdens)+.5));
        if (x >= hlim)                  /* right region */
                return(viewray(org, vd, &rightview,
                                (double)(x-hlim)/(2*sampdens)+.5,
                                (double)y/(2*sampdens)+.5));
                                        /* central region */
                                /* avoid over-counting of poles */
        t = random() % vsize;
        if (t*t >= (long)vsize*vsize - (long)y*y)
                return(-1);
        if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0)
                return(-1);
        spinvector(vd, vd, ourview.vup, h_theta(x));
        return(0);
}


memerr(s)                       /* malloc failure */
char    *s;
{
        fprintf(stderr, "%s: out of memory for %s\n", progname, s);
        exit(1);
}


printsources()                  /* print out glare sources */
{
        register struct source  *sp;

        printf("BEGIN glare source\n");
        for (sp = donelist; sp != NULL; sp = sp->next)
                printf("\t%f %f %f\t%f\t%f\n",
                                sp->dir[0], sp->dir[1], sp->dir[2],
                                sp->dom, sp->brt);
        printf("END glare source\n");
}


printillum()                    /* print out indirect illuminances */
{
        register int    i;

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