--- ray/src/util/rmtxop.c	2019/08/12 01:20:26	2.13
+++ ray/src/util/rmtxop.c	2023/12/02 00:42:21	2.27
@@ -1,34 +1,32 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rmtxop.c,v 2.13 2019/08/12 01:20:26 greg Exp $";
+static const char RCSid[] = "$Id: rmtxop.c,v 2.27 2023/12/02 00:42:21 greg Exp $";
 #endif
 /*
  * General component matrix operations.
  */
 
-#include <stdio.h>
-#include <stdlib.h>
 #include <errno.h>
 #include "rtio.h"
 #include "resolu.h"
 #include "rmatrix.h"
 #include "platform.h"
 
-#define MAXCOMP		16		/* #components we support */
+#define MAXCOMP		MAXCSAMP	/* #components we support */
 
-static const char	stdin_name[] = "<stdin>";
-
-/* unary matrix operation(s) */
+/* 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 */
+	short		nsf;			/* number of scalars */
 	short		transpose;		/* do transpose? */
 } RUNARYOP;
 
-/* matrix input source and requested operation(s) */
+/* Matrix input source and requested operation(s) */
 typedef struct {
 	const char	*inspec;		/* input specification */
+	RMPref		rmp;			/* matrix preference */
 	RUNARYOP	preop;			/* unary operation(s) */
 	RMATRIX		*mtx;			/* original matrix if loaded */
 	int		binop;			/* binary op with next (or 0) */
@@ -40,34 +38,272 @@ int	verbose = 0;			/* verbose reporting? */
 static int
 loadmatrix(ROPMAT *rop)
 {
-	if (rop->mtx != NULL)
+	if (rop->mtx != NULL)		/* already loaded? */
 		return(0);
 
-	rop->mtx = rmx_load(rop->inspec == stdin_name ?
-				(const char *)NULL : rop->inspec);
-	if (rop->mtx == NULL) {
-		fputs(rop->inspec, stderr);
-		fputs(": cannot load matrix\n", stderr);
+	rop->mtx = rmx_load(rop->inspec, rop->rmp);
+
+	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);
 	}
-	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)
+{
+	const int	nc = rop->mtx->ncomp;
+	const float *	wlp = rop->mtx->wlpart;
+	int		i;
+
+	for (i = nc; i--; ) {
+		int	nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
+		int	nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
+		COLOR	crgb;
+		spec_rgb(crgb, nmStart, nmEnd);
+		rop->preop.cmat[r*nc+i] = crgb[p];
+	}
+}
+
+/* Compute conversion row from spectrum to one channel of XYZ */
+static void
+xyzrow(ROPMAT *rop, int r, int p)
+{
+	const int	nc = rop->mtx->ncomp;
+	const float *	wlp = rop->mtx->wlpart;
+	int		i;
+
+	for (i = nc; i--; ) {
+		int	nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
+		int	nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
+		COLOR	cxyz;
+		spec_cie(cxyz, nmStart, nmEnd);
+		rop->preop.cmat[r*nc+i] = cxyz[p];
+	}
+}
+
+/* Use the spectral sensitivity function to compute matrix coefficients */
+static void
+sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
+{
+	const int	nc = rop->mtx->ncomp;
+	int		i;
+
+	for (i = nc; i--; ) {
+		SCOLOR	sclr;
+		memset(sclr, 0, sizeof(COLORV)*nc);
+		sclr[i] = 1.f;
+		rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
+	}
+}
+
+/* Check/set symbolic transform */
+static int
+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",
+				rop->inspec, rop->preop.csym);
+		return(-1);
+	}
+	rop->preop.clen = strlen(rop->preop.csym) * nc;
+	if (rop->preop.clen > MAXCOMP*MAXCOMP) {
+		fprintf(stderr, "%s: -c '%s' results in too many components\n",
+				rop->inspec, rop->preop.csym);
+		return(-1);
+	}
+	for (j = 0; rop->preop.csym[j]; j++) {
+		int	comp = 0;
+		switch (rop->preop.csym[j]) {
+		case 'B':
+			++comp;
+			/* fall through */
+		case 'G':
+			++comp;
+			/* fall through */
+		case 'R':
+			if (dt == DTxyze) {
+				for (i = 3; i--; )
+					rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
+							xyz2rgbmat[comp][i];
+			} else if (nc == 3)
+				rop->preop.cmat[j*nc+comp] = 1.;
+			else
+				rgbrow(rop, j, comp);
+			break;
+		case 'Z':
+			++comp;
+			/* fall through */
+		case 'Y':
+			++comp;
+			/* fall through */
+		case 'X':
+			if (dt == DTxyze) {
+				rop->preop.cmat[j*nc+comp] = 1.;
+			} else if (nc == 3) {
+				for (i = 3; i--; )
+					rop->preop.cmat[j*nc+i] =
+							rgb2xyzmat[comp][i];
+			} else if (comp == CIEY)
+				sensrow(rop, j, scolor2photopic);
+			else
+				xyzrow(rop, j, comp);
+
+			for (i = nc*(dt != DTxyze); i--; )
+				rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
+			break;
+		case 'S':		/* scotopic (il)luminance */
+			sensrow(rop, j, scolor2scotopic);
+			for (i = nc; i--; )
+				rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
+			break;
+		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]);
+			return(-1);
+		}
+	}
+					/* return recommended output type */
+	if (!strcmp(rop->preop.csym, "XYZ")) {
+		if (dt <= DTspec)
+			return(DTxyze);
+	} else if (!strcmp(rop->preop.csym, "RGB")) {
+		if (dt <= DTspec)
+			return(DTrgbe);
+	}
+	if ((nc > 3) & (dt <= DTspec))
+		return(DTfloat);	/* probably not actual spectrum */
+	return(0);
+}
+
 /* Get matrix and perform unary operations */
 static RMATRIX *
 loadop(ROPMAT *rop)
 {
+	int	outtype = 0;
 	RMATRIX	*mres;
-	int	i;
+	int	i, j;
 
 	if (loadmatrix(rop) < 0)		/* make sure we're loaded */
 		return(NULL);
 
-	if (rop->preop.nsf > 0) {		/* apply scalar(s) */
-		if (rop->preop.clen > 0) {
-			fputs("Options -s and -c are exclusive\n", stderr);
+	if (rop->preop.csym &&			/* symbolic transform? */
+			(outtype = checksymbolic(rop)) < 0)
+		goto failure;
+	if (rop->preop.clen > 0) {		/* apply component transform? */
+		if (rop->preop.clen % rop->mtx->ncomp) {
+			fprintf(stderr, "%s: -c must have N x %d coefficients\n",
+					rop->inspec, rop->mtx->ncomp);
 			goto failure;
 		}
+		if (rop->preop.nsf > 0) {	/* scale transform, first */
+			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) {
+				fprintf(stderr, "%s: -s must have one or %d factors\n",
+						rop->inspec, 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];
+			}
+		}
+		mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
+					rop->preop.cmat);
+		if (mres == NULL) {
+			fprintf(stderr, "%s: matrix transform failed\n",
+						rop->inspec);
+			goto failure;
+		}
+		if (verbose)
+			fprintf(stderr, "%s: applied %d x %d transform%s\n",
+					rop->inspec, mres->ncomp,
+					rop->mtx->ncomp,
+					rop->preop.nsf ? " (* scalar)" : "");
+		rop->preop.nsf = 0;		/* now folded in */
+		if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
+			outtype = DTfloat;	/* probably not actual spectrum */
+		rmx_free(rop->mtx);
+		rop->mtx = mres;
+	}
+	if (rop->preop.nsf > 0) {		/* apply scalar(s)? */
 		if (rop->preop.nsf == 1) {
 			for (i = rop->mtx->ncomp; --i; )
 				rop->preop.sca[i] = rop->preop.sca[0];
@@ -89,27 +325,7 @@ loadop(ROPMAT *rop)
 			fputs(" )\n", stderr);
 		}
 	}
-	if (rop->preop.clen > 0) {	/* apply transform */
-		if (rop->preop.clen % rop->mtx->ncomp) {
-			fprintf(stderr, "%s: -c must have N x %d coefficients\n",
-					rop->inspec, rop->mtx->ncomp);
-			goto failure;
-		}
-		mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
-					rop->preop.cmat);
-		if (mres == NULL) {
-			fprintf(stderr, "%s: matrix transform failed\n",
-						rop->inspec);
-			goto failure;
-		}
-		if (verbose)
-			fprintf(stderr, "%s: applied %d x %d transform\n",
-					rop->inspec, mres->ncomp,
-					rop->mtx->ncomp);
-		rmx_free(rop->mtx);
-		rop->mtx = mres;
-	}
-	if (rop->preop.transpose) {	/* transpose matrix? */
+	if (rop->preop.transpose) {		/* transpose matrix? */
 		mres = rmx_transpose(rop->mtx);
 		if (mres == NULL) {
 			fputs(rop->inspec, stderr);
@@ -125,6 +341,8 @@ loadop(ROPMAT *rop)
 	}
 	mres = rop->mtx;
 	rop->mtx = NULL;
+	if (outtype)
+		mres->dtype = outtype;
 	return(mres);
 failure:
 	rmx_free(rop->mtx);
@@ -140,19 +358,22 @@ binaryop(const char *inspec, RMATRIX *mleft, int op, R
 
 	if ((mleft == NULL) | (mright == NULL))
 		return(NULL);
-
 	switch (op) {
 	case '.':			/* concatenate */
-		mres = rmx_multiply(mleft, mright);
+		if (mleft->ncomp != mright->ncomp) {
+			fputs(inspec, stderr);
+			fputs(": # components do not match\n", stderr);
+		} else if (mleft->ncols != mright->nrows) {
+			fputs(inspec, stderr);
+			fputs(": mismatched dimensions\n",
+					stderr);
+		} else
+			mres = rmx_multiply(mleft, mright);
 		rmx_free(mleft);
 		rmx_free(mright);
 		if (mres == NULL) {
 			fputs(inspec, stderr);
-			if (mleft->ncols != mright->nrows)
-				fputs(": mismatched dimensions for multiply\n",
-						stderr);
-			else
-				fputs(": concatenation failed\n", stderr);
+			fputs(": concatenation failed\n", stderr);
 			return(NULL);
 		}
 		if (verbose) {
@@ -229,12 +450,17 @@ op_right2left(ROPMAT *mop)
 	int	rpos = 0;
 					/* find end of list */
 	while (mop[rpos].binop)
-		rpos++;
+		if (mop[rpos++].binop != '.') {
+			fputs(
+		"Right-to-left evaluation only for matrix multiplication!\n",
+					stderr);
+			return(NULL);
+		}
 	mright = loadop(mop+rpos);
 	while (rpos-- > 0) {
 		if (mright == NULL)
 			break;
-		mright = binaryop(mop[rpos].inspec,
+		mright = binaryop(mop[rpos+1].inspec,
 				loadop(mop+rpos), mop[rpos].binop, mright);
 	}
 	return(mright);
@@ -284,15 +510,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)
@@ -304,25 +533,25 @@ 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 (nmats >= nall)
-			mop = grow_moparray(mop, nall += 2);
 		if (argv[i][0] && !argv[i][1] &&
 				strchr(".+*/", argv[i][0]) != NULL) {
-			if (mop[nmats].inspec == NULL || mop[nmats].binop) {
-				fprintf(stderr, "%s: missing matrix argument\n",
-						argv[0]);
+			if (!nmats || mop[nmats-1].binop) {
+				fprintf(stderr,
+			"%s: missing matrix argument before '%c' operation\n",
+						argv[0], argv[i][0]);
 				return(1);
 			}
-			mop[nmats++].binop = argv[i][0];
+			mop[nmats-1].binop = argv[i][0];
 		} else if (argv[i][0] != '-' || !argv[i][1]) {
 			if (argv[i][0] == '-') {
 				if (stdin_used++) {
@@ -334,6 +563,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++;
@@ -341,7 +572,7 @@ main(int argc, char *argv[])
 			int	n = argc-1 - i;
 			switch (argv[i][1]) {	/* get option */
 			case 'v':
-				verbose = !verbose;
+				verbose++;
 				break;
 			case 'f':
 				switch (argv[i][2]) {
@@ -369,42 +600,84 @@ main(int argc, char *argv[])
 				i += mop[nmats].preop.nsf =
 					get_factors(mop[nmats].preop.sca,
 							n, argv+i+1);
+				if (mop[nmats].preop.nsf <= 0) {
+					fprintf(stderr, "%s: -s missing arguments\n",
+							argv[0]);
+					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 && !isflt(argv[i+1])) {
+					mop[nmats].preop.csym = argv[++i];
+					mop[nmats].preop.clen = 0;
+					break;
+				}
 				if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
 				i += mop[nmats].preop.clen =
 					get_factors(mop[nmats].preop.cmat,
 							n, argv+i+1);
+				if (mop[nmats].preop.clen <= 0) {
+					fprintf(stderr, "%s: -c missing arguments\n",
+							argv[0]);
+					goto userr;
+				}
+				mop[nmats].preop.csym = NULL;
 				break;
+			case 'r':
+				if (argv[i][2] == 'f')
+					mop[nmats].rmp = RMPreflF;
+				else if (argv[i][2] == 'b')
+					mop[nmats].rmp = RMPreflB;
+				else
+					goto userr;
+				break;
 			default:
 				fprintf(stderr, "%s: unknown operation '%s'\n",
 						argv[0], argv[i]);
 				goto userr;
 			}
 		}
+		if (nmats >= nall)
+			mop = resize_moparr(mop, nall += 2);
 	}
 	if (mop[0].inspec == NULL)	/* nothing to do? */
 		goto userr;
+	if (mop[nmats-1].binop) {
+		fprintf(stderr,
+			"%s: missing matrix argument after '%c' operation\n",
+				argv[0], mop[nmats-1].binop);
+		return(1);
+	}
 					/* favor quicker concatenation */
-	mres = prefer_right2left(mop) ? op_right2left(mop) : op_left2right(mop);
-	if (!mres)
+	mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
+						: op_left2right(mop);
+	if (mop[nmats].mtx == NULL)
 		return(1);
-					/* write result to stdout */
-	if (outfmt == DTfromHeader)
-		outfmt = mres->dtype;
-	if (outfmt != DTascii)
-		SET_FILE_BINARY(stdout);
-	newheader("RADIANCE", stdout);
-	printargs(argc, argv, stdout);
-	if (!rmx_write(mres, outfmt, stdout)) {
-		fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
+					/* apply trailing unary operations */
+	mop[nmats].inspec = "trailing_ops";
+	mres = loadop(mop+nmats);
+	if (mres == NULL)
 		return(1);
+	if (outfmt == DTfromHeader)	/* check data type */
+		outfmt = mres->dtype;
+	if (outfmt == DTrgbe) {
+		if (mres->ncomp > 3)
+			outfmt = DTspec;
+		else if (mres->dtype == DTxyze)
+			outfmt = DTxyze;
 	}
-	/* rmx_free(mres); free(mop); */
-	return(0);
+	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 ..] m1 [.+*/] .. > mres\n",
+	"Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
 			argv[0]);
 	return(1);
 }