src/cv/bsdfrbf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrbf.c,v 2.11 2013/10/19 00:11:50 greg Exp $";
#endif
/*
 * Radial basis function representation for BSDF data.
 *
 *      G. Ward
 */

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