--- ray/src/util/rmatrix.c	2014/05/31 05:02:37	2.1
+++ ray/src/util/rmatrix.c	2014/08/01 23:37:24	2.5
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rmatrix.c,v 2.1 2014/05/31 05:02:37 greg Exp $";
+static const char RCSid[] = "$Id: rmatrix.c,v 2.5 2014/08/01 23:37:24 greg Exp $";
 #endif
 /*
  * General matrix operations.
@@ -12,9 +12,13 @@ static const char RCSid[] = "$Id: rmatrix.c,v 2.1 2014
 #include "resolu.h"
 #include "rmatrix.h"
 
+#define	MAX_INFO	16000
+
 typedef struct {
 	int	nrows, ncols, ncomp;
 	int	dtype;
+	int	info_len;
+	char	info[MAX_INFO];
 } DMINFO;
 
 /* Allocate a nr x nc matrix with n components */
@@ -30,14 +34,32 @@ rmx_alloc(int nr, int nc, int n)
 	if (dnew == NULL)
 		return(NULL);
 	dnew->nrows = nr; dnew->ncols = nc; dnew->ncomp = n;
+	dnew->info = NULL;
 	return(dnew);
 }
 
+/* Append header information associated with matrix data */
+int
+rmx_addinfo(RMATRIX *rm, const char *info)
+{
+	if (!info || !*info)
+		return(0);
+	if (!rm->info)
+		rm->info = (char *)malloc(strlen(info)+1);
+	else
+		rm->info = (char *)realloc(rm->info,
+				strlen(rm->info)+strlen(info)+1);
+	if (!rm->info)
+		return(0);
+	strcat(rm->info, info);
+	return(1);
+}
+
 static int
 get_dminfo(char *s, void *p)
 {
 	DMINFO	*ip = (DMINFO *)p;
-	char	fmt[32];
+	char	fmt[64];
 	int	i;
 
 	if (!strncmp(s, "NCOMP=", 6)) {
@@ -52,8 +74,20 @@ get_dminfo(char *s, void *p)
 		ip->ncols = atoi(s+6);
 		return(0);
 	}
-	if (!formatval(fmt, s))
+	if (!formatval(fmt, s)) {
+		if (headidval(fmt, s))
+			return(0);
+		while (*s) {
+			if (ip->info_len == MAX_INFO-2 &&
+					ip->info[MAX_INFO-3] != '\n')
+				ip->info[ip->info_len++] = '\n';
+			if (ip->info_len >= MAX_INFO-1)
+				break;
+			ip->info[ip->info_len++] = *s++;
+		}
+		ip->info[ip->info_len] = '\0';
 		return(0);
+	}
 	for (i = 1; i < DTend; i++)
 		if (!strcmp(fmt, cm_fmt_id[i])) {
 			ip->dtype = i;
@@ -173,22 +207,31 @@ rmx_load(const char *fname)
 #endif
 	dinfo.nrows = dinfo.ncols = dinfo.ncomp = 0;
 	dinfo.dtype = DTascii;
-	if (getheader(fp, &get_dminfo, &dinfo) < 0) {
+	dinfo.info_len = 0;
+	if (getheader(fp, get_dminfo, &dinfo) < 0) {
 		fclose(fp);
 		return(NULL);
 	}
-	if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze)) {
+	if ((dinfo.nrows <= 0) | (dinfo.ncols <= 0)) {
 		if (!fscnresolu(&dinfo.ncols, &dinfo.nrows, fp)) {
 			fclose(fp);
 			return(NULL);
 		}
-		dinfo.ncomp = 3;
+		if (dinfo.ncomp <= 0)
+			dinfo.ncomp = 3;
+		else if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze) &&
+				dinfo.ncomp != 3) {
+			fclose(fp);
+			return(NULL);
+		}
 	}
 	dnew = rmx_alloc(dinfo.nrows, dinfo.ncols, dinfo.ncomp);
 	if (dnew == NULL) {
 		fclose(fp);
 		return(NULL);
 	}
+	if (dinfo.info_len)
+		rmx_addinfo(dnew, dinfo.info);
 	switch (dinfo.dtype) {
 	case DTascii:
 		if (!rmx_load_ascii(dnew, fp))
@@ -309,6 +352,8 @@ rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
 	if ((rm == NULL) | (fp == NULL))
 		return(0);
 						/* complete header */
+	if (rm->info)
+		fputs(rm->info, fp);
 	if ((dtype != DTrgbe) & (dtype != DTxyze)) {
 		fprintf(fp, "NROWS=%d\n", rm->nrows);
 		fprintf(fp, "NCOLS=%d\n", rm->ncols);
@@ -377,34 +422,33 @@ rmx_copy(const RMATRIX *rm)
 	dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
 	if (dnew == NULL)
 		return(NULL);
+	rmx_addinfo(dnew, rm->info);
 	memcpy(dnew->mtx, rm->mtx,
 		sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
 	return(dnew);
 }
 
-/* Swap rows and columns in the given matrix in situ */
-int
-rmx_transpose(RMATRIX *rm)
+/* Allocate and assign transposed matrix */
+RMATRIX *
+rmx_transpose(const RMATRIX *rm)
 {
-	RMATRIX	dswap;
+	RMATRIX	*dnew;
 	int	i, j, k;
 
 	if (rm == NULL)
 		return(0);
-	dswap.nrows = rm->ncols;
-	dswap.ncols = rm->nrows;
-	dswap.ncomp = rm->ncomp;
-	for (i = 1; i < rm->nrows; i++)
-	    for (j = 0; j < i; j++)
-		for (k = rm->ncomp; k--; ) {
-			double	*opp = rm->mtx + rmx_indx(&dswap,j,i,k);
-			double	d = *opp;
-			*opp = rmx_lval(rm,i,j,k);
-			rmx_lval(rm,i,j,k) = d;
-		}
-	rm->nrows = dswap.nrows;
-	rm->ncols = dswap.ncols;
-	return(1);
+	dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
+	if (dnew == NULL)
+		return(NULL);
+	if (rm->info) {
+		rmx_addinfo(dnew, rm->info);
+		rmx_addinfo(dnew, "Transposed rows and columns\n");
+	}
+	for (i = dnew->nrows; i--; )
+	    for (j = dnew->ncols; j--; )
+		for (k = dnew->ncomp; k--; )
+			rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
+	return(dnew);
 }
 
 /* Multiply (concatenate) two matrices and allocate the result */