src/cv/bsdf2rad.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf2rad.c,v 2.11 2014/03/30 17:28:04 greg Exp $";
#endif
/*
 *  Plot 3-D BSDF output based on scattering interpolant or XML representation
 */

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "rtprocess.h"
#include "bsdfrep.h"

const float     colarr[6][3] = {
                .7, 1., .7,
                1., .7, .7,
                .7, .7, 1.,
                1., .5, 1.,
                1., 1., .5,
                .5, 1., 1.
        };

#ifdef _WIN32
char    validf[] = "-e \"valid(s,t)=X`SYS(s,t)^2+Y`SYS(s,t)^2+Z`SYS(s,t)^2-1e-7\"";
#else
char    validf[] = "-e 'valid(s,t)=X`SYS(s,t)^2+Y`SYS(s,t)^2+Z`SYS(s,t)^2-1e-7'";
#endif

char    *progname;

/* Produce a Radiance model plotting the indicated incident direction(s) */
int
main(int argc, char *argv[])
{
        int     showPeaks = 0;
        int     doTrans = 0;
        int     inpXML = -1;
        RBFNODE *rbf = NULL;
        FILE    *fp;
        char    buf[128];
        SDData  myBSDF;
        double  bsdf, min_log;
        FVECT   idir, odir;
        int     i, j, n;
                                                /* check arguments */
        progname = argv[0];
        if (argc > 1 && !strcmp(argv[1], "-p")) {
                ++showPeaks;
                ++argv; --argc;
        }
        if (argc > 1 && !strcmp(argv[1], "-t")) {
                ++doTrans;
                ++argv; --argc;
        }
        if (argc >= 4 && (n = strlen(argv[1])-4) > 0) {
                if (!strcasecmp(argv[1]+n, ".xml"))
                        inpXML = 1;
                else if (!strcasecmp(argv[1]+n, ".sir"))
                        inpXML = 0;
        }
        if (inpXML < 0) {
                fprintf(stderr, "Usage: %s [-p] bsdf.sir theta1 phi1 .. > output.rad\n", progname);
                fprintf(stderr, "   Or: %s [-t] bsdf.xml theta1 phi1 .. > output.rad\n", progname);
                return(1);
        }
                                                /* load input */
        if (inpXML) {
                SDclearBSDF(&myBSDF, argv[1]);
                if (SDreportError(SDloadFile(&myBSDF, argv[1]), stderr))
                        return(1);
                bsdf_min = 1./M_PI;
                if (myBSDF.rf != NULL && myBSDF.rLambFront.cieY < bsdf_min*M_PI)
                        bsdf_min = myBSDF.rLambFront.cieY/M_PI;
                if (myBSDF.rb != NULL && myBSDF.rLambBack.cieY < bsdf_min*M_PI)
                        bsdf_min = myBSDF.rLambBack.cieY/M_PI;
                if ((myBSDF.tf != NULL) | (myBSDF.tb != NULL) &&
                                myBSDF.tLamb.cieY < bsdf_min*M_PI)
                        bsdf_min = myBSDF.tLamb.cieY/M_PI;
                if (doTrans && (myBSDF.tf == NULL) & (myBSDF.tb == NULL)) {
                        fprintf(stderr, "%s: no transmitted component in '%s'\n",
                                        progname, argv[1]);
                        return(1);
                }
        } else {
                fp = fopen(argv[1], "rb");
                if (fp == NULL) {
                        fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
                                        progname, argv[1]);
                        return(1);
                }
                if (!load_bsdf_rep(fp))
                        return(1);
                fclose(fp);
        }
#ifdef DEBUG
        fprintf(stderr, "Minimum BSDF set to %.4f\n", bsdf_min);
#endif
        min_log = log(bsdf_min*.5 + 1e-5);
                                                /* output BSDF rep. */
        for (n = 0; (n < 6) & (2*n+3 < argc); n++) {
                double  theta = (M_PI/180.)*atof(argv[2*n+2]);
                double  phi = (M_PI/180.)*atof(argv[2*n+3]);
                if (theta < -FTINY) {
                        fprintf(stderr, "%s: theta values must be positive\n",
                                        progname);
                        return(1);
                }
                if (inpXML) {
                        input_orient = (theta <= M_PI/2.) ? 1 : -1;
                        output_orient = doTrans ? -input_orient : input_orient;
                }
                idir[2] = sin(theta);
                idir[0] = idir[2] * cos(phi);
                idir[1] = idir[2] * sin(phi);
                idir[2] = input_orient * sqrt(1. - idir[2]*idir[2]);
#ifdef DEBUG
                fprintf(stderr, "Computing BSDF for incident direction (%.1f,%.1f)\n",
                                get_theta180(idir), get_phi360(idir));
#endif
                if (!inpXML)
                        rbf = advect_rbf(idir, 15000);
#ifdef DEBUG
                if (inpXML)
                        fprintf(stderr, "Hemispherical %s: %.3f\n",
                                (output_orient > 0 ? "reflection" : "transmission"),
                                SDdirectHemi(idir, SDsampSp|SDsampDf |
                                                (output_orient > 0 ?
                                                 SDsampR : SDsampT), &myBSDF));
                else if (rbf == NULL)
                        fputs("Empty RBF\n", stderr);
                else
                        fprintf(stderr, "Hemispherical %s: %.3f\n",
                                (output_orient > 0 ? "reflection" : "transmission"),
                                rbf->vtotal);
#endif
                printf("void trans tmat\n0\n0\n7 %f %f %f .04 .04 .9 1\n",
                                colarr[n][0], colarr[n][1], colarr[n][2]);
                if (showPeaks && rbf != NULL) {
                        printf("void plastic pmat\n0\n0\n5 %f %f %f .04 .08\n",
                                1.-colarr[n][0], 1.-colarr[n][1], 1.-colarr[n][2]);
                        for (i = 0; i < rbf->nrbf; i++) {
                                ovec_from_pos(odir, rbf->rbfa[i].gx, rbf->rbfa[i].gy);
                                bsdf = eval_rbfrep(rbf, odir);
                                bsdf = log(bsdf + 1e-5) - min_log;
                                printf("pmat sphere p%d\n0\n0\n4 %f %f %f %f\n",
                                        i+1, odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf,
                                                .007*bsdf);
                        }
                }
                fflush(stdout);
                sprintf(buf, "gensurf tmat bsdf%d - - - %d %d %s", n+1,
                                                GRIDRES-1, GRIDRES-1, validf);
                fp = popen(buf, "w");
                if (fp == NULL) {
                        fprintf(stderr, "%s: cannot open '| %s'\n", progname, buf);
                        return(1);
                }
                for (i = 0; i < GRIDRES; i++)
                    for (j = 0; j < GRIDRES; j++) {
                        ovec_from_pos(odir, i, j);
                        if (inpXML) {
                                SDValue sval;
                                if (SDreportError(SDevalBSDF(&sval, odir,
                                                        idir, &myBSDF), stderr))
                                        return(1);
                                bsdf = sval.cieY;
                        } else
                                bsdf = eval_rbfrep(rbf, odir);
                        bsdf = log(bsdf + 1e-5) - min_log;
                        fprintf(fp, "%.8e %.8e %.8e\n",
                                        odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf);
                    }
                if (rbf != NULL)
                        free(rbf);
                if (pclose(fp))
                        return(1);
        }
        return(0);
}
Revision:	2.12
Committed:	Thu Aug 21 10:33:48 2014 UTC (9 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.11:	+3 -4 lines
Log Message:	Added grazing angle extrapolation to BSDF interpolation
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.12	static const char RCSid[] = "$Id: bsdf2rad.c,v 2.11 2014/03/30 17:28:04 greg Exp $";
3	greg	2.1	#endif
4			/*
5	greg	2.4	* Plot 3-D BSDF output based on scattering interpolant or XML representation
6	greg	2.1	*/
7
8			#define _USE_MATH_DEFINES
9			#include <stdio.h>
10	greg	2.4	#include <string.h>
11	greg	2.1	#include <stdlib.h>
12			#include <math.h>
13	greg	2.5	#include "rtprocess.h"
14	greg	2.1	#include "bsdfrep.h"
15
16			const float colarr[6][3] = {
17			.7, 1., .7,
18			1., .7, .7,
19			.7, .7, 1.,
20			1., .5, 1.,
21			1., 1., .5,
22			.5, 1., 1.
23			};
24
25	greg	2.11	#ifdef _WIN32
26	greg	2.10	char validf[] = "-e \"valid(s,t)=X`SYS(s,t)^2+Y`SYS(s,t)^2+Z`SYS(s,t)^2-1e-7\"";
27	greg	2.11	#else
28			char validf[] = "-e 'valid(s,t)=X`SYS(s,t)^2+Y`SYS(s,t)^2+Z`SYS(s,t)^2-1e-7'";
29			#endif
30	greg	2.10
31	greg	2.1	char *progname;
32
33			/* Produce a Radiance model plotting the indicated incident direction(s) */
34			int
35			main(int argc, char *argv[])
36			{
37	greg	2.2	int showPeaks = 0;
38	greg	2.4	int doTrans = 0;
39			int inpXML = -1;
40			RBFNODE *rbf = NULL;
41			FILE *fp;
42	greg	2.1	char buf[128];
43	greg	2.4	SDData myBSDF;
44	greg	2.1	double bsdf, min_log;
45	greg	2.4	FVECT idir, odir;
46	greg	2.1	int i, j, n;
47	greg	2.4	/* check arguments */
48	greg	2.1	progname = argv[0];
49	greg	2.2	if (argc > 1 && !strcmp(argv[1], "-p")) {
50			++showPeaks;
51			++argv; --argc;
52			}
53	greg	2.4	if (argc > 1 && !strcmp(argv[1], "-t")) {
54			++doTrans;
55			++argv; --argc;
56			}
57			if (argc >= 4 && (n = strlen(argv[1])-4) > 0) {
58			if (!strcasecmp(argv[1]+n, ".xml"))
59			inpXML = 1;
60			else if (!strcasecmp(argv[1]+n, ".sir"))
61			inpXML = 0;
62			}
63			if (inpXML < 0) {
64	greg	2.2	fprintf(stderr, "Usage: %s [-p] bsdf.sir theta1 phi1 .. > output.rad\n", progname);
65	greg	2.4	fprintf(stderr, " Or: %s [-t] bsdf.xml theta1 phi1 .. > output.rad\n", progname);
66	greg	2.1	return(1);
67			}
68			/* load input */
69	greg	2.4	if (inpXML) {
70			SDclearBSDF(&myBSDF, argv[1]);
71			if (SDreportError(SDloadFile(&myBSDF, argv[1]), stderr))
72			return(1);
73			bsdf_min = 1./M_PI;
74			if (myBSDF.rf != NULL && myBSDF.rLambFront.cieY < bsdf_min*M_PI)
75			bsdf_min = myBSDF.rLambFront.cieY/M_PI;
76			if (myBSDF.rb != NULL && myBSDF.rLambBack.cieY < bsdf_min*M_PI)
77			bsdf_min = myBSDF.rLambBack.cieY/M_PI;
78			if ((myBSDF.tf != NULL) \| (myBSDF.tb != NULL) &&
79			myBSDF.tLamb.cieY < bsdf_min*M_PI)
80			bsdf_min = myBSDF.tLamb.cieY/M_PI;
81			if (doTrans && (myBSDF.tf == NULL) & (myBSDF.tb == NULL)) {
82			fprintf(stderr, "%s: no transmitted component in '%s'\n",
83			progname, argv[1]);
84			return(1);
85			}
86			} else {
87			fp = fopen(argv[1], "rb");
88			if (fp == NULL) {
89			fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
90			progname, argv[1]);
91			return(1);
92			}
93			if (!load_bsdf_rep(fp))
94			return(1);
95			fclose(fp);
96	greg	2.1	}
97	greg	2.4	#ifdef DEBUG
98			fprintf(stderr, "Minimum BSDF set to %.4f\n", bsdf_min);
99			#endif
100	greg	2.6	min_log = log(bsdf_min*.5 + 1e-5);
101	greg	2.2	/* output BSDF rep. */
102	greg	2.1	for (n = 0; (n < 6) & (2*n+3 < argc); n++) {
103	greg	2.9	double theta = (M_PI/180.)atof(argv[2n+2]);
104			double phi = (M_PI/180.)atof(argv[2n+3]);
105			if (theta < -FTINY) {
106			fprintf(stderr, "%s: theta values must be positive\n",
107			progname);
108			return(1);
109			}
110	greg	2.4	if (inpXML) {
111	greg	2.9	input_orient = (theta <= M_PI/2.) ? 1 : -1;
112	greg	2.4	output_orient = doTrans ? -input_orient : input_orient;
113			}
114	greg	2.9	idir[2] = sin(theta);
115			idir[0] = idir[2] * cos(phi);
116			idir[1] = idir[2] * sin(phi);
117	greg	2.4	idir[2] = input_orient * sqrt(1. - idir[2]*idir[2]);
118	greg	2.3	#ifdef DEBUG
119	greg	2.4	fprintf(stderr, "Computing BSDF for incident direction (%.1f,%.1f)\n",
120			get_theta180(idir), get_phi360(idir));
121	greg	2.3	#endif
122	greg	2.4	if (!inpXML)
123			rbf = advect_rbf(idir, 15000);
124	greg	2.3	#ifdef DEBUG
125	greg	2.4	if (inpXML)
126			fprintf(stderr, "Hemispherical %s: %.3f\n",
127			(output_orient > 0 ? "reflection" : "transmission"),
128			SDdirectHemi(idir, SDsampSp\|SDsampDf \|
129			(output_orient > 0 ?
130			SDsampR : SDsampT), &myBSDF));
131			else if (rbf == NULL)
132			fputs("Empty RBF\n", stderr);
133	greg	2.1	else
134	greg	2.4	fprintf(stderr, "Hemispherical %s: %.3f\n",
135			(output_orient > 0 ? "reflection" : "transmission"),
136			rbf->vtotal);
137	greg	2.3	#endif
138	greg	2.2	printf("void trans tmat\n0\n0\n7 %f %f %f .04 .04 .9 1\n",
139			colarr[n][0], colarr[n][1], colarr[n][2]);
140			if (showPeaks && rbf != NULL) {
141			printf("void plastic pmat\n0\n0\n5 %f %f %f .04 .08\n",
142			1.-colarr[n][0], 1.-colarr[n][1], 1.-colarr[n][2]);
143			for (i = 0; i < rbf->nrbf; i++) {
144	greg	2.4	ovec_from_pos(odir, rbf->rbfa[i].gx, rbf->rbfa[i].gy);
145	greg	2.12	bsdf = eval_rbfrep(rbf, odir);
146	greg	2.7	bsdf = log(bsdf + 1e-5) - min_log;
147	greg	2.2	printf("pmat sphere p%d\n0\n0\n4 %f %f %f %f\n",
148	greg	2.4	i+1, odir[0]bsdf, odir[1]bsdf, odir[2]*bsdf,
149	greg	2.2	.007*bsdf);
150			}
151			}
152			fflush(stdout);
153	greg	2.10	sprintf(buf, "gensurf tmat bsdf%d - - - %d %d %s", n+1,
154			GRIDRES-1, GRIDRES-1, validf);
155	greg	2.2	fp = popen(buf, "w");
156			if (fp == NULL) {
157			fprintf(stderr, "%s: cannot open '\| %s'\n", progname, buf);
158			return(1);
159			}
160	greg	2.1	for (i = 0; i < GRIDRES; i++)
161			for (j = 0; j < GRIDRES; j++) {
162	greg	2.4	ovec_from_pos(odir, i, j);
163			if (inpXML) {
164			SDValue sval;
165			if (SDreportError(SDevalBSDF(&sval, odir,
166			idir, &myBSDF), stderr))
167			return(1);
168			bsdf = sval.cieY;
169			} else
170	greg	2.12	bsdf = eval_rbfrep(rbf, odir);
171	greg	2.7	bsdf = log(bsdf + 1e-5) - min_log;
172	greg	2.1	fprintf(fp, "%.8e %.8e %.8e\n",
173	greg	2.4	odir[0]bsdf, odir[1]bsdf, odir[2]*bsdf);
174	greg	2.1	}
175			if (rbf != NULL)
176			free(rbf);
177			if (pclose(fp))
178			return(1);
179			}
180			return(0);
181			}