--- ray/src/util/rmtxop.c	2023/11/28 16:36:50	2.23
+++ ray/src/util/rmtxop.c	2023/12/03 02:28:33	2.28
@@ -1,13 +1,11 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rmtxop.c,v 2.23 2023/11/28 16:36:50 greg Exp $";
+static const char RCSid[] = "$Id: rmtxop.c,v 2.28 2023/12/03 02:28:33 greg Exp $";
 #endif
 /*
  * General component matrix operations.
  */
 
-#include <stdlib.h>
 #include <errno.h>
-#include <ctype.h>
 #include "rtio.h"
 #include "resolu.h"
 #include "rmatrix.h"
@@ -17,12 +15,12 @@ static const char RCSid[] = "$Id: rmtxop.c,v 2.23 2023
 
 /* Unary matrix operation(s) */
 typedef struct {
-	double		sca[MAXCOMP];		/* scalar coefficients */
 	double		cmat[MAXCOMP*MAXCOMP];	/* component transformation */
-	short		nsf;			/* number of scalars */
+	double		sca[MAXCOMP];		/* scalar coefficients */
+	const char	*csym;			/* symbolic coefs or file */
 	short		clen;			/* number of coefficients */
-	char		csym[11];		/* symbolic coefficients */
-	char		transpose;		/* do transpose? */
+	short		nsf;			/* number of scalars */
+	short		transpose;		/* do transpose? */
 } RUNARYOP;
 
 /* Matrix input source and requested operation(s) */
@@ -48,6 +46,66 @@ loadmatrix(ROPMAT *rop)
 	return(!rop->mtx ? -1 : 1);
 }
 
+static int	checksymbolic(ROPMAT *rop);
+
+/* Check/set transform based on a reference input file */
+static int
+checkreffile(ROPMAT *rop)
+{
+	static const char	*curRF = NULL;
+	static RMATRIX		refm;
+	const int		nc = rop->mtx->ncomp;
+	int			i;
+
+	if (!curRF || strcmp(rop->preop.csym, curRF)) {
+		FILE	*fp = fopen(rop->preop.csym, "rb");
+		if (!rmx_load_header(&refm, fp)) {
+			fprintf(stderr, "%s: cannot read info header\n",
+					rop->preop.csym);
+			curRF = NULL;
+			if (fp) fclose(fp);
+			return(-1);
+		}
+		fclose(fp);
+		curRF = rop->preop.csym;
+	}
+	if ((refm.ncomp == 3) & (refm.dtype != DTspec)) {
+		rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
+		return(checksymbolic(rop));
+	}
+	if (refm.ncomp == 2) {
+		fprintf(stderr, "%s: cannot convert to 2 components\n",
+				curRF);
+		return(-1);
+	}
+	if (refm.ncomp == 1) {
+		rop->preop.csym = "Y";		/* XXX big assumption */
+		return(checksymbolic(rop));
+	}
+	if (refm.ncomp == nc &&
+			!memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
+		return(0);			/* nothing to do */
+
+	if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
+		fprintf(stderr, "%s: cannot resample from %d to %d components\n",
+				curRF, nc, refm.ncomp);
+		return(-1);
+	}
+	rop->preop.clen = refm.ncomp * nc;	/* compute spec to ref */
+
+	for (i = 0; i < nc; i++) {
+		SCOLOR	scstim, scresp;
+		int	j;
+		memset(scstim, 0, sizeof(COLORV)*nc);
+		scstim[i] = 1.f;
+		convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
+				scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
+		for (j = refm.ncomp; j-- > 0; )
+			rop->preop.cmat[j*nc + i] = scresp[j];
+	}
+	return(0);
+}
+
 /* Compute conversion row from spectrum to one channel of RGB */
 static void
 rgbrow(ROPMAT *rop, int r, int p)
@@ -91,8 +149,8 @@ sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, in
 
 	for (i = nc; i--; ) {
 		SCOLOR	sclr;
-		scolorblack(sclr);
-		sclr[i] = 1;
+		memset(sclr, 0, sizeof(COLORV)*nc);
+		sclr[i] = 1.f;
 		rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
 	}
 }
@@ -104,6 +162,9 @@ checksymbolic(ROPMAT *rop)
 	const int	nc = rop->mtx->ncomp;
 	const int	dt = rop->mtx->dtype;
 	int		i, j;
+					/* check suffix => reference file */
+	if (strchr(rop->preop.csym, '.') > rop->preop.csym)
+		return(checkreffile(rop));
 
 	if (nc < 3) {
 		fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
@@ -156,16 +217,20 @@ checksymbolic(ROPMAT *rop)
 			for (i = nc*(dt != DTxyze); i--; )
 				rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
 			break;
-		case 'S':
+		case 'S':		/* scotopic (il)luminance */
 			sensrow(rop, j, scolor2scotopic);
 			for (i = nc; i--; )
 				rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
 			break;
-		case 'M':
+		case 'M':		/* melanopic (il)luminance */
 			sensrow(rop, j, scolor2melanopic);
 			for (i = nc; i--; )
 				rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
 			break;
+		case 'A':		/* average component */
+			for (i = nc; i--; )
+				rop->preop.cmat[j*nc+i] = 1./(double)nc;
+			break;
 		default:
 			fprintf(stderr, "%s: -c '%c' unsupported\n",
 				rop->inspec, rop->preop.csym[j]);
@@ -196,7 +261,7 @@ loadop(ROPMAT *rop)
 	if (loadmatrix(rop) < 0)		/* make sure we're loaded */
 		return(NULL);
 
-	if (rop->preop.csym[0] &&		/* symbolic transform? */
+	if (rop->preop.csym &&			/* symbolic transform? */
 			(outtype = checksymbolic(rop)) < 0)
 		goto failure;
 	if (rop->preop.clen > 0) {		/* apply component transform? */
@@ -209,15 +274,16 @@ loadop(ROPMAT *rop)
 			if (rop->preop.nsf == 1) {
 				for (i = rop->preop.clen; i--; )
 					rop->preop.cmat[i] *= rop->preop.sca[0];
- 			} else if (rop->preop.nsf != rop->mtx->ncomp) {
+ 			} else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
 				fprintf(stderr, "%s: -s must have one or %d factors\n",
-						rop->inspec, rop->mtx->ncomp);
+						rop->inspec,
+						rop->preop.clen/rop->mtx->ncomp);
 				goto failure;
 			} else {
-				for (j = rop->preop.clen/rop->preop.nsf; j--; )
-					for (i = rop->preop.nsf; i--; )
-						rop->preop.cmat[j*rop->preop.nsf+i] *=
-								rop->preop.sca[i];
+				for (i = rop->preop.nsf; i--; )
+					for (j = rop->mtx->ncomp; j--; )
+						rop->preop.cmat[i*rop->mtx->ncomp+j]
+								*= rop->preop.sca[i];
 			}
 		}
 		mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
@@ -232,7 +298,7 @@ loadop(ROPMAT *rop)
 					rop->inspec, mres->ncomp,
 					rop->mtx->ncomp,
 					rop->preop.nsf ? " (* scalar)" : "");
-		rop->preop.nsf = 0;
+		rop->preop.nsf = 0;		/* now folded in */
 		if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
 			outtype = DTfloat;	/* probably not actual spectrum */
 		rmx_free(rop->mtx);
@@ -445,15 +511,18 @@ get_factors(double da[], int n, char *av[])
 }
 
 static ROPMAT *
-grow_moparray(ROPMAT *mop, int n2alloc)
+resize_moparr(ROPMAT *mop, int n2alloc)
 {
 	int	nmats = 0;
+	int	i;
 
 	while (mop[nmats++].binop)
 		;
+	for (i = nmats; i > n2alloc; i--)
+		rmx_free(mop[i].mtx);
 	mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
 	if (mop == NULL) {
-		fputs("Out of memory in grow_moparray()\n", stderr);
+		fputs("Out of memory in resize_moparr()\n", stderr);
 		exit(1);
 	}
 	if (n2alloc > nmats)
@@ -465,13 +534,14 @@ grow_moparray(ROPMAT *mop, int n2alloc)
 int
 main(int argc, char *argv[])
 {
-	int	outfmt = DTfromHeader;
-	int	nall = 2;
-	ROPMAT	*mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
-	int	nmats = 0;
-	RMATRIX	*mres = NULL;
-	int	stdin_used = 0;
-	int	i;
+	int		outfmt = DTfromHeader;
+	const char	*defCsym = NULL;
+	int		nall = 2;
+	ROPMAT		*mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
+	int		nmats = 0;
+	RMATRIX		*mres = NULL;
+	int		stdin_used = 0;
+	int		i;
 					/* get options and arguments */
 	for (i = 1; i < argc; i++) {
 		if (argv[i][0] && !argv[i][1] &&
@@ -494,6 +564,8 @@ main(int argc, char *argv[])
 				mop[nmats].inspec = stdin_name;
 			} else
 				mop[nmats].inspec = argv[i];
+			if (!mop[nmats].preop.csym)
+				mop[nmats].preop.csym = defCsym;
 			if (nmats > 0 && !mop[nmats-1].binop)
 				mop[nmats-1].binop = '.';
 			nmats++;
@@ -535,11 +607,15 @@ main(int argc, char *argv[])
 					goto userr;
 				}
 				break;
+			case 'C':
+				if (!n || isflt(argv[i+1]))
+					goto userr;
+				defCsym = mop[nmats].preop.csym = argv[++i];
+				mop[nmats].preop.clen = 0;
+				break;
 			case 'c':
-				if (n && isupper(argv[i+1][0])) {
-					strlcpy(mop[nmats].preop.csym,
-						argv[++i],
-						sizeof(mop[0].preop.csym));
+				if (n && !isflt(argv[i+1])) {
+					mop[nmats].preop.csym = argv[++i];
 					mop[nmats].preop.clen = 0;
 					break;
 				}
@@ -552,7 +628,7 @@ main(int argc, char *argv[])
 							argv[0]);
 					goto userr;
 				}
-				mop[nmats].preop.csym[0] = '\0';
+				mop[nmats].preop.csym = NULL;
 				break;
 			case 'r':
 				if (argv[i][2] == 'f')
@@ -569,7 +645,7 @@ main(int argc, char *argv[])
 			}
 		}
 		if (nmats >= nall)
-			mop = grow_moparray(mop, nall += 2);
+			mop = resize_moparr(mop, nall += 2);
 	}
 	if (mop[0].inspec == NULL)	/* nothing to do? */
 		goto userr;
@@ -597,15 +673,12 @@ main(int argc, char *argv[])
 		else if (mres->dtype == DTxyze)
 			outfmt = DTxyze;
 	}
-	if (outfmt != DTascii)		/* write result to stdout */
-		SET_FILE_BINARY(stdout);
-	newheader("RADIANCE", stdout);
+	newheader("RADIANCE", stdout);	/* write result to stdout */
 	printargs(argc, argv, stdout);
-
 	return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
 userr:
 	fprintf(stderr,
-	"Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
+	"Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
 			argv[0]);
 	return(1);
 }