src/cv/bsdfrep.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfrep.c,v 2.23 2014/03/24 06:07:46 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"
                                /* name and manufacturer if known */
char                    bsdf_name[256];
char                    bsdf_manuf[256];
                                /* 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 */
unsigned long           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 = (cos(phi) < 1.-FTINY)*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 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);
}

/* Compute volume associated with Gaussian lobe */
double
rbf_volume(const RBFVAL *rbfp)
{
        double  rad = R2ANG(rbfp->crad);
        FVECT   odir;
        double  elev, integ;
                                /* infinite integral approximation */
        integ = (2.*M_PI) * rbfp->peak * rad*rad;
                                /* check if we're near horizon */
        ovec_from_pos(odir, rbfp->gx, rbfp->gy);
        elev = output_orient*odir[2];
                                /* apply cut-off correction if > 1% */
        if (elev < 2.8*rad) {
                /* elev = asin(elev);   /* this is so crude, anyway... */
                integ *= 1. - .5*exp(-.5*elev*elev/(rad*rad));
        }
        return(integ);
}

/* Evaluate RBF for DSF at the given normalized outgoing direction */
double
eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
{
        const double    rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
        double          minval = bsdf_min*output_orient*outvec[2];
        int             pos[2];
        double          res = 0;
        const RBFVAL    *rbfp;
        FVECT           odir;
        double          rad2;
        int             n;
                                /* check for wrong side */
        if (outvec[2] > 0 ^ output_orient > 0)
                return(.0);
                                /* use minimum if no information avail. */
        if (rp == NULL)
                return(minval);
                                /* optimization for fast lobe culling */
        pos_from_vec(pos, outvec);
                                /* sum radial basis function */
        rbfp = rp->rbfa;
        for (n = rp->nrbf; n--; rbfp++) {
                int     d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
                                (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
                rad2 = R2ANG(rbfp->crad);
                rad2 *= rad2;
                if (d2 > rad2*rfact2)
                        continue;
                ovec_from_pos(odir, rbfp->gx, rbfp->gy);
                res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
        }
        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 at (%.1f,%.1f)\n",
                                        progname, get_theta180(newrbf->invec),
                                        get_phi360(newrbf->invec));
                        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));
}

/* Advect and allocate new RBF along edge (internal call) */
RBFNODE *
e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
{
        double          cthresh = FTINY;
        RBFNODE         *rbf;
        int             n, i, j;
        double          t, full_dist;
                                                /* get relative position */
        t = Acos(DOT(invec, mig->rbfv[0]->invec));
        if (t < M_PI/grid_res) {                /* near first DSF */
                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
                rbf = (RBFNODE *)malloc(n);
                if (rbf == NULL)
                        goto memerr;
                memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
                rbf->next = NULL; rbf->ejl = NULL;
                return(rbf);
        }
        full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
                rbf = (RBFNODE *)malloc(n);
                if (rbf == NULL)
                        goto memerr;
                memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
                rbf->next = NULL; rbf->ejl = NULL;
                return(rbf);
        }
        t /= full_dist;
tryagain:
        n = 0;                                  /* count migrating particles */
        for (i = 0; i < mtx_nrows(mig); i++)
            for (j = 0; j < mtx_ncols(mig); j++)
                n += (mtx_coef(mig,i,j) > cthresh);
                                                /* are we over our limit? */
        if ((lobe_lim > 0) & (n > lobe_lim)) {
                cthresh = cthresh*2. + 10.*FTINY;
                goto tryagain;
        }
#ifdef DEBUG
        fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
                        mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
#endif
        rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
        if (rbf == NULL)
                goto memerr;
        rbf->next = NULL; rbf->ejl = NULL;
        VCOPY(rbf->invec, invec);
        rbf->nrbf = n;
        rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
        n = 0;                                  /* advect RBF lobes */
        for (i = 0; i < mtx_nrows(mig); i++) {
            const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
            const float         peak0 = rbf0i->peak;
            const double        rad0 = R2ANG(rbf0i->crad);
            FVECT               v0;
            float               mv;
            ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
            for (j = 0; j < mtx_ncols(mig); j++)
                if ((mv = mtx_coef(mig,i,j)) > cthresh) {
                        const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
                        double          rad2;
                        FVECT           v;
                        int             pos[2];
                        rad2 = R2ANG(rbf1j->crad);
                        rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
                                                rad0*rad0/rad2;
                        rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
                        ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
                        geodesic(v, v0, v, t, GEOD_REL);
                        pos_from_vec(pos, v);
                        rbf->rbfa[n].gx = pos[0];
                        rbf->rbfa[n].gy = pos[1];
                        ++n;
                }
        }
        rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
        return(rbf);
memerr:
        fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
        exit(1);
        return(NULL);   /* pro forma return */
}

/* 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);
        }
        bsdf_name[0] = '\0';
        bsdf_manuf[0] = '\0';
        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 */
        if (bsdf_name[0])
                fprintf(ofp, "NAME=%s\n", bsdf_name);
        if (bsdf_manuf[0])
                fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
        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, "NAME=", 5)) {
                strcpy(bsdf_name, s+5);
                bsdf_name[strlen(bsdf_name)-1] = '\0';
        }
        if (!strncmp(s, "MANUFACT=", 9)) {
                strcpy(bsdf_manuf, s+9);
                bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
        }
        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 |
                        (grid_res < 16) | (grid_res > 256)) {
                fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
                                progname);
                return(0);
        }
        memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
        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;
                *newrbf = rbfh;
                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.24
Committed:	Tue Mar 25 14:55:35 2014 UTC (10 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R2
Changes since 2.23:	+2 -2 lines
Log Message:	Made check against non-rotation more robust
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.24	static const char RCSid[] = "$Id: bsdfrep.c,v 2.23 2014/03/24 06:07:46 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.19	/* name and manufacturer if known */
18			char bsdf_name[256];
19			char bsdf_manuf[256];
20	greg	2.5	/* active grid resolution */
21			int grid_res = GRIDRES;
22
23	greg	2.4	/* coverage/symmetry using INP_QUAD? flags */
24	greg	2.1	int inp_coverage = 0;
25			/* all incident angles in-plane so far? */
26			int single_plane_incident = -1;
27
28			/* input/output orientations */
29			int input_orient = 0;
30			int output_orient = 0;
31
32	greg	2.12	/* BSDF histogram */
33	greg	2.15	unsigned long bsdf_hist[HISTLEN];
34	greg	2.12
35			/* BSDF value for boundary regions */
36			double bsdf_min = 0;
37
38	greg	2.1	/* processed incident DSF measurements */
39			RBFNODE *dsf_list = NULL;
40
41			/* RBF-linking matrices (edges) */
42			MIGRATION *mig_list = NULL;
43
44			/* current input direction */
45			double theta_in_deg, phi_in_deg;
46
47			/* Register new input direction */
48			int
49			new_input_direction(double new_theta, double new_phi)
50			{
51			if (!input_orient) /* check input orientation */
52			input_orient = 1 - 2*(new_theta > 90.);
53			else if (input_orient > 0 ^ new_theta < 90.) {
54			fprintf(stderr,
55			"%s: Cannot handle input angles on both sides of surface\n",
56			progname);
57			return(0);
58			}
59			/* normalize angle ranges */
60			while (new_theta < -180.)
61			new_theta += 360.;
62			while (new_theta > 180.)
63			new_theta -= 360.;
64			if (new_theta < 0) {
65			new_theta = -new_theta;
66			new_phi += 180.;
67			}
68	greg	2.6	if ((theta_in_deg = new_theta) < 1.0)
69			return(1); /* don't rely on phi near normal */
70	greg	2.1	while (new_phi < 0)
71			new_phi += 360.;
72			while (new_phi >= 360.)
73			new_phi -= 360.;
74			if (single_plane_incident > 0) /* check input coverage */
75			single_plane_incident = (round(new_phi) == round(phi_in_deg));
76			else if (single_plane_incident < 0)
77			single_plane_incident = 1;
78			phi_in_deg = new_phi;
79			if ((1. < new_phi) & (new_phi < 89.))
80			inp_coverage \|= INP_QUAD1;
81			else if ((91. < new_phi) & (new_phi < 179.))
82			inp_coverage \|= INP_QUAD2;
83			else if ((181. < new_phi) & (new_phi < 269.))
84			inp_coverage \|= INP_QUAD3;
85			else if ((271. < new_phi) & (new_phi < 359.))
86			inp_coverage \|= INP_QUAD4;
87			return(1);
88			}
89
90			/* Apply symmetry to the given vector based on distribution */
91			int
92			use_symmetry(FVECT vec)
93			{
94	greg	2.11	const double phi = get_phi360(vec);
95	greg	2.1
96			switch (inp_coverage) {
97			case INP_QUAD1\|INP_QUAD2\|INP_QUAD3\|INP_QUAD4:
98			break;
99			case INP_QUAD1\|INP_QUAD2:
100			if ((-FTINY > phi) \| (phi > 180.+FTINY))
101			goto mir_y;
102			break;
103			case INP_QUAD2\|INP_QUAD3:
104			if ((90.-FTINY > phi) \| (phi > 270.+FTINY))
105			goto mir_x;
106			break;
107			case INP_QUAD3\|INP_QUAD4:
108			if ((180.-FTINY > phi) \| (phi > 360.+FTINY))
109			goto mir_y;
110			break;
111			case INP_QUAD4\|INP_QUAD1:
112			if ((270.-FTINY > phi) & (phi > 90.+FTINY))
113			goto mir_x;
114			break;
115			case INP_QUAD1:
116			if ((-FTINY > phi) \| (phi > 90.+FTINY))
117			switch ((int)(phi*(1./90.))) {
118			case 1: goto mir_x;
119			case 2: goto mir_xy;
120			case 3: goto mir_y;
121			}
122			break;
123			case INP_QUAD2:
124			if ((90.-FTINY > phi) \| (phi > 180.+FTINY))
125			switch ((int)(phi*(1./90.))) {
126			case 0: goto mir_x;
127			case 2: goto mir_y;
128			case 3: goto mir_xy;
129			}
130			break;
131			case INP_QUAD3:
132			if ((180.-FTINY > phi) \| (phi > 270.+FTINY))
133			switch ((int)(phi*(1./90.))) {
134			case 0: goto mir_xy;
135			case 1: goto mir_y;
136			case 3: goto mir_x;
137			}
138			break;
139			case INP_QUAD4:
140			if ((270.-FTINY > phi) \| (phi > 360.+FTINY))
141			switch ((int)(phi*(1./90.))) {
142			case 0: goto mir_y;
143			case 1: goto mir_xy;
144			case 2: goto mir_x;
145			}
146			break;
147			default:
148			fprintf(stderr, "%s: Illegal input coverage (%d)\n",
149			progname, inp_coverage);
150			exit(1);
151			}
152			return(0); /* in range */
153			mir_x:
154			vec[0] = -vec[0];
155			return(MIRROR_X);
156			mir_y:
157			vec[1] = -vec[1];
158			return(MIRROR_Y);
159			mir_xy:
160			vec[0] = -vec[0];
161			vec[1] = -vec[1];
162			return(MIRROR_X\|MIRROR_Y);
163			}
164
165			/* Reverse symmetry based on what was done before */
166			void
167			rev_symmetry(FVECT vec, int sym)
168			{
169			if (sym & MIRROR_X)
170			vec[0] = -vec[0];
171			if (sym & MIRROR_Y)
172			vec[1] = -vec[1];
173			}
174
175			/* Reverse symmetry for an RBF distribution */
176			void
177			rev_rbf_symmetry(RBFNODE *rbf, int sym)
178			{
179			int n;
180
181			rev_symmetry(rbf->invec, sym);
182			if (sym & MIRROR_X)
183			for (n = rbf->nrbf; n-- > 0; )
184	greg	2.5	rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
185	greg	2.1	if (sym & MIRROR_Y)
186			for (n = rbf->nrbf; n-- > 0; )
187	greg	2.5	rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
188	greg	2.1	}
189
190	greg	2.6	/* Rotate RBF to correspond to given incident vector */
191			void
192			rotate_rbf(RBFNODE *rbf, const FVECT invec)
193			{
194			static const FVECT vnorm = {.0, .0, 1.};
195			const double phi = atan2(invec[1],invec[0]) -
196			atan2(rbf->invec[1],rbf->invec[0]);
197			FVECT outvec;
198			int pos[2];
199			int n;
200	greg	2.8
201	greg	2.24	for (n = (cos(phi) < 1.-FTINY)*rbf->nrbf; n-- > 0; ) {
202	greg	2.6	ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
203			spinvector(outvec, outvec, vnorm, phi);
204			pos_from_vec(pos, outvec);
205			rbf->rbfa[n].gx = pos[0];
206			rbf->rbfa[n].gy = pos[1];
207			}
208			VCOPY(rbf->invec, invec);
209			}
210
211	greg	2.1	/* Compute outgoing vector from grid position */
212			void
213			ovec_from_pos(FVECT vec, int xpos, int ypos)
214			{
215			double uv[2];
216			double r2;
217
218	greg	2.8	SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
219	greg	2.1	/* uniform hemispherical projection */
220			r2 = uv[0]uv[0] + uv[1]uv[1];
221			vec[0] = vec[1] = sqrt(2. - r2);
222			vec[0] *= uv[0];
223			vec[1] *= uv[1];
224			vec[2] = output_orient*(1. - r2);
225			}
226
227			/* Compute grid position from normalized input/output vector */
228			void
229			pos_from_vec(int pos[2], const FVECT vec)
230			{
231			double sq[2]; /* uniform hemispherical projection */
232			double norm = 1./sqrt(1. + fabs(vec[2]));
233
234			SDdisk2square(sq, vec[0]norm, vec[1]norm);
235
236	greg	2.5	pos[0] = (int)(sq[0]*grid_res);
237			pos[1] = (int)(sq[1]*grid_res);
238	greg	2.1	}
239
240	greg	2.14	/* Compute volume associated with Gaussian lobe */
241			double
242			rbf_volume(const RBFVAL *rbfp)
243			{
244			double rad = R2ANG(rbfp->crad);
245			FVECT odir;
246			double elev, integ;
247			/* infinite integral approximation */
248			integ = (2.M_PI) rbfp->peak * rad*rad;
249			/* check if we're near horizon */
250			ovec_from_pos(odir, rbfp->gx, rbfp->gy);
251			elev = output_orient*odir[2];
252			/* apply cut-off correction if > 1% */
253			if (elev < 2.8*rad) {
254			/* elev = asin(elev); /* this is so crude, anyway... */
255			integ = 1. - .5exp(-.5elevelev/(rad*rad));
256			}
257			return(integ);
258			}
259
260	greg	2.1	/* Evaluate RBF for DSF at the given normalized outgoing direction */
261			double
262			eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
263			{
264	greg	2.17	const double rfact2 = (38./M_PI/M_PI)(grid_resgrid_res);
265	greg	2.13	double minval = bsdf_minoutput_orientoutvec[2];
266	greg	2.16	int pos[2];
267	greg	2.8	double res = 0;
268	greg	2.1	const RBFVAL *rbfp;
269			FVECT odir;
270	greg	2.16	double rad2;
271	greg	2.1	int n;
272	greg	2.14	/* check for wrong side */
273			if (outvec[2] > 0 ^ output_orient > 0)
274			return(.0);
275	greg	2.12	/* use minimum if no information avail. */
276	greg	2.14	if (rp == NULL)
277	greg	2.13	return(minval);
278	greg	2.16	/* optimization for fast lobe culling */
279			pos_from_vec(pos, outvec);
280	greg	2.14	/* sum radial basis function */
281	greg	2.1	rbfp = rp->rbfa;
282			for (n = rp->nrbf; n--; rbfp++) {
283	greg	2.16	int d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
284			(pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
285			rad2 = R2ANG(rbfp->crad);
286			rad2 *= rad2;
287	greg	2.17	if (d2 > rad2*rfact2)
288	greg	2.16	continue;
289	greg	2.1	ovec_from_pos(odir, rbfp->gx, rbfp->gy);
290	greg	2.16	res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
291	greg	2.1	}
292	greg	2.13	if (res < minval) /* never return less than minval */
293			return(minval);
294	greg	2.1	return(res);
295			}
296
297			/* Insert a new directional scattering function in our global list */
298			int
299			insert_dsf(RBFNODE *newrbf)
300			{
301			RBFNODE rbf, rbf_last;
302			int pos;
303			/* check for redundant meas. */
304			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
305			if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
306			fprintf(stderr,
307	greg	2.22	"%s: Duplicate incident measurement ignored at (%.1f,%.1f)\n",
308			progname, get_theta180(newrbf->invec),
309			get_phi360(newrbf->invec));
310	greg	2.1	free(newrbf);
311			return(-1);
312			}
313			/* keep in ascending theta order */
314			for (rbf_last = NULL, rbf = dsf_list; rbf != NULL;
315			rbf_last = rbf, rbf = rbf->next)
316			if (single_plane_incident && input_orient*rbf->invec[2] <
317			input_orient*newrbf->invec[2])
318			break;
319			if (rbf_last == NULL) { /* insert new node in list */
320			newrbf->ord = 0;
321			newrbf->next = dsf_list;
322			dsf_list = newrbf;
323			} else {
324			newrbf->ord = rbf_last->ord + 1;
325			newrbf->next = rbf;
326			rbf_last->next = newrbf;
327			}
328			rbf_last = newrbf;
329			while (rbf != NULL) { /* update ordinal positions */
330			rbf->ord = rbf_last->ord + 1;
331			rbf_last = rbf;
332			rbf = rbf->next;
333			}
334			return(newrbf->ord);
335			}
336
337			/* Get the DSF indicated by its ordinal position */
338			RBFNODE *
339			get_dsf(int ord)
340			{
341			RBFNODE *rbf;
342
343			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
344	greg	2.3	if (rbf->ord == ord)
345	greg	2.1	return(rbf);
346			return(NULL);
347			}
348
349			/* Get triangle surface orientation (unnormalized) */
350			void
351			tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
352			{
353			FVECT v2minus1, v3minus2;
354
355			VSUB(v2minus1, v2, v1);
356			VSUB(v3minus2, v3, v2);
357			VCROSS(vres, v2minus1, v3minus2);
358			}
359
360			/* Determine if vertex order is reversed (inward normal) */
361			int
362			is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
363			{
364			FVECT tor;
365
366			tri_orient(tor, v1, v2, v3);
367
368			return(DOT(tor, v2) < 0.);
369			}
370
371			/* Find vertices completing triangles on either side of the given edge */
372			int
373			get_triangles(RBFNODE rbfv[2], const MIGRATION mig)
374			{
375	greg	2.4	const MIGRATION ej1, ej2;
376	greg	2.1	RBFNODE *tv;
377
378			rbfv[0] = rbfv[1] = NULL;
379			if (mig == NULL)
380			return(0);
381	greg	2.4	for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
382			ej1 = nextedge(mig->rbfv[0],ej1)) {
383			if (ej1 == mig)
384	greg	2.1	continue;
385	greg	2.4	tv = opp_rbf(mig->rbfv[0],ej1);
386	greg	2.1	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
387			if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
388			rbfv[is_rev_tri(mig->rbfv[0]->invec,
389			mig->rbfv[1]->invec,
390			tv->invec)] = tv;
391			break;
392			}
393			}
394			return((rbfv[0] != NULL) + (rbfv[1] != NULL));
395			}
396
397	greg	2.20	/* Advect and allocate new RBF along edge (internal call) */
398			RBFNODE *
399			e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
400			{
401			double cthresh = FTINY;
402			RBFNODE *rbf;
403			int n, i, j;
404			double t, full_dist;
405			/* get relative position */
406			t = Acos(DOT(invec, mig->rbfv[0]->invec));
407			if (t < M_PI/grid_res) { /* near first DSF */
408			n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
409			rbf = (RBFNODE *)malloc(n);
410			if (rbf == NULL)
411			goto memerr;
412			memcpy(rbf, mig->rbfv[0], n); /* just duplicate */
413			rbf->next = NULL; rbf->ejl = NULL;
414			return(rbf);
415			}
416			full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
417			if (t > full_dist-M_PI/grid_res) { /* near second DSF */
418			n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
419			rbf = (RBFNODE *)malloc(n);
420			if (rbf == NULL)
421			goto memerr;
422			memcpy(rbf, mig->rbfv[1], n); /* just duplicate */
423			rbf->next = NULL; rbf->ejl = NULL;
424			return(rbf);
425			}
426			t /= full_dist;
427			tryagain:
428			n = 0; /* count migrating particles */
429			for (i = 0; i < mtx_nrows(mig); i++)
430			for (j = 0; j < mtx_ncols(mig); j++)
431			n += (mtx_coef(mig,i,j) > cthresh);
432			/* are we over our limit? */
433			if ((lobe_lim > 0) & (n > lobe_lim)) {
434			cthresh = cthresh2. + 10.FTINY;
435			goto tryagain;
436			}
437			#ifdef DEBUG
438			fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
439			mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
440			#endif
441			rbf = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(n-1));
442			if (rbf == NULL)
443			goto memerr;
444			rbf->next = NULL; rbf->ejl = NULL;
445			VCOPY(rbf->invec, invec);
446			rbf->nrbf = n;
447			rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
448			n = 0; /* advect RBF lobes */
449			for (i = 0; i < mtx_nrows(mig); i++) {
450			const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i];
451			const float peak0 = rbf0i->peak;
452			const double rad0 = R2ANG(rbf0i->crad);
453			FVECT v0;
454			float mv;
455			ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
456			for (j = 0; j < mtx_ncols(mig); j++)
457			if ((mv = mtx_coef(mig,i,j)) > cthresh) {
458			const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j];
459			double rad2;
460			FVECT v;
461			int pos[2];
462			rad2 = R2ANG(rbf1j->crad);
463			rad2 = rad0rad0(1.-t) + rad2rad2t;
464			rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
465			rad0*rad0/rad2;
466			rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
467			ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
468			geodesic(v, v0, v, t, GEOD_REL);
469			pos_from_vec(pos, v);
470			rbf->rbfa[n].gx = pos[0];
471			rbf->rbfa[n].gy = pos[1];
472			++n;
473			}
474			}
475			rbf->vtotal = mig->rbfv[0]->vtotal; / turn ratio into actual */
476			return(rbf);
477			memerr:
478			fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
479			exit(1);
480			return(NULL); /* pro forma return */
481			}
482
483	greg	2.4	/* Clear our BSDF representation and free memory */
484			void
485			clear_bsdf_rep(void)
486			{
487			while (mig_list != NULL) {
488			MIGRATION *mig = mig_list;
489			mig_list = mig->next;
490			free(mig);
491			}
492			while (dsf_list != NULL) {
493			RBFNODE *rbf = dsf_list;
494			dsf_list = rbf->next;
495			free(rbf);
496			}
497	greg	2.19	bsdf_name[0] = '\0';
498			bsdf_manuf[0] = '\0';
499	greg	2.4	inp_coverage = 0;
500			single_plane_incident = -1;
501			input_orient = output_orient = 0;
502	greg	2.5	grid_res = GRIDRES;
503	greg	2.4	}
504
505	greg	2.1	/* Write our BSDF mesh interpolant out to the given binary stream */
506			void
507			save_bsdf_rep(FILE *ofp)
508			{
509			RBFNODE *rbf;
510			MIGRATION *mig;
511			int i, n;
512			/* finish header */
513	greg	2.19	if (bsdf_name[0])
514			fprintf(ofp, "NAME=%s\n", bsdf_name);
515			if (bsdf_manuf[0])
516			fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
517	greg	2.2	fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
518			fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
519	greg	2.5	fprintf(ofp, "GRIDRES=%d\n", grid_res);
520	greg	2.12	fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
521	greg	2.1	fputformat(BSDFREP_FMT, ofp);
522			fputc('\n', ofp);
523			/* write each DSF */
524			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
525			putint(rbf->ord, 4, ofp);
526			putflt(rbf->invec[0], ofp);
527			putflt(rbf->invec[1], ofp);
528			putflt(rbf->invec[2], ofp);
529			putflt(rbf->vtotal, ofp);
530			putint(rbf->nrbf, 4, ofp);
531			for (i = 0; i < rbf->nrbf; i++) {
532			putflt(rbf->rbfa[i].peak, ofp);
533			putint(rbf->rbfa[i].crad, 2, ofp);
534			putint(rbf->rbfa[i].gx, 1, ofp);
535			putint(rbf->rbfa[i].gy, 1, ofp);
536			}
537			}
538			putint(-1, 4, ofp); /* terminator */
539			/* write each migration matrix */
540	greg	2.2	for (mig = mig_list; mig != NULL; mig = mig->next) {
541			int zerocnt = 0;
542	greg	2.1	putint(mig->rbfv[0]->ord, 4, ofp);
543			putint(mig->rbfv[1]->ord, 4, ofp);
544	greg	2.2	/* write out as sparse data */
545	greg	2.1	n = mtx_nrows(mig) * mtx_ncols(mig);
546	greg	2.2	for (i = 0; i < n; i++) {
547	greg	2.3	if (zerocnt == 0xff) {
548			putint(0xff, 1, ofp); zerocnt = 0;
549	greg	2.2	}
550			if (mig->mtx[i] != 0) {
551			putint(zerocnt, 1, ofp); zerocnt = 0;
552			putflt(mig->mtx[i], ofp);
553			} else
554			++zerocnt;
555			}
556			putint(zerocnt, 1, ofp);
557	greg	2.1	}
558			putint(-1, 4, ofp); /* terminator */
559			putint(-1, 4, ofp);
560			if (fflush(ofp) == EOF) {
561			fprintf(stderr, "%s: error writing BSDF interpolant\n",
562			progname);
563			exit(1);
564			}
565			}
566
567	greg	2.2	/* Check header line for critical information */
568			static int
569			headline(char s, void p)
570			{
571			char fmt[32];
572
573	greg	2.19	if (!strncmp(s, "NAME=", 5)) {
574			strcpy(bsdf_name, s+5);
575			bsdf_name[strlen(bsdf_name)-1] = '\0';
576			}
577			if (!strncmp(s, "MANUFACT=", 9)) {
578			strcpy(bsdf_manuf, s+9);
579			bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
580			}
581	greg	2.2	if (!strncmp(s, "SYMMETRY=", 9)) {
582			inp_coverage = atoi(s+9);
583			single_plane_incident = !inp_coverage;
584			return(0);
585			}
586			if (!strncmp(s, "IO_SIDES=", 9)) {
587			sscanf(s+9, "%d %d", &input_orient, &output_orient);
588			return(0);
589			}
590	greg	2.5	if (!strncmp(s, "GRIDRES=", 8)) {
591			sscanf(s+8, "%d", &grid_res);
592			return(0);
593			}
594	greg	2.12	if (!strncmp(s, "BSDFMIN=", 8)) {
595			sscanf(s+8, "%lf", &bsdf_min);
596			return(0);
597			}
598	greg	2.2	if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
599			return(-1);
600			return(0);
601			}
602
603	greg	2.1	/* Read a BSDF mesh interpolant from the given binary stream */
604			int
605			load_bsdf_rep(FILE *ifp)
606			{
607			RBFNODE rbfh;
608			int from_ord, to_ord;
609			int i;
610	greg	2.4
611			clear_bsdf_rep();
612	greg	2.5	if (ifp == NULL)
613			return(0);
614	greg	2.21	if (getheader(ifp, headline, NULL) < 0 \|\| (single_plane_incident < 0) \|
615	greg	2.23	!input_orient \| !output_orient \|
616			(grid_res < 16) \| (grid_res > 256)) {
617	greg	2.1	fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
618			progname);
619			return(0);
620			}
621	greg	2.18	memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
622	greg	2.1	while ((rbfh.ord = getint(4, ifp)) >= 0) {
623			RBFNODE *newrbf;
624
625			rbfh.invec[0] = getflt(ifp);
626			rbfh.invec[1] = getflt(ifp);
627			rbfh.invec[2] = getflt(ifp);
628	greg	2.9	if (normalize(rbfh.invec) == 0) {
629			fprintf(stderr, "%s: zero incident vector\n", progname);
630			return(0);
631			}
632	greg	2.3	rbfh.vtotal = getflt(ifp);
633	greg	2.1	rbfh.nrbf = getint(4, ifp);
634			newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
635			sizeof(RBFVAL)*(rbfh.nrbf-1));
636			if (newrbf == NULL)
637			goto memerr;
638	greg	2.18	*newrbf = rbfh;
639	greg	2.1	for (i = 0; i < rbfh.nrbf; i++) {
640			newrbf->rbfa[i].peak = getflt(ifp);
641			newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
642			newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
643			newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
644			}
645			if (feof(ifp))
646			goto badEOF;
647			/* insert in global list */
648			if (insert_dsf(newrbf) != rbfh.ord) {
649			fprintf(stderr, "%s: error adding DSF\n", progname);
650			return(0);
651			}
652			}
653			/* read each migration matrix */
654			while ((from_ord = getint(4, ifp)) >= 0 &&
655			(to_ord = getint(4, ifp)) >= 0) {
656			RBFNODE *from_rbf = get_dsf(from_ord);
657			RBFNODE *to_rbf = get_dsf(to_ord);
658			MIGRATION *newmig;
659			int n;
660
661			if ((from_rbf == NULL) \| (to_rbf == NULL)) {
662			fprintf(stderr,
663			"%s: bad DSF reference in migration edge\n",
664			progname);
665			return(0);
666			}
667			n = from_rbf->nrbf * to_rbf->nrbf;
668			newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
669			sizeof(float)*(n-1));
670			if (newmig == NULL)
671			goto memerr;
672			newmig->rbfv[0] = from_rbf;
673			newmig->rbfv[1] = to_rbf;
674	greg	2.2	memset(newmig->mtx, 0, sizeof(float)*n);
675			for (i = 0; ; ) { /* read sparse data */
676			int zc = getint(1, ifp) & 0xff;
677			if ((i += zc) >= n)
678			break;
679	greg	2.3	if (zc == 0xff)
680			continue;
681	greg	2.2	newmig->mtx[i++] = getflt(ifp);
682			}
683	greg	2.1	if (feof(ifp))
684			goto badEOF;
685			/* insert in edge lists */
686			newmig->enxt[0] = from_rbf->ejl;
687			from_rbf->ejl = newmig;
688			newmig->enxt[1] = to_rbf->ejl;
689			to_rbf->ejl = newmig;
690			/* push onto global list */
691			newmig->next = mig_list;
692			mig_list = newmig;
693			}
694			return(1); /* success! */
695			memerr:
696			fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
697			exit(1);
698			badEOF:
699			fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
700			return(0);
701			}