src/cv/bsdfrep.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrep.c,v 2.12 2013/03/20 01:00:22 greg Exp $";
#endif
/*
 * Support BSDF representation as radial basis functions.
 *
 *      G. Ward
 */

#define _USE_MATH_DEFINES
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "rtio.h"
#include "resolu.h"
#include "bsdfrep.h"
                                /* active grid resolution */
int                     grid_res = GRIDRES;

                                /* coverage/symmetry using INP_QUAD? flags */
int                     inp_coverage = 0;
                                /* all incident angles in-plane so far? */
int                     single_plane_incident = -1;

                                /* input/output orientations */
int                     input_orient = 0;
int                     output_orient = 0;

                                /* BSDF histogram */
int                     bsdf_hist[HISTLEN];

                                /* BSDF value for boundary regions */
double                  bsdf_min = 0;

                                /* processed incident DSF measurements */
RBFNODE                 *dsf_list = NULL;

                                /* RBF-linking matrices (edges) */
MIGRATION               *mig_list = NULL;

                                /* current input direction */
double                  theta_in_deg, phi_in_deg;

/* Register new input direction */
int
new_input_direction(double new_theta, double new_phi)
{
        if (!input_orient)              /* check input orientation */
                input_orient = 1 - 2*(new_theta > 90.);
        else if (input_orient > 0 ^ new_theta < 90.) {
                fprintf(stderr,
                "%s: Cannot handle input angles on both sides of surface\n",
                                progname);
                return(0);
        }
                                        /* normalize angle ranges */
        while (new_theta < -180.)
                new_theta += 360.;
        while (new_theta > 180.)
                new_theta -= 360.;
        if (new_theta < 0) {
                new_theta = -new_theta;
                new_phi += 180.;
        }
        if ((theta_in_deg = new_theta) < 1.0)
                return(1);              /* don't rely on phi near normal */
        while (new_phi < 0)
                new_phi += 360.;
        while (new_phi >= 360.)
                new_phi -= 360.;
        if (single_plane_incident > 0)  /* check input coverage */
                single_plane_incident = (round(new_phi) == round(phi_in_deg));
        else if (single_plane_incident < 0)
                single_plane_incident = 1;
        phi_in_deg = new_phi;
        if ((1. < new_phi) & (new_phi < 89.))
                inp_coverage |= INP_QUAD1;
        else if ((91. < new_phi) & (new_phi < 179.))
                inp_coverage |= INP_QUAD2;
        else if ((181. < new_phi) & (new_phi < 269.))
                inp_coverage |= INP_QUAD3;
        else if ((271. < new_phi) & (new_phi < 359.))
                inp_coverage |= INP_QUAD4;
        return(1);
}

/* Apply symmetry to the given vector based on distribution */
int
use_symmetry(FVECT vec)
{
        const double    phi = get_phi360(vec);

        switch (inp_coverage) {
        case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
                break;
        case INP_QUAD1|INP_QUAD2:
                if ((-FTINY > phi) | (phi > 180.+FTINY))
                        goto mir_y;
                break;
        case INP_QUAD2|INP_QUAD3:
                if ((90.-FTINY > phi) | (phi > 270.+FTINY))
                        goto mir_x;
                break;
        case INP_QUAD3|INP_QUAD4:
                if ((180.-FTINY > phi) | (phi > 360.+FTINY))
                        goto mir_y;
                break;
        case INP_QUAD4|INP_QUAD1:
                if ((270.-FTINY > phi) & (phi > 90.+FTINY))
                        goto mir_x;
                break;
        case INP_QUAD1:
                if ((-FTINY > phi) | (phi > 90.+FTINY))
                        switch ((int)(phi*(1./90.))) {
                        case 1: goto mir_x;
                        case 2: goto mir_xy;
                        case 3: goto mir_y;
                        }
                break;
        case INP_QUAD2:
                if ((90.-FTINY > phi) | (phi > 180.+FTINY))
                        switch ((int)(phi*(1./90.))) {
                        case 0: goto mir_x;
                        case 2: goto mir_y;
                        case 3: goto mir_xy;
                        }
                break;
        case INP_QUAD3:
                if ((180.-FTINY > phi) | (phi > 270.+FTINY))
                        switch ((int)(phi*(1./90.))) {
                        case 0: goto mir_xy;
                        case 1: goto mir_y;
                        case 3: goto mir_x;
                        }
                break;
        case INP_QUAD4:
                if ((270.-FTINY > phi) | (phi > 360.+FTINY))
                        switch ((int)(phi*(1./90.))) {
                        case 0: goto mir_y;
                        case 1: goto mir_xy;
                        case 2: goto mir_x;
                        }
                break;
        default:
                fprintf(stderr, "%s: Illegal input coverage (%d)\n",
                                        progname, inp_coverage);
                exit(1);
        }
        return(0);              /* in range */
mir_x:
        vec[0] = -vec[0];
        return(MIRROR_X);
mir_y:
        vec[1] = -vec[1];
        return(MIRROR_Y);
mir_xy:
        vec[0] = -vec[0];
        vec[1] = -vec[1];
        return(MIRROR_X|MIRROR_Y);
}

/* Reverse symmetry based on what was done before */
void
rev_symmetry(FVECT vec, int sym)
{
        if (sym & MIRROR_X)
                vec[0] = -vec[0];
        if (sym & MIRROR_Y)
                vec[1] = -vec[1];
}

/* Reverse symmetry for an RBF distribution */
void
rev_rbf_symmetry(RBFNODE *rbf, int sym)
{
        int     n;

        rev_symmetry(rbf->invec, sym);
        if (sym & MIRROR_X)
                for (n = rbf->nrbf; n-- > 0; )
                        rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
        if (sym & MIRROR_Y)
                for (n = rbf->nrbf; n-- > 0; )
                        rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
}

/* Rotate RBF to correspond to given incident vector */
void
rotate_rbf(RBFNODE *rbf, const FVECT invec)
{
        static const FVECT      vnorm = {.0, .0, 1.};
        const double            phi = atan2(invec[1],invec[0]) -
                                        atan2(rbf->invec[1],rbf->invec[0]);
        FVECT                   outvec;
        int                     pos[2];
        int                     n;

        for (n = ((-.01 > phi) | (phi > .01))*rbf->nrbf; n-- > 0; ) {
                ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
                spinvector(outvec, outvec, vnorm, phi);
                pos_from_vec(pos, outvec);
                rbf->rbfa[n].gx = pos[0];
                rbf->rbfa[n].gy = pos[1];
        }
        VCOPY(rbf->invec, invec);
}

/* Compute volume associated with Gaussian lobe */
double
rbf_volume(const RBFVAL *rbfp)
{
        double  rad = R2ANG(rbfp->crad);

        return((2.*M_PI) * rbfp->peak * rad*rad);
}

/* Compute outgoing vector from grid position */
void
ovec_from_pos(FVECT vec, int xpos, int ypos)
{
        double  uv[2];
        double  r2;
        
        SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
                                /* uniform hemispherical projection */
        r2 = uv[0]*uv[0] + uv[1]*uv[1];
        vec[0] = vec[1] = sqrt(2. - r2);
        vec[0] *= uv[0];
        vec[1] *= uv[1];
        vec[2] = output_orient*(1. - r2);
}

/* Compute grid position from normalized input/output vector */
void
pos_from_vec(int pos[2], const FVECT vec)
{
        double  sq[2];          /* uniform hemispherical projection */
        double  norm = 1./sqrt(1. + fabs(vec[2]));

        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);

        pos[0] = (int)(sq[0]*grid_res);
        pos[1] = (int)(sq[1]*grid_res);
}

/* Evaluate RBF for DSF at the given normalized outgoing direction */
double
eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
{
        double          minval = bsdf_min*output_orient*outvec[2];
        double          res = 0;
        const RBFVAL    *rbfp;
        FVECT           odir;
        double          sig2;
        int             n;
                                /* use minimum if no information avail. */
        if (rp == NULL) {
                if (outvec[2] > 0 ^ output_orient > 0)
                        return(.0);
                return(minval);
        }
        rbfp = rp->rbfa;
        for (n = rp->nrbf; n--; rbfp++) {
                ovec_from_pos(odir, rbfp->gx, rbfp->gy);
                sig2 = R2ANG(rbfp->crad);
                sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
                if (sig2 > -19.)
                        res += rbfp->peak * exp(sig2);
        }
        if (res < minval)       /* never return less than minval */
                return(minval);
        return(res);
}

/* Insert a new directional scattering function in our global list */
int
insert_dsf(RBFNODE *newrbf)
{
        RBFNODE         *rbf, *rbf_last;
        int             pos;
                                        /* check for redundant meas. */
        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
                if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
                        fprintf(stderr,
                                "%s: Duplicate incident measurement (ignored)\n",
                                        progname);
                        free(newrbf);
                        return(-1);
                }
                                        /* keep in ascending theta order */
        for (rbf_last = NULL, rbf = dsf_list; rbf != NULL;
                                        rbf_last = rbf, rbf = rbf->next)
                if (single_plane_incident && input_orient*rbf->invec[2] <
                                                input_orient*newrbf->invec[2])
                        break;
        if (rbf_last == NULL) {         /* insert new node in list */
                newrbf->ord = 0;
                newrbf->next = dsf_list;
                dsf_list = newrbf;
        } else {
                newrbf->ord = rbf_last->ord + 1;
                newrbf->next = rbf;
                rbf_last->next = newrbf;
        }
        rbf_last = newrbf;
        while (rbf != NULL) {           /* update ordinal positions */
                rbf->ord = rbf_last->ord + 1;
                rbf_last = rbf;
                rbf = rbf->next;
        }
        return(newrbf->ord);
}

/* Get the DSF indicated by its ordinal position */
RBFNODE *
get_dsf(int ord)
{
        RBFNODE         *rbf;

        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
                if (rbf->ord == ord)
                        return(rbf);
        return(NULL);
}

/* Get triangle surface orientation (unnormalized) */
void
tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
{
        FVECT   v2minus1, v3minus2;

        VSUB(v2minus1, v2, v1);
        VSUB(v3minus2, v3, v2);
        VCROSS(vres, v2minus1, v3minus2);
}

/* Determine if vertex order is reversed (inward normal) */
int
is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
{
        FVECT   tor;

        tri_orient(tor, v1, v2, v3);

        return(DOT(tor, v2) < 0.);
}

/* Find vertices completing triangles on either side of the given edge */
int
get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
{
        const MIGRATION *ej1, *ej2;
        RBFNODE         *tv;

        rbfv[0] = rbfv[1] = NULL;
        if (mig == NULL)
                return(0);
        for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
                                ej1 = nextedge(mig->rbfv[0],ej1)) {
                if (ej1 == mig)
                        continue;
                tv = opp_rbf(mig->rbfv[0],ej1);
                for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
                        if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
                                rbfv[is_rev_tri(mig->rbfv[0]->invec,
                                                mig->rbfv[1]->invec,
                                                tv->invec)] = tv;
                                break;
                        }
        }
        return((rbfv[0] != NULL) + (rbfv[1] != NULL));
}

/* Clear our BSDF representation and free memory */
void
clear_bsdf_rep(void)
{
        while (mig_list != NULL) {
                MIGRATION       *mig = mig_list;
                mig_list = mig->next;
                free(mig);
        }
        while (dsf_list != NULL) {
                RBFNODE         *rbf = dsf_list;
                dsf_list = rbf->next;
                free(rbf);
        }
        inp_coverage = 0;
        single_plane_incident = -1;
        input_orient = output_orient = 0;
        grid_res = GRIDRES;
}

/* Write our BSDF mesh interpolant out to the given binary stream */
void
save_bsdf_rep(FILE *ofp)
{
        RBFNODE         *rbf;
        MIGRATION       *mig;
        int             i, n;
                                        /* finish header */
        fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
        fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
        fprintf(ofp, "GRIDRES=%d\n", grid_res);
        fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
        fputformat(BSDFREP_FMT, ofp);
        fputc('\n', ofp);
                                        /* write each DSF */
        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
                putint(rbf->ord, 4, ofp);
                putflt(rbf->invec[0], ofp);
                putflt(rbf->invec[1], ofp);
                putflt(rbf->invec[2], ofp);
                putflt(rbf->vtotal, ofp);
                putint(rbf->nrbf, 4, ofp);
                for (i = 0; i < rbf->nrbf; i++) {
                        putflt(rbf->rbfa[i].peak, ofp);
                        putint(rbf->rbfa[i].crad, 2, ofp);
                        putint(rbf->rbfa[i].gx, 1, ofp);
                        putint(rbf->rbfa[i].gy, 1, ofp);
                }
        }
        putint(-1, 4, ofp);             /* terminator */
                                        /* write each migration matrix */
        for (mig = mig_list; mig != NULL; mig = mig->next) {
                int     zerocnt = 0;
                putint(mig->rbfv[0]->ord, 4, ofp);
                putint(mig->rbfv[1]->ord, 4, ofp);
                                        /* write out as sparse data */
                n = mtx_nrows(mig) * mtx_ncols(mig);
                for (i = 0; i < n; i++) {
                        if (zerocnt == 0xff) {
                                putint(0xff, 1, ofp); zerocnt = 0;
                        }
                        if (mig->mtx[i] != 0) {
                                putint(zerocnt, 1, ofp); zerocnt = 0;
                                putflt(mig->mtx[i], ofp);
                        } else
                                ++zerocnt;
                }
                putint(zerocnt, 1, ofp);
        }
        putint(-1, 4, ofp);             /* terminator */
        putint(-1, 4, ofp);
        if (fflush(ofp) == EOF) {
                fprintf(stderr, "%s: error writing BSDF interpolant\n",
                                progname);
                exit(1);
        }
}

/* Check header line for critical information */
static int
headline(char *s, void *p)
{
        char    fmt[32];

        if (!strncmp(s, "SYMMETRY=", 9)) {
                inp_coverage = atoi(s+9);
                single_plane_incident = !inp_coverage;
                return(0);
        }
        if (!strncmp(s, "IO_SIDES=", 9)) {
                sscanf(s+9, "%d %d", &input_orient, &output_orient);
                return(0);
        }
        if (!strncmp(s, "GRIDRES=", 8)) {
                sscanf(s+8, "%d", &grid_res);
                return(0);
        }
        if (!strncmp(s, "BSDFMIN=", 8)) {
                sscanf(s+8, "%lf", &bsdf_min);
                return(0);
        }
        if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
                return(-1);
        return(0);
}

/* Read a BSDF mesh interpolant from the given binary stream */
int
load_bsdf_rep(FILE *ifp)
{
        RBFNODE         rbfh;
        int             from_ord, to_ord;
        int             i;

        clear_bsdf_rep();
        if (ifp == NULL)
                return(0);
        if (getheader(ifp, headline, NULL) < 0 || single_plane_incident < 0 |
                        !input_orient | !output_orient) {
                fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
                                progname);
                return(0);
        }
        rbfh.next = NULL;               /* read each DSF */
        rbfh.ejl = NULL;
        while ((rbfh.ord = getint(4, ifp)) >= 0) {
                RBFNODE         *newrbf;

                rbfh.invec[0] = getflt(ifp);
                rbfh.invec[1] = getflt(ifp);
                rbfh.invec[2] = getflt(ifp);
                if (normalize(rbfh.invec) == 0) {
                        fprintf(stderr, "%s: zero incident vector\n", progname);
                        return(0);
                }
                rbfh.vtotal = getflt(ifp);
                rbfh.nrbf = getint(4, ifp);
                newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
                                        sizeof(RBFVAL)*(rbfh.nrbf-1));
                if (newrbf == NULL)
                        goto memerr;
                memcpy(newrbf, &rbfh, sizeof(RBFNODE)-sizeof(RBFVAL));
                for (i = 0; i < rbfh.nrbf; i++) {
                        newrbf->rbfa[i].peak = getflt(ifp);
                        newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
                        newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
                        newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
                }
                if (feof(ifp))
                        goto badEOF;
                                        /* insert in global list */
                if (insert_dsf(newrbf) != rbfh.ord) {
                        fprintf(stderr, "%s: error adding DSF\n", progname);
                        return(0);
                }
        }
                                        /* read each migration matrix */
        while ((from_ord = getint(4, ifp)) >= 0 &&
                        (to_ord = getint(4, ifp)) >= 0) {
                RBFNODE         *from_rbf = get_dsf(from_ord);
                RBFNODE         *to_rbf = get_dsf(to_ord);
                MIGRATION       *newmig;
                int             n;

                if ((from_rbf == NULL) | (to_rbf == NULL)) {
                        fprintf(stderr,
                                "%s: bad DSF reference in migration edge\n",
                                        progname);
                        return(0);
                }
                n = from_rbf->nrbf * to_rbf->nrbf;
                newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
                                                sizeof(float)*(n-1));
                if (newmig == NULL)
                        goto memerr;
                newmig->rbfv[0] = from_rbf;
                newmig->rbfv[1] = to_rbf;
                memset(newmig->mtx, 0, sizeof(float)*n);
                for (i = 0; ; ) {       /* read sparse data */
                        int     zc = getint(1, ifp) & 0xff;
                        if ((i += zc) >= n)
                                break;
                        if (zc == 0xff)
                                continue;
                        newmig->mtx[i++] = getflt(ifp);
                }
                if (feof(ifp))
                        goto badEOF;
                                        /* insert in edge lists */
                newmig->enxt[0] = from_rbf->ejl;
                from_rbf->ejl = newmig;
                newmig->enxt[1] = to_rbf->ejl;
                to_rbf->ejl = newmig;
                                        /* push onto global list */
                newmig->next = mig_list;
                mig_list = newmig;
        }
        return(1);                      /* success! */
memerr:
        fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
        exit(1);
badEOF:
        fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
        return(0);
}
Revision:	2.13
Committed:	Sun Mar 24 17:22:23 2013 UTC (11 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.12:	+5 -2 lines
Log Message:	Make sure not to output less than diffuse minimum
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.13	static const char RCSid[] = "$Id: bsdfrep.c,v 2.12 2013/03/20 01:00:22 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Support BSDF representation as radial basis functions.
6			*
7			* G. Ward
8			*/
9
10			#define _USE_MATH_DEFINES
11			#include <stdlib.h>
12	greg	2.2	#include <string.h>
13	greg	2.1	#include <math.h>
14			#include "rtio.h"
15			#include "resolu.h"
16			#include "bsdfrep.h"
17	greg	2.5	/* active grid resolution */
18			int grid_res = GRIDRES;
19
20	greg	2.4	/* coverage/symmetry using INP_QUAD? flags */
21	greg	2.1	int inp_coverage = 0;
22			/* all incident angles in-plane so far? */
23			int single_plane_incident = -1;
24
25			/* input/output orientations */
26			int input_orient = 0;
27			int output_orient = 0;
28
29	greg	2.12	/* BSDF histogram */
30			int bsdf_hist[HISTLEN];
31
32			/* BSDF value for boundary regions */
33			double bsdf_min = 0;
34
35	greg	2.1	/* processed incident DSF measurements */
36			RBFNODE *dsf_list = NULL;
37
38			/* RBF-linking matrices (edges) */
39			MIGRATION *mig_list = NULL;
40
41			/* current input direction */
42			double theta_in_deg, phi_in_deg;
43
44			/* Register new input direction */
45			int
46			new_input_direction(double new_theta, double new_phi)
47			{
48			if (!input_orient) /* check input orientation */
49			input_orient = 1 - 2*(new_theta > 90.);
50			else if (input_orient > 0 ^ new_theta < 90.) {
51			fprintf(stderr,
52			"%s: Cannot handle input angles on both sides of surface\n",
53			progname);
54			return(0);
55			}
56			/* normalize angle ranges */
57			while (new_theta < -180.)
58			new_theta += 360.;
59			while (new_theta > 180.)
60			new_theta -= 360.;
61			if (new_theta < 0) {
62			new_theta = -new_theta;
63			new_phi += 180.;
64			}
65	greg	2.6	if ((theta_in_deg = new_theta) < 1.0)
66			return(1); /* don't rely on phi near normal */
67	greg	2.1	while (new_phi < 0)
68			new_phi += 360.;
69			while (new_phi >= 360.)
70			new_phi -= 360.;
71			if (single_plane_incident > 0) /* check input coverage */
72			single_plane_incident = (round(new_phi) == round(phi_in_deg));
73			else if (single_plane_incident < 0)
74			single_plane_incident = 1;
75			phi_in_deg = new_phi;
76			if ((1. < new_phi) & (new_phi < 89.))
77			inp_coverage \|= INP_QUAD1;
78			else if ((91. < new_phi) & (new_phi < 179.))
79			inp_coverage \|= INP_QUAD2;
80			else if ((181. < new_phi) & (new_phi < 269.))
81			inp_coverage \|= INP_QUAD3;
82			else if ((271. < new_phi) & (new_phi < 359.))
83			inp_coverage \|= INP_QUAD4;
84			return(1);
85			}
86
87			/* Apply symmetry to the given vector based on distribution */
88			int
89			use_symmetry(FVECT vec)
90			{
91	greg	2.11	const double phi = get_phi360(vec);
92	greg	2.1
93			switch (inp_coverage) {
94			case INP_QUAD1\|INP_QUAD2\|INP_QUAD3\|INP_QUAD4:
95			break;
96			case INP_QUAD1\|INP_QUAD2:
97			if ((-FTINY > phi) \| (phi > 180.+FTINY))
98			goto mir_y;
99			break;
100			case INP_QUAD2\|INP_QUAD3:
101			if ((90.-FTINY > phi) \| (phi > 270.+FTINY))
102			goto mir_x;
103			break;
104			case INP_QUAD3\|INP_QUAD4:
105			if ((180.-FTINY > phi) \| (phi > 360.+FTINY))
106			goto mir_y;
107			break;
108			case INP_QUAD4\|INP_QUAD1:
109			if ((270.-FTINY > phi) & (phi > 90.+FTINY))
110			goto mir_x;
111			break;
112			case INP_QUAD1:
113			if ((-FTINY > phi) \| (phi > 90.+FTINY))
114			switch ((int)(phi*(1./90.))) {
115			case 1: goto mir_x;
116			case 2: goto mir_xy;
117			case 3: goto mir_y;
118			}
119			break;
120			case INP_QUAD2:
121			if ((90.-FTINY > phi) \| (phi > 180.+FTINY))
122			switch ((int)(phi*(1./90.))) {
123			case 0: goto mir_x;
124			case 2: goto mir_y;
125			case 3: goto mir_xy;
126			}
127			break;
128			case INP_QUAD3:
129			if ((180.-FTINY > phi) \| (phi > 270.+FTINY))
130			switch ((int)(phi*(1./90.))) {
131			case 0: goto mir_xy;
132			case 1: goto mir_y;
133			case 3: goto mir_x;
134			}
135			break;
136			case INP_QUAD4:
137			if ((270.-FTINY > phi) \| (phi > 360.+FTINY))
138			switch ((int)(phi*(1./90.))) {
139			case 0: goto mir_y;
140			case 1: goto mir_xy;
141			case 2: goto mir_x;
142			}
143			break;
144			default:
145			fprintf(stderr, "%s: Illegal input coverage (%d)\n",
146			progname, inp_coverage);
147			exit(1);
148			}
149			return(0); /* in range */
150			mir_x:
151			vec[0] = -vec[0];
152			return(MIRROR_X);
153			mir_y:
154			vec[1] = -vec[1];
155			return(MIRROR_Y);
156			mir_xy:
157			vec[0] = -vec[0];
158			vec[1] = -vec[1];
159			return(MIRROR_X\|MIRROR_Y);
160			}
161
162			/* Reverse symmetry based on what was done before */
163			void
164			rev_symmetry(FVECT vec, int sym)
165			{
166			if (sym & MIRROR_X)
167			vec[0] = -vec[0];
168			if (sym & MIRROR_Y)
169			vec[1] = -vec[1];
170			}
171
172			/* Reverse symmetry for an RBF distribution */
173			void
174			rev_rbf_symmetry(RBFNODE *rbf, int sym)
175			{
176			int n;
177
178			rev_symmetry(rbf->invec, sym);
179			if (sym & MIRROR_X)
180			for (n = rbf->nrbf; n-- > 0; )
181	greg	2.5	rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
182	greg	2.1	if (sym & MIRROR_Y)
183			for (n = rbf->nrbf; n-- > 0; )
184	greg	2.5	rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
185	greg	2.1	}
186
187	greg	2.6	/* Rotate RBF to correspond to given incident vector */
188			void
189			rotate_rbf(RBFNODE *rbf, const FVECT invec)
190			{
191			static const FVECT vnorm = {.0, .0, 1.};
192			const double phi = atan2(invec[1],invec[0]) -
193			atan2(rbf->invec[1],rbf->invec[0]);
194			FVECT outvec;
195			int pos[2];
196			int n;
197	greg	2.8
198	greg	2.10	for (n = ((-.01 > phi) \| (phi > .01))*rbf->nrbf; n-- > 0; ) {
199	greg	2.6	ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
200			spinvector(outvec, outvec, vnorm, phi);
201			pos_from_vec(pos, outvec);
202			rbf->rbfa[n].gx = pos[0];
203			rbf->rbfa[n].gy = pos[1];
204			}
205			VCOPY(rbf->invec, invec);
206			}
207
208	greg	2.1	/* Compute volume associated with Gaussian lobe */
209			double
210			rbf_volume(const RBFVAL *rbfp)
211			{
212			double rad = R2ANG(rbfp->crad);
213
214			return((2.M_PI) rbfp->peak * rad*rad);
215			}
216
217			/* Compute outgoing vector from grid position */
218			void
219			ovec_from_pos(FVECT vec, int xpos, int ypos)
220			{
221			double uv[2];
222			double r2;
223
224	greg	2.8	SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
225	greg	2.1	/* uniform hemispherical projection */
226			r2 = uv[0]uv[0] + uv[1]uv[1];
227			vec[0] = vec[1] = sqrt(2. - r2);
228			vec[0] *= uv[0];
229			vec[1] *= uv[1];
230			vec[2] = output_orient*(1. - r2);
231			}
232
233			/* Compute grid position from normalized input/output vector */
234			void
235			pos_from_vec(int pos[2], const FVECT vec)
236			{
237			double sq[2]; /* uniform hemispherical projection */
238			double norm = 1./sqrt(1. + fabs(vec[2]));
239
240			SDdisk2square(sq, vec[0]norm, vec[1]norm);
241
242	greg	2.5	pos[0] = (int)(sq[0]*grid_res);
243			pos[1] = (int)(sq[1]*grid_res);
244	greg	2.1	}
245
246			/* Evaluate RBF for DSF at the given normalized outgoing direction */
247			double
248			eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
249			{
250	greg	2.13	double minval = bsdf_minoutput_orientoutvec[2];
251	greg	2.8	double res = 0;
252	greg	2.1	const RBFVAL *rbfp;
253			FVECT odir;
254			double sig2;
255			int n;
256	greg	2.12	/* use minimum if no information avail. */
257			if (rp == NULL) {
258			if (outvec[2] > 0 ^ output_orient > 0)
259			return(.0);
260	greg	2.13	return(minval);
261	greg	2.12	}
262	greg	2.1	rbfp = rp->rbfa;
263			for (n = rp->nrbf; n--; rbfp++) {
264			ovec_from_pos(odir, rbfp->gx, rbfp->gy);
265			sig2 = R2ANG(rbfp->crad);
266			sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
267			if (sig2 > -19.)
268			res += rbfp->peak * exp(sig2);
269			}
270	greg	2.13	if (res < minval) /* never return less than minval */
271			return(minval);
272	greg	2.1	return(res);
273			}
274
275			/* Insert a new directional scattering function in our global list */
276			int
277			insert_dsf(RBFNODE *newrbf)
278			{
279			RBFNODE rbf, rbf_last;
280			int pos;
281			/* check for redundant meas. */
282			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
283			if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
284			fprintf(stderr,
285			"%s: Duplicate incident measurement (ignored)\n",
286			progname);
287			free(newrbf);
288			return(-1);
289			}
290			/* keep in ascending theta order */
291			for (rbf_last = NULL, rbf = dsf_list; rbf != NULL;
292			rbf_last = rbf, rbf = rbf->next)
293			if (single_plane_incident && input_orient*rbf->invec[2] <
294			input_orient*newrbf->invec[2])
295			break;
296			if (rbf_last == NULL) { /* insert new node in list */
297			newrbf->ord = 0;
298			newrbf->next = dsf_list;
299			dsf_list = newrbf;
300			} else {
301			newrbf->ord = rbf_last->ord + 1;
302			newrbf->next = rbf;
303			rbf_last->next = newrbf;
304			}
305			rbf_last = newrbf;
306			while (rbf != NULL) { /* update ordinal positions */
307			rbf->ord = rbf_last->ord + 1;
308			rbf_last = rbf;
309			rbf = rbf->next;
310			}
311			return(newrbf->ord);
312			}
313
314			/* Get the DSF indicated by its ordinal position */
315			RBFNODE *
316			get_dsf(int ord)
317			{
318			RBFNODE *rbf;
319
320			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
321	greg	2.3	if (rbf->ord == ord)
322	greg	2.1	return(rbf);
323			return(NULL);
324			}
325
326			/* Get triangle surface orientation (unnormalized) */
327			void
328			tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
329			{
330			FVECT v2minus1, v3minus2;
331
332			VSUB(v2minus1, v2, v1);
333			VSUB(v3minus2, v3, v2);
334			VCROSS(vres, v2minus1, v3minus2);
335			}
336
337			/* Determine if vertex order is reversed (inward normal) */
338			int
339			is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
340			{
341			FVECT tor;
342
343			tri_orient(tor, v1, v2, v3);
344
345			return(DOT(tor, v2) < 0.);
346			}
347
348			/* Find vertices completing triangles on either side of the given edge */
349			int
350			get_triangles(RBFNODE rbfv[2], const MIGRATION mig)
351			{
352	greg	2.4	const MIGRATION ej1, ej2;
353	greg	2.1	RBFNODE *tv;
354
355			rbfv[0] = rbfv[1] = NULL;
356			if (mig == NULL)
357			return(0);
358	greg	2.4	for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
359			ej1 = nextedge(mig->rbfv[0],ej1)) {
360			if (ej1 == mig)
361	greg	2.1	continue;
362	greg	2.4	tv = opp_rbf(mig->rbfv[0],ej1);
363	greg	2.1	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
364			if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
365			rbfv[is_rev_tri(mig->rbfv[0]->invec,
366			mig->rbfv[1]->invec,
367			tv->invec)] = tv;
368			break;
369			}
370			}
371			return((rbfv[0] != NULL) + (rbfv[1] != NULL));
372			}
373
374	greg	2.4	/* Clear our BSDF representation and free memory */
375			void
376			clear_bsdf_rep(void)
377			{
378			while (mig_list != NULL) {
379			MIGRATION *mig = mig_list;
380			mig_list = mig->next;
381			free(mig);
382			}
383			while (dsf_list != NULL) {
384			RBFNODE *rbf = dsf_list;
385			dsf_list = rbf->next;
386			free(rbf);
387			}
388			inp_coverage = 0;
389			single_plane_incident = -1;
390			input_orient = output_orient = 0;
391	greg	2.5	grid_res = GRIDRES;
392	greg	2.4	}
393
394	greg	2.1	/* Write our BSDF mesh interpolant out to the given binary stream */
395			void
396			save_bsdf_rep(FILE *ofp)
397			{
398			RBFNODE *rbf;
399			MIGRATION *mig;
400			int i, n;
401			/* finish header */
402	greg	2.2	fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
403			fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
404	greg	2.5	fprintf(ofp, "GRIDRES=%d\n", grid_res);
405	greg	2.12	fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
406	greg	2.1	fputformat(BSDFREP_FMT, ofp);
407			fputc('\n', ofp);
408			/* write each DSF */
409			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
410			putint(rbf->ord, 4, ofp);
411			putflt(rbf->invec[0], ofp);
412			putflt(rbf->invec[1], ofp);
413			putflt(rbf->invec[2], ofp);
414			putflt(rbf->vtotal, ofp);
415			putint(rbf->nrbf, 4, ofp);
416			for (i = 0; i < rbf->nrbf; i++) {
417			putflt(rbf->rbfa[i].peak, ofp);
418			putint(rbf->rbfa[i].crad, 2, ofp);
419			putint(rbf->rbfa[i].gx, 1, ofp);
420			putint(rbf->rbfa[i].gy, 1, ofp);
421			}
422			}
423			putint(-1, 4, ofp); /* terminator */
424			/* write each migration matrix */
425	greg	2.2	for (mig = mig_list; mig != NULL; mig = mig->next) {
426			int zerocnt = 0;
427	greg	2.1	putint(mig->rbfv[0]->ord, 4, ofp);
428			putint(mig->rbfv[1]->ord, 4, ofp);
429	greg	2.2	/* write out as sparse data */
430	greg	2.1	n = mtx_nrows(mig) * mtx_ncols(mig);
431	greg	2.2	for (i = 0; i < n; i++) {
432	greg	2.3	if (zerocnt == 0xff) {
433			putint(0xff, 1, ofp); zerocnt = 0;
434	greg	2.2	}
435			if (mig->mtx[i] != 0) {
436			putint(zerocnt, 1, ofp); zerocnt = 0;
437			putflt(mig->mtx[i], ofp);
438			} else
439			++zerocnt;
440			}
441			putint(zerocnt, 1, ofp);
442	greg	2.1	}
443			putint(-1, 4, ofp); /* terminator */
444			putint(-1, 4, ofp);
445			if (fflush(ofp) == EOF) {
446			fprintf(stderr, "%s: error writing BSDF interpolant\n",
447			progname);
448			exit(1);
449			}
450			}
451
452	greg	2.2	/* Check header line for critical information */
453			static int
454			headline(char s, void p)
455			{
456			char fmt[32];
457
458			if (!strncmp(s, "SYMMETRY=", 9)) {
459			inp_coverage = atoi(s+9);
460			single_plane_incident = !inp_coverage;
461			return(0);
462			}
463			if (!strncmp(s, "IO_SIDES=", 9)) {
464			sscanf(s+9, "%d %d", &input_orient, &output_orient);
465			return(0);
466			}
467	greg	2.5	if (!strncmp(s, "GRIDRES=", 8)) {
468			sscanf(s+8, "%d", &grid_res);
469			return(0);
470			}
471	greg	2.12	if (!strncmp(s, "BSDFMIN=", 8)) {
472			sscanf(s+8, "%lf", &bsdf_min);
473			return(0);
474			}
475	greg	2.2	if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
476			return(-1);
477			return(0);
478			}
479
480	greg	2.1	/* Read a BSDF mesh interpolant from the given binary stream */
481			int
482			load_bsdf_rep(FILE *ifp)
483			{
484			RBFNODE rbfh;
485			int from_ord, to_ord;
486			int i;
487	greg	2.4
488			clear_bsdf_rep();
489	greg	2.5	if (ifp == NULL)
490			return(0);
491	greg	2.2	if (getheader(ifp, headline, NULL) < 0 \|\| single_plane_incident < 0 \|
492			!input_orient \| !output_orient) {
493	greg	2.1	fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
494			progname);
495			return(0);
496			}
497			rbfh.next = NULL; /* read each DSF */
498			rbfh.ejl = NULL;
499			while ((rbfh.ord = getint(4, ifp)) >= 0) {
500			RBFNODE *newrbf;
501
502			rbfh.invec[0] = getflt(ifp);
503			rbfh.invec[1] = getflt(ifp);
504			rbfh.invec[2] = getflt(ifp);
505	greg	2.9	if (normalize(rbfh.invec) == 0) {
506			fprintf(stderr, "%s: zero incident vector\n", progname);
507			return(0);
508			}
509	greg	2.3	rbfh.vtotal = getflt(ifp);
510	greg	2.1	rbfh.nrbf = getint(4, ifp);
511			newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
512			sizeof(RBFVAL)*(rbfh.nrbf-1));
513			if (newrbf == NULL)
514			goto memerr;
515	greg	2.8	memcpy(newrbf, &rbfh, sizeof(RBFNODE)-sizeof(RBFVAL));
516	greg	2.1	for (i = 0; i < rbfh.nrbf; i++) {
517			newrbf->rbfa[i].peak = getflt(ifp);
518			newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
519			newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
520			newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
521			}
522			if (feof(ifp))
523			goto badEOF;
524			/* insert in global list */
525			if (insert_dsf(newrbf) != rbfh.ord) {
526			fprintf(stderr, "%s: error adding DSF\n", progname);
527			return(0);
528			}
529			}
530			/* read each migration matrix */
531			while ((from_ord = getint(4, ifp)) >= 0 &&
532			(to_ord = getint(4, ifp)) >= 0) {
533			RBFNODE *from_rbf = get_dsf(from_ord);
534			RBFNODE *to_rbf = get_dsf(to_ord);
535			MIGRATION *newmig;
536			int n;
537
538			if ((from_rbf == NULL) \| (to_rbf == NULL)) {
539			fprintf(stderr,
540			"%s: bad DSF reference in migration edge\n",
541			progname);
542			return(0);
543			}
544			n = from_rbf->nrbf * to_rbf->nrbf;
545			newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
546			sizeof(float)*(n-1));
547			if (newmig == NULL)
548			goto memerr;
549			newmig->rbfv[0] = from_rbf;
550			newmig->rbfv[1] = to_rbf;
551	greg	2.2	memset(newmig->mtx, 0, sizeof(float)*n);
552			for (i = 0; ; ) { /* read sparse data */
553			int zc = getint(1, ifp) & 0xff;
554			if ((i += zc) >= n)
555			break;
556	greg	2.3	if (zc == 0xff)
557			continue;
558	greg	2.2	newmig->mtx[i++] = getflt(ifp);
559			}
560	greg	2.1	if (feof(ifp))
561			goto badEOF;
562			/* insert in edge lists */
563			newmig->enxt[0] = from_rbf->ejl;
564			from_rbf->ejl = newmig;
565			newmig->enxt[1] = to_rbf->ejl;
566			to_rbf->ejl = newmig;
567			/* push onto global list */
568			newmig->next = mig_list;
569			mig_list = newmig;
570			}
571			return(1); /* success! */
572			memerr:
573			fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
574			exit(1);
575			badEOF:
576			fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
577			return(0);
578			}