--- ray/src/util/findglare.c	1991/04/05 14:26:47	1.15
+++ ray/src/util/findglare.c	2014/07/29 21:36:08	2.14
@@ -1,9 +1,6 @@
-/* Copyright (c) 1991 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: findglare.c,v 2.14 2014/07/29 21:36:08 greg Exp $";
 #endif
-
 /*
  * Find glare sources in a scene or image.
  *
@@ -12,12 +9,12 @@ 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"};
-int	rtargc = 4;
+char	*rtargv[64] = {"rtrace", "-h-", "-ov", "-fff", "-ld-", "-i-", "-I-"};
+int	rtargc = 7;
 
 VIEW	ourview = STDVIEW;		/* our view */
 VIEW	pictview = STDVIEW;		/* picture view */
@@ -36,18 +33,27 @@ ANGLE	glarang[180] = {AEND};		/* glare calculation ang
 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 */
 
+static int angcmp(const void *ap1, const void *ap2);
+static void init(void);
+static void cleanup(void);
+static void printsources(void);
+static void printillum(void);
 
-main(argc, argv)
-int	argc;
-char	*argv[];
+
+
+int
+main(
+	int	argc,
+	char	*argv[]
+)
 {
+	int	combine = 1;
 	int	gotview = 0;
 	int	rval, i;
 	char	*err;
@@ -55,6 +61,14 @@ char	*argv[];
 	progname = argv[0];
 					/* process options */
 	for (i = 1; i < argc && argv[i][0] == '-'; i++) {
+						/* expand arguments */
+		while ((rval = expandarg(&argc, &argv, i)) > 0)
+			;
+		if (rval < 0) {
+			fprintf(stderr, "%s: cannot expand '%s'\n",
+					argv[0], argv[i]);
+			exit(1);
+		}
 		rval = getviewopt(&ourview, argc-i, argv+i);
 		if (rval >= 0) {
 			i += rval;
@@ -75,7 +89,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",
@@ -101,11 +115,29 @@ char	*argv[];
 		case 'p':
 			picture = argv[++i];
 			break;
+		case 'c':
+			combine = !combine;
+			break;
 		case 'd':
+			rtargv[rtargc++] = argv[i];
+			if (argv[i][2] != 'v')
+				rtargv[rtargc++] = argv[++i];
+			break;
 		case 'l':
+			if (argv[i][2] == 'd')
+				break;
+			/* FALL THROUGH */
+		case 's':
+		case 'P':
+		case 'n':
 			rtargv[rtargc++] = argv[i];
 			rtargv[rtargc++] = argv[++i];
 			break;
+		case 'w':
+		case 'u':
+		case 'b':
+			rtargv[rtargc++] = argv[i];
+			break;
 		case 'a':
 			rtargv[rtargc++] = argv[i];
 			if (argv[i][2] == 'v') {
@@ -114,6 +146,14 @@ char	*argv[];
 			}
 			rtargv[rtargc++] = argv[++i];
 			break;
+		case 'm':
+			rtargv[rtargc++] = argv[i];
+			if ((argv[i][2] == 'e') | (argv[i][2] == 'a')) {
+				rtargv[rtargc++] = argv[++i];
+				rtargv[rtargc++] = argv[++i];
+			}
+			rtargv[rtargc++] = argv[++i];
+			break;
 		default:
 			goto userr;
 		}
@@ -127,7 +167,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);
@@ -154,7 +194,7 @@ char	*argv[];
 					progname);
 			exit(1);
 		}
-		copystruct(&ourview, &pictview);
+		ourview = pictview;
 	} else if (picture != NULL && !VEQ(ourview.vp, pictview.vp)) {
 		fprintf(stderr, "%s: picture must have same viewpoint\n",
 				progname);
@@ -179,12 +219,17 @@ char	*argv[];
 	if (threshold <= FTINY)
 		comp_thresh();			/* compute glare threshold */
 	analyze();				/* analyze view */
+	if (combine)
+		absorb_specks();		/* eliminate tiny sources */
 	cleanup();				/* tidy up */
 						/* print header */
+	newheader("RADIANCE", stdout);
 	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);
@@ -196,8 +241,28 @@ userr:
 }
 
 
-init()				/* initialize global variables */
+static int
+angcmp(		/* compare two angles */
+	const void	*ap1,
+	const void	*ap2
+)
 {
+	register int	a1, a2;
+
+	a1 = *(ANGLE *)ap1;
+	a2 = *(ANGLE *)ap2;
+	if (a1 == a2) {
+		fprintf(stderr, "%s: duplicate glare angle (%d)\n",
+				progname, a1);
+		exit(1);
+	}
+	return(a1-a2);
+}
+
+
+static void
+init(void)				/* initialize global variables */
+{
 	double	d;
 	register int	i;
 
@@ -207,8 +272,9 @@ init()				/* initialize global variables */
 						/* 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",
+	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);
 	}
@@ -219,16 +285,15 @@ init()				/* initialize global variables */
 		maxtheta = (PI/180.)*glarang[nglarangs-1];
 	else
 		maxtheta = 0.0;
-	hlim = sampdens*maxtheta;
-	hsize = hlim + sampdens - 1;
+	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);
+	leftview = ourview;
+	rightview = ourview;
 	spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta);
 	spinvector(rightview.vdir, ourview.vdir, ourview.vup, -maxtheta);
 	setview(&leftview);
@@ -272,64 +337,89 @@ init()				/* initialize global variables */
 }
 
 
-cleanup()				/* close files, wait for children */
+static void
+cleanup(void)				/* 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 %ld%% of samples\n",
-				progname, 100L*npixmiss/npixinvw);
+		fprintf(stderr, "%s: warning -- missing %d%% of samples\n",
+				progname, (int)(100L*npixmiss/npixinvw));
 }
 
 
-compdir(vd, x, y)			/* compute direction for x,y */
-FVECT	vd;
-int	x, y;
+extern int
+compdir(			/* compute direction for x,y */
+	FVECT	vd,
+	int	x,
+	int	y
+)
 {
-	static int	cury = 10000;
-	static double	err, cmpval;
-	long	t;
+	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,
-				(double)(x+hlim)/(2*sampdens)+.5,
+				(double)(x+hl)/(2*sampdens)+.5,
 				(double)y/(2*sampdens)+.5));
-	if (x >= hlim)			/* right region */
+	}
+	if (x >= hl) {			/* right region */
+		if (x >= hl+sampdens)
+			return(-1);
 		return(viewray(org, vd, &rightview,
-				(double)(x-hlim)/(2*sampdens)+.5,
+				(double)(x-hl)/(2*sampdens)+.5,
 				(double)y/(2*sampdens)+.5));
-					/* central region */
-				/* avoid over-counting of poles */
-	if (cury != y) {
-		err = 0.0;
-		cmpval = sqrt(1.0 - (double)((long)y*y)/((long)vsize*vsize));
-		cury = y;
 	}
-	err += cmpval;
-	if (err <= 0.5)
-		return(-1);
-	err -= 1.0;
+					/* 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);
 }
 
 
-memerr(s)			/* malloc failure */
-char	*s;
+extern double
+pixsize(		/* return the solid angle of pixel at (x,y) */
+	int	x,
+	int	y
+)
 {
+	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)((long)xo*xo + (long)y*y)/((long)sampdens*sampdens);
+	if (disc <= FTINY*FTINY)
+		return(0.);
+	return(1./(sampdens*sampdens*sqrt(disc)));
+}
+
+
+extern void
+memerr(			/* malloc failure */
+	char	*s
+)
+{
 	fprintf(stderr, "%s: out of memory for %s\n", progname, s);
 	exit(1);
 }
 
 
-printsources()			/* print out glare sources */
+static void
+printsources(void)			/* print out glare sources */
 {
 	register struct source	*sp;
 
@@ -342,7 +432,8 @@ printsources()			/* print out glare sources */
 }
 
 
-printillum()			/* print out indirect illuminances */
+static void
+printillum(void)			/* print out indirect illuminances */
 {
 	register int	i;