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.2 by greg, Fri Jun 19 06:49:42 2009 UTC vs.
Revision 2.24 by greg, Fri Jan 11 05:07:47 2013 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 160 | Line 164 | cm_load(const char *fname, int nrows, int ncols, int d
164                  int     maxrow = (nrows > 0 ? nrows : 32000);
165                  int     r, c;
166                  for (r = 0; r < maxrow; r++) {
167 <                    if (r >= cm->nrows)         /* need more space? */
167 >                    if (r >= cm->nrows)                 /* need more space? */
168                          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 <                                if ((nrows <= 0) & (r > 0) & (c == 0)) {
172 >                                if ((nrows <= 0) & (r > 0) & !c) {
173                                          cm = cm_resize(cm, maxrow=r);
174                                          break;
175                                  } else
# Line 193 | Line 197 | cm_load(const char *fname, int nrows, int ncols, int d
197                                          if (nread == cm->nrows*cm->ncols)
198                                                          /* need more space? */
199                                                  cm = cm_resize(cm, 2*cm->nrows);
200 <                                        else if (nread % cm->ncols == 0)
200 >                                        else if (nread && !(nread % cm->ncols))
201                                                          /* seem to be  done */
202                                                  cm = cm_resize(cm, nread/cm->ncols);
203                                          else            /* ended mid-row */
# 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 > /* Extract a column vector from a matrix */
256   static CMATRIX *
257 + cm_column(const CMATRIX *cm, int c)
258 + {
259 +        CMATRIX *cvr;
260 +        int     dr;
261 +
262 +        if ((c < 0) | (c >= cm->ncols))
263 +                return(NULL);
264 +        cvr = cm_alloc(cm->nrows, 1);
265 +        if (cvr == NULL)
266 +                return(NULL);
267 +        for (dr = 0; dr < cm->nrows; dr++) {
268 +                const COLORV    *sp = cm_lval(cm,dr,c);
269 +                COLORV          *dp = cv_lval(cvr,dr);
270 +                dp[0] = sp[0];
271 +                dp[1] = sp[1];
272 +                dp[2] = sp[2];
273 +        }
274 +        return(cvr);
275 + }
276 +
277 + /* Scale a matrix by a single value */
278 + static CMATRIX *
279 + cm_scale(const CMATRIX *cm1, const COLOR sca)
280 + {
281 +        CMATRIX *cmr;
282 +        int     dr, dc;
283 +
284 +        cmr = cm_alloc(cm1->nrows, cm1->ncols);
285 +        if (cmr == NULL)
286 +                return(NULL);
287 +        for (dr = 0; dr < cmr->nrows; dr++)
288 +            for (dc = 0; dc < cmr->ncols; dc++) {
289 +                const COLORV    *sp = cm_lval(cm1,dr,dc);
290 +                COLORV          *dp = cm_lval(cmr,dr,dc);
291 +                dp[0] = sp[0] * sca[0];
292 +                dp[1] = sp[1] * sca[1];
293 +                dp[2] = sp[2] * sca[2];
294 +            }
295 +        return(cmr);
296 + }
297 +
298 + /* Multiply two matrices (or a matrix and a vector) and allocate the result */
299 + static CMATRIX *
300   cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
301   {
302          CMATRIX *cmr;
# Line 289 | Line 336 | cm_print(const CMATRIX *cm, FILE *fp)
336          }
337   }
338  
339 < /* convert a BSDF to our matrix representation */
339 > /* Convert a BSDF to our matrix representation */
340   static CMATRIX *
341 < cm_bsdf(const struct BSDF_data *bsdf)
341 > cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
342   {
343          CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
344 <        COLORV  *mp = cm->cmem;
344 >        int     nbadohm = 0;
345 >        int     nneg = 0;
346          int     r, c;
347 <        
348 <        for (r = 0; r < cm->nrows; r++)
349 <                for (c = 0; c < cm->ncols; c++, mp += 3)
350 <                        mp[0] = mp[1] = mp[2] = BSDF_value(bsdf,c,r);
347 >                                        /* loop over incident angles */
348 >        for (c = 0; c < cm->ncols; c++) {
349 >                const double    dom = mBSDF_incohm(bsdf,c);
350 >                                        /* projected solid angle */
351 >                nbadohm += (dom <= 0);
352 >
353 >                for (r = 0; r < cm->nrows; r++) {
354 >                        float   f = mBSDF_value(bsdf,c,r);
355 >                        COLORV  *mp = cm_lval(cm,r,c);
356 >                                        /* check BSDF value */
357 >                        if ((f <= 0) | (dom <= 0)) {
358 >                                nneg += (f < -FTINY);
359 >                                f = .0f;
360 >                        }
361 >                        copycolor(mp, specCol);
362 >                        scalecolor(mp, f);
363 >                        addcolor(mp, bsdfLamb);
364 >                        scalecolor(mp, dom);
365 >                }
366 >        }
367 >        if (nneg | nbadohm) {
368 >                sprintf(errmsg,
369 >                    "BTDF has %d negatives and %d bad incoming solid angles",
370 >                                nneg, nbadohm);
371 >                error(WARNING, errmsg);
372 >        }
373          return(cm);
374   }
375  
376 < /* Sum together a set of images and write result to stdout */
376 > /* Convert between input and output indices for reciprocity */
377   static int
378 < sum_images(const char *fspec, const CMATRIX *cv)
378 > recip_out_from_in(const SDMat *bsdf, int in_recip)
379   {
380 +        FVECT   v;
381 +
382 +        if (!mBSDF_incvec(v, bsdf, in_recip+.5))
383 +                return(in_recip);               /* XXX should be error! */
384 +        v[2] = -v[2];
385 +        return(mBSDF_outndx(bsdf, v));
386 + }
387 +
388 + /* Convert between output and input indices for reciprocity */
389 + static int
390 + recip_in_from_out(const SDMat *bsdf, int out_recip)
391 + {
392 +        FVECT   v;
393 +
394 +        if (!mBSDF_outvec(v, bsdf, out_recip+.5))
395 +                return(out_recip);              /* XXX should be error! */
396 +        v[2] = -v[2];
397 +        return(mBSDF_incndx(bsdf, v));
398 + }
399 +
400 + /* Convert a BSDF to our matrix representation, applying reciprocity */
401 + static CMATRIX *
402 + cm_bsdf_recip(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
403 + {
404 +        CMATRIX *cm = cm_alloc(bsdf->ninc, bsdf->nout);
405 +        int     nbadohm = 0;
406 +        int     nneg = 0;
407 +        int     r, c;
408 +                                        /* loop over incident angles */
409 +        for (c = 0; c < cm->ncols; c++) {
410 +                const int       ro = recip_out_from_in(bsdf,c);
411 +                const double    dom = mBSDF_outohm(bsdf,ro);
412 +                                        /* projected solid angle */
413 +                nbadohm += (dom <= 0);
414 +
415 +                for (r = 0; r < cm->nrows; r++) {
416 +                        const int       ri = recip_in_from_out(bsdf,r);
417 +                        float           f = mBSDF_value(bsdf,ri,ro);
418 +                        COLORV          *mp = cm_lval(cm,r,c);
419 +                                        /* check BSDF value */
420 +                        if ((f <= 0) | (dom <= 0)) {
421 +                                nneg += (f < -FTINY);
422 +                                f = .0f;
423 +                        }
424 +                        copycolor(mp, specCol);
425 +                        scalecolor(mp, f);
426 +                        addcolor(mp, bsdfLamb);
427 +                        scalecolor(mp, dom);
428 +                }
429 +        }
430 +        if (nneg | nbadohm) {
431 +                sprintf(errmsg,
432 +                    "BTDF has %d negatives and %d bad incoming solid angles",
433 +                                nneg, nbadohm);
434 +                error(WARNING, errmsg);
435 +        }
436 +        return(cm);
437 + }
438 +
439 + /* Load and convert a matrix BSDF from the given XML file */
440 + static CMATRIX *
441 + cm_loadBSDF(char *fname, COLOR cLamb)
442 + {
443 +        CMATRIX         *Tmat;
444 +        char            *fpath;
445 +        int             recip;
446 +        SDError         ec;
447 +        SDData          myBSDF;
448 +        SDSpectralDF    *tdf;
449 +        COLOR           bsdfLamb, specCol;
450 +                                        /* find path to BSDF file */
451 +        fpath = getpath(fname, getrlibpath(), R_OK);
452 +        if (fpath == NULL) {
453 +                sprintf(errmsg, "cannot find BSDF file '%s'", fname);
454 +                error(USER, errmsg);
455 +        }
456 +        SDclearBSDF(&myBSDF, fname);    /* load XML and check type */
457 +        ec = SDloadFile(&myBSDF, fpath);
458 +        if (ec)
459 +                error(USER, transSDError(ec));
460 +        ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb);
461 +        recip = (myBSDF.tb == NULL);
462 +        tdf = recip ? myBSDF.tf : myBSDF.tb;
463 +        if (tdf == NULL) {              /* no non-Lambertian transmission? */
464 +                if (cLamb != NULL)
465 +                        copycolor(cLamb, bsdfLamb);
466 +                SDfreeBSDF(&myBSDF);
467 +                return(NULL);
468 +        }
469 +        if (tdf->ncomp != 1 || tdf->comp[0].func != &SDhandleMtx) {
470 +                sprintf(errmsg, "unsupported BSDF '%s'", fpath);
471 +                error(USER, errmsg);
472 +        }
473 +                                        /* convert BTDF to matrix */
474 +        ccy2rgb(&tdf->comp[0].cspec[0], 1., specCol);
475 +        Tmat = recip ? cm_bsdf_recip(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist)
476 +                        : cm_bsdf(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist);
477 +        if (cLamb != NULL)              /* Lambertian is included */
478 +                setcolor(cLamb, .0, .0, .0);
479 +                                        /* free BSDF and return */
480 +        SDfreeBSDF(&myBSDF);
481 +        return(Tmat);
482 + }
483 +
484 + /* Sum together a set of images and write result to fout */
485 + static int
486 + sum_images(const char *fspec, const CMATRIX *cv, FILE *fout)
487 + {
488          int     myDT = DTfromHeader;
489 <        CMATRIX *pmat;
490 <        COLOR   *scanline;
491 <        int     myXR, myYR;
489 >        COLOR   *scanline = NULL;
490 >        CMATRIX *pmat = NULL;
491 >        int     myXR=0, myYR=0;
492          int     i, y;
493  
494          if (cv->ncols != 1)
# Line 321 | Line 499 | sum_images(const char *fspec, const CMATRIX *cv)
499                  FILE            *fp;
500                  int             dt, xr, yr;
501                  COLORV          *psp;
502 +                                                        /* check for zero */
503 +                if ((scv[RED] == 0) & (scv[GRN] == 0) & (scv[BLU] == 0) &&
504 +                                (myDT != DTfromHeader) | (i < cv->nrows-1))
505 +                        continue;
506                                                          /* open next picture */
507                  sprintf(fname, fspec, i);
508                  if ((fp = fopen(fname, "r")) == NULL) {
# Line 343 | Line 525 | sum_images(const char *fspec, const CMATRIX *cv)
525                          pmat = cm_alloc(myYR, myXR);
526                          memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
527                                                          /* finish header */
528 <                        fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
529 <                        fputc('\n', stdout);
530 <                        fprtresolu(myXR, myYR, stdout);
531 <                        fflush(stdout);
528 >                        fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout);
529 >                        fputc('\n', fout);
530 >                        fprtresolu(myXR, myYR, fout);
531 >                        fflush(fout);
532                  } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
533                          sprintf(errmsg, "picture '%s' format/size mismatch",
534                                          fname);
# Line 371 | Line 553 | sum_images(const char *fspec, const CMATRIX *cv)
553          free(scanline);
554                                                          /* write scanlines */
555          for (y = 0; y < myYR; y++)
556 <                if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
556 >                if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
557                          return(0);
558          cm_free(pmat);                                  /* all done */
559 <        return(fflush(stdout) == 0);
559 >        return(fflush(fout) == 0);
560   }
561  
562 < /* check to see if a string contains a %d specification */
563 < int
564 < hasDecimalSpec(const char *s)
562 > /* check to see if a string contains a %d or %o specification */
563 > static int
564 > hasNumberFormat(const char *s)
565   {
566 <        while (*s && *s != '%')
385 <                s++;
386 <        if (!*s)
566 >        if (s == NULL)
567                  return(0);
388        do
389                ++s;
390        while (isdigit(*s));
568  
569 <        return(*s == 'd');
569 >        while (*s) {
570 >                while (*s != '%')
571 >                        if (!*s++)
572 >                                return(0);
573 >                if (*++s == '%') {              /* ignore "%%" */
574 >                        ++s;
575 >                        continue;
576 >                }
577 >                while (isdigit(*s))             /* field length */
578 >                        ++s;
579 >                                                /* field we'll use? */
580 >                if ((*s == 'd') | (*s == 'i') | (*s == 'o') |
581 >                                        (*s == 'x') | (*s == 'X'))
582 >                        return(1);
583 >        }
584 >        return(0);                              /* didn't find one */
585   }
586  
587   int
588   main(int argc, char *argv[])
589   {
590 <        CMATRIX                 *tvec, *Dmat, *Tmat, *ivec, *cvec;
591 <        struct BSDF_data        *btdf;
590 >        int             nsteps = 1;
591 >        char            *ofspec = NULL;
592 >        FILE            *ofp = stdout;
593 >        CMATRIX         *cmtx;          /* component vector/matrix result */
594 >        char            fnbuf[256];
595 >        int             a, i;
596  
597          progname = argv[0];
598 +                                        /* get options */
599 +        for (a = 1; a < argc-1 && argv[a][0] == '-'; a++)
600 +                switch (argv[0][1]) {
601 +                case 'n':
602 +                        nsteps = atoi(argv[++a]);
603 +                        if (nsteps <= 0)
604 +                                goto userr;
605 +                        break;
606 +                case 'o':
607 +                        ofspec = argv[++a];
608 +                        break;
609 +                default:
610 +                        goto userr;
611 +                }
612 +        if ((argc-a < 1) | (argc-a > 4))
613 +                goto userr;
614  
615 <        if ((argc < 4) | (argc > 5)) {
616 <                fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
617 <                                progname);
618 <                return(1);
615 >        if (argc-a > 2) {                       /* VTDs expression */
616 >                CMATRIX *smtx, *Dmat, *Tmat, *imtx;
617 >                COLOR   tLamb;
618 >                                                /* get sky vector/matrix */
619 >                smtx = cm_load(argv[a+3], 0, nsteps, DTascii);
620 >                                                /* load BSDF */
621 >                Tmat = cm_loadBSDF(argv[a+1], tLamb);
622 >                                                /* load Daylight matrix */
623 >                Dmat = cm_load(argv[a+2], Tmat==NULL ? 0 : Tmat->ncols,
624 >                                        smtx->nrows, DTfromHeader);
625 >                                                /* multiply vector through */
626 >                imtx = cm_multiply(Dmat, smtx);
627 >                cm_free(Dmat); cm_free(smtx);
628 >                if (Tmat == NULL) {             /* diffuse only */
629 >                        cmtx = cm_scale(imtx, tLamb);
630 >                } else {                        /* else apply BTDF matrix */
631 >                        cmtx = cm_multiply(Tmat, imtx);
632 >                        cm_free(Tmat);
633 >                }
634 >                cm_free(imtx);
635 >        } else {                                /* sky vector/matrix only */
636 >                cmtx = cm_load(argv[a+1], 0, nsteps, DTascii);
637          }
638 <        tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
639 <        Dmat = cm_load(argv[3], 0, tvec->nrows, DTfromHeader);
640 <        btdf = load_BSDF(argv[2]);
641 <        if (btdf == NULL)
412 <                return(1);
413 <        if (btdf->ninc != Dmat->nrows) {
414 <                sprintf(errmsg, "Incoming BTDF dir (%d) mismatch to D (%d)",
415 <                                btdf->ninc, Dmat->nrows);
416 <                error(USER, errmsg);
638 >                                                /* prepare output stream */
639 >        if ((ofspec != NULL) & (nsteps == 1) && hasNumberFormat(ofspec)) {
640 >                sprintf(fnbuf, ofspec, 1);
641 >                ofspec = fnbuf;
642          }
643 <                                                /* multiply vector through */
644 <        ivec = cm_multiply(Dmat, tvec);
645 <        cm_free(Dmat); cm_free(tvec);
646 <        Tmat = cm_bsdf(btdf);                   /* convert BTDF to matrix */
422 <        free_BSDF(btdf);
423 <        cvec = cm_multiply(Tmat, ivec);         /* cvec = component vector */
424 <        cm_free(Tmat); cm_free(ivec);
425 <        if (hasDecimalSpec(argv[1])) {          /* generating image */
426 <                SET_FILE_BINARY(stdout);
427 <                newheader("RADIANCE", stdout);
428 <                printargs(argc, argv, stdout);
429 <                fputnow(stdout);
430 <                if (!sum_images(argv[1], cvec))
643 >        if (ofspec != NULL && !hasNumberFormat(ofspec)) {
644 >                if ((ofp = fopen(ofspec, "w")) == NULL) {
645 >                        fprintf(stderr, "%s: cannot open '%s' for output\n",
646 >                                        progname, ofspec);
647                          return(1);
648 <        } else {                                /* generating vector */
649 <                CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
650 <                CMATRIX *rvec = cm_multiply(Vmat, cvec);
648 >                }
649 >                ofspec = NULL;                  /* only need to open once */
650 >        }
651 >        if (hasNumberFormat(argv[a])) {         /* generating image(s) */
652 >                if (ofspec == NULL) {
653 >                        SET_FILE_BINARY(ofp);
654 >                        newheader("RADIANCE", ofp);
655 >                        printargs(argc, argv, ofp);
656 >                        fputnow(ofp);
657 >                }
658 >                if (nsteps > 1)                 /* multiple output frames? */
659 >                        for (i = 0; i < nsteps; i++) {
660 >                                CMATRIX *cvec = cm_column(cmtx, i);
661 >                                if (ofspec != NULL) {
662 >                                        sprintf(fnbuf, ofspec, i+1);
663 >                                        if ((ofp = fopen(fnbuf, "wb")) == NULL) {
664 >                                                fprintf(stderr,
665 >                                        "%s: cannot open '%s' for output\n",
666 >                                                        progname, fnbuf);
667 >                                                return(1);
668 >                                        }
669 >                                        newheader("RADIANCE", ofp);
670 >                                        printargs(argc, argv, ofp);
671 >                                        fputnow(ofp);
672 >                                }
673 >                                fprintf(ofp, "FRAME=%d\n", i+1);
674 >                                if (!sum_images(argv[a], cvec, ofp))
675 >                                        return(1);
676 >                                if (ofspec != NULL) {
677 >                                        if (fclose(ofp) == EOF) {
678 >                                                fprintf(stderr,
679 >                                                "%s: error writing to '%s'\n",
680 >                                                        progname, fnbuf);
681 >                                                return(1);
682 >                                        }
683 >                                        ofp = stdout;
684 >                                }
685 >                                cm_free(cvec);
686 >                        }
687 >                else if (!sum_images(argv[a], cmtx, ofp))
688 >                        return(1);
689 >        } else {                                /* generating vector/matrix */
690 >                CMATRIX *Vmat = cm_load(argv[a], 0, cmtx->nrows, DTfromHeader);
691 >                CMATRIX *rmtx = cm_multiply(Vmat, cmtx);
692                  cm_free(Vmat);
693 <                cm_print(rvec, stdout);
694 <                cm_free(rvec);
693 >                if (ofspec != NULL)             /* multiple vector files? */
694 >                        for (i = 0; i < nsteps; i++) {
695 >                                CMATRIX *rvec = cm_column(rmtx, i);
696 >                                sprintf(fnbuf, ofspec, i+1);
697 >                                if ((ofp = fopen(fnbuf, "w")) == NULL) {
698 >                                        fprintf(stderr,
699 >                                        "%s: cannot open '%s' for output\n",
700 >                                                        progname, fnbuf);
701 >                                        return(1);
702 >                                }
703 >                                cm_print(rvec, ofp);
704 >                                if (fclose(ofp) == EOF) {
705 >                                        fprintf(stderr,
706 >                                                "%s: error writing to '%s'\n",
707 >                                                        progname, fnbuf);
708 >                                        return(1);
709 >                                }
710 >                                ofp = stdout;
711 >                                cm_free(rvec);
712 >                        }
713 >                else
714 >                        cm_print(rmtx, ofp);
715 >                cm_free(rmtx);
716          }
717 <        cm_free(cvec);                          /* final clean-up */
717 >        if (fflush(ofp) == EOF) {               /* final clean-up */
718 >                fprintf(stderr, "%s: write error on output\n", progname);
719 >                return(1);
720 >        }
721 >        cm_free(cmtx);
722          return(0);
723 + userr:
724 +        fprintf(stderr, "Usage: %s [-n nsteps][-o ospec] DCspec [skyf]\n",
725 +                                progname);
726 +        fprintf(stderr, "   or: %s [-n nsteps][-o ospec] Vspec Tbsdf.xml Dmat.dat [skyf]\n",
727 +                                progname);
728 +        return(1);
729   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines