src/cv/bsdfrbf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrbf.c,v 2.12 2013/10/21 18:33:15 greg Exp $";
#endif
/*
 * Radial basis function representation for BSDF data.
 *
 *      G. Ward
 */

/****************************************************************
1) Collect samples into a grid using the Shirley-Chiu
        angular mapping from a hemisphere to a square.

2) Compute an adaptive quadtree by subdividing the grid so that
        each leaf node has at least one sample up to as many
        samples as fit nicely on a plane to within a certain
        MSE tolerance.

3) Place one Gaussian lobe at each leaf node in the quadtree,
        sizing it to have a radius equal to the leaf size and
        a volume equal to the energy in that node.
*****************************************************************/

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

#ifndef RSCA
#define RSCA            2.2             /* radius scaling factor (empirical) */
#endif
#ifndef SMOOTH_MSE
#define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
#endif
#ifndef SMOOTH_MSER
#define SMOOTH_MSER     0.07            /* acceptable relative MSE */
#endif
#define MAX_RAD         (GRIDRES/8)     /* maximum lobe radius */

#define RBFALLOCB       10              /* RBF allocation block size */

                                /* our loaded grid for this incident angle */
GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];

/* Start new DSF input grid */
void
new_bsdf_data(double new_theta, double new_phi)
{
        if (!new_input_direction(new_theta, new_phi))
                exit(1);
        memset(dsf_grid, 0, sizeof(dsf_grid));
}

/* Add BSDF data point */
void
add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
{
        FVECT   ovec;
        int     pos[2];

        if (!output_orient)             /* check output orientation */
                output_orient = 1 - 2*(theta_out > 90.);
        else if (output_orient > 0 ^ theta_out < 90.) {
                fputs("Cannot handle output angles on both sides of surface\n",
                                stderr);
                exit(1);
        }
        ovec[2] = sin((M_PI/180.)*theta_out);
        ovec[0] = cos((M_PI/180.)*phi_out) * ovec[2];
        ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
        ovec[2] = sqrt(1. - ovec[2]*ovec[2]);

        if (val <= 0)                   /* truncate to zero */
                val = 0;
        else if (!isDSF)
                val *= ovec[2];         /* convert from BSDF to DSF */

                                        /* update BSDF histogram */
        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
                ++bsdf_hist[histndx(val/ovec[2])];

        pos_from_vec(pos, ovec);

        dsf_grid[pos[0]][pos[1]].vsum += val;
        dsf_grid[pos[0]][pos[1]].nval++;
}

/* Compute minimum BSDF from histogram (does not clear) */
static void
comp_bsdf_min()
{
        int     cnt;
        int     i, target;

        cnt = 0;
        for (i = HISTLEN; i--; )
                cnt += bsdf_hist[i];
        if (!cnt) {                             /* shouldn't happen */
                bsdf_min = 0;
                return;
        }
        target = cnt/100;                       /* ignore bottom 1% */
        cnt = 0;
        for (i = 0; cnt <= target; i++)
                cnt += bsdf_hist[i];
        bsdf_min = histval(i-1);
}

/* Determine if the given region is empty of grid samples */
static int
empty_region(int x0, int x1, int y0, int y1)
{
        int     x, y;

        for (x = x0; x < x1; x++)
            for (y = y0; y < y1; y++)
                if (dsf_grid[x][y].nval)
                        return(0);
        return(1);
}

/* Determine if the given region is smooth enough to be a single lobe */
static int
smooth_region(int x0, int x1, int y0, int y1)
{
        RREAL   rMtx[3][3];
        FVECT   xvec;
        double  A, B, C, nvs, sqerr;
        int     x, y, n;
                                        /* compute planar regression */
        memset(rMtx, 0, sizeof(rMtx));
        memset(xvec, 0, sizeof(xvec));
        for (x = x0; x < x1; x++)
            for (y = y0; y < y1; y++)
                if ((n = dsf_grid[x][y].nval) > 0) {
                        double  z = dsf_grid[x][y].vsum;
                        rMtx[0][0] += x*x*(double)n;
                        rMtx[0][1] += x*y*(double)n;
                        rMtx[0][2] += x*(double)n;
                        rMtx[1][1] += y*y*(double)n;
                        rMtx[1][2] += y*(double)n;
                        rMtx[2][2] += (double)n;
                        xvec[0] += x*z;
                        xvec[1] += y*z;
                        xvec[2] += z;
                }
        rMtx[1][0] = rMtx[0][1];
        rMtx[2][1] = rMtx[1][2];
        nvs = rMtx[2][2];
        if (SDinvXform(rMtx, rMtx) != SDEnone)
                return(0);
        A = DOT(rMtx[0], xvec);
        B = DOT(rMtx[1], xvec);
        C = DOT(rMtx[2], xvec);
        sqerr = 0.0;                    /* compute mean squared error */
        for (x = x0; x < x1; x++)
            for (y = y0; y < y1; y++)
                if ((n = dsf_grid[x][y].nval) > 0) {
                        double  d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
                        sqerr += n*d*d;
                }
        if (sqerr <= nvs*SMOOTH_MSE)    /* below absolute MSE threshold? */
                return(1);
                                        /* OR below relative MSE threshold? */
        return(sqerr*nvs <= xvec[2]*xvec[2]*SMOOTH_MSER);
}

/* Create new lobe based on integrated samples in region */
static void
create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
{
        double  vtot = 0.0;
        int     nv = 0;
        double  rad;
        int     x, y;
                                        /* compute average for region */
        for (x = x0; x < x1; x++)
            for (y = y0; y < y1; y++) {
                vtot += dsf_grid[x][y].vsum;
                nv += dsf_grid[x][y].nval;
            }
        if (!nv) {
                fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
                                progname);
                exit(1);
        }
                                        /* peak value based on integral */
        vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
        rad = (RSCA/(double)GRIDRES)*(x1-x0);
        rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
        rvp->crad = ANG2R(rad);
        rvp->gx = (x0+x1)>>1;
        rvp->gy = (y0+y1)>>1;
}

/* Recursive function to build radial basis function representation */
static int
build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, int y0, int y1)
{
        int     xmid = (x0+x1)>>1;
        int     ymid = (y0+y1)>>1;
        int     branched[4];
        int     nadded, nleaves;
                                        /* need to make this a leaf? */
        if (empty_region(x0, xmid, y0, ymid) ||
                        empty_region(xmid, x1, y0, ymid) ||
                        empty_region(x0, xmid, ymid, y1) ||
                        empty_region(xmid, x1, ymid, y1))
                return(0);
                                        /* add children (branches+leaves) */
        if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
                return(-1);
        if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
                return(-1);
        if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
                return(-1);
        if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
                return(-1);
        nadded = branched[0] + branched[1] + branched[2] + branched[3];
        nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
        if (!nleaves)                   /* nothing but branches? */
                return(nadded);
                                        /* combine 4 leaves into 1? */
        if (nleaves == 4 && x1-x0 < MAX_RAD && smooth_region(x0, x1, y0, y1))
                return(0);
                                        /* need more array space? */
        if ((*np+nleaves-1)>>RBFALLOCB != (*np-1)>>RBFALLOCB) {
                *arp = (RBFVAL *)realloc(*arp,
                                sizeof(RBFVAL)*(*np+nleaves-1+(1<<RBFALLOCB)));
                if (*arp == NULL)
                        return(-1);
        }
                                        /* create lobes for leaves */
        if (!branched[0])
                create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
        if (!branched[1])
                create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
        if (!branched[2])
                create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
        if (!branched[3])
                create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
        nadded += nleaves;
        return(nadded);
}

/* Count up filled nodes and build RBF representation from current grid */ 
RBFNODE *
make_rbfrep()
{
        RBFNODE *newnode;
        RBFVAL  *rbfarr;
        int     nn;
                                /* compute minimum BSDF */
        comp_bsdf_min();
                                /* create RBF node list */
        rbfarr = NULL; nn = 0;
        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
                goto memerr;
                                /* (re)allocate RBF array */
        newnode = (RBFNODE *)realloc(rbfarr,
                        sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
        if (newnode == NULL)
                goto memerr;
                                /* copy computed lobes into RBF node */
        memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
        newnode->ord = -1;
        newnode->next = NULL;
        newnode->ejl = NULL;
        newnode->invec[2] = sin((M_PI/180.)*theta_in_deg);
        newnode->invec[0] = cos((M_PI/180.)*phi_in_deg)*newnode->invec[2];
        newnode->invec[1] = sin((M_PI/180.)*phi_in_deg)*newnode->invec[2];
        newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
        newnode->vtotal = .0;
        newnode->nrbf = nn;
                                /* compute sum for normalization */
        while (nn-- > 0)
                newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
#ifdef DEBUG
        fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
        fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
                        get_theta180(newnode->invec), get_phi360(newnode->invec),
                        newnode->vtotal);
#endif
        insert_dsf(newnode);

        return(newnode);
memerr:
        fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
        exit(1);
}
Revision:	2.13
Committed:	Mon Oct 21 21:48:42 2013 UTC (10 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.12:	+22 -8 lines
Log Message:	Fixed potential bug in smoothness criterion
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.13	static const char RCSid[] = "$Id: bsdfrbf.c,v 2.12 2013/10/21 18:33:15 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Radial basis function representation for BSDF data.
6			*
7			* G. Ward
8			*/
9
10	greg	2.13	/****************************************************************
11			1) Collect samples into a grid using the Shirley-Chiu
12			angular mapping from a hemisphere to a square.
13
14			2) Compute an adaptive quadtree by subdividing the grid so that
15			each leaf node has at least one sample up to as many
16			samples as fit nicely on a plane to within a certain
17			MSE tolerance.
18
19			3) Place one Gaussian lobe at each leaf node in the quadtree,
20			sizing it to have a radius equal to the leaf size and
21			a volume equal to the energy in that node.
22			*****************************************************************/
23
24	greg	2.1	#define _USE_MATH_DEFINES
25			#include <stdio.h>
26			#include <stdlib.h>
27			#include <string.h>
28			#include <math.h>
29			#include "bsdfrep.h"
30
31	greg	2.12	#ifndef RSCA
32			#define RSCA 2.2 /* radius scaling factor (empirical) */
33	greg	2.9	#endif
34	greg	2.12	#ifndef SMOOTH_MSE
35			#define SMOOTH_MSE 5e-5 /* acceptable mean squared error */
36	greg	2.1	#endif
37	greg	2.12	#ifndef SMOOTH_MSER
38			#define SMOOTH_MSER 0.07 /* acceptable relative MSE */
39	greg	2.7	#endif
40	greg	2.12	#define MAX_RAD (GRIDRES/8) /* maximum lobe radius */
41
42			#define RBFALLOCB 10 /* RBF allocation block size */
43
44	greg	2.1	/* our loaded grid for this incident angle */
45			GRIDVAL dsf_grid[GRIDRES][GRIDRES];
46
47			/* Start new DSF input grid */
48			void
49			new_bsdf_data(double new_theta, double new_phi)
50			{
51			if (!new_input_direction(new_theta, new_phi))
52			exit(1);
53			memset(dsf_grid, 0, sizeof(dsf_grid));
54			}
55
56			/* Add BSDF data point */
57			void
58			add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
59			{
60			FVECT ovec;
61			int pos[2];
62
63			if (!output_orient) /* check output orientation */
64			output_orient = 1 - 2*(theta_out > 90.);
65			else if (output_orient > 0 ^ theta_out < 90.) {
66			fputs("Cannot handle output angles on both sides of surface\n",
67			stderr);
68			exit(1);
69			}
70			ovec[2] = sin((M_PI/180.)*theta_out);
71			ovec[0] = cos((M_PI/180.)phi_out) ovec[2];
72			ovec[1] = sin((M_PI/180.)phi_out) ovec[2];
73			ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
74
75	greg	2.8	if (val <= 0) /* truncate to zero */
76			val = 0;
77			else if (!isDSF)
78	greg	2.1	val = ovec[2]; / convert from BSDF to DSF */
79
80	greg	2.4	/* update BSDF histogram */
81			if (val < BSDF2BIGovec[2] && val > BSDF2SMLovec[2])
82			++bsdf_hist[histndx(val/ovec[2])];
83
84	greg	2.1	pos_from_vec(pos, ovec);
85
86			dsf_grid[pos[0]][pos[1]].vsum += val;
87			dsf_grid[pos[0]][pos[1]].nval++;
88			}
89
90	greg	2.11	/* Compute minimum BSDF from histogram (does not clear) */
91	greg	2.5	static void
92			comp_bsdf_min()
93			{
94			int cnt;
95			int i, target;
96
97			cnt = 0;
98			for (i = HISTLEN; i--; )
99			cnt += bsdf_hist[i];
100			if (!cnt) { /* shouldn't happen */
101			bsdf_min = 0;
102			return;
103			}
104			target = cnt/100; /* ignore bottom 1% */
105			cnt = 0;
106			for (i = 0; cnt <= target; i++)
107			cnt += bsdf_hist[i];
108			bsdf_min = histval(i-1);
109			}
110
111	greg	2.12	/* Determine if the given region is empty of grid samples */
112	greg	2.6	static int
113	greg	2.12	empty_region(int x0, int x1, int y0, int y1)
114	greg	2.6	{
115	greg	2.12	int x, y;
116
117			for (x = x0; x < x1; x++)
118			for (y = y0; y < y1; y++)
119			if (dsf_grid[x][y].nval)
120			return(0);
121			return(1);
122			}
123
124			/* Determine if the given region is smooth enough to be a single lobe */
125			static int
126			smooth_region(int x0, int x1, int y0, int y1)
127			{
128			RREAL rMtx[3][3];
129			FVECT xvec;
130			double A, B, C, nvs, sqerr;
131			int x, y, n;
132			/* compute planar regression */
133			memset(rMtx, 0, sizeof(rMtx));
134			memset(xvec, 0, sizeof(xvec));
135			for (x = x0; x < x1; x++)
136			for (y = y0; y < y1; y++)
137			if ((n = dsf_grid[x][y].nval) > 0) {
138			double z = dsf_grid[x][y].vsum;
139	greg	2.13	rMtx[0][0] += xx(double)n;
140			rMtx[0][1] += xy(double)n;
141			rMtx[0][2] += x*(double)n;
142			rMtx[1][1] += yy(double)n;
143			rMtx[1][2] += y*(double)n;
144			rMtx[2][2] += (double)n;
145	greg	2.12	xvec[0] += x*z;
146			xvec[1] += y*z;
147			xvec[2] += z;
148			}
149			rMtx[1][0] = rMtx[0][1];
150			rMtx[2][1] = rMtx[1][2];
151			nvs = rMtx[2][2];
152			if (SDinvXform(rMtx, rMtx) != SDEnone)
153			return(0);
154			A = DOT(rMtx[0], xvec);
155			B = DOT(rMtx[1], xvec);
156			C = DOT(rMtx[2], xvec);
157			sqerr = 0.0; /* compute mean squared error */
158			for (x = x0; x < x1; x++)
159			for (y = y0; y < y1; y++)
160			if ((n = dsf_grid[x][y].nval) > 0) {
161			double d = Ax + By + C - dsf_grid[x][y].vsum/n;
162			sqerr += ndd;
163	greg	2.6	}
164	greg	2.12	if (sqerr <= nvsSMOOTH_MSE) / below absolute MSE threshold? */
165			return(1);
166	greg	2.13	/* OR below relative MSE threshold? */
167	greg	2.12	return(sqerrnvs <= xvec[2]xvec[2]*SMOOTH_MSER);
168			}
169
170			/* Create new lobe based on integrated samples in region */
171			static void
172			create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
173			{
174			double vtot = 0.0;
175			int nv = 0;
176			double rad;
177			int x, y;
178			/* compute average for region */
179			for (x = x0; x < x1; x++)
180			for (y = y0; y < y1; y++) {
181			vtot += dsf_grid[x][y].vsum;
182			nv += dsf_grid[x][y].nval;
183			}
184			if (!nv) {
185			fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
186			progname);
187			exit(1);
188	greg	2.6	}
189	greg	2.12	/* peak value based on integral */
190			vtot = (x1-x0)(y1-y0)(2.M_PI/GRIDRES/GRIDRES)/(double)nv;
191			rad = (RSCA/(double)GRIDRES)*(x1-x0);
192			rvp->peak = vtot / ((2.M_PI) rad*rad);
193			rvp->crad = ANG2R(rad);
194			rvp->gx = (x0+x1)>>1;
195			rvp->gy = (y0+y1)>>1;
196	greg	2.6	}
197
198	greg	2.12	/* Recursive function to build radial basis function representation */
199	greg	2.6	static int
200	greg	2.12	build_rbfrep(RBFVAL *arp, int np, int x0, int x1, int y0, int y1)
201	greg	2.6	{
202	greg	2.12	int xmid = (x0+x1)>>1;
203			int ymid = (y0+y1)>>1;
204			int branched[4];
205			int nadded, nleaves;
206			/* need to make this a leaf? */
207			if (empty_region(x0, xmid, y0, ymid) \|\|
208			empty_region(xmid, x1, y0, ymid) \|\|
209			empty_region(x0, xmid, ymid, y1) \|\|
210			empty_region(xmid, x1, ymid, y1))
211			return(0);
212			/* add children (branches+leaves) */
213			if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
214			return(-1);
215			if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
216			return(-1);
217			if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
218			return(-1);
219			if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
220			return(-1);
221			nadded = branched[0] + branched[1] + branched[2] + branched[3];
222			nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
223			if (!nleaves) /* nothing but branches? */
224			return(nadded);
225			/* combine 4 leaves into 1? */
226			if (nleaves == 4 && x1-x0 < MAX_RAD && smooth_region(x0, x1, y0, y1))
227			return(0);
228			/* need more array space? */
229			if ((np+nleaves-1)>>RBFALLOCB != (np-1)>>RBFALLOCB) {
230			arp = (RBFVAL )realloc(*arp,
231			sizeof(RBFVAL)(np+nleaves-1+(1<<RBFALLOCB)));
232			if (*arp == NULL)
233			return(-1);
234	greg	2.6	}
235	greg	2.12	/* create lobes for leaves */
236			if (!branched[0])
237			create_lobe(arp+(np)++, x0, xmid, y0, ymid);
238			if (!branched[1])
239			create_lobe(arp+(np)++, xmid, x1, y0, ymid);
240			if (!branched[2])
241			create_lobe(arp+(np)++, x0, xmid, ymid, y1);
242			if (!branched[3])
243			create_lobe(arp+(np)++, xmid, x1, ymid, y1);
244			nadded += nleaves;
245			return(nadded);
246	greg	2.6	}
247
248	greg	2.1	/* Count up filled nodes and build RBF representation from current grid */
249			RBFNODE *
250	greg	2.12	make_rbfrep()
251	greg	2.1	{
252	greg	2.12	RBFNODE *newnode;
253			RBFVAL *rbfarr;
254	greg	2.1	int nn;
255	greg	2.5	/* compute minimum BSDF */
256			comp_bsdf_min();
257	greg	2.12	/* create RBF node list */
258			rbfarr = NULL; nn = 0;
259			if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
260			goto memerr;
261			/* (re)allocate RBF array */
262			newnode = (RBFNODE *)realloc(rbfarr,
263			sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
264	greg	2.2	if (newnode == NULL)
265			goto memerr;
266	greg	2.12	/* copy computed lobes into RBF node */
267			memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
268	greg	2.1	newnode->ord = -1;
269			newnode->next = NULL;
270			newnode->ejl = NULL;
271			newnode->invec[2] = sin((M_PI/180.)*theta_in_deg);
272			newnode->invec[0] = cos((M_PI/180.)phi_in_deg)newnode->invec[2];
273			newnode->invec[1] = sin((M_PI/180.)phi_in_deg)newnode->invec[2];
274			newnode->invec[2] = input_orientsqrt(1. - newnode->invec[2]newnode->invec[2]);
275	greg	2.12	newnode->vtotal = .0;
276	greg	2.1	newnode->nrbf = nn;
277	greg	2.12	/* compute sum for normalization */
278			while (nn-- > 0)
279			newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
280	greg	2.3	#ifdef DEBUG
281	greg	2.12	fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
282	greg	2.3	fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
283			get_theta180(newnode->invec), get_phi360(newnode->invec),
284			newnode->vtotal);
285			#endif
286	greg	2.1	insert_dsf(newnode);
287
288			return(newnode);
289	greg	2.2	memerr:
290			fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
291			exit(1);
292	greg	2.1	}