--- ray/src/util/rmatrix.c 2022/03/05 01:45:21 2.53 +++ ray/src/util/rmatrix.c 2023/11/28 16:19:31 2.64 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: rmatrix.c,v 2.53 2022/03/05 01:45:21 greg Exp $"; +static const char RCSid[] = "$Id: rmatrix.c,v 2.64 2023/11/28 16:19:31 greg Exp $"; #endif /* * General matrix operations. @@ -16,7 +16,7 @@ static const char RCSid[] = "$Id: rmatrix.c,v 2.53 202 #include #endif -static char rmx_mismatch_warn[] = "WARNING: data type mismatch\n"; +static const char rmx_mismatch_warn[] = "WARNING: data type mismatch\n"; #define array_size(rm) (sizeof(double)*(rm)->nrows*(rm)->ncols*(rm)->ncomp) #define mapped_size(rm) ((char *)(rm)->mtx + array_size(rm) - (char *)(rm)->mapped) @@ -25,13 +25,19 @@ static char rmx_mismatch_warn[] = "WARNING: data type RMATRIX * rmx_new(int nr, int nc, int n) { - RMATRIX *dnew = (RMATRIX *)calloc(1, sizeof(RMATRIX)); + RMATRIX *dnew; + if (n <= 0) + return(NULL); + + dnew = (RMATRIX *)calloc(1, sizeof(RMATRIX)); if (dnew) { dnew->dtype = DTdouble; dnew->nrows = nr; dnew->ncols = nc; dnew->ncomp = n; + setcolor(dnew->cexp, 1.f, 1.f, 1.f); + memcpy(dnew->wlpart, WLPART, sizeof(dnew->wlpart)); } return(dnew); } @@ -43,6 +49,8 @@ rmx_prepare(RMATRIX *rm) if (!rm) return(0); if (rm->mtx) return(1); + if ((rm->nrows <= 0) | (rm->ncols <= 0) | (rm->ncomp <= 0)) + return(0); rm->mtx = (double *)malloc(array_size(rm)); return(rm->mtx != NULL); } @@ -80,8 +88,8 @@ rmx_free(RMATRIX *rm) int rmx_newtype(int dtyp1, int dtyp2) { - if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) | - (dtyp2==DTxyze) | (dtyp2==DTrgbe) + if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) | (dtyp1==DTspec) | + (dtyp2==DTxyze) | (dtyp2==DTrgbe) | (dtyp2==DTspec) && dtyp1 != dtyp2) return(0); if (dtyp1 < dtyp2) @@ -118,10 +126,10 @@ get_dminfo(char *s, void *p) char fmt[MAXFMTLEN]; int i; - if (headidval(fmt, s)) + if (headidval(NULL, s)) return(0); - if (!strncmp(s, "NCOMP=", 6)) { - ip->ncomp = atoi(s+6); + if (isncomp(s)) { + ip->ncomp = ncompval(s); return(0); } if (!strncmp(s, "NROWS=", 6)) { @@ -147,6 +155,10 @@ get_dminfo(char *s, void *p) multcolor(ip->cexp, ctmp); return(0); } + if (iswlsplit(s)) { + wlsplitval(ip->wlpart, s); + return(0); + } if (!formatval(fmt, s)) { rmx_addinfo(ip, s); return(0); @@ -206,7 +218,7 @@ rmx_load_double(RMATRIX *rm, FILE *fp) { int i; #ifdef MAP_FILE - long pos; /* map memory to file if possible */ + long pos; /* map memory for file > 1MB if possible */ if (!rm->swapin && array_size(rm) >= 1L<<20 && (pos = ftell(fp)) >= 0 && !(pos % sizeof(double))) { rm->mapped = mmap(NULL, array_size(rm)+pos, PROT_READ|PROT_WRITE, @@ -214,7 +226,7 @@ rmx_load_double(RMATRIX *rm, FILE *fp) if (rm->mapped != MAP_FAILED) { rm->mtx = (double *)rm->mapped + pos/sizeof(double); return(1); - } + } /* else fall back on reading into memory */ rm->mapped = NULL; } #endif @@ -238,10 +250,12 @@ rmx_load_rgbe(RMATRIX *rm, FILE *fp) if (!scan) return(0); - if (!rmx_prepare(rm)) + if (!rmx_prepare(rm)) { + free(scan); return(0); + } for (i = 0; i < rm->nrows; i++) { - double *dp = rmx_lval(rm,i,j); + double *dp = rmx_lval(rm,i,0); if (freadscan(scan, rm->ncols, fp) < 0) { free(scan); return(0); @@ -256,23 +270,122 @@ rmx_load_rgbe(RMATRIX *rm, FILE *fp) return(1); } +static int +rmx_load_spec(RMATRIX *rm, FILE *fp) +{ + uby8 *scan; + SCOLOR scol; + int i, j, k; + + if (rm->ncomp < 3) + return(0); + scan = (uby8 *)malloc((rm->ncomp+1)*rm->ncols); + if (!scan) + return(0); + if (!rmx_prepare(rm)) { + free(scan); + return(0); + } + for (i = 0; i < rm->nrows; i++) { + double *dp = rmx_lval(rm,i,0); + if (freadscolrs(scan, rm->ncomp, rm->ncols, fp) < 0) { + free(scan); + return(0); + } + for (j = 0; j < rm->ncols; j++) { + scolr2scolor(scol, scan+j*(rm->ncomp+1), rm->ncomp); + for (k = 0; k < rm->ncomp; k++) + *dp++ = scol[k]; + } + } + free(scan); + return(1); +} + +/* Read matrix header from input stream (cannot be XML) */ +int +rmx_load_header(RMATRIX *rm, FILE *fp) +{ + if (!rm | !fp) + return(0); + if (rm->info) { /* clear state */ + free(rm->info); + rm->info = NULL; + } + if (rm->mtx) { /* ...and data */ +#ifdef MAP_FILE + if (rm->mapped) { + munmap(rm->mapped, mapped_size(rm)); + rm->mapped = NULL; + } else +#endif + free(rm->mtx); + rm->mtx = NULL; + } + if (rm->nrows | rm->ncols | !rm->dtype) { + rm->nrows = rm->ncols = 0; + rm->ncomp = 3; + setcolor(rm->cexp, 1.f, 1.f, 1.f); + memcpy(rm->wlpart, WLPART, sizeof(rm->wlpart)); + rm->swapin = 0; + } + SET_FILE_BINARY(fp); + rm->dtype = DTascii; /* assumed w/o FORMAT */ + if (getheader(fp, get_dminfo, rm) < 0) { + fputs("Unrecognized matrix format\n", stderr); + return(0); + } + /* resolution string? */ + if ((rm->nrows <= 0) | (rm->ncols <= 0)) { + if (!fscnresolu(&rm->ncols, &rm->nrows, fp)) + return(0); + if ((rm->dtype == DTrgbe) | (rm->dtype == DTxyze) && + rm->ncomp != 3) + return(0); + } + return(1); +} + +/* Allocate & load post-header data from stream given type set in rm->dtype */ +int +rmx_load_data(RMATRIX *rm, FILE *fp) +{ + switch (rm->dtype) { + case DTascii: + SET_FILE_TEXT(fp); + return(rmx_load_ascii(rm, fp)); + case DTfloat: + return(rmx_load_float(rm, fp)); + case DTdouble: + return(rmx_load_double(rm, fp)); + case DTrgbe: + case DTxyze: + return(rmx_load_rgbe(rm, fp)); + case DTspec: + return(rmx_load_spec(rm, fp)); + default: + fputs("Unsupported data type in rmx_loaddata()\n", stderr); + } + return(0); +} + /* Load matrix from supported file type */ RMATRIX * rmx_load(const char *inspec, RMPref rmp) { FILE *fp; RMATRIX *dnew; + int ok; if (!inspec) inspec = stdin_name; else if (!*inspec) return(NULL); - if (inspec == stdin_name) { /* reading from stdin? */ + if (inspec == stdin_name) /* reading from stdin? */ fp = stdin; - } else if (inspec[0] == '!') { - if (!(fp = popen(inspec+1, "r"))) - return(NULL); - } else { + else if (inspec[0] == '!') + fp = popen(inspec+1, "r"); + else { const char *sp = inspec; /* check suffix */ while (*sp) ++sp; @@ -286,73 +399,26 @@ rmx_load(const char *inspec, RMPref rmp) dnew = rmx_from_cmatrix(cm); cm_free(cm); dnew->dtype = DTascii; - return(dnew); - } - /* else open it ourselves */ - if (!(fp = fopen(inspec, "r"))) - return(NULL); + return(dnew); /* return here */ + } /* else open it ourselves */ + fp = fopen(inspec, "r"); } - SET_FILE_BINARY(fp); + if (!fp) + return(NULL); #ifdef getc_unlocked flockfile(fp); #endif - if (!(dnew = rmx_new(0,0,3))) { - fclose(fp); + /* load header info */ + if (!rmx_load_header(dnew = rmx_new(0,0,3), fp)) { + fprintf(stderr, "Bad header in: %s\n", inspec); + if (inspec[0] == '!') pclose(fp); + else fclose(fp); + rmx_free(dnew); return(NULL); } - dnew->dtype = DTascii; /* assumed w/o FORMAT */ - dnew->cexp[0] = dnew->cexp[1] = dnew->cexp[2] = 1.f; - if (getheader(fp, get_dminfo, dnew) < 0) { - fclose(fp); - return(NULL); - } - if ((dnew->nrows <= 0) | (dnew->ncols <= 0)) { - if (!fscnresolu(&dnew->ncols, &dnew->nrows, fp)) { - fclose(fp); - return(NULL); - } - if ((dnew->dtype == DTrgbe) | (dnew->dtype == DTxyze) && - dnew->ncomp != 3) { - fclose(fp); - return(NULL); - } - } - switch (dnew->dtype) { - case DTascii: - SET_FILE_TEXT(fp); - if (!rmx_load_ascii(dnew, fp)) - goto loaderr; - dnew->dtype = DTascii; /* should leave double? */ - break; - case DTfloat: - if (!rmx_load_float(dnew, fp)) - goto loaderr; - dnew->dtype = DTfloat; - break; - case DTdouble: - if (!rmx_load_double(dnew, fp)) - goto loaderr; - dnew->dtype = DTdouble; - break; - case DTrgbe: - case DTxyze: - if (!rmx_load_rgbe(dnew, fp)) - goto loaderr; - /* undo exposure? */ - if ((dnew->cexp[0] != 1.f) | (dnew->cexp[1] != 1.f) | - (dnew->cexp[2] != 1.f)) { - double cmlt[3]; - cmlt[0] = 1./dnew->cexp[0]; - cmlt[1] = 1./dnew->cexp[1]; - cmlt[2] = 1./dnew->cexp[2]; - rmx_scale(dnew, cmlt); - } - dnew->swapin = 0; - break; - default: - goto loaderr; - } - if (fp != stdin) { + ok = rmx_load_data(dnew, fp); /* allocate & load data */ + + if (fp != stdin) { /* close input stream */ if (inspec[0] == '!') pclose(fp); else @@ -362,27 +428,44 @@ rmx_load(const char *inspec, RMPref rmp) else funlockfile(fp); #endif + if (!ok) { /* load failure? */ + fprintf(stderr, "Error loading data from: %s\n", inspec); + rmx_free(dnew); + return(NULL); + } + /* undo exposure? */ + if ((dnew->cexp[0] != 1.f) | + (dnew->cexp[1] != 1.f) | (dnew->cexp[2] != 1.f)) { + double cmlt[MAXCSAMP]; + int i; + cmlt[0] = 1./dnew->cexp[0]; + cmlt[1] = 1./dnew->cexp[1]; + cmlt[2] = 1./dnew->cexp[2]; + if (dnew->ncomp > MAXCSAMP) { + fprintf(stderr, "Excess spectral components in: %s\n", + inspec); + rmx_free(dnew); + return(NULL); + } + for (i = dnew->ncomp; i-- > 3; ) + cmlt[i] = cmlt[1]; + rmx_scale(dnew, cmlt); + setcolor(dnew->cexp, 1.f, 1.f, 1.f); + } return(dnew); -loaderr: /* should report error? */ - if (inspec[0] == '!') - pclose(fp); - else - fclose(fp); - rmx_free(dnew); - return(NULL); } static int rmx_write_ascii(const RMATRIX *rm, FILE *fp) { const char *fmt = (rm->dtype == DTfloat) ? " %.7e" : - (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) ? " %.3e" : - " %.15e" ; + (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) | + (rm->dtype == DTspec) ? " %.3e" : " %.15e" ; int i, j, k; for (i = 0; i < rm->nrows; i++) { for (j = 0; j < rm->ncols; j++) { - const double *dp = rmx_lval(rm,i,j); + const double *dp = rmx_val(rm,i,j); for (k = 0; k < rm->ncomp; k++) fprintf(fp, fmt, dp[k]); fputc('\t', fp); @@ -404,7 +487,7 @@ rmx_write_float(const RMATRIX *rm, FILE *fp) } for (i = 0; i < rm->nrows; i++) for (j = 0; j < rm->ncols; j++) { - const double *dp = rmx_lval(rm,i,j); + const double *dp = rmx_val(rm,i,j); for (k = rm->ncomp; k--; ) val[k] = (float)dp[k]; if (putbinary(val, sizeof(float), rm->ncomp, fp) != rm->ncomp) @@ -419,7 +502,7 @@ rmx_write_double(const RMATRIX *rm, FILE *fp) int i; for (i = 0; i < rm->nrows; i++) - if (putbinary(rmx_lval(rm,i,0), sizeof(double)*rm->ncomp, + if (putbinary(rmx_val(rm,i,0), sizeof(double)*rm->ncomp, rm->ncols, fp) != rm->ncols) return(0); return(1); @@ -435,8 +518,11 @@ rmx_write_rgbe(const RMATRIX *rm, FILE *fp) return(0); for (i = 0; i < rm->nrows; i++) { for (j = rm->ncols; j--; ) { - const double *dp = rmx_lval(rm,i,j); - setcolr(scan[j], dp[0], dp[1], dp[2]); + const double *dp = rmx_val(rm,i,j); + if (rm->ncomp == 1) + setcolr(scan[j], dp[0], dp[0], dp[0]); + else + setcolr(scan[j], dp[0], dp[1], dp[2]); } if (fwritecolrs(scan, rm->ncols, fp) < 0) { free(scan); @@ -447,6 +533,32 @@ rmx_write_rgbe(const RMATRIX *rm, FILE *fp) return(1); } +static int +rmx_write_spec(const RMATRIX *rm, FILE *fp) +{ + uby8 *scan = (uby8 *)malloc((rm->ncomp+1)*rm->ncols); + int ok = 1; + SCOLOR scol; + int i, j, k; + + if (!scan) + return(0); + for (i = 0; i < rm->nrows; i++) { + for (j = rm->ncols; j--; ) { + const double *dp = rmx_val(rm,i,j); + for (k = rm->ncomp; k--; ) + scol[k] = dp[k]; + scolor2scolr(scan+j*(rm->ncomp+1), scol, rm->ncomp); + } + if (fwritescolrs(scan, rm->ncomp, rm->ncols, fp) < 0) { + ok = 0; + break; + } + } + free(scan); + return(ok); +} + /* Check if CIE XYZ primaries were specified */ static int findCIEprims(const char *info) @@ -475,8 +587,7 @@ rmx_write(const RMATRIX *rm, int dtype, FILE *fp) #ifdef getc_unlocked flockfile(fp); #endif - /* complete header */ - if (rm->info) + if (rm->info) /* complete header */ fputs(rm->info, fp); if (dtype == DTfromHeader) dtype = rm->dtype; @@ -485,26 +596,28 @@ rmx_write(const RMATRIX *rm, int dtype, FILE *fp) dtype = DTxyze; else if ((dtype == DTxyze) & (rm->dtype == DTrgbe)) dtype = DTrgbe; + /* write exposure? */ + if (rm->ncomp == 3 && (rm->cexp[RED] != rm->cexp[GRN]) | + (rm->cexp[GRN] != rm->cexp[BLU])) + fputcolcor(rm->cexp, fp); + else if (rm->cexp[GRN] != 1.f) + fputexpos(rm->cexp[GRN], fp); if ((dtype != DTrgbe) & (dtype != DTxyze)) { - fprintf(fp, "NROWS=%d\n", rm->nrows); - fprintf(fp, "NCOLS=%d\n", rm->ncols); - fprintf(fp, "NCOMP=%d\n", rm->ncomp); - } else if (rm->ncomp != 3) { /* wrong # components? */ - CMATRIX *cm; /* convert & write */ - if (rm->ncomp != 1) /* only convert grayscale */ - return(0); - if (!(cm = cm_from_rmatrix(rm))) - return(0); - fputformat(cm_fmt_id[dtype], fp); - fputc('\n', fp); - ok = cm_write(cm, dtype, fp); - cm_free(cm); - return(ok); - } + if (dtype != DTspec) { + fprintf(fp, "NROWS=%d\n", rm->nrows); + fprintf(fp, "NCOLS=%d\n", rm->ncols); + } else if (rm->ncomp < 3) + return(0); /* bad # components */ + fputncomp(rm->ncomp, fp); + if (dtype == DTspec || (rm->ncomp > 3 && + memcmp(rm->wlpart, WLPART, sizeof(WLPART)))) + fputwlsplit(rm->wlpart, fp); + } else if ((rm->ncomp != 3) & (rm->ncomp != 1)) + return(0); /* wrong # components */ if ((dtype == DTfloat) | (dtype == DTdouble)) fputendian(fp); /* important to record */ fputformat(cm_fmt_id[dtype], fp); - fputc('\n', fp); + fputc('\n', fp); /* end of header */ switch (dtype) { /* write data */ case DTascii: ok = rmx_write_ascii(rm, fp); @@ -520,6 +633,10 @@ rmx_write(const RMATRIX *rm, int dtype, FILE *fp) fprtresolu(rm->ncols, rm->nrows, fp); ok = rmx_write_rgbe(rm, fp); break; + case DTspec: + fprtresolu(rm->ncols, rm->nrows, fp); + ok = rmx_write_spec(rm, fp); + break; default: return(0); } @@ -527,6 +644,7 @@ rmx_write(const RMATRIX *rm, int dtype, FILE *fp) #ifdef getc_unlocked funlockfile(fp); #endif + if (!ok) fputs("Error writing matrix\n", stderr); return(ok); } @@ -561,6 +679,8 @@ rmx_copy(const RMATRIX *rm) return(NULL); rmx_addinfo(dnew, rm->info); dnew->dtype = rm->dtype; + copycolor(dnew->cexp, rm->cexp); + memcpy(dnew->wlpart, rm->wlpart, sizeof(dnew->wlpart)); memcpy(dnew->mtx, rm->mtx, array_size(dnew)); return(dnew); } @@ -570,7 +690,7 @@ RMATRIX * rmx_transpose(const RMATRIX *rm) { RMATRIX *dnew; - int i, j, k; + int i, j; if (!rm) return(0); @@ -590,9 +710,11 @@ rmx_transpose(const RMATRIX *rm) rmx_addinfo(dnew, "Transposed rows and columns\n"); } dnew->dtype = rm->dtype; - for (i = dnew->nrows; i--; ) - for (j = dnew->ncols; j--; ) - memcpy(rmx_lval(dnew,i,j), rmx_lval(rm,i,j), + copycolor(dnew->cexp, rm->cexp); + memcpy(dnew->wlpart, rm->wlpart, sizeof(dnew->wlpart)); + for (j = dnew->ncols; j--; ) + for (i = dnew->nrows; i--; ) + memcpy(rmx_lval(dnew,i,j), rmx_val(rm,j,i), sizeof(double)*dnew->ncomp); return(dnew); } @@ -619,7 +741,7 @@ rmx_multiply(const RMATRIX *m1, const RMATRIX *m2) for (k = mres->ncomp; k--; ) { double d = 0; for (h = m1->ncols; h--; ) - d += rmx_lval(m1,i,h)[k] * rmx_lval(m2,h,j)[k]; + d += rmx_val(m1,i,h)[k] * rmx_val(m2,h,j)[k]; rmx_lval(mres,i,j)[k] = d; } return(mres); @@ -646,7 +768,7 @@ rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide if (divide) { double d; if (m2->ncomp == 1) { - d = rmx_lval(m2,i,j)[0]; + d = rmx_val(m2,i,j)[0]; if (d == 0) { ++zeroDivides; for (k = m1->ncomp; k--; ) @@ -658,7 +780,7 @@ rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide } } else for (k = m1->ncomp; k--; ) { - d = rmx_lval(m2,i,j)[k]; + d = rmx_val(m2,i,j)[k]; if (d == 0) { ++zeroDivides; rmx_lval(m1,i,j)[k] = 0; @@ -667,12 +789,12 @@ rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide } } else { if (m2->ncomp == 1) { - const double d = rmx_lval(m2,i,j)[0]; + const double d = rmx_val(m2,i,j)[0]; for (k = m1->ncomp; k--; ) rmx_lval(m1,i,j)[k] *= d; } else for (k = m1->ncomp; k--; ) - rmx_lval(m1,i,j)[k] *= rmx_lval(m2,i,j)[k]; + rmx_lval(m1,i,j)[k] *= rmx_val(m2,i,j)[k]; } if (zeroDivides) { rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n"); @@ -708,7 +830,7 @@ rmx_sum(RMATRIX *msum, const RMATRIX *madd, const doub rmx_addinfo(msum, rmx_mismatch_warn); for (i = msum->nrows; i--; ) for (j = msum->ncols; j--; ) { - const double *da = rmx_lval(madd,i,j); + const double *da = rmx_val(madd,i,j); double *ds = rmx_lval(msum,i,j); for (k = msum->ncomp; k--; ) ds[k] += sf[k] * da[k]; @@ -734,6 +856,7 @@ rmx_scale(RMATRIX *rm, const double sf[]) } if (rm->info) rmx_addinfo(rm, "Applied scalar\n"); + /* XXX: should record as exposure for COLR and SCOLR types? */ return(1); } @@ -759,7 +882,7 @@ rmx_transform(const RMATRIX *msrc, int n, const double dnew->dtype = msrc->dtype; for (i = dnew->nrows; i--; ) for (j = dnew->ncols; j--; ) { - const double *ds = rmx_lval(msrc,i,j); + const double *ds = rmx_val(msrc,i,j); for (kd = dnew->ncomp; kd--; ) { double d = 0; for (ks = msrc->ncomp; ks--; ) @@ -801,19 +924,30 @@ cm_from_rmatrix(const RMATRIX *rm) int i, j; CMATRIX *cnew; - if (!rm || !rm->mtx | ((rm->ncomp != 3) & (rm->ncomp != 1))) + if (!rm || !rm->mtx | (rm->ncomp == 2)) return(NULL); cnew = cm_alloc(rm->nrows, rm->ncols); if (!cnew) return(NULL); for (i = cnew->nrows; i--; ) for (j = cnew->ncols; j--; ) { - const double *dp = rmx_lval(rm,i,j); + const double *dp = rmx_val(rm,i,j); COLORV *cv = cm_lval(cnew,i,j); - if (rm->ncomp == 1) - setcolor(cv, dp[0], dp[0], dp[0]); - else + switch (rm->ncomp) { + case 3: setcolor(cv, dp[0], dp[1], dp[2]); + break; + case 1: + setcolor(cv, dp[0], dp[0], dp[0]); + break; + default: { + SCOLOR scol; + int k; + for (k = rm->ncomp; k--; ) + scol[k] = dp[k]; + scolor2color(cv, scol, rm->ncomp, rm->wlpart); + } break; + } } return(cnew); }