ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
(Generate patch)

Comparing ray/src/util/dctimestep.c (file contents):
Revision 2.9 by greg, Sun Jun 21 21:53:10 2009 UTC vs.
Revision 2.22 by greg, Tue Oct 23 22:24:09 2012 UTC

# Line 10 | Line 10 | static const char RCSid[] = "$Id$";
10   #include <ctype.h>
11   #include "standard.h"
12   #include "platform.h"
13 + #include "paths.h"
14   #include "color.h"
15   #include "resolu.h"
16   #include "bsdf.h"
17 + #include "bsdf_m.h"
18  
19   char    *progname;                      /* global argv[0] */
20  
# Line 38 | Line 40 | cm_alloc(int nrows, int ncols)
40          CMATRIX *cm;
41  
42          if ((nrows <= 0) | (ncols <= 0))
43 <                return(NULL);
43 >                error(USER, "attempt to create empty matrix");
44          cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
45                                  3*sizeof(COLORV)*(nrows*ncols - 1));
46          if (cm == NULL)
# Line 108 | Line 110 | cm_load(const char *fname, int nrows, int ncols, int d
110          CMATRIX *cm;
111          FILE    *fp = stdin;
112  
113 +        if (ncols <= 0)
114 +                error(USER, "Non-positive number of columns");
115          if (fname == NULL)
116                  fname = "<stdin>";
117          else if ((fp = fopen(fname, "r")) == NULL) {
# Line 248 | Line 252 | not_handled:
252          return(NULL);   /* gratis return */
253   }
254  
255 < /* Multiply two matrices (or a matrix and a vector) and allocate the result*/
255 > /* 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 + static CMATRIX *
278   cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
279   {
280          CMATRIX *cmr;
# Line 289 | Line 314 | cm_print(const CMATRIX *cm, FILE *fp)
314          }
315   }
316  
317 < /* convert a BSDF to our matrix representation */
317 > /* Convert a BSDF to our matrix representation */
318   static CMATRIX *
319 < cm_bsdf(const struct BSDF_data *bsdf)
319 > cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
320   {
321          CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
322          int     nbadohm = 0;
323          int     nneg = 0;
324          int     r, c;
325 <        
325 >                                        /* loop over incident angles */
326          for (c = 0; c < cm->ncols; c++) {
327 <                float   dom = getBSDF_incohm(bsdf,c);
328 <                FVECT   v;
329 <                
305 <                if (dom <= .0) {
306 <                        nbadohm++;
307 <                        continue;
308 <                }
309 <                if (!getBSDF_incvec(v,bsdf,c) || v[2] > FTINY)
310 <                        error(USER, "illegal incoming BTDF direction");
311 <                dom *= -v[2];
327 >                const double    dom = mBSDF_incohm(bsdf,c);
328 >                                        /* projected solid angle */
329 >                nbadohm += (dom <= 0);
330  
331                  for (r = 0; r < cm->nrows; r++) {
332 <                        float   f = BSDF_value(bsdf,c,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 +                        }
339 +                        copycolor(mp, specCol);
340 +                        scalecolor(mp, f);
341 +                        addcolor(mp, bsdfLamb);
342 +                        scalecolor(mp, dom);
343 +                }
344 +        }
345 +        if (nneg | nbadohm) {
346 +                sprintf(errmsg,
347 +                    "BTDF has %d negatives and %d bad incoming solid angles",
348 +                                nneg, nbadohm);
349 +                error(WARNING, errmsg);
350 +        }
351 +        return(cm);
352 + }
353  
354 <                        if (f <= .0) {
354 > /* 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 <                        mp[0] = mp[1] = mp[2] = f * dom;
402 >                        copycolor(mp, specCol);
403 >                        scalecolor(mp, f);
404 >                        addcolor(mp, bsdfLamb);
405 >                        scalecolor(mp, dom);
406                  }
407          }
408 <        if (nneg || nbadohm) {
408 >        if (nneg | nbadohm) {
409                  sprintf(errmsg,
410                      "BTDF has %d negatives and %d bad incoming solid angles",
411                                  nneg, nbadohm);
# Line 330 | Line 414 | cm_bsdf(const struct BSDF_data *bsdf)
414          return(cm);
415   }
416  
417 + /* Load and convert a matrix BSDF from the given XML file */
418 + static CMATRIX *
419 + cm_loadBSDF(char *fname, COLOR cLamb)
420 + {
421 +        CMATRIX         *Tmat;
422 +        char            *fpath;
423 +        int             recip;
424 +        SDError         ec;
425 +        SDData          myBSDF;
426 +        SDSpectralDF    *tdf;
427 +        COLOR           bsdfLamb, specCol;
428 +                                        /* 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 +        SDclearBSDF(&myBSDF, fname);    /* load XML and check type */
435 +        ec = SDloadFile(&myBSDF, fpath);
436 +        if (ec)
437 +                error(USER, transSDError(ec));
438 +        ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb);
439 +        recip = (myBSDF.tf == NULL);
440 +        tdf = recip ? myBSDF.tb : myBSDF.tf;
441 +        if (tdf == NULL) {              /* no non-Lambertian transmission? */
442 +                if (cLamb != NULL)
443 +                        copycolor(cLamb, bsdfLamb);
444 +                SDfreeBSDF(&myBSDF);
445 +                return(NULL);
446 +        }
447 +        if (tdf->ncomp != 1 || tdf->comp[0].func != &SDhandleMtx) {
448 +                sprintf(errmsg, "unsupported BSDF '%s'", fpath);
449 +                error(USER, errmsg);
450 +        }
451 +                                        /* convert BTDF to matrix */
452 +        ccy2rgb(&tdf->comp[0].cspec[0], 1., specCol);
453 +        Tmat = recip ? cm_bsdf_recip(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist)
454 +                        : cm_bsdf(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist);
455 +        if (cLamb != NULL)              /* Lambertian is included */
456 +                setcolor(cLamb, .0, .0, .0);
457 +                                        /* free BSDF and return */
458 +        SDfreeBSDF(&myBSDF);
459 +        return(Tmat);
460 + }
461 +
462   /* Sum together a set of images and write result to stdout */
463   static int
464 < sum_images(const char *fspec, const CMATRIX *cv)
464 > sum_images(const char *fspec, const CMATRIX *cv, FILE *fout)
465   {
466          int     myDT = DTfromHeader;
467          CMATRIX *pmat;
# Line 348 | Line 477 | sum_images(const char *fspec, const CMATRIX *cv)
477                  FILE            *fp;
478                  int             dt, xr, yr;
479                  COLORV          *psp;
480 +                                                        /* check for zero */
481 +                if ((scv[RED] == 0) & (scv[GRN] == 0) & (scv[BLU] == 0) &&
482 +                                (myDT != DTfromHeader) | (i < cv->nrows-1))
483 +                        continue;
484                                                          /* open next picture */
485                  sprintf(fname, fspec, i);
486                  if ((fp = fopen(fname, "r")) == NULL) {
# Line 370 | Line 503 | sum_images(const char *fspec, const CMATRIX *cv)
503                          pmat = cm_alloc(myYR, myXR);
504                          memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
505                                                          /* finish header */
506 <                        fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
507 <                        fputc('\n', stdout);
508 <                        fprtresolu(myXR, myYR, stdout);
509 <                        fflush(stdout);
506 >                        fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout);
507 >                        fputc('\n', fout);
508 >                        fprtresolu(myXR, myYR, fout);
509 >                        fflush(fout);
510                  } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
511                          sprintf(errmsg, "picture '%s' format/size mismatch",
512                                          fname);
# Line 398 | Line 531 | sum_images(const char *fspec, const CMATRIX *cv)
531          free(scanline);
532                                                          /* write scanlines */
533          for (y = 0; y < myYR; y++)
534 <                if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
534 >                if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
535                          return(0);
536          cm_free(pmat);                                  /* all done */
537 <        return(fflush(stdout) == 0);
537 >        return(fflush(fout) == 0);
538   }
539  
540 < /* check to see if a string contains a %d specification */
541 < int
542 < hasDecimalSpec(const char *s)
540 > /* check to see if a string contains a %d or %o specification */
541 > static int
542 > hasNumberFormat(const char *s)
543   {
544 <        while (*s && *s != '%')
545 <                s++;
546 <        if (!*s)
547 <                return(0);
548 <        do
549 <                ++s;
550 <        while (isdigit(*s));
551 <
552 <        return(*s == 'd');
544 >        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   }
561  
562   int
563   main(int argc, char *argv[])
564   {
565 <        CMATRIX                 *tvec, *Dmat, *Tmat, *ivec, *cvec;
426 <        struct BSDF_data        *btdf;
565 >        CMATRIX                 *cvec;          /* component vector result */
566  
567          progname = argv[0];
568  
569 <        if ((argc < 4) | (argc > 5)) {
570 <                fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
571 <                                progname);
569 >        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                  return(1);
574          }
575 <                                                /* load Tregenza vector */
576 <        tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
577 <                                                /* load BTDF */
578 <        btdf = load_BSDF(argv[2]);
579 <        if (btdf == NULL)
580 <                return(1);
575 >
576 >        if (argc > 3) {                         /* VTDs expression */
577 >                CMATRIX *svec, *Dmat, *Tmat, *ivec;
578 >                COLOR   tLamb;
579 >                                                /* get sky vector */
580 >                svec = cm_load(argv[4], 0, 1, DTascii);
581 >                                                /* load BSDF */
582 >                Tmat = cm_loadBSDF(argv[2], tLamb);
583                                                  /* load Daylight matrix */
584 <        Dmat = cm_load(argv[3], btdf->ninc, tvec->nrows, DTfromHeader);
584 >                Dmat = cm_load(argv[3], Tmat==NULL ? 0 : Tmat->ncols,
585 >                                        svec->nrows, DTfromHeader);
586                                                  /* multiply vector through */
587 <        ivec = cm_multiply(Dmat, tvec);
588 <        cm_free(Dmat); cm_free(tvec);
589 <        Tmat = cm_bsdf(btdf);                   /* convert BTDF to matrix */
590 <        free_BSDF(btdf);
591 <        cvec = cm_multiply(Tmat, ivec);         /* cvec = component vector */
592 <        cm_free(Tmat); cm_free(ivec);
593 <        if (hasDecimalSpec(argv[1])) {          /* generating image */
587 >                ivec = cm_multiply(Dmat, svec);
588 >                cm_free(Dmat); cm_free(svec);
589 >                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 >        } else {                                /* else just use sky vector */
597 >                cvec = cm_load(argv[2], 0, 1, DTascii);
598 >        }
599 >
600 >        if (hasNumberFormat(argv[1])) {         /* generating image */
601                  SET_FILE_BINARY(stdout);
602                  newheader("RADIANCE", stdout);
603                  printargs(argc, argv, stdout);
604                  fputnow(stdout);
605 <                if (!sum_images(argv[1], cvec))
605 >                if (!sum_images(argv[1], cvec, stdout))
606                          return(1);
607          } else {                                /* generating vector */
608                  CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines