src/cv/bsdfrbf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrbf.c,v 2.2 2012/11/13 04:23:38 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.7             /* radius scaling factor (empirical) */
#endif
                                /* 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 (!isDSF)
                val *= ovec[2];         /* convert from BSDF to DSF */

        pos_from_vec(pos, ovec);

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

/* Compute radii for non-empty bins */
/* (distance to furthest empty bin for which non-empty bin is the closest) */
static void
compute_radii(void)
{
        unsigned int    fill_grid[GRIDRES][GRIDRES];
        unsigned short  fill_cnt[GRIDRES][GRIDRES];
        FVECT           ovec0, ovec1;
        double          ang2, lastang2;
        int             r, i, j, jn, ii, jj, inear, jnear;

        r = GRIDRES/2;                          /* proceed in zig-zag */
        for (i = 0; i < GRIDRES; i++)
            for (jn = 0; jn < GRIDRES; jn++) {
                j = (i&1) ? jn : GRIDRES-1-jn;
                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
                        continue;
                ovec_from_pos(ovec0, i, j);
                inear = jnear = -1;             /* find nearest non-empty */
                lastang2 = M_PI*M_PI;
                for (ii = i-r; ii <= i+r; ii++) {
                    if (ii < 0) continue;
                    if (ii >= GRIDRES) break;
                    for (jj = j-r; jj <= j+r; jj++) {
                        if (jj < 0) continue;
                        if (jj >= GRIDRES) break;
                        if (!dsf_grid[ii][jj].nval)
                                continue;
                        ovec_from_pos(ovec1, ii, jj);
                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
                        if (ang2 >= lastang2)
                                continue;
                        lastang2 = ang2;
                        inear = ii; jnear = jj;
                    }
                }
                if (inear < 0) {
                        fprintf(stderr,
                                "%s: Could not find non-empty neighbor!\n",
                                        progname);
                        exit(1);
                }
                ang2 = sqrt(lastang2);
                r = ANG2R(ang2);                /* record if > previous */
                if (r > dsf_grid[inear][jnear].crad)
                        dsf_grid[inear][jnear].crad = r;
                                                /* next search radius */
                r = ang2*(2.*GRIDRES/M_PI) + 3;
            }
                                                /* blur radii over hemisphere */
        memset(fill_grid, 0, sizeof(fill_grid));
        memset(fill_cnt, 0, sizeof(fill_cnt));
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++) {
                if (!dsf_grid[i][j].crad)
                        continue;               /* missing distance */
                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
                for (ii = i-r; ii <= i+r; ii++) {
                    if (ii < 0) continue;
                    if (ii >= GRIDRES) break;
                    for (jj = j-r; jj <= j+r; jj++) {
                        if (jj < 0) continue;
                        if (jj >= GRIDRES) break;
                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
                                continue;
                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
                        fill_cnt[ii][jj]++;
                    }
                }
            }
                                                /* copy back blurred radii */
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++)
                if (fill_cnt[i][j])
                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
}

/* Cull points for more uniform distribution, leave all nval 0 or 1 */
static void
cull_values(void)
{
        FVECT   ovec0, ovec1;
        double  maxang, maxang2;
        int     i, j, ii, jj, r;
                                                /* simple greedy algorithm */
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++) {
                if (!dsf_grid[i][j].nval)
                        continue;
                if (!dsf_grid[i][j].crad)
                        continue;               /* shouldn't happen */
                ovec_from_pos(ovec0, i, j);
                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
                if (maxang > ovec0[2])          /* clamp near horizon */
                        maxang = ovec0[2];
                r = maxang*(2.*GRIDRES/M_PI) + 1;
                maxang2 = maxang*maxang;
                for (ii = i-r; ii <= i+r; ii++) {
                    if (ii < 0) continue;
                    if (ii >= GRIDRES) break;
                    for (jj = j-r; jj <= j+r; jj++) {
                        if (jj < 0) continue;
                        if (jj >= GRIDRES) break;
                        if (!dsf_grid[ii][jj].nval)
                                continue;
                        if ((ii == i) & (jj == j))
                                continue;       /* don't get self-absorbed */
                        ovec_from_pos(ovec1, ii, jj);
                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
                                continue;
                                                /* absorb sum */
                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
                                                /* keep value, though */
                        dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
                        dsf_grid[ii][jj].nval = 0;
                    }
                }
            }
                                                /* final averaging pass */
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++)
                if (dsf_grid[i][j].nval > 1) {
                        dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
                        dsf_grid[i][j].nval = 1;
                }
}

/* Count up filled nodes and build RBF representation from current grid */ 
RBFNODE *
make_rbfrep(void)
{
        int     niter = 16;
        double  lastVar, thisVar = 100.;
        int     nn;
        RBFNODE *newnode;
        RBFVAL  *itera;
        int     i, j;
                                /* compute RBF radii */
        compute_radii();
                                /* coagulate lobes */
        cull_values();
        nn = 0;                 /* count selected bins */
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++)
                nn += dsf_grid[i][j].nval;
                                /* allocate RBF array */
        newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
        if (newnode == NULL)
                goto memerr;
        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;
        nn = 0;                 /* fill RBF array */
        for (i = 0; i < GRIDRES; i++)
            for (j = 0; j < GRIDRES; j++)
                if (dsf_grid[i][j].nval) {
                        newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
                        newnode->rbfa[nn].gx = i;
                        newnode->rbfa[nn].gy = j;
                        ++nn;
                }
                                /* iterate to improve interpolation accuracy */
        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
        if (itera == NULL)
                goto memerr;
        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
        do {
                double  dsum = 0, dsum2 = 0;
                nn = 0;
                for (i = 0; i < GRIDRES; i++)
                    for (j = 0; j < GRIDRES; j++)
                        if (dsf_grid[i][j].nval) {
                                FVECT   odir;
                                double  corr;
                                ovec_from_pos(odir, i, j);
                                itera[nn++].peak *= corr =
                                        dsf_grid[i][j].vsum /
                                                eval_rbfrep(newnode, odir);
                                dsum += 1. - corr;
                                dsum2 += (1.-corr)*(1.-corr);
                        }
                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
                lastVar = thisVar;
                thisVar = dsum2/(double)nn;
#ifdef DEBUG
                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
                                        100.*dsum/(double)nn,
                                        100.*sqrt(thisVar));
#endif
        } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);

        free(itera);
        nn = 0;                 /* compute sum for normalization */
        while (nn < newnode->nrbf)
                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
#ifdef DEBUG
        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.3
Committed:	Thu Nov 22 06:07:17 2012 UTC (11 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.2:	+6 -2 lines
Log Message:	Bug fix in geodesic() and other minor improvements
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.3	static const char RCSid[] = "$Id: bsdfrbf.c,v 2.2 2012/11/13 04:23:38 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			#ifndef RSCA
18			#define RSCA 2.7 /* radius scaling factor (empirical) */
19			#endif
20			/* our loaded grid for this incident angle */
21			GRIDVAL dsf_grid[GRIDRES][GRIDRES];
22
23			/* Start new DSF input grid */
24			void
25			new_bsdf_data(double new_theta, double new_phi)
26			{
27			if (!new_input_direction(new_theta, new_phi))
28			exit(1);
29			memset(dsf_grid, 0, sizeof(dsf_grid));
30			}
31
32			/* Add BSDF data point */
33			void
34			add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
35			{
36			FVECT ovec;
37			int pos[2];
38
39			if (!output_orient) /* check output orientation */
40			output_orient = 1 - 2*(theta_out > 90.);
41			else if (output_orient > 0 ^ theta_out < 90.) {
42			fputs("Cannot handle output angles on both sides of surface\n",
43			stderr);
44			exit(1);
45			}
46			ovec[2] = sin((M_PI/180.)*theta_out);
47			ovec[0] = cos((M_PI/180.)phi_out) ovec[2];
48			ovec[1] = sin((M_PI/180.)phi_out) ovec[2];
49			ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
50
51			if (!isDSF)
52			val = ovec[2]; / convert from BSDF to DSF */
53
54			pos_from_vec(pos, ovec);
55
56			dsf_grid[pos[0]][pos[1]].vsum += val;
57			dsf_grid[pos[0]][pos[1]].nval++;
58			}
59
60			/* Compute radii for non-empty bins */
61			/* (distance to furthest empty bin for which non-empty bin is the closest) */
62			static void
63			compute_radii(void)
64			{
65			unsigned int fill_grid[GRIDRES][GRIDRES];
66			unsigned short fill_cnt[GRIDRES][GRIDRES];
67			FVECT ovec0, ovec1;
68			double ang2, lastang2;
69			int r, i, j, jn, ii, jj, inear, jnear;
70
71			r = GRIDRES/2; /* proceed in zig-zag */
72			for (i = 0; i < GRIDRES; i++)
73			for (jn = 0; jn < GRIDRES; jn++) {
74			j = (i&1) ? jn : GRIDRES-1-jn;
75			if (dsf_grid[i][j].nval) /* find empty grid pos. */
76			continue;
77			ovec_from_pos(ovec0, i, j);
78			inear = jnear = -1; /* find nearest non-empty */
79			lastang2 = M_PI*M_PI;
80			for (ii = i-r; ii <= i+r; ii++) {
81			if (ii < 0) continue;
82			if (ii >= GRIDRES) break;
83			for (jj = j-r; jj <= j+r; jj++) {
84			if (jj < 0) continue;
85			if (jj >= GRIDRES) break;
86			if (!dsf_grid[ii][jj].nval)
87			continue;
88			ovec_from_pos(ovec1, ii, jj);
89			ang2 = 2. - 2.*DOT(ovec0,ovec1);
90			if (ang2 >= lastang2)
91			continue;
92			lastang2 = ang2;
93			inear = ii; jnear = jj;
94			}
95			}
96			if (inear < 0) {
97			fprintf(stderr,
98			"%s: Could not find non-empty neighbor!\n",
99			progname);
100			exit(1);
101			}
102			ang2 = sqrt(lastang2);
103			r = ANG2R(ang2); /* record if > previous */
104			if (r > dsf_grid[inear][jnear].crad)
105			dsf_grid[inear][jnear].crad = r;
106			/* next search radius */
107			r = ang2(2.GRIDRES/M_PI) + 3;
108			}
109			/* blur radii over hemisphere */
110			memset(fill_grid, 0, sizeof(fill_grid));
111			memset(fill_cnt, 0, sizeof(fill_cnt));
112			for (i = 0; i < GRIDRES; i++)
113			for (j = 0; j < GRIDRES; j++) {
114			if (!dsf_grid[i][j].crad)
115			continue; /* missing distance */
116			r = R2ANG(dsf_grid[i][j].crad)(2.RSCA*GRIDRES/M_PI);
117			for (ii = i-r; ii <= i+r; ii++) {
118			if (ii < 0) continue;
119			if (ii >= GRIDRES) break;
120			for (jj = j-r; jj <= j+r; jj++) {
121			if (jj < 0) continue;
122			if (jj >= GRIDRES) break;
123			if ((ii-i)(ii-i) + (jj-j)(jj-j) > r*r)
124			continue;
125			fill_grid[ii][jj] += dsf_grid[i][j].crad;
126			fill_cnt[ii][jj]++;
127			}
128			}
129			}
130			/* copy back blurred radii */
131			for (i = 0; i < GRIDRES; i++)
132			for (j = 0; j < GRIDRES; j++)
133			if (fill_cnt[i][j])
134			dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
135			}
136
137			/* Cull points for more uniform distribution, leave all nval 0 or 1 */
138			static void
139			cull_values(void)
140			{
141			FVECT ovec0, ovec1;
142			double maxang, maxang2;
143			int i, j, ii, jj, r;
144			/* simple greedy algorithm */
145			for (i = 0; i < GRIDRES; i++)
146			for (j = 0; j < GRIDRES; j++) {
147			if (!dsf_grid[i][j].nval)
148			continue;
149			if (!dsf_grid[i][j].crad)
150			continue; /* shouldn't happen */
151			ovec_from_pos(ovec0, i, j);
152			maxang = 2.*R2ANG(dsf_grid[i][j].crad);
153			if (maxang > ovec0[2]) /* clamp near horizon */
154			maxang = ovec0[2];
155			r = maxang(2.GRIDRES/M_PI) + 1;
156			maxang2 = maxang*maxang;
157			for (ii = i-r; ii <= i+r; ii++) {
158			if (ii < 0) continue;
159			if (ii >= GRIDRES) break;
160			for (jj = j-r; jj <= j+r; jj++) {
161			if (jj < 0) continue;
162			if (jj >= GRIDRES) break;
163			if (!dsf_grid[ii][jj].nval)
164			continue;
165			if ((ii == i) & (jj == j))
166			continue; /* don't get self-absorbed */
167			ovec_from_pos(ovec1, ii, jj);
168			if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
169			continue;
170			/* absorb sum */
171			dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
172			dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
173			/* keep value, though */
174			dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
175			dsf_grid[ii][jj].nval = 0;
176			}
177			}
178			}
179			/* final averaging pass */
180			for (i = 0; i < GRIDRES; i++)
181			for (j = 0; j < GRIDRES; j++)
182			if (dsf_grid[i][j].nval > 1) {
183			dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
184			dsf_grid[i][j].nval = 1;
185			}
186			}
187
188			/* Count up filled nodes and build RBF representation from current grid */
189			RBFNODE *
190			make_rbfrep(void)
191			{
192			int niter = 16;
193			double lastVar, thisVar = 100.;
194			int nn;
195			RBFNODE *newnode;
196	greg	2.2	RBFVAL *itera;
197	greg	2.1	int i, j;
198			/* compute RBF radii */
199			compute_radii();
200			/* coagulate lobes */
201			cull_values();
202			nn = 0; /* count selected bins */
203			for (i = 0; i < GRIDRES; i++)
204			for (j = 0; j < GRIDRES; j++)
205			nn += dsf_grid[i][j].nval;
206			/* allocate RBF array */
207			newnode = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(nn-1));
208	greg	2.2	if (newnode == NULL)
209			goto memerr;
210	greg	2.1	newnode->ord = -1;
211			newnode->next = NULL;
212			newnode->ejl = NULL;
213			newnode->invec[2] = sin((M_PI/180.)*theta_in_deg);
214			newnode->invec[0] = cos((M_PI/180.)phi_in_deg)newnode->invec[2];
215			newnode->invec[1] = sin((M_PI/180.)phi_in_deg)newnode->invec[2];
216			newnode->invec[2] = input_orientsqrt(1. - newnode->invec[2]newnode->invec[2]);
217			newnode->vtotal = 0;
218			newnode->nrbf = nn;
219			nn = 0; /* fill RBF array */
220			for (i = 0; i < GRIDRES; i++)
221			for (j = 0; j < GRIDRES; j++)
222			if (dsf_grid[i][j].nval) {
223			newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
224			newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
225			newnode->rbfa[nn].gx = i;
226			newnode->rbfa[nn].gy = j;
227			++nn;
228			}
229			/* iterate to improve interpolation accuracy */
230	greg	2.2	itera = (RBFVAL )malloc(sizeof(RBFVAL)newnode->nrbf);
231			if (itera == NULL)
232			goto memerr;
233			memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
234	greg	2.1	do {
235			double dsum = 0, dsum2 = 0;
236			nn = 0;
237			for (i = 0; i < GRIDRES; i++)
238			for (j = 0; j < GRIDRES; j++)
239			if (dsf_grid[i][j].nval) {
240			FVECT odir;
241			double corr;
242			ovec_from_pos(odir, i, j);
243	greg	2.2	itera[nn++].peak *= corr =
244	greg	2.1	dsf_grid[i][j].vsum /
245			eval_rbfrep(newnode, odir);
246	greg	2.2	dsum += 1. - corr;
247			dsum2 += (1.-corr)*(1.-corr);
248	greg	2.1	}
249	greg	2.2	memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
250	greg	2.1	lastVar = thisVar;
251			thisVar = dsum2/(double)nn;
252			#ifdef DEBUG
253			fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
254			100.*dsum/(double)nn,
255			100.*sqrt(thisVar));
256			#endif
257			} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
258
259	greg	2.2	free(itera);
260	greg	2.1	nn = 0; /* compute sum for normalization */
261			while (nn < newnode->nrbf)
262			newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
263	greg	2.3	#ifdef DEBUG
264			fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
265			get_theta180(newnode->invec), get_phi360(newnode->invec),
266			newnode->vtotal);
267			#endif
268	greg	2.1	insert_dsf(newnode);
269
270			return(newnode);
271	greg	2.2	memerr:
272			fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
273			exit(1);
274	greg	2.1	}