--- ray/src/util/rtcontrib.c	2007/11/17 01:13:51	1.45
+++ ray/src/util/rtcontrib.c	2010/09/26 15:44:00	1.57
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rtcontrib.c,v 1.45 2007/11/17 01:13:51 greg Exp $";
+static const char RCSid[] = "$Id: rtcontrib.c,v 1.57 2010/09/26 15:44:00 greg Exp $";
 #endif
 /*
  * Gather rtrace output to compute contributions from particular sources
@@ -14,6 +14,7 @@ static const char RCSid[] = "$Id: rtcontrib.c,v 1.45 2
 #include  "color.h"
 #include  "resolu.h"
 #include  "lookup.h"
+#include  "random.h"
 #include  "calcomp.h"
 
 #ifndef	MAXMODLIST
@@ -99,7 +100,7 @@ struct rtproc {
 
 					/* rtrace command and defaults */
 char		*rtargv[256+2*MAXMODLIST] = { "rtrace",
-				"-dj", ".5", "-dr", "3",
+				"-dj", ".9", "-dr", "3",
 				"-ab", "1", "-ad", "350", };
 
 int  rtargc = 9;
@@ -131,11 +132,13 @@ int		inpfmt = 'a';		/* input format */
 int		outfmt = 'a';		/* output format */
 
 int		header = 1;		/* output header? */
-int		accumulate = 0;		/* accumulate ray values? */
 int		force_open = 0;		/* truncate existing output? */
+int		recover = 0;		/* recover previous output? */
+int		accumulate = 1;		/* input rays per output record */
 int		xres = 0;		/* horiz. output resolution */
 int		yres = 0;		/* vert. output resolution */
 
+int		account;		/* current accumulation count */
 unsigned long	raysleft;		/* number of rays left to trace */
 long		waitflush;		/* how long until next flush */
 
@@ -165,7 +168,7 @@ void process_queue(void);
 
 void put_contrib(const DCOLOR cnt, FILE *fout);
 void add_contrib(const char *modn);
-void done_contrib(void);
+void done_contrib(int navg);
 
 /* return number of open rtrace processes */
 static int
@@ -220,7 +223,6 @@ main(int argc, char *argv[])
 {
 	int	contrib = 0;
 	int	nprocs = 1;
-	int	recover = 0;
 	char	*curout = NULL;
 	char	*binval = NULL;
 	int	bincnt = 0;
@@ -245,7 +247,7 @@ main(int argc, char *argv[])
 		while ((j = expandarg(&argc, &argv, i)) > 0)
 			;
 		if (j < 0) {
-			fprintf(stderr, "%s: cannot expand '%s'",
+			fprintf(stderr, "%s: cannot expand '%s'\n",
 					argv[0], argv[i]);
 			exit(1);
 		}
@@ -274,26 +276,13 @@ main(int argc, char *argv[])
 					continue;
 				}
 				break;
-			case 'c':		/* accumulate ray values */
-				switch (argv[i][2]) {
-				case '\0':
-					accumulate = !accumulate;
-					continue;
-				case '+': case '1':
-				case 'T': case 't':
-				case 'Y': case 'y':
-					accumulate = 1;
-					continue;
-				case '-': case '0':
-				case 'F': case 'f':
-				case 'N': case 'n':
-					accumulate = 0;
-					continue;
-				}
-				break;
+			case 'c':		/* input rays per output */
+				if (argv[i][2] || i >= argc-2) break;
+				accumulate = atoi(argv[++i]);
+				continue;
 			case 'r':		/* recover output */
 				if (argv[i][2]) break;
-				recover++;
+				recover = 1;
 				continue;
 			case 'h':		/* output header? */
 				switch (argv[i][2]) {
@@ -328,7 +317,7 @@ main(int argc, char *argv[])
 					continue;
 				}
 				if (argv[i][2] == 'o') {
-					force_open++;
+					force_open = 1;
 					continue;
 				}
 				setformat(argv[i]+2);
@@ -352,7 +341,7 @@ main(int argc, char *argv[])
 			case 'b':		/* bin expression/count */
 				if (i >= argc-2) break;
 				if (argv[i][2] == 'n') {
-					bincnt = atoi(argv[++i]);
+					bincnt = (int)(eval(argv[++i]) + .5);
 					continue;
 				}
 				if (argv[i][2]) break;
@@ -380,7 +369,7 @@ main(int argc, char *argv[])
 			}
 		rtargv[rtargc++] = argv[i];	/* assume rtrace option */
 	}
-	if (accumulate)		/* no output flushing for single record */
+	if (accumulate <= 0)	/* no output flushing for single record */
 		xres = yres = 0;
 				/* set global argument list */
 	gargc = argc; gargv = argv;
@@ -390,15 +379,16 @@ main(int argc, char *argv[])
 	rtargv[rtargc++] = contrib ? RTCONTRIB : RTCOEFF;
 				/* just asking for defaults? */
 	if (!strcmp(argv[i], "-defaults")) {
-		char	sxres[16], syres[16];
+		char	nps[8], sxres[16], syres[16];
 		char	*rtpath;
-		printf("-n  %-2d\t\t\t\t# number of processes\n", nprocs);
+		printf("-c %-5d\t\t\t# accumulated rays per record\n",
+				accumulate);
 		printf("-V%c\t\t\t\t# output %s\n", contrib ? '+' : '-',
 				contrib ? "contributions" : "coefficients");
-		printf("-c%c\t\t\t\t# %s\n", accumulate ? '+' : '-',
-				accumulate ? "accumulate ray values" :
-					"one output record per ray");
 		fflush(stdout);			/* report OUR options */
+		rtargv[rtargc++] = "-n";
+		sprintf(nps, "%d", nprocs);
+		rtargv[rtargc++] = nps;
 		rtargv[rtargc++] = header ? "-h+" : "-h-";
 		sprintf(fmt, "-f%c%c", inpfmt, outfmt);
 		rtargv[rtargc++] = fmt;
@@ -437,14 +427,11 @@ main(int argc, char *argv[])
 		error(USER, "missing octree argument");
 	rtargv[rtargc++] = octree = argv[i];
 	rtargv[rtargc] = NULL;
-				/* start rtrace */
+				/* start rtrace & recover if requested */
 	init(nprocs);
-	if (recover)		/* perform recovery if requested */
-		if ((force_open = accumulate))
-			reload_output();
-		else
-			recover_output(stdin);
-	trace_contribs(stdin);	/* compute contributions */
+				/* compute contributions */
+	trace_contribs(stdin);
+				/* clean up */
 	quit(0);
 }
 
@@ -557,7 +544,16 @@ init(int np)
 			raysleft = yres;
 	} else
 		raysleft = 0;
-	waitflush = xres;
+	if ((account = accumulate) > 0)
+		raysleft *= accumulate;
+	waitflush = (yres > 0) & (xres > 1) ? 0 : xres;
+	if (!recover)
+		return;
+					/* recover previous values */
+	if (accumulate <= 0)
+		reload_output();
+	else
+		recover_output(stdin);
 }
 
 /* grow modifier to accommodate more bins */
@@ -587,7 +583,7 @@ addmodifier(char *modn, char *outf, char *binv, int bi
 		error(USER, errmsg);
 	}
 	if (nmods >= MAXMODLIST)
-		error(USER, "too many modifiers");
+		error(INTERNAL, "too many modifiers");
 	modname[nmods++] = modn;	/* XXX assumes static string */
 	lep->key = modn;		/* XXX assumes static string */
 	mp = (MODCONT *)malloc(sizeof(MODCONT));
@@ -595,15 +591,20 @@ addmodifier(char *modn, char *outf, char *binv, int bi
 		error(SYSTEM, "out of memory in addmodifier");
 	mp->outspec = outf;		/* XXX assumes static string */
 	mp->modname = modn;		/* XXX assumes static string */
-	if (binv != NULL)
-		mp->binv = eparse(binv);
-	else
-		mp->binv = eparse("0");
-	mp->nbins = 1;
+	if (binv == NULL)
+		binv = "0";		/* use single bin if unspecified */
+	mp->binv = eparse(binv);
+	if (mp->binv->type == NUM) {	/* check value if constant */
+		bincnt = (int)(evalue(mp->binv) + 1.5);
+		if (bincnt != 1) {
+			sprintf(errmsg, "illegal non-zero constant for bin (%s)",
+					binv);
+			error(USER, errmsg);
+		}
+	}
+	mp->nbins = 1;			/* initialize results holder */
 	setcolor(mp->cbin[0], 0., 0., 0.);
-	if (mp->binv->type == NUM)	/* assume one bin if constant */
-		bincnt = 1;
-	else if (bincnt > 1)
+	if (bincnt > 1)
 		mp = growmodifier(mp, bincnt);
 	lep->data = (char *)mp;
 					/* allocate output streams */
@@ -666,6 +667,10 @@ ofname(char *oname, const char *ospec, const char *mna
 				mnp = cp;
 				break;
 			case 'd':
+			case 'i':
+			case 'o':
+			case 'x':
+			case 'X':
 				if (bnp != NULL)
 					return -1;
 				bnp = cp;
@@ -732,8 +737,6 @@ printresolu(FILE *fout, int xr, int yr)
 {
 	if ((xr > 0) & (yr > 0))	/* resolution string */
 		fprtresolu(xr, yr, fout);
-	if (xres > 0)			/* global flush flag */
-		fflush(fout);
 }
 
 /* Get output stream pointer (open and write header if new and noopen==0) */
@@ -754,6 +757,8 @@ getostream(const char *ospec, const char *mname, int b
 			if (header)
 				printheader(stdout, NULL);
 			printresolu(stdout, xres, yres);
+			if (waitflush > 0)
+				fflush(stdout);
 			stdos.xr = xres; stdos.yr = yres;
 			using_stdout = 1;
 		}
@@ -779,7 +784,8 @@ getostream(const char *ospec, const char *mname, int b
 		sop->ofp = NULL;		/* open iff noopen==0 */
 		sop->xr = xres; sop->yr = yres;
 		lep->data = (char *)sop;
-		if (!force_open && access(oname, F_OK) == 0) {
+		if (!sop->outpipe & !force_open & !recover &&
+				access(oname, F_OK) == 0) {
 			errno = EEXIST;		/* file exists */
 			goto openerr;
 		}
@@ -808,9 +814,12 @@ getostream(const char *ospec, const char *mname, int b
 			*cp = '\0';
 			printheader(sop->ofp, info);
 		}
+		if (accumulate > 0) {		/* global resolution */
+			sop->xr = xres; sop->yr = yres;
+		}
 		printresolu(sop->ofp, sop->xr, sop->yr);
 						/* play catch-up */
-		for (i = accumulate ? 0 : lastdone; i--; ) {
+		for (i = accumulate > 0 ? lastdone/accumulate : 0; i--; ) {
 			int	j = sop->reclen;
 			if (j <= 0) j = 1;
 			while (j--)
@@ -818,7 +827,7 @@ getostream(const char *ospec, const char *mname, int b
 			if (outfmt == 'a')
 				putc('\n', sop->ofp);
 		}
-		if (xres > 0)
+		if (waitflush > 0)
 			fflush(sop->ofp);
 	}
 	sop->reclen += noopen;			/* add to length if noopen */
@@ -935,9 +944,7 @@ put_contrib(const DCOLOR cnt, FILE *fout)
 		fprintf(fout, "%.6e\t%.6e\t%.6e\t", cnt[0], cnt[1], cnt[2]);
 		break;
 	case 'f':
-		fv[0] = cnt[0];
-		fv[1] = cnt[1];
-		fv[2] = cnt[2];
+		copycolor(fv, cnt);
 		fwrite(fv, sizeof(float), 3, fout);
 		break;
 	case 'd':
@@ -950,18 +957,27 @@ put_contrib(const DCOLOR cnt, FILE *fout)
 	default:
 		error(INTERNAL, "botched output format");
 	}
+	if (waitflush < 0 && frandom() < 0.001)
+		fflush(fout);			/* staggers writes */
 }
 
 /* output ray tallies and clear for next accumulation */
 void
-done_contrib(void)
+done_contrib(int navg)
 {
+	double		sf = 1.;
 	int		i, j;
 	MODCONT		*mp;
 	STREAMOUT	*sop;
+						/* set average scaling */
+	if (navg > 1)
+		sf = 1. / (double)navg;
 						/* output modifiers in order */
 	for (i = 0; i < nmods; i++) {
 		mp = (MODCONT *)lu_find(&modconttab,modname[i])->data;
+		if (navg > 1)			/* average scaling */
+			for (j = mp->nbins; j--; )
+				scalecolor(mp->cbin[j], sf);
 		sop = getostream(mp->outspec, mp->modname, 0,0);
 		put_contrib(mp->cbin[0], sop->ofp);
 		if (mp->nbins > 3 &&		/* minor optimization */
@@ -980,7 +996,7 @@ done_contrib(void)
 	if (using_stdout & (outfmt == 'a'))
 		putc('\n', stdout);
 	if (!waitflush) {
-		waitflush = xres;
+		waitflush = (yres > 0) & (xres > 1) ? 0 : xres;
 		if (using_stdout)
 			fflush(stdout);
 	}
@@ -1045,8 +1061,9 @@ process_queue(void)
 			cp += sizeof(float)*9; n -= sizeof(float)*9;
 			add_contrib(modname);
 		}
-		if (!accumulate)
-			done_contrib();	/* sum up contributions & output */
+					/* time to produce record? */
+		if (account > 0 && !--account)
+			done_contrib(account = accumulate);
 		lastdone = rtp->raynum;
 		if (rtp->buf != NULL)	/* free up buffer space */
 			free(rtp->buf);
@@ -1133,12 +1150,19 @@ get_rproc(void)
 void
 trace_contribs(FILE *fin)
 {
+	static int	ignore_warning_given = 0;
 	char		inpbuf[128];
 	int		iblen;
 	struct rtproc	*rtp;
 						/* loop over input */
 	while ((iblen = getinp(inpbuf, fin)) >= 0) {
-		if (!iblen ||			/* need reset? */
+		if (!iblen && accumulate != 1) {
+			if (!ignore_warning_given++)
+				error(WARNING,
+				"dummy ray(s) ignored during accumulation\n");
+			continue;
+		}
+		if (!iblen ||			/* need flush/reset? */
 				queue_length() > 10*nrtprocs() ||
 				lastray+1 < lastray) {
 			while (wait_rproc() != NULL)
@@ -1151,9 +1175,9 @@ trace_contribs(FILE *fin)
 		if (iblen) {			/* trace ray if valid */
 			writebuf(rtp->pd.w, inpbuf, iblen);
 		} else {			/* else bypass dummy ray */
-			queue_raytree(rtp);	/* empty tree */
-			if ((yres <= 0) | (waitflush > 1))
-				waitflush = 1;	/* flush after this */
+			queue_raytree(rtp);	/* queue empty ray/record */
+			if ((yres <= 0) | (xres <= 0))
+				waitflush = 1;	/* flush right after */
 		}
 		process_queue();		/* catch up with results */
 		if (raysleft && !--raysleft)
@@ -1161,27 +1185,57 @@ trace_contribs(FILE *fin)
 	}
 	while (wait_rproc() != NULL)		/* process outstanding rays */
 		process_queue();
-	if (accumulate)
-		done_contrib();			/* output tallies */
+	if (accumulate <= 0)
+		done_contrib(0);		/* output tallies */
+	else if (account < accumulate) {
+		error(WARNING, "partial accumulation in final record");
+		done_contrib(accumulate - account);
+	}
 	if (raysleft)
 		error(USER, "unexpected EOF on input");
 	lu_done(&ofiletab);			/* close output files */
 }
 
-/* seek on the given output file */
+/* get ray contribution from previous file */
 static int
-myseeko(const LUENT *e, void *p)
+get_contrib(DCOLOR cnt, FILE *finp)
 {
+	COLOR	fv;
+	COLR	cv;
+
+	switch (outfmt) {
+	case 'a':
+		return(fscanf(finp,"%lf %lf %lf",&cnt[0],&cnt[1],&cnt[2]) == 3);
+	case 'f':
+		if (fread(fv, sizeof(fv[0]), 3, finp) != 3)
+			return(0);
+		copycolor(cnt, fv);
+		return(1);
+	case 'd':
+		return(fread(cnt, sizeof(cnt[0]), 3, finp) == 3);
+	case 'c':
+		if (fread(cv, sizeof(cv), 1, finp) != 1)
+			return(0);
+		colr_color(fv, cv);
+		copycolor(cnt, fv);
+		return(1);
+	default:
+		error(INTERNAL, "botched output format");
+	}
+	return(0);	/* pro forma return */
+}
+
+/* close output file opened for input */
+static int
+myclose(const LUENT *e, void *p)
+{
 	STREAMOUT	*sop = (STREAMOUT *)e->data;
-	off_t		nbytes = *(off_t *)p;
 	
-	if (sop->reclen > 1)
-		nbytes = nbytes * sop->reclen;
-	if (fseeko(sop->ofp, nbytes, SEEK_CUR) < 0) {
-		sprintf(errmsg, "seek error on file '%s'", e->key);
-		error(SYSTEM, errmsg);
-	}
-	return 0;
+	if (sop->ofp == NULL)
+		return(0);
+	fclose(sop->ofp);
+	sop->ofp = NULL;
+	return(0);
 }
 
 /* load previously accumulated values */
@@ -1195,10 +1249,10 @@ reload_output(void)
 	char		*fmode = "rb";
 	char		*outvfmt;
 	LUENT		*ment, *oent;
+	int		xr, yr;
 	STREAMOUT	sout;
+	DCOLOR		rgbv;
 
-	error(INTERNAL, "-r option not yet implemented with -c");
-
 	switch (outfmt) {
 	case 'a':
 		outvfmt = "ascii";
@@ -1217,7 +1271,7 @@ reload_output(void)
 		error(INTERNAL, "botched output format");
 		return;
 	}
-						/* check modifier outputs */
+						/* reload modifier values */
 	for (i = 0; i < nmods; i++) {
 		ment = lu_find(&modconttab,modname[i]);
 		mp = (MODCONT *)ment->data;
@@ -1225,7 +1279,7 @@ reload_output(void)
 			error(USER, "cannot reload from stdout");
 		if (mp->outspec[0] == '!')
 			error(USER, "cannot reload from command");
-		for (j = 0; ; j++) {		/* check each bin's file */
+		for (j = 0; ; j++) {		/* load each modifier bin */
 			ofl = ofname(oname, mp->outspec, mp->modname, j);
 			if (ofl < 0)
 				error(USER, "bad output file specification");
@@ -1234,38 +1288,77 @@ reload_output(void)
 				sout = *(STREAMOUT *)oent->data;
 			} else {
 				sout.reclen = 0;
+				sout.outpipe = 0;
+				sout.xr = xres; sout.yr = yres;
 				sout.ofp = NULL;
 			}
-			if (sout.ofp != NULL) {	/* already open? */
-				if (ofl & OF_BIN)
-					continue;
-				break;
+			if (sout.ofp == NULL) {	/* open output as input */
+				sout.ofp = fopen(oname, fmode);
+				if (sout.ofp == NULL) {
+					if (j)
+						break;	/* assume end of modifier */
+					sprintf(errmsg, "missing reload file '%s'",
+							oname);
+					error(WARNING, errmsg);
+					break;
+				}
+				if (header && checkheader(sout.ofp, outvfmt, NULL) != 1) {
+					sprintf(errmsg, "format mismatch for '%s'",
+							oname);
+					error(USER, errmsg);
+				}
+				if ((sout.xr > 0) & (sout.yr > 0) &&
+						(!fscnresolu(&xr, &yr, sout.ofp) ||
+							(xr != sout.xr) |
+							(yr != sout.yr))) {
+					sprintf(errmsg, "resolution mismatch for '%s'",
+							oname);
+					error(USER, errmsg);
+				}
 			}
-						/* open output */
-			sout.ofp = fopen(oname, fmode);
-			if (sout.ofp == NULL) {
-				if (j)
-					break;	/* assume end of modifier */
-				sprintf(errmsg, "missing reload file '%s'",
-						oname);
-				error(WARNING, errmsg);
+							/* read in RGB value */
+			if (!get_contrib(rgbv, sout.ofp)) {
+				if (!j) {
+					fclose(sout.ofp);
+					break;		/* ignore empty file */
+				}
+				if (j < mp->nbins) {
+					sprintf(errmsg, "missing data in '%s'",
+							oname);
+					error(USER, errmsg);
+				}
 				break;
 			}
-			if (oent->key == NULL)	/* new entry */
+			if (j >= mp->nbins)		/* grow modifier size */
+				ment->data = (char *)(mp = growmodifier(mp, j+1));
+			copycolor(mp->cbin[j], rgbv);
+			if (oent->key == NULL)		/* new file entry */
 				oent->key = strcpy((char *)
 						malloc(strlen(oname)+1), oname);
 			if (oent->data == NULL)
 				oent->data = (char *)malloc(sizeof(STREAMOUT));
 			*(STREAMOUT *)oent->data = sout;
-			if (!(ofl & OF_BIN))
-				break;		/* no bin separation */
 		}
-		if (j > mp->nbins)		/* reallocate modifier bins */
-			ment->data = (char *)(mp = growmodifier(mp, j));
 	}
-	lu_done(&ofiletab);			/* close all files */
+	lu_doall(&ofiletab, myclose, NULL);	/* close all files */
 }
 
+/* seek on the given output file */
+static int
+myseeko(const LUENT *e, void *p)
+{
+	STREAMOUT	*sop = (STREAMOUT *)e->data;
+	off_t		nbytes = *(off_t *)p;
+	
+	if (sop->reclen > 1)
+		nbytes = nbytes * sop->reclen;
+	if (fseeko(sop->ofp, nbytes, SEEK_CUR) < 0) {
+		sprintf(errmsg, "seek error on file '%s'", e->key);
+		error(SYSTEM, errmsg);
+	}
+	return 0;
+}
+
 /* recover output if possible */
 void
 recover_output(FILE *fin)
@@ -1319,6 +1412,7 @@ recover_output(FILE *fin)
 				sout = *(STREAMOUT *)oent->data;
 			} else {
 				sout.reclen = 0;
+				sout.outpipe = 0;
 				sout.ofp = NULL;
 			}
 			if (sout.ofp != NULL) {	/* already open? */
@@ -1359,10 +1453,11 @@ recover_output(FILE *fin)
 						oname);
 				error(USER, errmsg);
 			}
-			if ((xres > 0) & (yres > 0) &&
+			sout.xr = xres; sout.yr = yres;
+			if ((sout.xr > 0) & (sout.yr > 0) &&
 					(!fscnresolu(&xr, &yr, sout.ofp) ||
-						xr != xres ||
-						yr != yres)) {
+						(xr != sout.xr) |
+						(yr != sout.yr))) {
 				sprintf(errmsg, "resolution mismatch for '%s'",
 						oname);
 				error(USER, errmsg);
@@ -1391,7 +1486,7 @@ recover_output(FILE *fin)
 		error(WARNING, "no output files to recover");
 		return;
 	}
-	if (raysleft && lastout >= raysleft) {
+	if (raysleft && lastout >= raysleft/accumulate) {
 		error(WARNING, "output appears to be complete");
 		/* XXX should read & discard input? */
 		quit(0);
@@ -1403,7 +1498,7 @@ recover_output(FILE *fin)
 	for (nvals = 0; nvals < lastout; nvals++)
 		if (getinp(oname, fin) < 0)
 			error(USER, "unexpected EOF on input");
-	lastray = lastdone = (unsigned long)lastout;
+	lastray = lastdone = (unsigned long)lastout * accumulate;
 	if (raysleft)
 		raysleft -= lastray;
 }