| 1 | greg | 2.1 | #ifndef lint | 
| 2 | greg | 2.22 | static const char RCSid[] = "$Id: dctimestep.c,v 2.21 2012/10/16 15:53:04 greg Exp $"; | 
| 3 | greg | 2.1 | #endif | 
| 4 |  |  | /* | 
| 5 |  |  | * Compute time-step result using Daylight Coefficient method. | 
| 6 |  |  | * | 
| 7 |  |  | *      G. Ward | 
| 8 |  |  | */ | 
| 9 |  |  |  | 
| 10 |  |  | #include <ctype.h> | 
| 11 |  |  | #include "standard.h" | 
| 12 |  |  | #include "platform.h" | 
| 13 | greg | 2.15 | #include "paths.h" | 
| 14 | greg | 2.1 | #include "color.h" | 
| 15 |  |  | #include "resolu.h" | 
| 16 |  |  | #include "bsdf.h" | 
| 17 | greg | 2.15 | #include "bsdf_m.h" | 
| 18 | greg | 2.1 |  | 
| 19 |  |  | char    *progname;                      /* global argv[0] */ | 
| 20 |  |  |  | 
| 21 | greg | 2.2 | /* Data types for file loading */ | 
| 22 | greg | 2.1 | enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze}; | 
| 23 |  |  |  | 
| 24 |  |  | /* A color coefficient matrix -- vectors have ncols==1 */ | 
| 25 |  |  | typedef struct { | 
| 26 |  |  | int     nrows, ncols; | 
| 27 |  |  | COLORV  cmem[3];                /* extends struct */ | 
| 28 |  |  | } CMATRIX; | 
| 29 |  |  |  | 
| 30 |  |  | #define COLSPEC (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf") | 
| 31 |  |  |  | 
| 32 |  |  | #define cm_lval(cm,r,c) ((cm)->cmem + 3*((r)*(cm)->ncols + (c))) | 
| 33 |  |  |  | 
| 34 |  |  | #define cv_lval(cm,i)   ((cm)->cmem + 3*(i)) | 
| 35 |  |  |  | 
| 36 |  |  | /* Allocate a color coefficient matrix */ | 
| 37 |  |  | static CMATRIX * | 
| 38 |  |  | cm_alloc(int nrows, int ncols) | 
| 39 |  |  | { | 
| 40 |  |  | CMATRIX *cm; | 
| 41 |  |  |  | 
| 42 |  |  | if ((nrows <= 0) | (ncols <= 0)) | 
| 43 | greg | 2.10 | error(USER, "attempt to create empty matrix"); | 
| 44 | greg | 2.1 | cm = (CMATRIX *)malloc(sizeof(CMATRIX) + | 
| 45 |  |  | 3*sizeof(COLORV)*(nrows*ncols - 1)); | 
| 46 |  |  | if (cm == NULL) | 
| 47 |  |  | error(SYSTEM, "out of memory in cm_alloc()"); | 
| 48 |  |  | cm->nrows = nrows; | 
| 49 |  |  | cm->ncols = ncols; | 
| 50 |  |  | return(cm); | 
| 51 |  |  | } | 
| 52 |  |  |  | 
| 53 |  |  | #define cm_free(cm)     free(cm) | 
| 54 |  |  |  | 
| 55 |  |  | /* Resize color coefficient matrix */ | 
| 56 |  |  | static CMATRIX * | 
| 57 |  |  | cm_resize(CMATRIX *cm, int nrows) | 
| 58 |  |  | { | 
| 59 |  |  | if (nrows == cm->nrows) | 
| 60 |  |  | return(cm); | 
| 61 |  |  | if (nrows <= 0) { | 
| 62 |  |  | cm_free(cm); | 
| 63 |  |  | return(NULL); | 
| 64 |  |  | } | 
| 65 |  |  | cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) + | 
| 66 |  |  | 3*sizeof(COLORV)*(nrows*cm->ncols - 1)); | 
| 67 |  |  | if (cm == NULL) | 
| 68 |  |  | error(SYSTEM, "out of memory in cm_resize()"); | 
| 69 |  |  | cm->nrows = nrows; | 
| 70 |  |  | return(cm); | 
| 71 |  |  | } | 
| 72 |  |  |  | 
| 73 |  |  | /* Load header to obtain data type */ | 
| 74 |  |  | static int | 
| 75 |  |  | getDT(char *s, void *p) | 
| 76 |  |  | { | 
| 77 |  |  | char    fmt[32]; | 
| 78 |  |  |  | 
| 79 |  |  | if (formatval(fmt, s)) { | 
| 80 |  |  | if (!strcmp(fmt, "ascii")) | 
| 81 |  |  | *((int *)p) = DTascii; | 
| 82 |  |  | else if (!strcmp(fmt, "float")) | 
| 83 |  |  | *((int *)p) = DTfloat; | 
| 84 |  |  | else if (!strcmp(fmt, "double")) | 
| 85 |  |  | *((int *)p) = DTdouble; | 
| 86 |  |  | else if (!strcmp(fmt, COLRFMT)) | 
| 87 |  |  | *((int *)p) = DTrgbe; | 
| 88 |  |  | else if (!strcmp(fmt, CIEFMT)) | 
| 89 |  |  | *((int *)p) = DTxyze; | 
| 90 |  |  | } | 
| 91 |  |  | return(0); | 
| 92 |  |  | } | 
| 93 |  |  |  | 
| 94 |  |  | static int | 
| 95 |  |  | getDTfromHeader(FILE *fp) | 
| 96 |  |  | { | 
| 97 |  |  | int     dt = DTfromHeader; | 
| 98 |  |  |  | 
| 99 |  |  | if (getheader(fp, getDT, &dt) < 0) | 
| 100 |  |  | error(SYSTEM, "header read error"); | 
| 101 |  |  | if (dt == DTfromHeader) | 
| 102 |  |  | error(USER, "missing data format in header"); | 
| 103 |  |  | return(dt); | 
| 104 |  |  | } | 
| 105 |  |  |  | 
| 106 |  |  | /* Allocate and load a matrix from the given file (or stdin if NULL) */ | 
| 107 |  |  | static CMATRIX * | 
| 108 |  |  | cm_load(const char *fname, int nrows, int ncols, int dtype) | 
| 109 |  |  | { | 
| 110 |  |  | CMATRIX *cm; | 
| 111 |  |  | FILE    *fp = stdin; | 
| 112 |  |  |  | 
| 113 | greg | 2.12 | if (ncols <= 0) | 
| 114 |  |  | error(USER, "Non-positive number of columns"); | 
| 115 | greg | 2.1 | if (fname == NULL) | 
| 116 |  |  | fname = "<stdin>"; | 
| 117 |  |  | else if ((fp = fopen(fname, "r")) == NULL) { | 
| 118 |  |  | sprintf(errmsg, "cannot open file '%s'", fname); | 
| 119 |  |  | error(SYSTEM, errmsg); | 
| 120 |  |  | } | 
| 121 |  |  | if (dtype != DTascii) | 
| 122 |  |  | SET_FILE_BINARY(fp); | 
| 123 |  |  | if (dtype == DTfromHeader) | 
| 124 |  |  | dtype = getDTfromHeader(fp); | 
| 125 |  |  | switch (dtype) { | 
| 126 |  |  | case DTascii: | 
| 127 |  |  | case DTfloat: | 
| 128 |  |  | case DTdouble: | 
| 129 |  |  | break; | 
| 130 |  |  | default: | 
| 131 |  |  | error(USER, "unexpected data type in cm_load()"); | 
| 132 |  |  | } | 
| 133 |  |  | if (nrows <= 0) {                       /* don't know length? */ | 
| 134 |  |  | int     guessrows = 147;        /* usually big enough */ | 
| 135 |  |  | if ((dtype != DTascii) & (fp != stdin)) { | 
| 136 |  |  | long    startpos = ftell(fp); | 
| 137 |  |  | if (fseek(fp, 0L, SEEK_END) == 0) { | 
| 138 |  |  | long    endpos = ftell(fp); | 
| 139 |  |  | long    elemsiz = 3*(dtype==DTfloat ? | 
| 140 |  |  | sizeof(float) : sizeof(double)); | 
| 141 | greg | 2.2 |  | 
| 142 |  |  | if ((endpos - startpos) % (ncols*elemsiz)) { | 
| 143 |  |  | sprintf(errmsg, | 
| 144 |  |  | "improper length for binary file '%s'", | 
| 145 |  |  | fname); | 
| 146 |  |  | error(USER, errmsg); | 
| 147 |  |  | } | 
| 148 |  |  | guessrows = (endpos - startpos)/(ncols*elemsiz); | 
| 149 | greg | 2.1 | if (fseek(fp, startpos, SEEK_SET) < 0) { | 
| 150 |  |  | sprintf(errmsg, | 
| 151 |  |  | "fseek() error on file '%s'", | 
| 152 |  |  | fname); | 
| 153 |  |  | error(SYSTEM, errmsg); | 
| 154 |  |  | } | 
| 155 | greg | 2.2 | nrows = guessrows;      /* we're confident */ | 
| 156 | greg | 2.1 | } | 
| 157 |  |  | } | 
| 158 |  |  | cm = cm_alloc(guessrows, ncols); | 
| 159 |  |  | } else | 
| 160 |  |  | cm = cm_alloc(nrows, ncols); | 
| 161 |  |  | if (cm == NULL) | 
| 162 |  |  | return(NULL); | 
| 163 |  |  | if (dtype == DTascii) {                         /* read text file */ | 
| 164 |  |  | int     maxrow = (nrows > 0 ? nrows : 32000); | 
| 165 |  |  | int     r, c; | 
| 166 |  |  | for (r = 0; r < maxrow; r++) { | 
| 167 | greg | 2.6 | if (r >= cm->nrows)                 /* need more space? */ | 
| 168 | greg | 2.1 | cm = cm_resize(cm, 2*cm->nrows); | 
| 169 |  |  | for (c = 0; c < ncols; c++) { | 
| 170 |  |  | COLORV  *cv = cm_lval(cm,r,c); | 
| 171 |  |  | if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3) | 
| 172 | greg | 2.6 | if ((nrows <= 0) & (r > 0) & !c) { | 
| 173 | greg | 2.1 | cm = cm_resize(cm, maxrow=r); | 
| 174 |  |  | break; | 
| 175 |  |  | } else | 
| 176 |  |  | goto EOFerror; | 
| 177 |  |  | } | 
| 178 |  |  | } | 
| 179 |  |  | while ((c = getc(fp)) != EOF) | 
| 180 |  |  | if (!isspace(c)) { | 
| 181 |  |  | sprintf(errmsg, | 
| 182 |  |  | "unexpected data at end of ascii file %s", | 
| 183 |  |  | fname); | 
| 184 |  |  | error(WARNING, errmsg); | 
| 185 |  |  | break; | 
| 186 |  |  | } | 
| 187 |  |  | } else {                                        /* read binary file */ | 
| 188 |  |  | if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) : | 
| 189 |  |  | sizeof(double))) { | 
| 190 |  |  | int     nread = 0; | 
| 191 |  |  | do {                            /* read all we can */ | 
| 192 |  |  | nread += fread(cm->cmem + 3*nread, | 
| 193 |  |  | 3*sizeof(COLORV), | 
| 194 |  |  | cm->nrows*cm->ncols - nread, | 
| 195 |  |  | fp); | 
| 196 |  |  | if (nrows <= 0) {       /* unknown length */ | 
| 197 |  |  | if (nread == cm->nrows*cm->ncols) | 
| 198 |  |  | /* need more space? */ | 
| 199 |  |  | cm = cm_resize(cm, 2*cm->nrows); | 
| 200 | greg | 2.6 | else if (nread && !(nread % cm->ncols)) | 
| 201 | greg | 2.1 | /* seem to be  done */ | 
| 202 |  |  | cm = cm_resize(cm, nread/cm->ncols); | 
| 203 |  |  | else            /* ended mid-row */ | 
| 204 |  |  | goto EOFerror; | 
| 205 |  |  | } else if (nread < cm->nrows*cm->ncols) | 
| 206 |  |  | goto EOFerror; | 
| 207 |  |  | } while (nread < cm->nrows*cm->ncols); | 
| 208 |  |  |  | 
| 209 |  |  | } else if (dtype == DTdouble) { | 
| 210 |  |  | double  dc[3];                  /* load from double */ | 
| 211 |  |  | COLORV  *cvp = cm->cmem; | 
| 212 |  |  | int     n = nrows*ncols; | 
| 213 |  |  |  | 
| 214 |  |  | if (n <= 0) | 
| 215 |  |  | goto not_handled; | 
| 216 |  |  | while (n--) { | 
| 217 |  |  | if (fread(dc, sizeof(double), 3, fp) != 3) | 
| 218 |  |  | goto EOFerror; | 
| 219 |  |  | copycolor(cvp, dc); | 
| 220 |  |  | cvp += 3; | 
| 221 |  |  | } | 
| 222 |  |  | } else /* dtype == DTfloat */ { | 
| 223 |  |  | float   fc[3];                  /* load from float */ | 
| 224 |  |  | COLORV  *cvp = cm->cmem; | 
| 225 |  |  | int     n = nrows*ncols; | 
| 226 |  |  |  | 
| 227 |  |  | if (n <= 0) | 
| 228 |  |  | goto not_handled; | 
| 229 |  |  | while (n--) { | 
| 230 |  |  | if (fread(fc, sizeof(float), 3, fp) != 3) | 
| 231 |  |  | goto EOFerror; | 
| 232 |  |  | copycolor(cvp, fc); | 
| 233 |  |  | cvp += 3; | 
| 234 |  |  | } | 
| 235 |  |  | } | 
| 236 |  |  | if (getc(fp) != EOF) { | 
| 237 |  |  | sprintf(errmsg, | 
| 238 |  |  | "unexpected data at end of binary file %s", | 
| 239 |  |  | fname); | 
| 240 |  |  | error(WARNING, errmsg); | 
| 241 |  |  | } | 
| 242 |  |  | } | 
| 243 |  |  | if (fp != stdin) | 
| 244 |  |  | fclose(fp); | 
| 245 |  |  | return(cm); | 
| 246 |  |  | EOFerror: | 
| 247 |  |  | sprintf(errmsg, "unexpected EOF reading %s", | 
| 248 |  |  | fname); | 
| 249 |  |  | error(USER, errmsg); | 
| 250 |  |  | not_handled: | 
| 251 |  |  | error(INTERNAL, "unhandled data size or length in cm_load()"); | 
| 252 |  |  | return(NULL);   /* gratis return */ | 
| 253 |  |  | } | 
| 254 |  |  |  | 
| 255 | greg | 2.20 | /* Scale a matrix by a single value */ | 
| 256 |  |  | static CMATRIX * | 
| 257 |  |  | cm_scale(const CMATRIX *cm1, const COLOR sca) | 
| 258 |  |  | { | 
| 259 |  |  | CMATRIX *cmr; | 
| 260 |  |  | int     dr, dc; | 
| 261 |  |  |  | 
| 262 |  |  | cmr = cm_alloc(cm1->nrows, cm1->ncols); | 
| 263 |  |  | if (cmr == NULL) | 
| 264 |  |  | return(NULL); | 
| 265 |  |  | for (dr = 0; dr < cmr->nrows; dr++) | 
| 266 |  |  | for (dc = 0; dc < cmr->ncols; dc++) { | 
| 267 |  |  | const COLORV    *sp = cm_lval(cm1,dr,dc); | 
| 268 |  |  | COLORV          *dp = cm_lval(cmr,dr,dc); | 
| 269 |  |  | dp[0] = sp[0] * sca[0]; | 
| 270 |  |  | dp[1] = sp[1] * sca[1]; | 
| 271 |  |  | dp[2] = sp[2] * sca[2]; | 
| 272 |  |  | } | 
| 273 |  |  | return(cmr); | 
| 274 |  |  | } | 
| 275 |  |  |  | 
| 276 |  |  | /* Multiply two matrices (or a matrix and a vector) and allocate the result */ | 
| 277 | greg | 2.1 | static CMATRIX * | 
| 278 |  |  | cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2) | 
| 279 |  |  | { | 
| 280 |  |  | CMATRIX *cmr; | 
| 281 |  |  | int     dr, dc, i; | 
| 282 |  |  |  | 
| 283 |  |  | if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows)) | 
| 284 |  |  | error(INTERNAL, "matrix dimension mismatch in cm_multiply()"); | 
| 285 |  |  | cmr = cm_alloc(cm1->nrows, cm2->ncols); | 
| 286 |  |  | if (cmr == NULL) | 
| 287 |  |  | return(NULL); | 
| 288 |  |  | for (dr = 0; dr < cmr->nrows; dr++) | 
| 289 |  |  | for (dc = 0; dc < cmr->ncols; dc++) { | 
| 290 |  |  | COLORV  *dp = cm_lval(cmr,dr,dc); | 
| 291 |  |  | dp[0] = dp[1] = dp[2] = 0; | 
| 292 |  |  | for (i = 0; i < cm1->ncols; i++) { | 
| 293 |  |  | const COLORV        *cp1 = cm_lval(cm1,dr,i); | 
| 294 |  |  | const COLORV        *cp2 = cm_lval(cm2,i,dc); | 
| 295 |  |  | dp[0] += cp1[0] * cp2[0]; | 
| 296 |  |  | dp[1] += cp1[1] * cp2[1]; | 
| 297 |  |  | dp[2] += cp1[2] * cp2[2]; | 
| 298 |  |  | } | 
| 299 |  |  | } | 
| 300 |  |  | return(cmr); | 
| 301 |  |  | } | 
| 302 |  |  |  | 
| 303 |  |  | /* print out matrix as ASCII text -- no header */ | 
| 304 |  |  | static void | 
| 305 |  |  | cm_print(const CMATRIX *cm, FILE *fp) | 
| 306 |  |  | { | 
| 307 |  |  | int             r, c; | 
| 308 |  |  | const COLORV    *mp = cm->cmem; | 
| 309 |  |  |  | 
| 310 |  |  | for (r = 0; r < cm->nrows; r++) { | 
| 311 |  |  | for (c = 0; c < cm->ncols; c++, mp += 3) | 
| 312 |  |  | fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]); | 
| 313 |  |  | fputc('\n', fp); | 
| 314 |  |  | } | 
| 315 |  |  | } | 
| 316 |  |  |  | 
| 317 | greg | 2.15 | /* Convert a BSDF to our matrix representation */ | 
| 318 | greg | 2.1 | static CMATRIX * | 
| 319 | greg | 2.15 | cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf) | 
| 320 | greg | 2.1 | { | 
| 321 | greg | 2.16 | CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc); | 
| 322 | greg | 2.3 | int     nbadohm = 0; | 
| 323 |  |  | int     nneg = 0; | 
| 324 | greg | 2.1 | int     r, c; | 
| 325 | greg | 2.22 | /* loop over incident angles */ | 
| 326 | greg | 2.16 | for (c = 0; c < cm->ncols; c++) { | 
| 327 | greg | 2.17 | const double    dom = mBSDF_incohm(bsdf,c); | 
| 328 |  |  | /* projected solid angle */ | 
| 329 | greg | 2.16 | nbadohm += (dom <= 0); | 
| 330 |  |  |  | 
| 331 |  |  | for (r = 0; r < cm->nrows; r++) { | 
| 332 |  |  | float   f = mBSDF_value(bsdf,c,r); | 
| 333 |  |  | COLORV  *mp = cm_lval(cm,r,c); | 
| 334 |  |  | /* check BSDF value */ | 
| 335 |  |  | if ((f <= 0) | (dom <= 0)) { | 
| 336 |  |  | nneg += (f < -FTINY); | 
| 337 |  |  | f = .0f; | 
| 338 | greg | 2.3 | } | 
| 339 | greg | 2.16 | copycolor(mp, specCol); | 
| 340 |  |  | scalecolor(mp, f); | 
| 341 |  |  | addcolor(mp, bsdfLamb); | 
| 342 |  |  | scalecolor(mp, dom); | 
| 343 | greg | 2.3 | } | 
| 344 | greg | 2.7 | } | 
| 345 | greg | 2.15 | if (nneg | nbadohm) { | 
| 346 | greg | 2.3 | sprintf(errmsg, | 
| 347 | greg | 2.4 | "BTDF has %d negatives and %d bad incoming solid angles", | 
| 348 |  |  | nneg, nbadohm); | 
| 349 | greg | 2.3 | error(WARNING, errmsg); | 
| 350 |  |  | } | 
| 351 | greg | 2.1 | return(cm); | 
| 352 |  |  | } | 
| 353 |  |  |  | 
| 354 | greg | 2.22 | /* Convert between input and output indices for reciprocity */ | 
| 355 |  |  | static int | 
| 356 |  |  | recip_out_from_in(const SDMat *bsdf, int in_recip) | 
| 357 |  |  | { | 
| 358 |  |  | FVECT   v; | 
| 359 |  |  |  | 
| 360 |  |  | if (!mBSDF_incvec(v, bsdf, in_recip+.5)) | 
| 361 |  |  | return(in_recip);               /* XXX should be error! */ | 
| 362 |  |  | v[2] = -v[2]; | 
| 363 |  |  | return(mBSDF_outndx(bsdf, v)); | 
| 364 |  |  | } | 
| 365 |  |  |  | 
| 366 |  |  | /* Convert between output and input indices for reciprocity */ | 
| 367 |  |  | static int | 
| 368 |  |  | recip_in_from_out(const SDMat *bsdf, int out_recip) | 
| 369 |  |  | { | 
| 370 |  |  | FVECT   v; | 
| 371 |  |  |  | 
| 372 |  |  | if (!mBSDF_outvec(v, bsdf, out_recip+.5)) | 
| 373 |  |  | return(out_recip);              /* XXX should be error! */ | 
| 374 |  |  | v[2] = -v[2]; | 
| 375 |  |  | return(mBSDF_incndx(bsdf, v)); | 
| 376 |  |  | } | 
| 377 |  |  |  | 
| 378 |  |  | /* Convert a BSDF to our matrix representation, applying reciprocity */ | 
| 379 |  |  | static CMATRIX * | 
| 380 |  |  | cm_bsdf_recip(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf) | 
| 381 |  |  | { | 
| 382 |  |  | CMATRIX *cm = cm_alloc(bsdf->ninc, bsdf->nout); | 
| 383 |  |  | int     nbadohm = 0; | 
| 384 |  |  | int     nneg = 0; | 
| 385 |  |  | int     r, c; | 
| 386 |  |  | /* loop over incident angles */ | 
| 387 |  |  | for (c = 0; c < cm->ncols; c++) { | 
| 388 |  |  | const int       ro = recip_out_from_in(bsdf,c); | 
| 389 |  |  | const double    dom = mBSDF_outohm(bsdf,ro); | 
| 390 |  |  | /* projected solid angle */ | 
| 391 |  |  | nbadohm += (dom <= 0); | 
| 392 |  |  |  | 
| 393 |  |  | for (r = 0; r < cm->nrows; r++) { | 
| 394 |  |  | const int       ri = recip_in_from_out(bsdf,r); | 
| 395 |  |  | float           f = mBSDF_value(bsdf,ri,ro); | 
| 396 |  |  | COLORV          *mp = cm_lval(cm,r,c); | 
| 397 |  |  | /* check BSDF value */ | 
| 398 |  |  | if ((f <= 0) | (dom <= 0)) { | 
| 399 |  |  | nneg += (f < -FTINY); | 
| 400 |  |  | f = .0f; | 
| 401 |  |  | } | 
| 402 |  |  | copycolor(mp, specCol); | 
| 403 |  |  | scalecolor(mp, f); | 
| 404 |  |  | addcolor(mp, bsdfLamb); | 
| 405 |  |  | scalecolor(mp, dom); | 
| 406 |  |  | } | 
| 407 |  |  | } | 
| 408 |  |  | if (nneg | nbadohm) { | 
| 409 |  |  | sprintf(errmsg, | 
| 410 |  |  | "BTDF has %d negatives and %d bad incoming solid angles", | 
| 411 |  |  | nneg, nbadohm); | 
| 412 |  |  | error(WARNING, errmsg); | 
| 413 |  |  | } | 
| 414 |  |  | return(cm); | 
| 415 |  |  | } | 
| 416 |  |  |  | 
| 417 | greg | 2.15 | /* Load and convert a matrix BSDF from the given XML file */ | 
| 418 |  |  | static CMATRIX * | 
| 419 | greg | 2.20 | cm_loadBSDF(char *fname, COLOR cLamb) | 
| 420 | greg | 2.15 | { | 
| 421 | greg | 2.21 | CMATRIX         *Tmat; | 
| 422 |  |  | char            *fpath; | 
| 423 | greg | 2.22 | int             recip; | 
| 424 | greg | 2.21 | SDError         ec; | 
| 425 |  |  | SDData          myBSDF; | 
| 426 |  |  | SDSpectralDF    *tdf; | 
| 427 |  |  | COLOR           bsdfLamb, specCol; | 
| 428 | greg | 2.15 | /* find path to BSDF file */ | 
| 429 |  |  | fpath = getpath(fname, getrlibpath(), R_OK); | 
| 430 |  |  | if (fpath == NULL) { | 
| 431 |  |  | sprintf(errmsg, "cannot find BSDF file '%s'", fname); | 
| 432 |  |  | error(USER, errmsg); | 
| 433 |  |  | } | 
| 434 | greg | 2.16 | SDclearBSDF(&myBSDF, fname);    /* load XML and check type */ | 
| 435 | greg | 2.15 | ec = SDloadFile(&myBSDF, fpath); | 
| 436 |  |  | if (ec) | 
| 437 |  |  | error(USER, transSDError(ec)); | 
| 438 | greg | 2.20 | ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb); | 
| 439 | greg | 2.22 | recip = (myBSDF.tf == NULL); | 
| 440 |  |  | tdf = recip ? myBSDF.tb : myBSDF.tf; | 
| 441 | greg | 2.21 | if (tdf == NULL) {              /* no non-Lambertian transmission? */ | 
| 442 | greg | 2.20 | if (cLamb != NULL) | 
| 443 |  |  | copycolor(cLamb, bsdfLamb); | 
| 444 |  |  | SDfreeBSDF(&myBSDF); | 
| 445 |  |  | return(NULL); | 
| 446 |  |  | } | 
| 447 | greg | 2.21 | if (tdf->ncomp != 1 || tdf->comp[0].func != &SDhandleMtx) { | 
| 448 | greg | 2.15 | sprintf(errmsg, "unsupported BSDF '%s'", fpath); | 
| 449 |  |  | error(USER, errmsg); | 
| 450 |  |  | } | 
| 451 |  |  | /* convert BTDF to matrix */ | 
| 452 | greg | 2.21 | ccy2rgb(&tdf->comp[0].cspec[0], 1., specCol); | 
| 453 | greg | 2.22 | Tmat = recip ? cm_bsdf_recip(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist) | 
| 454 |  |  | : cm_bsdf(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist); | 
| 455 | greg | 2.20 | if (cLamb != NULL)              /* Lambertian is included */ | 
| 456 |  |  | setcolor(cLamb, .0, .0, .0); | 
| 457 | greg | 2.15 | /* free BSDF and return */ | 
| 458 |  |  | SDfreeBSDF(&myBSDF); | 
| 459 |  |  | return(Tmat); | 
| 460 |  |  | } | 
| 461 |  |  |  | 
| 462 | greg | 2.1 | /* Sum together a set of images and write result to stdout */ | 
| 463 |  |  | static int | 
| 464 | greg | 2.13 | sum_images(const char *fspec, const CMATRIX *cv, FILE *fout) | 
| 465 | greg | 2.1 | { | 
| 466 |  |  | int     myDT = DTfromHeader; | 
| 467 |  |  | CMATRIX *pmat; | 
| 468 |  |  | COLOR   *scanline; | 
| 469 |  |  | int     myXR, myYR; | 
| 470 |  |  | int     i, y; | 
| 471 |  |  |  | 
| 472 |  |  | if (cv->ncols != 1) | 
| 473 |  |  | error(INTERNAL, "expected vector in sum_images()"); | 
| 474 |  |  | for (i = 0; i < cv->nrows; i++) { | 
| 475 |  |  | const COLORV    *scv = cv_lval(cv,i); | 
| 476 |  |  | char            fname[1024]; | 
| 477 |  |  | FILE            *fp; | 
| 478 |  |  | int             dt, xr, yr; | 
| 479 |  |  | COLORV          *psp; | 
| 480 | greg | 2.14 | /* check for zero */ | 
| 481 | greg | 2.19 | if ((scv[RED] == 0) & (scv[GRN] == 0) & (scv[BLU] == 0) && | 
| 482 |  |  | (myDT != DTfromHeader) | (i < cv->nrows-1)) | 
| 483 | greg | 2.14 | continue; | 
| 484 | greg | 2.1 | /* open next picture */ | 
| 485 |  |  | sprintf(fname, fspec, i); | 
| 486 |  |  | if ((fp = fopen(fname, "r")) == NULL) { | 
| 487 |  |  | sprintf(errmsg, "cannot open picture '%s'", fname); | 
| 488 |  |  | error(SYSTEM, errmsg); | 
| 489 |  |  | } | 
| 490 |  |  | SET_FILE_BINARY(fp); | 
| 491 |  |  | dt = getDTfromHeader(fp); | 
| 492 |  |  | if ((dt != DTrgbe) & (dt != DTxyze) || | 
| 493 |  |  | !fscnresolu(&xr, &yr, fp)) { | 
| 494 |  |  | sprintf(errmsg, "file '%s' not a picture", fname); | 
| 495 |  |  | error(USER, errmsg); | 
| 496 |  |  | } | 
| 497 |  |  | if (myDT == DTfromHeader) {             /* on first one */ | 
| 498 |  |  | myDT = dt; | 
| 499 |  |  | myXR = xr; myYR = yr; | 
| 500 |  |  | scanline = (COLOR *)malloc(sizeof(COLOR)*myXR); | 
| 501 |  |  | if (scanline == NULL) | 
| 502 |  |  | error(SYSTEM, "out of memory in sum_images()"); | 
| 503 |  |  | pmat = cm_alloc(myYR, myXR); | 
| 504 |  |  | memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR); | 
| 505 | greg | 2.2 | /* finish header */ | 
| 506 | greg | 2.13 | fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout); | 
| 507 |  |  | fputc('\n', fout); | 
| 508 |  |  | fprtresolu(myXR, myYR, fout); | 
| 509 |  |  | fflush(fout); | 
| 510 | greg | 2.1 | } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) { | 
| 511 |  |  | sprintf(errmsg, "picture '%s' format/size mismatch", | 
| 512 |  |  | fname); | 
| 513 |  |  | error(USER, errmsg); | 
| 514 |  |  | } | 
| 515 |  |  | psp = pmat->cmem; | 
| 516 |  |  | for (y = 0; y < yr; y++) {              /* read it in */ | 
| 517 |  |  | int     x; | 
| 518 |  |  | if (freadscan(scanline, xr, fp) < 0) { | 
| 519 |  |  | sprintf(errmsg, "error reading picture '%s'", | 
| 520 |  |  | fname); | 
| 521 |  |  | error(SYSTEM, errmsg); | 
| 522 |  |  | } | 
| 523 |  |  | /* sum in scanline */ | 
| 524 |  |  | for (x = 0; x < xr; x++, psp += 3) { | 
| 525 |  |  | multcolor(scanline[x], scv); | 
| 526 |  |  | addcolor(psp, scanline[x]); | 
| 527 |  |  | } | 
| 528 |  |  | } | 
| 529 |  |  | fclose(fp);                             /* done this picture */ | 
| 530 |  |  | } | 
| 531 |  |  | free(scanline); | 
| 532 |  |  | /* write scanlines */ | 
| 533 |  |  | for (y = 0; y < myYR; y++) | 
| 534 | greg | 2.13 | if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0) | 
| 535 | greg | 2.1 | return(0); | 
| 536 |  |  | cm_free(pmat);                                  /* all done */ | 
| 537 | greg | 2.13 | return(fflush(fout) == 0); | 
| 538 | greg | 2.1 | } | 
| 539 |  |  |  | 
| 540 | greg | 2.11 | /* check to see if a string contains a %d or %o specification */ | 
| 541 |  |  | static int | 
| 542 |  |  | hasNumberFormat(const char *s) | 
| 543 | greg | 2.1 | { | 
| 544 | greg | 2.18 | while (*s) { | 
| 545 |  |  | while (*s != '%') | 
| 546 |  |  | if (!*s++) | 
| 547 |  |  | return(0); | 
| 548 |  |  | if (*++s == '%') {              /* ignore "%%" */ | 
| 549 |  |  | ++s; | 
| 550 |  |  | continue; | 
| 551 |  |  | } | 
| 552 |  |  | while (isdigit(*s))             /* field length */ | 
| 553 |  |  | ++s; | 
| 554 |  |  | /* field we'll use? */ | 
| 555 |  |  | if ((*s == 'd') | (*s == 'i') | (*s == 'o') | | 
| 556 |  |  | (*s == 'x') | (*s == 'X')) | 
| 557 |  |  | return(1); | 
| 558 |  |  | } | 
| 559 |  |  | return(0);                              /* didn't find one */ | 
| 560 | greg | 2.1 | } | 
| 561 |  |  |  | 
| 562 |  |  | int | 
| 563 |  |  | main(int argc, char *argv[]) | 
| 564 |  |  | { | 
| 565 | greg | 2.20 | CMATRIX                 *cvec;          /* component vector result */ | 
| 566 | greg | 2.1 |  | 
| 567 |  |  | progname = argv[0]; | 
| 568 |  |  |  | 
| 569 | greg | 2.12 | if ((argc < 2) | (argc > 5)) { | 
| 570 |  |  | fprintf(stderr, "Usage: %s DCspec [tregvec]\n", progname); | 
| 571 |  |  | fprintf(stderr, "   or: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n", | 
| 572 |  |  | progname); | 
| 573 | greg | 2.1 | return(1); | 
| 574 |  |  | } | 
| 575 | greg | 2.12 |  | 
| 576 | greg | 2.13 | if (argc > 3) {                         /* VTDs expression */ | 
| 577 | greg | 2.12 | CMATRIX *svec, *Dmat, *Tmat, *ivec; | 
| 578 | greg | 2.20 | COLOR   tLamb; | 
| 579 | greg | 2.12 | /* get sky vector */ | 
| 580 |  |  | svec = cm_load(argv[4], 0, 1, DTascii); | 
| 581 |  |  | /* load BSDF */ | 
| 582 | greg | 2.20 | Tmat = cm_loadBSDF(argv[2], tLamb); | 
| 583 | greg | 2.5 | /* load Daylight matrix */ | 
| 584 | greg | 2.20 | Dmat = cm_load(argv[3], Tmat==NULL ? 0 : Tmat->ncols, | 
| 585 |  |  | svec->nrows, DTfromHeader); | 
| 586 | greg | 2.1 | /* multiply vector through */ | 
| 587 | greg | 2.12 | ivec = cm_multiply(Dmat, svec); | 
| 588 |  |  | cm_free(Dmat); cm_free(svec); | 
| 589 | greg | 2.20 | if (Tmat == NULL) {             /* diffuse only */ | 
| 590 |  |  | cvec = cm_scale(ivec, tLamb); | 
| 591 |  |  | } else {                        /* else apply BTDF matrix */ | 
| 592 |  |  | cvec = cm_multiply(Tmat, ivec); | 
| 593 |  |  | cm_free(Tmat); | 
| 594 |  |  | } | 
| 595 |  |  | cm_free(ivec); | 
| 596 | greg | 2.12 | } else {                                /* else just use sky vector */ | 
| 597 |  |  | cvec = cm_load(argv[2], 0, 1, DTascii); | 
| 598 |  |  | } | 
| 599 |  |  |  | 
| 600 | greg | 2.11 | if (hasNumberFormat(argv[1])) {         /* generating image */ | 
| 601 | greg | 2.1 | SET_FILE_BINARY(stdout); | 
| 602 |  |  | newheader("RADIANCE", stdout); | 
| 603 |  |  | printargs(argc, argv, stdout); | 
| 604 |  |  | fputnow(stdout); | 
| 605 | greg | 2.12 | if (!sum_images(argv[1], cvec, stdout)) | 
| 606 | greg | 2.1 | return(1); | 
| 607 |  |  | } else {                                /* generating vector */ | 
| 608 |  |  | CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader); | 
| 609 |  |  | CMATRIX *rvec = cm_multiply(Vmat, cvec); | 
| 610 |  |  | cm_free(Vmat); | 
| 611 |  |  | cm_print(rvec, stdout); | 
| 612 |  |  | cm_free(rvec); | 
| 613 |  |  | } | 
| 614 |  |  | cm_free(cvec);                          /* final clean-up */ | 
| 615 |  |  | return(0); | 
| 616 |  |  | } |