--- ray/src/util/findglare.c	1991/03/18 12:15:38	1.1
+++ ray/src/util/findglare.c	1993/02/03 17:55:34	2.5
@@ -12,14 +12,14 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #include "glare.h"
 
-#define FEQ(a,b)	((a)-(b)<=FTINY&&(a)-(b)<=FTINY)
+#define FEQ(a,b)	((a)-(b)<=FTINY&&(b)-(a)<=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"};
+char	*rtargv[32] = {"rtrace", "-h-", "-ov", "-fff"};
 int	rtargc = 4;
 
-VIEW	ourview;			/* our view */
+VIEW	ourview = STDVIEW;		/* our view */
 VIEW	pictview = STDVIEW;		/* picture view */
 VIEW	leftview, rightview;		/* leftmost and rightmost views */
 
@@ -29,18 +29,25 @@ 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	combine = 1;
 	int	gotview = 0;
 	int	rval, i;
 	char	*err;
@@ -55,6 +62,12 @@ char	*argv[];
 			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++;
@@ -62,7 +75,7 @@ char	*argv[];
 			}
 			if (argv[i][2] != 'f')
 				goto userr;
-			rval = viewfile(argv[++i], &ourview, 0, 0);
+			rval = viewfile(argv[++i], &ourview, NULL);
 			if (rval < 0) {
 				fprintf(stderr,
 				"%s: cannot open view file \"%s\"\n",
@@ -88,11 +101,24 @@ char	*argv[];
 		case 'p':
 			picture = argv[++i];
 			break;
+		case 'c':
+			combine = !combine;
+			break;
 		case 'd':
+			if (argv[i][2] == 'v') {
+				rtargv[rtargc++] = argv[i];
+				break;
+			}
+			/* FALL THROUGH */
 		case 'l':
+		case 's':
+		case 'P':
 			rtargv[rtargc++] = argv[i];
 			rtargv[rtargc++] = argv[++i];
 			break;
+		case 'w':
+			rtargv[rtargc++] = argv[i];
+			break;
 		case 'a':
 			rtargv[rtargc++] = argv[i];
 			if (argv[i][2] == 'v') {
@@ -114,7 +140,7 @@ char	*argv[];
 	}
 						/* get view */
 	if (picture != NULL) {
-		rval = viewfile(picture, &pictview, 0, 0);
+		rval = viewfile(picture, &pictview, NULL);
 		if (rval < 0) {
 			fprintf(stderr, "%s: cannot open picture file \"%s\"\n",
 					progname, picture);
@@ -163,14 +189,19 @@ char	*argv[];
 		exit(1);
 	}
 	init();					/* initialize program */
-	comp_thresh();				/* compute glare threshold */
+	if (threshold <= FTINY)
+		comp_thresh();			/* compute glare threshold */
 	analyze();				/* analyze view */
+	if (combine)
+		absorb_specks();		/* eliminate tiny sources */
 	cleanup();				/* tidy up */
 						/* print header */
 	printargs(argc, argv, stdout);
 	fputs(VIEWSTR, stdout);
 	fprintview(&ourview, stdout);
-	printf("\n\n");
+	printf("\n");
+	fputformat("ascii", stdout);
+	printf("\n");
 	printsources();				/* print glare sources */
 	printillum();				/* print illuminances */
 	exit(0);
@@ -182,6 +213,23 @@ userr:
 }
 
 
+int
+angcmp(ap1, ap2)		/* compare two angles */
+ANGLE	*ap1, *ap2;
+{
+	register int	a1, a2;
+
+	a1 = *ap1;
+	a2 = *ap2;
+	if (a1 == a2) {
+		fprintf(stderr, "%s: duplicate glare angle (%d)\n",
+				progname, a1);
+		exit(1);
+	}
+	return(a1-a2);
+}
+
+
 init()				/* initialize global variables */
 {
 	double	d;
@@ -193,46 +241,50 @@ init()				/* initialize global variables */
 						/* set direction vectors */
 	for (i = 0; glarang[i] != AEND; i++)
 		;
-	if (glarang[0] <= 0 || glarang[i-1] >= 180) {
-		fprintf(stderr, "%s: glare angles must be between 1 and 179\n",
+	qsort(glarang, i, sizeof(ANGLE), angcmp);
+	if (i > 0 && (glarang[0] <= 0 || glarang[i-1] > 180)) {
+		fprintf(stderr, "%s: glare angles must be between 1 and 180\n",
 				progname);
 		exit(1);
 	}
 	nglarangs = i;
 	/* nglardirs = 2*nglarangs + 1; */
-	/* maxtheta = (PI/180.)*glarang[nglarangs-1]; */
-	/* vsize = SAMPDENS; */
-	hlim = SAMPDENS*maxtheta;
-	hsize = SAMPDENS + hlim;
-	if (hsize > (int)(PI*SAMPDENS))
-		hsize = PI*SAMPDENS;
+	/* vsize = sampdens - 1; */
+	if (nglarangs > 0)
+		maxtheta = (PI/180.)*glarang[nglarangs-1];
+	else
+		maxtheta = 0.0;
+	hsize = hlim(0) + 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);
+	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].lsin =
-	-(indirect[nglarangs].rsin = sin(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].lsin =
-		-(indirect[nglarangs+i+1].rsin = sin(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].lsin =
-		-(indirect[nglarangs-i-1].rsin = sin(d));
-		indirect[nglarangs-i-1].theta = -(PI/180.)*glarang[i];
-		indirect[nglarangs+i+1].theta = (PI/180.)*glarang[i];
+		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) {
@@ -257,11 +309,14 @@ init()				/* initialize global variables */
 cleanup()				/* close files, wait for children */
 {
 	if (verbose)
-		fprintf(stderr, "%s: cleaning up...\n", progname);
+		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 %d%% of samples\n",
+				progname, (int)(100L*npixmiss/npixinvw));
 }
 
 
@@ -269,47 +324,57 @@ compdir(vd, x, y)			/* compute direction for x,y */
 FVECT	vd;
 int	x, y;
 {
+	int	hl;
 	FVECT	org;			/* dummy variable */
 
-	if (x <= -hlim)			/* left region */
+	hl = hlim(y);
+	if (x <= -hl) {			/* left region */
+		if (x <= -hl-sampdens)
+			return(-1);
 		return(viewray(org, vd, &leftview,
-				(x+hlim)/(2.*SAMPDENS)+.5,
-				y/(2.*SAMPDENS)+.5));
-	if (x >= hlim)			/* right region */
+				(double)(x+hl)/(2*sampdens)+.5,
+				(double)y/(2*sampdens)+.5));
+	}
+	if (x >= hl) {			/* right region */
+		if (x >= hl+sampdens)
+			return(-1);
 		return(viewray(org, vd, &rightview,
-				(x-hlim)/(2.*SAMPDENS)+.5,
-				y/(2.*SAMPDENS)+.5));
-						/* central region */
-	if (viewray(org, vd, &ourview, .5, y/(2.*SAMPDENS)+.5) < 0)
+				(double)(x-hl)/(2*sampdens)+.5,
+				(double)y/(2*sampdens)+.5));
+	}
+					/* central region */
+	if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0)
 		return(-1);
-	spinvector(vd, vd, ourview.vup, h_theta(x));
+	spinvector(vd, vd, ourview.vup, h_theta(x,y));
 	return(0);
 }
 
 
-spinvector(vres, vorig, vnorm, theta)	/* rotate vector around normal */
-FVECT  vres, vorig, vnorm;
-double  theta;
+double
+pixsize(x, y)		/* return the solid angle of pixel at (x,y) */
+int	x, y;
 {
-	extern double  sin(), cos();
-	double  sint, cost, dotp;
-	FVECT  vperp;
-	register int  i;
-	
-	sint = sin(theta);
-	cost = cos(theta);
-	dotp = DOT(vorig, vnorm);
-	fcross(vperp, vnorm, vorig);
-	for (i = 0; i < 3; i++)
-		vres[i] = vnorm[i]*dotp*(1.-cost) +
-				vorig[i]*cost + vperp[i]*sint;
+	register int	hl, xo;
+	double	disc;
+
+	hl = hlim(y);
+	if (x < -hl)
+		xo = x+hl;
+	else if (x > hl)
+		xo = x-hl;
+	else
+		xo = 0;
+	disc = 1. - (double)(xo*xo + y*y)/(sampdens*sampdens);
+	if (disc <= FTINY)
+		return(0.);
+	return(1./(sampdens*sampdens*sqrt(disc)));
 }
 
 
 memerr(s)			/* malloc failure */
 char	*s;
 {
-	fprintf(stderr, "%s: out of memory for %s\n", s);
+	fprintf(stderr, "%s: out of memory for %s\n", progname, s);
 	exit(1);
 }
 
@@ -333,7 +398,8 @@ printillum()			/* print out indirect illuminances */
 
 	printf("BEGIN indirect illuminance\n");
 	for (i = 0; i < nglardirs; i++)
-		printf("\t%f\t%f\n", (180.0/PI)*indirect[i].theta,
-				PI * indirect[i].sum / (double)indirect[i].n);
-	printf("END indirect illuminances\n");
+		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");
 }