src/common/bsdf_t.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf_t.c,v 3.4 2011/02/19 01:48:59 greg Exp $";
#endif
/*
 *  bsdf_t.c
 *  
 *  Definitions for variable-resolution BSDF trees
 *
 *  Created by Greg Ward on 2/2/11.
 *
 */

#include "rtio.h"
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "ezxml.h"
#include "bsdf.h"
#include "bsdf_t.h"

/* Allocate a new scattering distribution node */
static SDNode *
SDnewNode(int nd, int lg)
{
        SDNode  *st;
        
        if (nd <= 0) {
                strcpy(SDerrorDetail, "Zero dimension BSDF node request");
                return NULL;
        }
        if (nd > SD_MAXDIM) {
                sprintf(SDerrorDetail, "Illegal BSDF dimension (%d > %d)",
                                nd, SD_MAXDIM);
                return NULL;
        }
        if (lg < 0) {
                st = (SDNode *)malloc(sizeof(SDNode) +
                                ((1<<nd) - 1)*sizeof(st->u.t[0]));
                if (st != NULL)
                        memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
        } else
                st = (SDNode *)malloc(sizeof(SDNode) +
                                ((1 << nd*lg) - 1)*sizeof(st->u.v[0]));
                
        if (st == NULL) {
                if (lg < 0)
                        sprintf(SDerrorDetail,
                                "Cannot allocate %d branch BSDF tree", nd);
                else
                        sprintf(SDerrorDetail,
                                "Cannot allocate %d BSDF leaves", 1 << nd*lg);
                return NULL;
        }
        st->ndim = nd;
        st->log2GR = lg;
        return st;
}

/* Free an SD tree */
static void
SDfreeTre(void *p)
{
        SDNode  *st = (SDNode *)p;
        int     i;

        if (st == NULL)
                return;
        for (i = (st->log2GR < 0) << st->ndim; i--; )
                SDfreeTre(st->u.t[i]);
        free((void *)st);
}

#if 0   /* Unused and untested routines */

/* Add up N-dimensional hypercube array values over the given box */
static double
SDiterSum(const float *va, int nd, int siz, const int *imin, const int *imax)
{
        double          sum = .0;
        unsigned        skipsiz = 1;
        int             i;
        
        for (i = nd; --i > 0; )
                skipsiz *= siz;
        if (skipsiz == 1)
                for (i = *imin; i < *imax; i++)
                        sum += va[i];
        else
                for (i = *imin; i < *imax; i++)
                        sum += SDiterSum(va + i*skipsiz,
                                        nd-1, siz, imin+1, imax+1);
        return sum;
}

/* Average BSDF leaves over an orthotope defined by the unit hypercube */
static double
SDavgBox(const SDNode *st, const double *bmin, const double *bmax)
{
        int             imin[SD_MAXDIM], imax[SD_MAXDIM];
        unsigned        n;
        int             i;

        if (!st)
                return .0;
                                        /* check box limits */
        for (i = st->ndim; i--; ) {
                if (bmin[i] >= 1.)
                        return .0;
                if (bmax[i] <= .0)
                        return .0;
                if (bmin[i] >= bmax[i])
                        return .0;
        }
        if (st->log2GR < 0) {           /* iterate on subtree */
                double          sum = .0, wsum = 1e-20;
                double          sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;

                for (n = 1 << st->ndim; n--; ) {
                        w = 1.;
                        for (i = st->ndim; i--; ) {
                                sbmin[i] = 2.*bmin[i];
                                sbmax[i] = 2.*bmax[i];
                                if (n & 1<<i) {
                                        sbmin[i] -= 1.;
                                        sbmax[i] -= 1.;
                                }
                                if (sbmin[i] < .0) sbmin[i] = .0;
                                if (sbmax[i] > 1.) sbmax[i] = 1.;
                                w *= sbmax[i] - sbmin[i];
                        }
                        if (w > 1e-10) {
                                sum += w * SDavgBox(st->u.t[n], sbmin, sbmax);
                                wsum += w;
                        }
                }
                return sum / wsum;
        }
        n = 1;                          /* iterate over leaves */
        for (i = st->ndim; i--; ) {
                imin[i] = (bmin[i] <= .0) ? 0 
                                : (int)((1 << st->log2GR)*bmin[i]);
                imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
                                : (int)((1 << st->log2GR)*bmax[i] + .999999);
                n *= imax[i] - imin[i];
        }
        if (!n)
                return .0;
        
        return SDiterSum(st->u.v, st->ndim, 1 << st->log2GR, imin, imax) /
                        (double)n;
}

#endif  /* 0 */

/* Look up tree value at the given grid position */
static float
SDlookupTre(const SDNode *st, const double *pos)
{
        double  spos[SD_MAXDIM];
        int     i, n, t;
                                        /* climb the tree */
        while (st->log2GR < 0) {
                n = 0;                  /* move to appropriate branch */
                for (i = st->ndim; i--; ) {
                        spos[i] = 2.*pos[i];
                        t = (spos[i] >= 1.);
                        n |= t<<i;
                        spos[i] -= (double)t;
                }
                st = st->u.t[n];        /* avoids tail recursion */
                pos = spos;
        }
        n = t = 0;                      /* find grid array index */
        for (i = st->ndim; i--; ) {
                n += (int)((1<<st->log2GR)*pos[i]) << t;
                t += st->log2GR;
        }
        return st->u.v[n];              /* XXX no interpolation */
}

/* Compute non-diffuse component for variable-resolution BSDF */
static int
SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
                                const FVECT inVec, const void *dist)
{
        const SDNode    *st = (const SDNode *)dist;
        double          gridPos[4];
                                        /* convert vector coordinates */
        if (st->ndim == 3) {            /* reduce for isotropic BSDF? */
                static const FVECT      zvec = {.0, .0, 1.};
                FVECT                   rOutVec;
                spinvector(rOutVec, outVec, zvec, -atan2(inVec[1],inVec[0]));
                SDdisk2square(gridPos, rOutVec[0], rOutVec[1]);
                gridPos[2] = .5 - .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
        } else if (st->ndim == 4) {     /* XXX no component identity checks */
                SDdisk2square(gridPos, outVec[0], outVec[1]);
                SDdisk2square(gridPos+2, -inVec[0], -inVec[1]);
        } else
                return 0;               /* should be internal error */
                                        /* get nearest BSDF value */
        coef[0] = SDlookupTre(st, gridPos);
        return 1;                       /* monochromatic for now */
}

/* Load a variable-resolution BSDF tree from an open XML file */
SDError
SDloadTre(SDData *sd, ezxml_t wtl)
{
        return SDEsupport;
}

/* Variable resolution BSDF methods */
SDFunc SDhandleTre = {
        &SDgetTreBSDF,
        NULL,
        NULL,
        NULL,
        &SDfreeTre,
};
Revision:	3.5
Committed:	Tue Apr 19 21:31:22 2011 UTC (13 years ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.4:	+61 -7 lines
Log Message:	Fixed interface for determining BSDF solid angles
#	User	Rev	Content
1	greg	3.2	#ifndef lint
2	greg	3.5	static const char RCSid[] = "$Id: bsdf_t.c,v 3.4 2011/02/19 01:48:59 greg Exp $";
3	greg	3.2	#endif
4	greg	3.1	/*
5			* bsdf_t.c
6			*
7			* Definitions for variable-resolution BSDF trees
8			*
9			* Created by Greg Ward on 2/2/11.
10			*
11			*/
12
13	greg	3.3	#include "rtio.h"
14	greg	3.1	#include <stdlib.h>
15	greg	3.3	#include <math.h>
16			#include <ctype.h>
17	greg	3.1	#include "ezxml.h"
18			#include "bsdf.h"
19			#include "bsdf_t.h"
20
21			/* Allocate a new scattering distribution node */
22			static SDNode *
23			SDnewNode(int nd, int lg)
24			{
25			SDNode *st;
26
27			if (nd <= 0) {
28			strcpy(SDerrorDetail, "Zero dimension BSDF node request");
29			return NULL;
30			}
31			if (nd > SD_MAXDIM) {
32			sprintf(SDerrorDetail, "Illegal BSDF dimension (%d > %d)",
33			nd, SD_MAXDIM);
34			return NULL;
35			}
36			if (lg < 0) {
37			st = (SDNode *)malloc(sizeof(SDNode) +
38			((1<<nd) - 1)*sizeof(st->u.t[0]));
39			if (st != NULL)
40	greg	3.5	memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
41	greg	3.1	} else
42			st = (SDNode *)malloc(sizeof(SDNode) +
43			((1 << ndlg) - 1)sizeof(st->u.v[0]));
44
45			if (st == NULL) {
46			if (lg < 0)
47			sprintf(SDerrorDetail,
48			"Cannot allocate %d branch BSDF tree", nd);
49			else
50			sprintf(SDerrorDetail,
51			"Cannot allocate %d BSDF leaves", 1 << nd*lg);
52			return NULL;
53			}
54			st->ndim = nd;
55			st->log2GR = lg;
56			return st;
57			}
58
59			/* Free an SD tree */
60			static void
61	greg	3.5	SDfreeTre(void *p)
62	greg	3.1	{
63			SDNode st = (SDNode )p;
64			int i;
65
66			if (st == NULL)
67			return;
68			for (i = (st->log2GR < 0) << st->ndim; i--; )
69	greg	3.5	SDfreeTre(st->u.t[i]);
70	greg	3.1	free((void *)st);
71			}
72
73	greg	3.5	#if 0 /* Unused and untested routines */
74
75	greg	3.1	/* Add up N-dimensional hypercube array values over the given box */
76			static double
77			SDiterSum(const float va, int nd, int siz, const int imin, const int *imax)
78			{
79			double sum = .0;
80			unsigned skipsiz = 1;
81			int i;
82
83			for (i = nd; --i > 0; )
84			skipsiz *= siz;
85			if (skipsiz == 1)
86			for (i = imin; i < imax; i++)
87			sum += va[i];
88			else
89			for (i = imin; i < imax; i++)
90			sum += SDiterSum(va + i*skipsiz,
91			nd-1, siz, imin+1, imax+1);
92			return sum;
93			}
94
95			/* Average BSDF leaves over an orthotope defined by the unit hypercube */
96			static double
97			SDavgBox(const SDNode st, const double bmin, const double *bmax)
98			{
99			int imin[SD_MAXDIM], imax[SD_MAXDIM];
100			unsigned n;
101			int i;
102
103			if (!st)
104			return .0;
105			/* check box limits */
106			for (i = st->ndim; i--; ) {
107			if (bmin[i] >= 1.)
108			return .0;
109			if (bmax[i] <= .0)
110			return .0;
111			if (bmin[i] >= bmax[i])
112			return .0;
113			}
114			if (st->log2GR < 0) { /* iterate on subtree */
115			double sum = .0, wsum = 1e-20;
116			double sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;
117
118			for (n = 1 << st->ndim; n--; ) {
119			w = 1.;
120			for (i = st->ndim; i--; ) {
121			sbmin[i] = 2.*bmin[i];
122			sbmax[i] = 2.*bmax[i];
123			if (n & 1<<i) {
124			sbmin[i] -= 1.;
125			sbmax[i] -= 1.;
126			}
127			if (sbmin[i] < .0) sbmin[i] = .0;
128			if (sbmax[i] > 1.) sbmax[i] = 1.;
129			w *= sbmax[i] - sbmin[i];
130			}
131			if (w > 1e-10) {
132			sum += w * SDavgBox(st->u.t[n], sbmin, sbmax);
133			wsum += w;
134			}
135			}
136			return sum / wsum;
137			}
138			n = 1; /* iterate over leaves */
139			for (i = st->ndim; i--; ) {
140			imin[i] = (bmin[i] <= .0) ? 0
141			: (int)((1 << st->log2GR)*bmin[i]);
142			imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
143			: (int)((1 << st->log2GR)*bmax[i] + .999999);
144			n *= imax[i] - imin[i];
145			}
146			if (!n)
147			return .0;
148
149			return SDiterSum(st->u.v, st->ndim, 1 << st->log2GR, imin, imax) /
150			(double)n;
151			}
152
153	greg	3.5	#endif /* 0 */
154
155			/* Look up tree value at the given grid position */
156			static float
157			SDlookupTre(const SDNode st, const double pos)
158			{
159			double spos[SD_MAXDIM];
160			int i, n, t;
161			/* climb the tree */
162			while (st->log2GR < 0) {
163			n = 0; /* move to appropriate branch */
164			for (i = st->ndim; i--; ) {
165			spos[i] = 2.*pos[i];
166			t = (spos[i] >= 1.);
167			n \|= t<<i;
168			spos[i] -= (double)t;
169			}
170			st = st->u.t[n]; /* avoids tail recursion */
171			pos = spos;
172			}
173			n = t = 0; /* find grid array index */
174			for (i = st->ndim; i--; ) {
175			n += (int)((1<<st->log2GR)*pos[i]) << t;
176			t += st->log2GR;
177			}
178			return st->u.v[n]; /* XXX no interpolation */
179			}
180
181			/* Compute non-diffuse component for variable-resolution BSDF */
182			static int
183			SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
184			const FVECT inVec, const void *dist)
185			{
186			const SDNode st = (const SDNode )dist;
187			double gridPos[4];
188			/* convert vector coordinates */
189			if (st->ndim == 3) { /* reduce for isotropic BSDF? */
190			static const FVECT zvec = {.0, .0, 1.};
191			FVECT rOutVec;
192			spinvector(rOutVec, outVec, zvec, -atan2(inVec[1],inVec[0]));
193			SDdisk2square(gridPos, rOutVec[0], rOutVec[1]);
194			gridPos[2] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
195			} else if (st->ndim == 4) { /* XXX no component identity checks */
196			SDdisk2square(gridPos, outVec[0], outVec[1]);
197			SDdisk2square(gridPos+2, -inVec[0], -inVec[1]);
198			} else
199			return 0; /* should be internal error */
200			/* get nearest BSDF value */
201			coef[0] = SDlookupTre(st, gridPos);
202			return 1; /* monochromatic for now */
203			}
204
205	greg	3.1	/* Load a variable-resolution BSDF tree from an open XML file */
206			SDError
207	greg	3.4	SDloadTre(SDData *sd, ezxml_t wtl)
208	greg	3.1	{
209			return SDEsupport;
210			}
211
212			/* Variable resolution BSDF methods */
213	greg	3.5	SDFunc SDhandleTre = {
214			&SDgetTreBSDF,
215	greg	3.1	NULL,
216			NULL,
217			NULL,
218	greg	3.5	&SDfreeTre,
219	greg	3.1	};