Back to index

radiance  4R0+20100331
dctimestep.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: dctimestep.c,v 2.11 2009/12/09 20:33:48 greg Exp $";
00003 #endif
00004 /*
00005  * Compute time-step result using Daylight Coefficient method.
00006  *
00007  *     G. Ward
00008  */
00009 
00010 #include <ctype.h>
00011 #include "standard.h"
00012 #include "platform.h"
00013 #include "color.h"
00014 #include "resolu.h"
00015 #include "bsdf.h"
00016 
00017 char   *progname;                  /* global argv[0] */
00018 
00019 /* Data types for file loading */
00020 enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze};
00021 
00022 /* A color coefficient matrix -- vectors have ncols==1 */
00023 typedef struct {
00024        int    nrows, ncols;
00025        COLORV cmem[3];             /* extends struct */
00026 } CMATRIX;
00027 
00028 #define COLSPEC      (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf")
00029 
00030 #define cm_lval(cm,r,c)     ((cm)->cmem + 3*((r)*(cm)->ncols + (c)))
00031 
00032 #define cv_lval(cm,i)       ((cm)->cmem + 3*(i))
00033 
00034 /* Allocate a color coefficient matrix */
00035 static CMATRIX *
00036 cm_alloc(int nrows, int ncols)
00037 {
00038        CMATRIX       *cm;
00039 
00040        if ((nrows <= 0) | (ncols <= 0))
00041               error(USER, "attempt to create empty matrix");
00042        cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
00043                             3*sizeof(COLORV)*(nrows*ncols - 1));
00044        if (cm == NULL)
00045               error(SYSTEM, "out of memory in cm_alloc()");
00046        cm->nrows = nrows;
00047        cm->ncols = ncols;
00048        return(cm);
00049 }
00050 
00051 #define cm_free(cm)  free(cm)
00052 
00053 /* Resize color coefficient matrix */
00054 static CMATRIX *
00055 cm_resize(CMATRIX *cm, int nrows)
00056 {
00057        if (nrows == cm->nrows)
00058               return(cm);
00059        if (nrows <= 0) {
00060               cm_free(cm);
00061               return(NULL);
00062        }
00063        cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
00064                      3*sizeof(COLORV)*(nrows*cm->ncols - 1));
00065        if (cm == NULL)
00066               error(SYSTEM, "out of memory in cm_resize()");
00067        cm->nrows = nrows;
00068        return(cm);
00069 }
00070 
00071 /* Load header to obtain data type */
00072 static int
00073 getDT(char *s, void *p)
00074 {
00075        char   fmt[32];
00076        
00077        if (formatval(fmt, s)) {
00078               if (!strcmp(fmt, "ascii"))
00079                      *((int *)p) = DTascii;
00080               else if (!strcmp(fmt, "float"))
00081                      *((int *)p) = DTfloat;
00082               else if (!strcmp(fmt, "double"))
00083                      *((int *)p) = DTdouble;
00084               else if (!strcmp(fmt, COLRFMT))
00085                      *((int *)p) = DTrgbe;
00086               else if (!strcmp(fmt, CIEFMT))
00087                      *((int *)p) = DTxyze;
00088        }
00089        return(0);
00090 }
00091 
00092 static int
00093 getDTfromHeader(FILE *fp)
00094 {
00095        int    dt = DTfromHeader;
00096        
00097        if (getheader(fp, getDT, &dt) < 0)
00098               error(SYSTEM, "header read error");
00099        if (dt == DTfromHeader)
00100               error(USER, "missing data format in header");
00101        return(dt);
00102 }
00103 
00104 /* Allocate and load a matrix from the given file (or stdin if NULL) */
00105 static CMATRIX *
00106 cm_load(const char *fname, int nrows, int ncols, int dtype)
00107 {
00108        CMATRIX       *cm;
00109        FILE   *fp = stdin;
00110 
00111        if (fname == NULL)
00112               fname = "<stdin>";
00113        else if ((fp = fopen(fname, "r")) == NULL) {
00114               sprintf(errmsg, "cannot open file '%s'", fname);
00115               error(SYSTEM, errmsg);
00116        }
00117        if (dtype != DTascii)
00118               SET_FILE_BINARY(fp);
00119        if (dtype == DTfromHeader)
00120               dtype = getDTfromHeader(fp);
00121        switch (dtype) {
00122        case DTascii:
00123        case DTfloat:
00124        case DTdouble:
00125               break;
00126        default:
00127               error(USER, "unexpected data type in cm_load()");
00128        }
00129        if (nrows <= 0) {                  /* don't know length? */
00130               int    guessrows = 147;     /* usually big enough */
00131               if ((dtype != DTascii) & (fp != stdin)) {
00132                      long   startpos = ftell(fp);
00133                      if (fseek(fp, 0L, SEEK_END) == 0) {
00134                             long   endpos = ftell(fp);
00135                             long   elemsiz = 3*(dtype==DTfloat ?
00136                                        sizeof(float) : sizeof(double));
00137 
00138                             if ((endpos - startpos) % (ncols*elemsiz)) {
00139                                    sprintf(errmsg,
00140                                    "improper length for binary file '%s'",
00141                                                  fname);
00142                                    error(USER, errmsg);
00143                             }
00144                             guessrows = (endpos - startpos)/(ncols*elemsiz);
00145                             if (fseek(fp, startpos, SEEK_SET) < 0) {
00146                                    sprintf(errmsg,
00147                                           "fseek() error on file '%s'",
00148                                                  fname);
00149                                    error(SYSTEM, errmsg);
00150                             }
00151                             nrows = guessrows;   /* we're confident */
00152                      }
00153               }
00154               cm = cm_alloc(guessrows, ncols);
00155        } else
00156               cm = cm_alloc(nrows, ncols);
00157        if (cm == NULL)
00158               return(NULL);
00159        if (dtype == DTascii) {                          /* read text file */
00160               int    maxrow = (nrows > 0 ? nrows : 32000);
00161               int    r, c;
00162               for (r = 0; r < maxrow; r++) {
00163                   if (r >= cm->nrows)                   /* need more space? */
00164                      cm = cm_resize(cm, 2*cm->nrows);
00165                   for (c = 0; c < ncols; c++) {
00166                       COLORV       *cv = cm_lval(cm,r,c);
00167                      if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
00168                             if ((nrows <= 0) & (r > 0) & !c) {
00169                                    cm = cm_resize(cm, maxrow=r);
00170                                    break;
00171                             } else
00172                                    goto EOFerror;
00173                   }
00174               }
00175               while ((c = getc(fp)) != EOF)
00176                      if (!isspace(c)) {
00177                             sprintf(errmsg,
00178                             "unexpected data at end of ascii file %s",
00179                                           fname);
00180                             error(WARNING, errmsg);
00181                             break;
00182                      }
00183        } else {                                  /* read binary file */
00184               if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
00185                                                  sizeof(double))) {
00186                      int    nread = 0;
00187                      do {                        /* read all we can */
00188                             nread += fread(cm->cmem + 3*nread,
00189                                           3*sizeof(COLORV),
00190                                           cm->nrows*cm->ncols - nread,
00191                                           fp);
00192                             if (nrows <= 0) {    /* unknown length */
00193                                    if (nread == cm->nrows*cm->ncols)
00194                                                  /* need more space? */
00195                                           cm = cm_resize(cm, 2*cm->nrows);
00196                                    else if (nread && !(nread % cm->ncols))
00197                                                  /* seem to be  done */
00198                                           cm = cm_resize(cm, nread/cm->ncols);
00199                                    else          /* ended mid-row */
00200                                           goto EOFerror;
00201                             } else if (nread < cm->nrows*cm->ncols)
00202                                    goto EOFerror;
00203                      } while (nread < cm->nrows*cm->ncols);
00204 
00205               } else if (dtype == DTdouble) {
00206                      double dc[3];               /* load from double */
00207                      COLORV *cvp = cm->cmem;
00208                      int    n = nrows*ncols;
00209 
00210                      if (n <= 0)
00211                             goto not_handled;
00212                      while (n--) {
00213                             if (fread(dc, sizeof(double), 3, fp) != 3)
00214                                    goto EOFerror;
00215                             copycolor(cvp, dc);
00216                             cvp += 3;
00217                      }
00218               } else /* dtype == DTfloat */ {
00219                      float  fc[3];               /* load from float */
00220                      COLORV *cvp = cm->cmem;
00221                      int    n = nrows*ncols;
00222 
00223                      if (n <= 0)
00224                             goto not_handled;
00225                      while (n--) {
00226                             if (fread(fc, sizeof(float), 3, fp) != 3)
00227                                    goto EOFerror;
00228                             copycolor(cvp, fc);
00229                             cvp += 3;
00230                      }
00231               }
00232               if (getc(fp) != EOF) {
00233                             sprintf(errmsg,
00234                             "unexpected data at end of binary file %s",
00235                                           fname);
00236                             error(WARNING, errmsg);
00237               }
00238        }
00239        if (fp != stdin)
00240               fclose(fp);
00241        return(cm);
00242 EOFerror:
00243        sprintf(errmsg, "unexpected EOF reading %s",
00244                      fname);
00245        error(USER, errmsg);
00246 not_handled:
00247        error(INTERNAL, "unhandled data size or length in cm_load()");
00248        return(NULL); /* gratis return */
00249 }
00250 
00251 /* Multiply two matrices (or a matrix and a vector) and allocate the result*/
00252 static CMATRIX *
00253 cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
00254 {
00255        CMATRIX       *cmr;
00256        int    dr, dc, i;
00257 
00258        if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
00259               error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
00260        cmr = cm_alloc(cm1->nrows, cm2->ncols);
00261        if (cmr == NULL)
00262               return(NULL);
00263        for (dr = 0; dr < cmr->nrows; dr++)
00264            for (dc = 0; dc < cmr->ncols; dc++) {
00265               COLORV *dp = cm_lval(cmr,dr,dc);
00266               dp[0] = dp[1] = dp[2] = 0;
00267               for (i = 0; i < cm1->ncols; i++) {
00268                   const COLORV     *cp1 = cm_lval(cm1,dr,i);
00269                   const COLORV     *cp2 = cm_lval(cm2,i,dc);
00270                   dp[0] += cp1[0] * cp2[0];
00271                   dp[1] += cp1[1] * cp2[1];
00272                   dp[2] += cp1[2] * cp2[2];
00273               }
00274            }
00275        return(cmr);
00276 }
00277 
00278 /* print out matrix as ASCII text -- no header */
00279 static void
00280 cm_print(const CMATRIX *cm, FILE *fp)
00281 {
00282        int           r, c;
00283        const COLORV  *mp = cm->cmem;
00284        
00285        for (r = 0; r < cm->nrows; r++) {
00286               for (c = 0; c < cm->ncols; c++, mp += 3)
00287                      fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
00288               fputc('\n', fp);
00289        }
00290 }
00291 
00292 /* convert a BSDF to our matrix representation */
00293 static CMATRIX *
00294 cm_bsdf(const struct BSDF_data *bsdf)
00295 {
00296        CMATRIX       *cm = cm_alloc(bsdf->nout, bsdf->ninc);
00297        int    nbadohm = 0;
00298        int    nneg = 0;
00299        int    r, c;
00300        
00301        for (c = 0; c < cm->ncols; c++) {
00302               float  dom = getBSDF_incohm(bsdf,c);
00303               FVECT  v;
00304               
00305               if (dom <= .0) {
00306                      nbadohm++;
00307                      continue;
00308               }
00309               if (!getBSDF_incvec(v,bsdf,c) || v[2] > FTINY)
00310                      error(USER, "illegal incoming BTDF direction");
00311               dom *= -v[2];
00312 
00313               for (r = 0; r < cm->nrows; r++) {
00314                      float  f = BSDF_value(bsdf,c,r);
00315                      COLORV *mp = cm_lval(cm,r,c);
00316 
00317                      if (f <= .0) {
00318                             nneg += (f < -FTINY);
00319                             f = .0f;
00320                      }
00321                      mp[0] = mp[1] = mp[2] = f * dom;
00322               }
00323        }
00324        if (nneg || nbadohm) {
00325               sprintf(errmsg,
00326                   "BTDF has %d negatives and %d bad incoming solid angles",
00327                             nneg, nbadohm);
00328               error(WARNING, errmsg);
00329        }
00330        return(cm);
00331 }
00332 
00333 /* Sum together a set of images and write result to stdout */
00334 static int
00335 sum_images(const char *fspec, const CMATRIX *cv)
00336 {
00337        int    myDT = DTfromHeader;
00338        CMATRIX       *pmat;
00339        COLOR  *scanline;
00340        int    myXR, myYR;
00341        int    i, y;
00342 
00343        if (cv->ncols != 1)
00344               error(INTERNAL, "expected vector in sum_images()");
00345        for (i = 0; i < cv->nrows; i++) {
00346               const COLORV  *scv = cv_lval(cv,i);
00347               char          fname[1024];
00348               FILE          *fp;
00349               int           dt, xr, yr;
00350               COLORV        *psp;
00351                                                  /* open next picture */
00352               sprintf(fname, fspec, i);
00353               if ((fp = fopen(fname, "r")) == NULL) {
00354                      sprintf(errmsg, "cannot open picture '%s'", fname);
00355                      error(SYSTEM, errmsg);
00356               }
00357               SET_FILE_BINARY(fp);
00358               dt = getDTfromHeader(fp);
00359               if ((dt != DTrgbe) & (dt != DTxyze) ||
00360                             !fscnresolu(&xr, &yr, fp)) {
00361                      sprintf(errmsg, "file '%s' not a picture", fname);
00362                      error(USER, errmsg);
00363               }
00364               if (myDT == DTfromHeader) {        /* on first one */
00365                      myDT = dt;
00366                      myXR = xr; myYR = yr;
00367                      scanline = (COLOR *)malloc(sizeof(COLOR)*myXR);
00368                      if (scanline == NULL)
00369                             error(SYSTEM, "out of memory in sum_images()");
00370                      pmat = cm_alloc(myYR, myXR);
00371                      memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
00372                                                  /* finish header */
00373                      fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
00374                      fputc('\n', stdout);
00375                      fprtresolu(myXR, myYR, stdout);
00376                      fflush(stdout);
00377               } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
00378                      sprintf(errmsg, "picture '%s' format/size mismatch",
00379                                    fname);
00380                      error(USER, errmsg);
00381               }
00382               psp = pmat->cmem;
00383               for (y = 0; y < yr; y++) {         /* read it in */
00384                      int    x;
00385                      if (freadscan(scanline, xr, fp) < 0) {
00386                             sprintf(errmsg, "error reading picture '%s'",
00387                                           fname);
00388                             error(SYSTEM, errmsg);
00389                      }
00390                                                  /* sum in scanline */
00391                      for (x = 0; x < xr; x++, psp += 3) {
00392                             multcolor(scanline[x], scv);
00393                             addcolor(psp, scanline[x]);
00394                      }
00395               }
00396               fclose(fp);                        /* done this picture */
00397        }
00398        free(scanline);
00399                                                  /* write scanlines */
00400        for (y = 0; y < myYR; y++)
00401               if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
00402                      return(0);
00403        cm_free(pmat);                                   /* all done */
00404        return(fflush(stdout) == 0);
00405 }
00406 
00407 /* check to see if a string contains a %d or %o specification */
00408 static int
00409 hasNumberFormat(const char *s)
00410 {
00411        while (*s && *s != '%')
00412               s++;
00413        if (!*s)
00414               return(0);
00415        do
00416               ++s;
00417        while (isdigit(*s));
00418 
00419        return(*s == 'd' | *s == 'i' | *s == 'o' | *s == 'x' | *s == 'X');
00420 }
00421 
00422 int
00423 main(int argc, char *argv[])
00424 {
00425        CMATRIX                     *tvec, *Dmat, *Tmat, *ivec, *cvec;
00426        struct BSDF_data     *btdf;
00427 
00428        progname = argv[0];
00429 
00430        if ((argc < 4) | (argc > 5)) {
00431               fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
00432                             progname);
00433               return(1);
00434        }
00435                                           /* load Tregenza vector */
00436        tvec = cm_load(argv[4], 0, 1, DTascii);   /* argv[4]==NULL iff argc==4 */
00437                                           /* load BTDF */
00438        btdf = load_BSDF(argv[2]);
00439        if (btdf == NULL)
00440               return(1);
00441                                           /* load Daylight matrix */
00442        Dmat = cm_load(argv[3], btdf->ninc, tvec->nrows, DTfromHeader);
00443                                           /* multiply vector through */
00444        ivec = cm_multiply(Dmat, tvec);
00445        cm_free(Dmat); cm_free(tvec);
00446        Tmat = cm_bsdf(btdf);                     /* convert BTDF to matrix */
00447        free_BSDF(btdf);
00448        cvec = cm_multiply(Tmat, ivec);           /* cvec = component vector */
00449        cm_free(Tmat); cm_free(ivec);
00450        if (hasNumberFormat(argv[1])) {           /* generating image */
00451               SET_FILE_BINARY(stdout);
00452               newheader("RADIANCE", stdout);
00453               printargs(argc, argv, stdout);
00454               fputnow(stdout);
00455               if (!sum_images(argv[1], cvec))
00456                      return(1);
00457        } else {                           /* generating vector */
00458               CMATRIX       *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
00459               CMATRIX       *rvec = cm_multiply(Vmat, cvec);
00460               cm_free(Vmat);
00461               cm_print(rvec, stdout);
00462               cm_free(rvec);
00463        }
00464        cm_free(cvec);                            /* final clean-up */
00465        return(0);
00466 }