src/cv/bsdfrbf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrbf.c,v 2.30 2016/12/06 23:39:54 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.0             /* radius scaling factor (empirical) */
#endif
#ifndef MAXSLOPE
#define MAXSLOPE        200.0           /* maximum slope for smooth region */
#endif
#ifndef SMOOTH_MSE
#define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
#endif
#ifndef SMOOTH_MSER
#define SMOOTH_MSER     0.01            /* acceptable relative MSE */
#endif
#define MAX_RAD         (GRIDRES/8)     /* maximum lobe radius */

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

                                        /* loaded grid or comparison DSFs */
GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
                                        /* allocated chrominance sums if any */
float                   (*spec_grid)[GRIDRES][GRIDRES];
int                     nspec_grid = 0;

/* Set up visible spectrum sampling */
void
set_spectral_samples(int nspec)
{
        if (rbf_colorimetry == RBCunknown) {
                if (nspec_grid > 0) {
                        free(spec_grid);
                        spec_grid = NULL;
                        nspec_grid = 0;
                }
                if (nspec == 1) {
                        rbf_colorimetry = RBCphotopic;
                        return;
                }
                if (nspec == 3) {
                        rbf_colorimetry = RBCtristimulus;
                        spec_grid = (float (*)[GRIDRES][GRIDRES])calloc(
                                        2*GRIDRES*GRIDRES, sizeof(float) );
                        if (spec_grid == NULL)
                                goto mem_error;
                        nspec_grid = 2;
                        return;
                }
                fprintf(stderr,
                        "%s: only 1 or 3 spectral samples currently supported\n",
                                progname);
                exit(1);
        }
        if (nspec != nspec_grid+1) {
                fprintf(stderr,
                        "%s: number of spectral samples cannot be changed\n",
                                progname);
                exit(1);
        }
        return;
mem_error:
        fprintf(stderr, "%s: out of memory in set_spectral_samples()\n",
                        progname);
        exit(1);
}

/* 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));
        if (nspec_grid > 0)
                memset(spec_grid, 0, sizeof(float)*GRIDRES*GRIDRES*nspec_grid);
}

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

        if (nspec_grid > 2) {
                fprintf(stderr, "%s: unsupported color space\n", progname);
                exit(1);
        }
        if (!output_orient)             /* check output orientation */
                output_orient = 1 - 2*(theta_out > 90.);
        else if (output_orient > 0 ^ theta_out < 90.) {
                fprintf(stderr,
                "%s: cannot handle output angles on both sides of surface\n",
                                progname);
                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]);
                                        /* BSDF to DSF correction */
        cfact = isDSF ? 1. : ovec[2];

        Yval = cfact * val[rbf_colorimetry==RBCtristimulus];
                                        /* update BSDF histogram */
        if (BSDF2SML*ovec[2] < Yval && Yval < BSDF2BIG*ovec[2])
                ++bsdf_hist[histndx(Yval/ovec[2])];

        pos_from_vec(pos, ovec);

        dsf_grid[pos[0]][pos[1]].sum.v += Yval;
        dsf_grid[pos[0]][pos[1]].sum.n++;
                                        /* add in X and Z values */
        if (rbf_colorimetry == RBCtristimulus) {
                spec_grid[0][pos[0]][pos[1]] += cfact * val[0];
                spec_grid[1][pos[0]][pos[1]] += cfact * val[2];
        }
}

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

        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].sum.n)
                        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].sum.n) > 0) {
                        double  z = dsf_grid[x][y].sum.v;
                        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][0] = rMtx[0][2];
        rMtx[2][1] = rMtx[1][2];
        nvs = rMtx[2][2];
        if (SDinvXform(rMtx, rMtx) != SDEnone)
                return(1);              /* colinear values */
        A = DOT(rMtx[0], xvec);
        B = DOT(rMtx[1], xvec);
        if (A*A + B*B > MAXSLOPE*MAXSLOPE)      /* too steep? */
                return(0);
        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].sum.n) > 0) {
                        double  d = A*x + B*y + C - dsf_grid[x][y].sum.v/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 int
create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
{
        double  vtot = 0.0;
        double  CIEXtot = 0.0, CIEZtot = 0.0;
        int     nv = 0;
        double  wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
        double  rad;
        int     x, y;
                                        /* compute average for region */
        for (x = x0; x < x1; x++)
            for (y = y0; y < y1; y++)
                if (dsf_grid[x][y].sum.n) {
                    const double        v = dsf_grid[x][y].sum.v;
                    const int           n = dsf_grid[x][y].sum.n;

                    if (v > 0) {
                        const double    wt = v / (double)n;
                        wxsum += wt * x;
                        wysum += wt * y;
                        wtsum += wt;
                    }
                    vtot += v;
                    nv += n;
                    if (rbf_colorimetry == RBCtristimulus) {
                        CIEXtot += spec_grid[0][x][y];
                        CIEZtot += spec_grid[1][x][y];
                    }
                }
        if (!nv) {
                fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
                                progname);
                exit(1);
        }
        if (vtot <= 0)                  /* only create positive lobes */
                return(0);
                                        /* assign color */
        if (rbf_colorimetry == RBCtristimulus) {
                const double    df = 1.0 / (CIEXtot + vtot + CIEZtot);
                C_COLOR         cclr;
                c_cset(&cclr, CIEXtot*df, vtot*df);
                rvp->chroma = c_encodeChroma(&cclr);
        } else
                rvp->chroma = c_dfchroma;
                                        /* 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);         /* put peak at centroid */
        rvp->gx = (int)(wxsum/wtsum + .5);
        rvp->gy = (int)(wysum/wtsum + .5);
        return(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)) {
                ++(*np); ++nadded;
        }
        if (!branched[1] && create_lobe(*arp+*np, xmid, x1, y0, ymid)) {
                ++(*np); ++nadded;
        }
        if (!branched[2] && create_lobe(*arp+*np, x0, xmid, ymid, y1)) {
                ++(*np); ++nadded;
        }
        if (!branched[3] && create_lobe(*arp+*np, xmid, x1, ymid, y1)) {
                ++(*np); ++nadded;
        }
        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) {
                if (nn)
                        goto memerr;
                fprintf(stderr,
                        "%s: warning - skipping bad incidence (%.1f,%.1f)\n",
                                progname, theta_in_deg, phi_in_deg);
                return(NULL);
        }
                                /* (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.31
Committed:	Wed Aug 1 17:00:22 2018 UTC (6 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2
Changes since 2.30:	+2 -2 lines
Log Message:	Reduced relative error tolerance value
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.31	static const char RCSid[] = "$Id: bsdfrbf.c,v 2.30 2016/12/06 23:39:54 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	greg	2.25	#define RSCA 2.0 /* radius scaling factor (empirical) */
33	greg	2.9	#endif
34	greg	2.29	#ifndef MAXSLOPE
35	greg	2.30	#define MAXSLOPE 200.0 /* maximum slope for smooth region */
36	greg	2.29	#endif
37	greg	2.12	#ifndef SMOOTH_MSE
38	greg	2.19	#define SMOOTH_MSE 5e-5 /* acceptable mean squared error */
39	greg	2.1	#endif
40	greg	2.12	#ifndef SMOOTH_MSER
41	greg	2.31	#define SMOOTH_MSER 0.01 /* acceptable relative MSE */
42	greg	2.7	#endif
43	greg	2.12	#define MAX_RAD (GRIDRES/8) /* maximum lobe radius */
44
45			#define RBFALLOCB 10 /* RBF allocation block size */
46
47	greg	2.25	/* loaded grid or comparison DSFs */
48	greg	2.1	GRIDVAL dsf_grid[GRIDRES][GRIDRES];
49	greg	2.27	/* allocated chrominance sums if any */
50			float (*spec_grid)[GRIDRES][GRIDRES];
51			int nspec_grid = 0;
52
53			/* Set up visible spectrum sampling */
54			void
55			set_spectral_samples(int nspec)
56			{
57			if (rbf_colorimetry == RBCunknown) {
58			if (nspec_grid > 0) {
59			free(spec_grid);
60			spec_grid = NULL;
61			nspec_grid = 0;
62			}
63			if (nspec == 1) {
64			rbf_colorimetry = RBCphotopic;
65			return;
66			}
67			if (nspec == 3) {
68			rbf_colorimetry = RBCtristimulus;
69			spec_grid = (float (*)[GRIDRES][GRIDRES])calloc(
70			2GRIDRESGRIDRES, sizeof(float) );
71			if (spec_grid == NULL)
72			goto mem_error;
73			nspec_grid = 2;
74			return;
75			}
76			fprintf(stderr,
77			"%s: only 1 or 3 spectral samples currently supported\n",
78			progname);
79			exit(1);
80			}
81			if (nspec != nspec_grid+1) {
82			fprintf(stderr,
83			"%s: number of spectral samples cannot be changed\n",
84			progname);
85			exit(1);
86			}
87			return;
88			mem_error:
89			fprintf(stderr, "%s: out of memory in set_spectral_samples()\n",
90			progname);
91			exit(1);
92			}
93	greg	2.1
94			/* Start new DSF input grid */
95			void
96			new_bsdf_data(double new_theta, double new_phi)
97			{
98			if (!new_input_direction(new_theta, new_phi))
99			exit(1);
100			memset(dsf_grid, 0, sizeof(dsf_grid));
101	greg	2.27	if (nspec_grid > 0)
102			memset(spec_grid, 0, sizeof(float)GRIDRESGRIDRES*nspec_grid);
103	greg	2.1	}
104
105			/* Add BSDF data point */
106			void
107	greg	2.27	add_bsdf_data(double theta_out, double phi_out, const double val[], int isDSF)
108	greg	2.1	{
109			FVECT ovec;
110	greg	2.27	double cfact, Yval;
111	greg	2.1	int pos[2];
112
113	greg	2.27	if (nspec_grid > 2) {
114			fprintf(stderr, "%s: unsupported color space\n", progname);
115			exit(1);
116			}
117	greg	2.1	if (!output_orient) /* check output orientation */
118			output_orient = 1 - 2*(theta_out > 90.);
119			else if (output_orient > 0 ^ theta_out < 90.) {
120	greg	2.27	fprintf(stderr,
121			"%s: cannot handle output angles on both sides of surface\n",
122			progname);
123	greg	2.1	exit(1);
124			}
125			ovec[2] = sin((M_PI/180.)*theta_out);
126			ovec[0] = cos((M_PI/180.)phi_out) ovec[2];
127			ovec[1] = sin((M_PI/180.)phi_out) ovec[2];
128			ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
129	greg	2.27	/* BSDF to DSF correction */
130			cfact = isDSF ? 1. : ovec[2];
131	greg	2.1
132	greg	2.27	Yval = cfact * val[rbf_colorimetry==RBCtristimulus];
133	greg	2.4	/* update BSDF histogram */
134	greg	2.27	if (BSDF2SMLovec[2] < Yval && Yval < BSDF2BIGovec[2])
135			++bsdf_hist[histndx(Yval/ovec[2])];
136	greg	2.4
137	greg	2.1	pos_from_vec(pos, ovec);
138
139	greg	2.27	dsf_grid[pos[0]][pos[1]].sum.v += Yval;
140	greg	2.20	dsf_grid[pos[0]][pos[1]].sum.n++;
141	greg	2.27	/* add in X and Z values */
142			if (rbf_colorimetry == RBCtristimulus) {
143	greg	2.28	spec_grid[0][pos[0]][pos[1]] += cfact * val[0];
144			spec_grid[1][pos[0]][pos[1]] += cfact * val[2];
145	greg	2.27	}
146	greg	2.1	}
147
148	greg	2.11	/* Compute minimum BSDF from histogram (does not clear) */
149	greg	2.5	static void
150			comp_bsdf_min()
151			{
152	greg	2.17	unsigned long cnt, target;
153			int i;
154	greg	2.5
155			cnt = 0;
156			for (i = HISTLEN; i--; )
157			cnt += bsdf_hist[i];
158			if (!cnt) { /* shouldn't happen */
159			bsdf_min = 0;
160			return;
161			}
162			target = cnt/100; /* ignore bottom 1% */
163			cnt = 0;
164			for (i = 0; cnt <= target; i++)
165			cnt += bsdf_hist[i];
166			bsdf_min = histval(i-1);
167			}
168
169	greg	2.12	/* Determine if the given region is empty of grid samples */
170	greg	2.6	static int
171	greg	2.12	empty_region(int x0, int x1, int y0, int y1)
172	greg	2.6	{
173	greg	2.12	int x, y;
174
175			for (x = x0; x < x1; x++)
176			for (y = y0; y < y1; y++)
177	greg	2.20	if (dsf_grid[x][y].sum.n)
178	greg	2.12	return(0);
179			return(1);
180			}
181
182			/* Determine if the given region is smooth enough to be a single lobe */
183			static int
184			smooth_region(int x0, int x1, int y0, int y1)
185			{
186			RREAL rMtx[3][3];
187			FVECT xvec;
188			double A, B, C, nvs, sqerr;
189			int x, y, n;
190			/* compute planar regression */
191			memset(rMtx, 0, sizeof(rMtx));
192			memset(xvec, 0, sizeof(xvec));
193			for (x = x0; x < x1; x++)
194			for (y = y0; y < y1; y++)
195	greg	2.20	if ((n = dsf_grid[x][y].sum.n) > 0) {
196			double z = dsf_grid[x][y].sum.v;
197	greg	2.13	rMtx[0][0] += xx(double)n;
198			rMtx[0][1] += xy(double)n;
199			rMtx[0][2] += x*(double)n;
200			rMtx[1][1] += yy(double)n;
201			rMtx[1][2] += y*(double)n;
202			rMtx[2][2] += (double)n;
203	greg	2.12	xvec[0] += x*z;
204			xvec[1] += y*z;
205			xvec[2] += z;
206			}
207			rMtx[1][0] = rMtx[0][1];
208	greg	2.15	rMtx[2][0] = rMtx[0][2];
209	greg	2.12	rMtx[2][1] = rMtx[1][2];
210			nvs = rMtx[2][2];
211			if (SDinvXform(rMtx, rMtx) != SDEnone)
212	greg	2.16	return(1); /* colinear values */
213	greg	2.12	A = DOT(rMtx[0], xvec);
214			B = DOT(rMtx[1], xvec);
215	greg	2.29	if (AA + BB > MAXSLOPEMAXSLOPE) / too steep? */
216			return(0);
217	greg	2.12	C = DOT(rMtx[2], xvec);
218			sqerr = 0.0; /* compute mean squared error */
219			for (x = x0; x < x1; x++)
220			for (y = y0; y < y1; y++)
221	greg	2.20	if ((n = dsf_grid[x][y].sum.n) > 0) {
222			double d = Ax + By + C - dsf_grid[x][y].sum.v/n;
223	greg	2.12	sqerr += ndd;
224	greg	2.6	}
225	greg	2.12	if (sqerr <= nvsSMOOTH_MSE) / below absolute MSE threshold? */
226			return(1);
227	greg	2.13	/* OR below relative MSE threshold? */
228	greg	2.12	return(sqerrnvs <= xvec[2]xvec[2]*SMOOTH_MSER);
229			}
230
231			/* Create new lobe based on integrated samples in region */
232	greg	2.21	static int
233	greg	2.12	create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
234			{
235			double vtot = 0.0;
236	greg	2.27	double CIEXtot = 0.0, CIEZtot = 0.0;
237	greg	2.12	int nv = 0;
238	greg	2.25	double wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
239	greg	2.12	double rad;
240			int x, y;
241			/* compute average for region */
242			for (x = x0; x < x1; x++)
243	greg	2.25	for (y = y0; y < y1; y++)
244			if (dsf_grid[x][y].sum.n) {
245			const double v = dsf_grid[x][y].sum.v;
246			const int n = dsf_grid[x][y].sum.n;
247
248			if (v > 0) {
249	greg	2.27	const double wt = v / (double)n;
250	greg	2.25	wxsum += wt * x;
251			wysum += wt * y;
252			wtsum += wt;
253			}
254			vtot += v;
255			nv += n;
256	greg	2.27	if (rbf_colorimetry == RBCtristimulus) {
257			CIEXtot += spec_grid[0][x][y];
258			CIEZtot += spec_grid[1][x][y];
259			}
260	greg	2.25	}
261	greg	2.12	if (!nv) {
262			fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
263			progname);
264			exit(1);
265	greg	2.6	}
266	greg	2.21	if (vtot <= 0) /* only create positive lobes */
267			return(0);
268	greg	2.27	/* assign color */
269			if (rbf_colorimetry == RBCtristimulus) {
270			const double df = 1.0 / (CIEXtot + vtot + CIEZtot);
271			C_COLOR cclr;
272			c_cset(&cclr, CIEXtotdf, vtotdf);
273			rvp->chroma = c_encodeChroma(&cclr);
274			} else
275			rvp->chroma = c_dfchroma;
276	greg	2.12	/* peak value based on integral */
277			vtot = (x1-x0)(y1-y0)(2.M_PI/GRIDRES/GRIDRES)/(double)nv;
278			rad = (RSCA/(double)GRIDRES)*(x1-x0);
279			rvp->peak = vtot / ((2.M_PI) rad*rad);
280	greg	2.25	rvp->crad = ANG2R(rad); /* put peak at centroid */
281			rvp->gx = (int)(wxsum/wtsum + .5);
282			rvp->gy = (int)(wysum/wtsum + .5);
283	greg	2.21	return(1);
284	greg	2.6	}
285
286	greg	2.12	/* Recursive function to build radial basis function representation */
287	greg	2.6	static int
288	greg	2.12	build_rbfrep(RBFVAL *arp, int np, int x0, int x1, int y0, int y1)
289	greg	2.6	{
290	greg	2.12	int xmid = (x0+x1)>>1;
291			int ymid = (y0+y1)>>1;
292			int branched[4];
293			int nadded, nleaves;
294			/* need to make this a leaf? */
295			if (empty_region(x0, xmid, y0, ymid) \|\|
296			empty_region(xmid, x1, y0, ymid) \|\|
297			empty_region(x0, xmid, ymid, y1) \|\|
298			empty_region(xmid, x1, ymid, y1))
299			return(0);
300			/* add children (branches+leaves) */
301			if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
302			return(-1);
303			if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
304			return(-1);
305			if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
306			return(-1);
307			if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
308			return(-1);
309			nadded = branched[0] + branched[1] + branched[2] + branched[3];
310			nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
311			if (!nleaves) /* nothing but branches? */
312			return(nadded);
313			/* combine 4 leaves into 1? */
314	greg	2.14	if ((nleaves == 4) & (x1-x0 <= MAX_RAD) &&
315			smooth_region(x0, x1, y0, y1))
316	greg	2.12	return(0);
317			/* need more array space? */
318			if ((np+nleaves-1)>>RBFALLOCB != (np-1)>>RBFALLOCB) {
319			arp = (RBFVAL )realloc(*arp,
320			sizeof(RBFVAL)(np+nleaves-1+(1<<RBFALLOCB)));
321			if (*arp == NULL)
322			return(-1);
323	greg	2.6	}
324	greg	2.12	/* create lobes for leaves */
325	greg	2.21	if (!branched[0] && create_lobe(arp+np, x0, xmid, y0, ymid)) {
326			++(*np); ++nadded;
327			}
328			if (!branched[1] && create_lobe(arp+np, xmid, x1, y0, ymid)) {
329			++(*np); ++nadded;
330			}
331			if (!branched[2] && create_lobe(arp+np, x0, xmid, ymid, y1)) {
332			++(*np); ++nadded;
333			}
334			if (!branched[3] && create_lobe(arp+np, xmid, x1, ymid, y1)) {
335			++(*np); ++nadded;
336			}
337	greg	2.12	return(nadded);
338	greg	2.6	}
339
340	greg	2.1	/* Count up filled nodes and build RBF representation from current grid */
341			RBFNODE *
342	greg	2.12	make_rbfrep()
343	greg	2.1	{
344	greg	2.12	RBFNODE *newnode;
345			RBFVAL *rbfarr;
346	greg	2.1	int nn;
347	greg	2.5	/* compute minimum BSDF */
348			comp_bsdf_min();
349	greg	2.12	/* create RBF node list */
350			rbfarr = NULL; nn = 0;
351	greg	2.22	if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0) {
352			if (nn)
353			goto memerr;
354	greg	2.24	fprintf(stderr,
355			"%s: warning - skipping bad incidence (%.1f,%.1f)\n",
356			progname, theta_in_deg, phi_in_deg);
357			return(NULL);
358	greg	2.22	}
359	greg	2.12	/* (re)allocate RBF array */
360			newnode = (RBFNODE *)realloc(rbfarr,
361			sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
362	greg	2.2	if (newnode == NULL)
363			goto memerr;
364	greg	2.12	/* copy computed lobes into RBF node */
365			memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
366	greg	2.1	newnode->ord = -1;
367			newnode->next = NULL;
368			newnode->ejl = NULL;
369			newnode->invec[2] = sin((M_PI/180.)*theta_in_deg);
370			newnode->invec[0] = cos((M_PI/180.)phi_in_deg)newnode->invec[2];
371			newnode->invec[1] = sin((M_PI/180.)phi_in_deg)newnode->invec[2];
372			newnode->invec[2] = input_orientsqrt(1. - newnode->invec[2]newnode->invec[2]);
373	greg	2.12	newnode->vtotal = .0;
374	greg	2.1	newnode->nrbf = nn;
375	greg	2.12	/* compute sum for normalization */
376			while (nn-- > 0)
377			newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
378	greg	2.3	#ifdef DEBUG
379	greg	2.12	fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
380	greg	2.3	fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
381			get_theta180(newnode->invec), get_phi360(newnode->invec),
382			newnode->vtotal);
383			#endif
384	greg	2.1	insert_dsf(newnode);
385			return(newnode);
386	greg	2.2	memerr:
387			fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
388			exit(1);
389	greg	2.1	}