src/cv/bsdfmesh.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfmesh.c,v 2.28 2014/03/26 22:29:08 greg Exp $";
#endif
/*
 * Create BSDF advection mesh from radial basis functions.
 *
 *      G. Ward
 */

#ifndef _WIN32
#include <unistd.h>
#include <sys/wait.h>
#include <sys/mman.h>
#endif
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "bsdfrep.h"

#ifndef NEIGH_FACT2
#define NEIGH_FACT2     0.1     /* empirical neighborhood distance weight */
#endif
                                /* number of processes to run */
int                     nprocs = 1;
                                /* number of children (-1 in child) */
static int              nchild = 0;

/* Create a new migration holder (sharing memory for multiprocessing) */
static MIGRATION *
new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
{
        size_t          memlen = sizeof(MIGRATION) +
                                sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
        MIGRATION       *newmig;
#ifdef _WIN32
        if (nprocs > 1)
                fprintf(stderr, "%s: warning - multiprocessing not supported\n",
                                progname);
        nprocs = 1;
        newmig = (MIGRATION *)malloc(memlen);
#else
        if (nprocs <= 1) {                      /* single process? */
                newmig = (MIGRATION *)malloc(memlen);
        } else {                                /* else need to share memory */
                newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
                                                MAP_ANON|MAP_SHARED, -1, 0);
                if ((void *)newmig == MAP_FAILED)
                        newmig = NULL;
        }
#endif
        if (newmig == NULL) {
                fprintf(stderr, "%s: cannot allocate new migration\n", progname);
                exit(1);
        }
        newmig->rbfv[0] = from_rbf;
        newmig->rbfv[1] = to_rbf;
                                                /* insert in edge lists */
        newmig->enxt[0] = from_rbf->ejl;
        from_rbf->ejl = newmig;
        newmig->enxt[1] = to_rbf->ejl;
        to_rbf->ejl = newmig;
        newmig->next = mig_list;                /* push onto global list */
        return(mig_list = newmig);
}

#ifdef _WIN32
#define await_children(n)       (void)(n)
#define run_subprocess()        0
#define end_subprocess()        (void)0
#else

/* Wait for the specified number of child processes to complete */
static void
await_children(int n)
{
        int     exit_status = 0;

        if (n > nchild)
                n = nchild;
        while (n-- > 0) {
                int     status;
                if (wait(&status) < 0) {
                        fprintf(stderr, "%s: missing child(ren)!\n", progname);
                        nchild = 0;
                        break;
                }
                --nchild;
                if (status) {                   /* something wrong */
                        if ((status = WEXITSTATUS(status)))
                                exit_status = status;
                        else
                                exit_status += !exit_status;
                        fprintf(stderr, "%s: subprocess died\n", progname);
                        n = nchild;             /* wait for the rest */
                }
        }
        if (exit_status)
                exit(exit_status);
}

/* Start child process if multiprocessing selected */
static pid_t
run_subprocess(void)
{
        int     status;
        pid_t   pid;

        if (nprocs <= 1)                        /* any children requested? */
                return(0);
        await_children(nchild + 1 - nprocs);    /* free up child process */
        if ((pid = fork())) {
                if (pid < 0) {
                        fprintf(stderr, "%s: cannot fork subprocess\n",
                                        progname);
                        await_children(nchild);
                        exit(1);
                }
                ++nchild;                       /* subprocess started */
                return(pid);
        }
        nchild = -1;
        return(0);                              /* put child to work */
}

/* If we are in subprocess, call exit */
#define end_subprocess()        if (nchild < 0) _exit(0); else

#endif  /* ! _WIN32 */

/* Compute normalized distribution scattering functions for comparison */
static void
compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
{
        const double    nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
        const double    nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
        int             x, y;
        FVECT           dv;

        for (x = GRIDRES; x--; )
            for (y = GRIDRES; y--; ) {
                ovec_from_pos(dv, x, y);        /* cube root (brightness) */
                dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
                dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
            }
}       

/* Compute neighborhood distance-squared (dissimilarity) */
static double
neighborhood_dist2(int x0, int y0, int x1, int y1)
{
        int     rad = GRIDRES>>5;
        double  sum2 = 0.;
        double  d;
        int     p[4];
        int     i, j;
                                                /* check radius */
        p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
        for (i = 4; i--; ) {
                if (p[i] < rad) rad = p[i];
                if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
        }
        for (i = -rad; i <= rad; i++)
            for (j = -rad; j <= rad; j++) {
                d = dsf_grid[x0+i][y0+j].val[0] -
                        dsf_grid[x1+i][y1+j].val[1];
                sum2 += d*d;
            }
        return(sum2 / (4*rad*(rad+1) + 1));
}

/* Compute distance between two RBF lobes */
double
lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
{
        FVECT   vfrom, vto;
        double  d, res;
                                        /* quadratic cost function */
        ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
        ovec_from_pos(vto, rbf2->gx, rbf2->gy);
        d = Acos(DOT(vfrom, vto));
        res = d*d;
        d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
        res += d*d;
                                        /* neighborhood difference */
        res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
                                                rbf2->gx, rbf2->gy );
        return(res);
}


/* Compute and insert migration along directed edge (may fork child) */
static MIGRATION *
create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
{
        MIGRATION       *newmig;
        int             i, j;
                                                /* check if exists already */
        for (newmig = from_rbf->ejl; newmig != NULL;
                        newmig = nextedge(from_rbf,newmig))
                if (newmig->rbfv[1] == to_rbf)
                        return(NULL);
                                                /* else allocate */
#ifdef DEBUG
        fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
                        get_theta180(from_rbf->invec),
                        get_phi360(from_rbf->invec));
        fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
                        get_theta180(to_rbf->invec),
                        get_phi360(to_rbf->invec), 
                        from_rbf->nrbf, to_rbf->nrbf);
#endif
        newmig = new_migration(from_rbf, to_rbf);
        if (run_subprocess())
                return(newmig);                 /* child continues */

                                                /* compute transport plan */
        compute_nDSFs(from_rbf, to_rbf);
        plan_transport(newmig);

        for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
            double      nf = rbf_volume(&from_rbf->rbfa[i]);
            if (nf <= FTINY) continue;
            nf = from_rbf->vtotal / nf;
            for (j = to_rbf->nrbf; j--; )
                mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
        }
        end_subprocess();                       /* exit here if subprocess */
        return(newmig);
}

/* Check if prospective vertex would create overlapping triangle */
static int
overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
{
        const MIGRATION *ej;
        RBFNODE         *vother[2];
        int             im_rev;
                                                /* find shared edge in mesh */
        for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
                const RBFNODE   *tv = opp_rbf(pv,ej);
                if (tv == bv0) {
                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
                                        ej->rbfv[1]->invec, bv1->invec);
                        break;
                }
                if (tv == bv1) {
                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
                                        ej->rbfv[1]->invec, bv0->invec);
                        break;
                }
        }
        if (!get_triangles(vother, ej))         /* triangle on same side? */
                return(0);
        return(vother[im_rev] != NULL);
}

/* Find convex hull vertex to complete triangle (oriented call) */
static RBFNODE *
find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
{
        FVECT   vmid, vejn, vp;
        RBFNODE *rbf, *rbfbest = NULL;
        double  dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;

        VSUB(vejn, rbf1->invec, rbf0->invec);
        VADD(vmid, rbf0->invec, rbf1->invec);
        if (normalize(vejn) == 0 || normalize(vmid) == 0)
                return(NULL);
                                                /* XXX exhaustive search */
        /* Find triangle with minimum rotation from perpendicular */
        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
                if ((rbf == rbf0) | (rbf == rbf1))
                        continue;
                tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
                if (DOT(vp, vmid) <= FTINY)
                        continue;               /* wrong orientation */
                area2 = .25*DOT(vp,vp);
                VSUB(vp, rbf->invec, vmid);
                dprod = -DOT(vp, vejn);
                VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
                dprod = DOT(vp, vmid) / VLEN(vp);
                if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
                        continue;               /* found better already */
                if (overlaps_tri(rbf0, rbf1, rbf))
                        continue;               /* overlaps another triangle */
                rbfbest = rbf;
                bestdprod = dprod;              /* new one to beat */
                bestarea2 = area2;
        }
        return(rbfbest);
}

/* Create new migration edge and grow mesh recursively around it */
static void
mesh_from_edge(MIGRATION *edge)
{
        MIGRATION       *ej0, *ej1;
        RBFNODE         *tvert[2];
        
        if (edge == NULL)
                return;
                                                /* triangle on either side? */
        get_triangles(tvert, edge);
        if (tvert[0] == NULL) {                 /* grow mesh on right */
                tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
                if (tvert[0] != NULL) {
                        if (tvert[0]->ord > edge->rbfv[0]->ord)
                                ej0 = create_migration(edge->rbfv[0], tvert[0]);
                        else
                                ej0 = create_migration(tvert[0], edge->rbfv[0]);
                        if (tvert[0]->ord > edge->rbfv[1]->ord)
                                ej1 = create_migration(edge->rbfv[1], tvert[0]);
                        else
                                ej1 = create_migration(tvert[0], edge->rbfv[1]);
                        mesh_from_edge(ej0);
                        mesh_from_edge(ej1);
                        return;
                }
        }
        if (tvert[1] == NULL) {                 /* grow mesh on left */
                tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
                if (tvert[1] != NULL) {
                        if (tvert[1]->ord > edge->rbfv[0]->ord)
                                ej0 = create_migration(edge->rbfv[0], tvert[1]);
                        else
                                ej0 = create_migration(tvert[1], edge->rbfv[0]);
                        if (tvert[1]->ord > edge->rbfv[1]->ord)
                                ej1 = create_migration(edge->rbfv[1], tvert[1]);
                        else
                                ej1 = create_migration(tvert[1], edge->rbfv[1]);
                        mesh_from_edge(ej0);
                        mesh_from_edge(ej1);
                }
        }
}

/* Add normal direction if missing */
static void
check_normal_incidence(void)
{
        static FVECT            norm_vec = {.0, .0, 1.};
        const int               saved_nprocs = nprocs;
        RBFNODE                 *near_rbf, *mir_rbf, *rbf;
        double                  bestd;
        int                     n;

        if (dsf_list == NULL)
                return;                         /* XXX should be error? */
        near_rbf = dsf_list;
        bestd = input_orient*near_rbf->invec[2];
        if (single_plane_incident) {            /* ordered plane incidence? */
                if (bestd >= 1.-2.*FTINY)
                        return;                 /* already have normal */
        } else {
                switch (inp_coverage) {
                case INP_QUAD1:
                case INP_QUAD2:
                case INP_QUAD3:
                case INP_QUAD4:
                        break;                  /* quadrilateral symmetry? */
                default:
                        return;                 /* else we can interpolate */
                }
                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
                        const double    d = input_orient*rbf->invec[2];
                        if (d >= 1.-2.*FTINY)
                                return;         /* seems we have normal */
                        if (d > bestd) {
                                near_rbf = rbf;
                                bestd = d;
                        }
                }
        }
        if (mig_list != NULL) {                 /* need to be called first */
                fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
                                progname);
                exit(1);
        }
#ifdef DEBUG
        fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
                        get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
#endif
                                                /* mirror nearest incidence */
        n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
        mir_rbf = (RBFNODE *)malloc(n);
        if (mir_rbf == NULL)
                goto memerr;
        memcpy(mir_rbf, near_rbf, n);
        mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
        mir_rbf->next = NULL;
        mir_rbf->ejl = NULL;
        rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
        nprocs = 1;                             /* compute migration matrix */
        if (create_migration(mir_rbf, near_rbf) == NULL)
                exit(1);                        /* XXX should never happen! */
        norm_vec[2] = input_orient;             /* interpolate normal dist. */
        rbf = e_advect_rbf(mig_list, norm_vec, 0);
        nprocs = saved_nprocs;                  /* final clean-up */
        free(mir_rbf);
        free(mig_list);
        mig_list = near_rbf->ejl = NULL;
        insert_dsf(rbf);                        /* insert interpolated normal */
        return;
memerr:
        fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
                                progname);
        exit(1);
}
        
/* Build our triangle mesh from recorded RBFs */
void
build_mesh(void)
{
        double          best2 = M_PI*M_PI;
        RBFNODE         *shrt_edj[2];
        RBFNODE         *rbf0, *rbf1;
                                                /* add normal if needed */
        check_normal_incidence();
                                                /* check if isotropic */
        if (single_plane_incident) {
                for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
                        if (rbf0->next != NULL)
                                create_migration(rbf0, rbf0->next);
                await_children(nchild);
                return;
        }
        shrt_edj[0] = shrt_edj[1] = NULL;       /* start w/ shortest edge */
        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
            for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
                double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
                if (dist2 < best2) {
                        shrt_edj[0] = rbf0;
                        shrt_edj[1] = rbf1;
                        best2 = dist2;
                }
        }
        if (shrt_edj[0] == NULL) {
                fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
                exit(1);
        }
                                                /* build mesh from this edge */
        if (shrt_edj[0]->ord < shrt_edj[1]->ord)
                mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
        else
                mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
                                                /* complete migrations */
        await_children(nchild);
}
Revision:	2.29
Committed:	Thu Mar 27 03:49:14 2014 UTC (10 years ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R2
Changes since 2.28:	+2 -2 lines
Log Message:	Minor fixes
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.29	static const char RCSid[] = "$Id: bsdfmesh.c,v 2.28 2014/03/26 22:29:08 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Create BSDF advection mesh from radial basis functions.
6			*
7			* G. Ward
8			*/
9
10			#ifndef _WIN32
11			#include <unistd.h>
12			#include <sys/wait.h>
13			#include <sys/mman.h>
14			#endif
15			#define _USE_MATH_DEFINES
16			#include <stdio.h>
17			#include <stdlib.h>
18			#include <string.h>
19			#include <math.h>
20			#include "bsdfrep.h"
21	greg	2.19
22			#ifndef NEIGH_FACT2
23	greg	2.21	#define NEIGH_FACT2 0.1 /* empirical neighborhood distance weight */
24	greg	2.19	#endif
25	greg	2.1	/* number of processes to run */
26			int nprocs = 1;
27			/* number of children (-1 in child) */
28			static int nchild = 0;
29
30	greg	2.2	/* Create a new migration holder (sharing memory for multiprocessing) */
31			static MIGRATION *
32			new_migration(RBFNODE from_rbf, RBFNODE to_rbf)
33			{
34			size_t memlen = sizeof(MIGRATION) +
35			sizeof(float)(from_rbf->nrbfto_rbf->nrbf - 1);
36			MIGRATION *newmig;
37			#ifdef _WIN32
38			if (nprocs > 1)
39			fprintf(stderr, "%s: warning - multiprocessing not supported\n",
40			progname);
41			nprocs = 1;
42			newmig = (MIGRATION *)malloc(memlen);
43			#else
44			if (nprocs <= 1) { /* single process? */
45			newmig = (MIGRATION *)malloc(memlen);
46			} else { /* else need to share memory */
47			newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ\|PROT_WRITE,
48			MAP_ANON\|MAP_SHARED, -1, 0);
49			if ((void *)newmig == MAP_FAILED)
50			newmig = NULL;
51			}
52			#endif
53			if (newmig == NULL) {
54			fprintf(stderr, "%s: cannot allocate new migration\n", progname);
55			exit(1);
56			}
57			newmig->rbfv[0] = from_rbf;
58			newmig->rbfv[1] = to_rbf;
59			/* insert in edge lists */
60			newmig->enxt[0] = from_rbf->ejl;
61			from_rbf->ejl = newmig;
62			newmig->enxt[1] = to_rbf->ejl;
63			to_rbf->ejl = newmig;
64			newmig->next = mig_list; /* push onto global list */
65			return(mig_list = newmig);
66			}
67
68			#ifdef _WIN32
69			#define await_children(n) (void)(n)
70			#define run_subprocess() 0
71			#define end_subprocess() (void)0
72			#else
73
74			/* Wait for the specified number of child processes to complete */
75			static void
76			await_children(int n)
77			{
78			int exit_status = 0;
79
80			if (n > nchild)
81			n = nchild;
82			while (n-- > 0) {
83			int status;
84			if (wait(&status) < 0) {
85			fprintf(stderr, "%s: missing child(ren)!\n", progname);
86			nchild = 0;
87			break;
88			}
89			--nchild;
90			if (status) { /* something wrong */
91			if ((status = WEXITSTATUS(status)))
92			exit_status = status;
93			else
94			exit_status += !exit_status;
95			fprintf(stderr, "%s: subprocess died\n", progname);
96			n = nchild; /* wait for the rest */
97			}
98			}
99			if (exit_status)
100			exit(exit_status);
101			}
102
103			/* Start child process if multiprocessing selected */
104			static pid_t
105			run_subprocess(void)
106			{
107			int status;
108			pid_t pid;
109
110			if (nprocs <= 1) /* any children requested? */
111			return(0);
112			await_children(nchild + 1 - nprocs); /* free up child process */
113			if ((pid = fork())) {
114			if (pid < 0) {
115			fprintf(stderr, "%s: cannot fork subprocess\n",
116			progname);
117	greg	2.6	await_children(nchild);
118	greg	2.2	exit(1);
119			}
120			++nchild; /* subprocess started */
121			return(pid);
122			}
123			nchild = -1;
124			return(0); /* put child to work */
125			}
126
127			/* If we are in subprocess, call exit */
128			#define end_subprocess() if (nchild < 0) _exit(0); else
129
130			#endif /* ! _WIN32 */
131
132	greg	2.19	/* Compute normalized distribution scattering functions for comparison */
133			static void
134			compute_nDSFs(const RBFNODE rbf0, const RBFNODE rbf1)
135			{
136			const double nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
137			const double nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
138			int x, y;
139			FVECT dv;
140
141			for (x = GRIDRES; x--; )
142			for (y = GRIDRES; y--; ) {
143	greg	2.20	ovec_from_pos(dv, x, y); /* cube root (brightness) */
144			dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
145			dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
146	greg	2.19	}
147			}
148
149			/* Compute neighborhood distance-squared (dissimilarity) */
150			static double
151			neighborhood_dist2(int x0, int y0, int x1, int y1)
152			{
153			int rad = GRIDRES>>5;
154			double sum2 = 0.;
155			double d;
156			int p[4];
157			int i, j;
158			/* check radius */
159			p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
160			for (i = 4; i--; ) {
161			if (p[i] < rad) rad = p[i];
162			if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
163			}
164			for (i = -rad; i <= rad; i++)
165			for (j = -rad; j <= rad; j++) {
166			d = dsf_grid[x0+i][y0+j].val[0] -
167			dsf_grid[x1+i][y1+j].val[1];
168			sum2 += d*d;
169			}
170			return(sum2 / (4rad(rad+1) + 1));
171			}
172
173	greg	2.27	/* Compute distance between two RBF lobes */
174			double
175			lobe_distance(RBFVAL rbf1, RBFVAL rbf2)
176			{
177			FVECT vfrom, vto;
178			double d, res;
179			/* quadratic cost function */
180			ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
181			ovec_from_pos(vto, rbf2->gx, rbf2->gy);
182			d = Acos(DOT(vfrom, vto));
183			res = d*d;
184			d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
185			res += d*d;
186			/* neighborhood difference */
187			res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
188			rbf2->gx, rbf2->gy );
189			return(res);
190	greg	2.1	}
191
192	greg	2.26
193	greg	2.1	/* Compute and insert migration along directed edge (may fork child) */
194			static MIGRATION *
195			create_migration(RBFNODE from_rbf, RBFNODE to_rbf)
196			{
197			MIGRATION *newmig;
198	greg	2.6	int i, j;
199	greg	2.1	/* check if exists already */
200			for (newmig = from_rbf->ejl; newmig != NULL;
201			newmig = nextedge(from_rbf,newmig))
202			if (newmig->rbfv[1] == to_rbf)
203			return(NULL);
204			/* else allocate */
205	greg	2.7	#ifdef DEBUG
206	greg	2.14	fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
207	greg	2.7	get_theta180(from_rbf->invec),
208			get_phi360(from_rbf->invec));
209	greg	2.14	fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
210	greg	2.7	get_theta180(to_rbf->invec),
211			get_phi360(to_rbf->invec),
212			from_rbf->nrbf, to_rbf->nrbf);
213			#endif
214	greg	2.1	newmig = new_migration(from_rbf, to_rbf);
215			if (run_subprocess())
216			return(newmig); /* child continues */
217	greg	2.27
218			/* compute transport plan */
219			compute_nDSFs(from_rbf, to_rbf);
220			plan_transport(newmig);
221	greg	2.6
222	greg	2.1	for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
223	greg	2.6	double nf = rbf_volume(&from_rbf->rbfa[i]);
224	greg	2.1	if (nf <= FTINY) continue;
225			nf = from_rbf->vtotal / nf;
226			for (j = to_rbf->nrbf; j--; )
227	greg	2.6	mtx_coef(newmig,i,j) = nf; / row now sums to 1.0 */
228	greg	2.1	}
229			end_subprocess(); /* exit here if subprocess */
230			return(newmig);
231			}
232
233			/* Check if prospective vertex would create overlapping triangle */
234			static int
235			overlaps_tri(const RBFNODE bv0, const RBFNODE bv1, const RBFNODE *pv)
236			{
237			const MIGRATION *ej;
238			RBFNODE *vother[2];
239			int im_rev;
240			/* find shared edge in mesh */
241			for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
242			const RBFNODE *tv = opp_rbf(pv,ej);
243			if (tv == bv0) {
244			im_rev = is_rev_tri(ej->rbfv[0]->invec,
245			ej->rbfv[1]->invec, bv1->invec);
246			break;
247			}
248			if (tv == bv1) {
249			im_rev = is_rev_tri(ej->rbfv[0]->invec,
250			ej->rbfv[1]->invec, bv0->invec);
251			break;
252			}
253			}
254			if (!get_triangles(vother, ej)) /* triangle on same side? */
255			return(0);
256			return(vother[im_rev] != NULL);
257			}
258
259	greg	2.14	/* Find convex hull vertex to complete triangle (oriented call) */
260	greg	2.1	static RBFNODE *
261			find_chull_vert(const RBFNODE rbf0, const RBFNODE rbf1)
262			{
263			FVECT vmid, vejn, vp;
264			RBFNODE rbf, rbfbest = NULL;
265			double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
266
267			VSUB(vejn, rbf1->invec, rbf0->invec);
268			VADD(vmid, rbf0->invec, rbf1->invec);
269			if (normalize(vejn) == 0 \|\| normalize(vmid) == 0)
270			return(NULL);
271			/* XXX exhaustive search */
272			/* Find triangle with minimum rotation from perpendicular */
273			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
274			if ((rbf == rbf0) \| (rbf == rbf1))
275			continue;
276			tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
277			if (DOT(vp, vmid) <= FTINY)
278			continue; /* wrong orientation */
279			area2 = .25*DOT(vp,vp);
280	greg	2.14	VSUB(vp, rbf->invec, vmid);
281	greg	2.1	dprod = -DOT(vp, vejn);
282			VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
283			dprod = DOT(vp, vmid) / VLEN(vp);
284			if (dprod <= bestdprod + FTINY(1 - 2(area2 < bestarea2)))
285			continue; /* found better already */
286			if (overlaps_tri(rbf0, rbf1, rbf))
287			continue; /* overlaps another triangle */
288			rbfbest = rbf;
289			bestdprod = dprod; /* new one to beat */
290			bestarea2 = area2;
291			}
292			return(rbfbest);
293			}
294
295			/* Create new migration edge and grow mesh recursively around it */
296			static void
297			mesh_from_edge(MIGRATION *edge)
298			{
299			MIGRATION ej0, ej1;
300			RBFNODE *tvert[2];
301
302			if (edge == NULL)
303			return;
304			/* triangle on either side? */
305			get_triangles(tvert, edge);
306			if (tvert[0] == NULL) { /* grow mesh on right */
307			tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
308			if (tvert[0] != NULL) {
309			if (tvert[0]->ord > edge->rbfv[0]->ord)
310			ej0 = create_migration(edge->rbfv[0], tvert[0]);
311			else
312			ej0 = create_migration(tvert[0], edge->rbfv[0]);
313			if (tvert[0]->ord > edge->rbfv[1]->ord)
314			ej1 = create_migration(edge->rbfv[1], tvert[0]);
315			else
316			ej1 = create_migration(tvert[0], edge->rbfv[1]);
317			mesh_from_edge(ej0);
318			mesh_from_edge(ej1);
319	greg	2.28	return;
320	greg	2.1	}
321	greg	2.28	}
322			if (tvert[1] == NULL) { /* grow mesh on left */
323	greg	2.1	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
324			if (tvert[1] != NULL) {
325			if (tvert[1]->ord > edge->rbfv[0]->ord)
326			ej0 = create_migration(edge->rbfv[0], tvert[1]);
327			else
328			ej0 = create_migration(tvert[1], edge->rbfv[0]);
329			if (tvert[1]->ord > edge->rbfv[1]->ord)
330			ej1 = create_migration(edge->rbfv[1], tvert[1]);
331			else
332			ej1 = create_migration(tvert[1], edge->rbfv[1]);
333			mesh_from_edge(ej0);
334			mesh_from_edge(ej1);
335			}
336			}
337			}
338	greg	2.15
339			/* Add normal direction if missing */
340			static void
341			check_normal_incidence(void)
342			{
343	greg	2.25	static FVECT norm_vec = {.0, .0, 1.};
344	greg	2.16	const int saved_nprocs = nprocs;
345			RBFNODE near_rbf, mir_rbf, *rbf;
346			double bestd;
347			int n;
348	greg	2.15
349			if (dsf_list == NULL)
350			return; /* XXX should be error? */
351			near_rbf = dsf_list;
352			bestd = input_orient*near_rbf->invec[2];
353			if (single_plane_incident) { /* ordered plane incidence? */
354			if (bestd >= 1.-2.*FTINY)
355			return; /* already have normal */
356			} else {
357			switch (inp_coverage) {
358			case INP_QUAD1:
359			case INP_QUAD2:
360			case INP_QUAD3:
361			case INP_QUAD4:
362			break; /* quadrilateral symmetry? */
363			default:
364			return; /* else we can interpolate */
365			}
366			for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
367			const double d = input_orient*rbf->invec[2];
368			if (d >= 1.-2.*FTINY)
369			return; /* seems we have normal */
370			if (d > bestd) {
371			near_rbf = rbf;
372			bestd = d;
373			}
374			}
375			}
376			if (mig_list != NULL) { /* need to be called first */
377			fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
378			progname);
379			exit(1);
380			}
381			#ifdef DEBUG
382			fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
383			get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
384			#endif
385			/* mirror nearest incidence */
386			n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
387			mir_rbf = (RBFNODE *)malloc(n);
388			if (mir_rbf == NULL)
389			goto memerr;
390			memcpy(mir_rbf, near_rbf, n);
391			mir_rbf->ord = near_rbf->ord - 1; /* not used, I think */
392			mir_rbf->next = NULL;
393	greg	2.22	mir_rbf->ejl = NULL;
394	greg	2.15	rev_rbf_symmetry(mir_rbf, MIRROR_X\|MIRROR_Y);
395			nprocs = 1; /* compute migration matrix */
396	greg	2.22	if (create_migration(mir_rbf, near_rbf) == NULL)
397	greg	2.15	exit(1); /* XXX should never happen! */
398	greg	2.25	norm_vec[2] = input_orient; /* interpolate normal dist. */
399	greg	2.29	rbf = e_advect_rbf(mig_list, norm_vec, 0);
400	greg	2.15	nprocs = saved_nprocs; /* final clean-up */
401			free(mir_rbf);
402			free(mig_list);
403			mig_list = near_rbf->ejl = NULL;
404			insert_dsf(rbf); /* insert interpolated normal */
405			return;
406			memerr:
407			fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
408			progname);
409			exit(1);
410			}
411	greg	2.1
412			/* Build our triangle mesh from recorded RBFs */
413			void
414			build_mesh(void)
415			{
416			double best2 = M_PI*M_PI;
417			RBFNODE *shrt_edj[2];
418			RBFNODE rbf0, rbf1;
419	greg	2.15	/* add normal if needed */
420			check_normal_incidence();
421	greg	2.1	/* check if isotropic */
422			if (single_plane_incident) {
423			for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
424			if (rbf0->next != NULL)
425			create_migration(rbf0, rbf0->next);
426			await_children(nchild);
427			return;
428			}
429			shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
430			for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
431			for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
432			double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
433			if (dist2 < best2) {
434			shrt_edj[0] = rbf0;
435			shrt_edj[1] = rbf1;
436			best2 = dist2;
437			}
438			}
439			if (shrt_edj[0] == NULL) {
440			fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
441			exit(1);
442			}
443			/* build mesh from this edge */
444			if (shrt_edj[0]->ord < shrt_edj[1]->ord)
445			mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
446			else
447			mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
448			/* complete migrations */
449			await_children(nchild);
450			}