src/cv/bsdfmesh.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfmesh.c,v 2.8 2013/03/20 01:00:22 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"
                                /* number of processes to run */
int                     nprocs = 1;
                                /* number of children (-1 in child) */
static int              nchild = 0;

typedef struct {
        int             nrows, ncols;   /* array size (matches migration) */
        float           *price;         /* migration prices */
        short           *sord;          /* sort for each row, low to high */
} PRICEMAT;                     /* sorted pricing matrix */

#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)

/* 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 */

/* Comparison routine needed for sorting price row */
static int
msrt_cmp(void *b, const void *p1, const void *p2)
{
        PRICEMAT        *pm = (PRICEMAT *)b;
        int             ri = ((const short *)p1 - pm->sord) / pm->ncols;
        float           c1 = pricerow(pm,ri)[*(const short *)p1];
        float           c2 = pricerow(pm,ri)[*(const short *)p2];

        if (c1 > c2) return(1);
        if (c1 < c2) return(-1);
        return(0);
}

/* Compute (and allocate) migration price matrix for optimization */
static void
price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
{
        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
        int     i, j;

        pm->nrows = from_rbf->nrbf;
        pm->ncols = to_rbf->nrbf;
        pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
        pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
        
        if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
                                progname);
                exit(1);
        }
        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);

        for (i = from_rbf->nrbf; i--; ) {
            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
            FVECT               vfrom;
            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
            for (j = to_rbf->nrbf; j--; ) {
                double          dprod = DOT(vfrom, vto[j]);
                pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
                psortrow(pm,i)[j] = j;
            }
            qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
        }
        free(vto);
}

/* Free price matrix */
static void
free_routes(PRICEMAT *pm)
{
        free(pm->price); pm->price = NULL;
        free(pm->sord); pm->sord = NULL;
}

/* Compute minimum (optimistic) cost for moving the given source material */
static double
min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
{
        double          total_cost = 0;
        int             j;

        if (amt2move <= FTINY)                  /* pre-emptive check */
                return(.0);
                                                /* move cheapest first */
        for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
                int     d = psortrow(pm,s)[j];
                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];

                total_cost += amt * pricerow(pm,s)[d];
                amt2move -= amt;
        }
        return(total_cost);
}

/* Take a step in migration by choosing optimal bucket to transfer */
static double
migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
{
        const double    maxamt = 1./(double)pm->ncols;
        const double    minamt = maxamt*5e-6;
        double          *src_cost;
        struct {
                int     s, d;   /* source and destination */
                double  price;  /* price estimate per amount moved */
                double  amt;    /* amount we can move */
        } cur, best;
        int             i;
                                                /* allocate cost array */
        src_cost = (double *)malloc(sizeof(double)*pm->nrows);
        if (src_cost == NULL) {
                fprintf(stderr, "%s: Out of memory in migration_step()\n",
                                progname);
                exit(1);
        }
        for (i = pm->nrows; i--; )              /* starting costs for diff. */
                src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);

                                                /* find best source & dest. */
        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
        for (cur.s = pm->nrows; cur.s--; ) {
            double      cost_others = 0;

            if (src_rem[cur.s] <= minamt)
                    continue;
                                                /* examine cheapest dest. */
            for (i = 0; i < pm->ncols; i++)
                if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
                        break;
            if (i >= pm->ncols)
                break;
            if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
                continue;                       /* no point checking further */
            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
                                src_rem[cur.s] : dst_rem[cur.d];
            if (cur.amt > maxamt) cur.amt = maxamt;
            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
            for (i = pm->nrows; i--; )
                if (i != cur.s)
                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
                                        - src_cost[i];
            dst_rem[cur.d] += cur.amt;          /* undo trial move */
            cur.price += cost_others/cur.amt;   /* adjust effective price */
            if (cur.price < best.price)         /* are we better than best? */
                    best = cur;
        }
        free(src_cost);                         /* finish up */

        if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
                return(.0);
                                                /* else make the actual move */
        mtx_coef(mig,best.s,best.d) += best.amt;
        src_rem[best.s] -= best.amt;
        dst_rem[best.d] -= best.amt;
        return(best.amt);
}

/* Compute and insert migration along directed edge (may fork child) */
static MIGRATION *
create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
{
        const double    end_thresh = 5e-6;
        PRICEMAT        pmtx;
        MIGRATION       *newmig;
        double          *src_rem, *dst_rem;
        double          total_rem = 1., move_amt;
        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) (%.0f,%.0f) ",
                        get_theta180(from_rbf->invec),
                        get_phi360(from_rbf->invec));
        fprintf(stderr, "to (%.0f,%.0f) 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 */
        price_routes(&pmtx, from_rbf, to_rbf);
        src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
        dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
        if ((src_rem == NULL) | (dst_rem == NULL)) {
                fprintf(stderr, "%s: Out of memory in create_migration()\n",
                                progname);
                exit(1);
        }
                                                /* starting quantities */
        memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
        for (i = from_rbf->nrbf; i--; )
                src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
        for (j = to_rbf->nrbf; j--; )
                dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;

        do {                                    /* move a bit at a time */
                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
                total_rem -= move_amt;
        } while ((total_rem > end_thresh) & (move_amt > 0));

        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 */
        free_routes(&pmtx);                     /* free working arrays */
        free(src_rem);
        free(dst_rem);
        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 context 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, rbf0->invec);
                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);
                }
        } else 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);
                }
        }
}
        
/* Build our triangle mesh from recorded RBFs */
void
build_mesh(void)
{
        double          best2 = M_PI*M_PI;
        RBFNODE         *shrt_edj[2];
        RBFNODE         *rbf0, *rbf1;
                                                /* 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.9
Committed:	Fri Jun 28 23:18:51 2013 UTC (12 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.8:	+1 -22 lines
Log Message:	Fixed calculation of bsdf_min
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.9	static const char RCSid[] = "$Id: bsdfmesh.c,v 2.8 2013/03/20 01:00:22 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			/* number of processes to run */
22			int nprocs = 1;
23			/* number of children (-1 in child) */
24			static int nchild = 0;
25
26	greg	2.3	typedef struct {
27			int nrows, ncols; /* array size (matches migration) */
28			float price; / migration prices */
29			short sord; / sort for each row, low to high */
30			} PRICEMAT; /* sorted pricing matrix */
31
32			#define pricerow(p,i) ((p)->price + (i)*(p)->ncols)
33			#define psortrow(p,i) ((p)->sord + (i)*(p)->ncols)
34
35	greg	2.2	/* Create a new migration holder (sharing memory for multiprocessing) */
36			static MIGRATION *
37			new_migration(RBFNODE from_rbf, RBFNODE to_rbf)
38			{
39			size_t memlen = sizeof(MIGRATION) +
40			sizeof(float)(from_rbf->nrbfto_rbf->nrbf - 1);
41			MIGRATION *newmig;
42			#ifdef _WIN32
43			if (nprocs > 1)
44			fprintf(stderr, "%s: warning - multiprocessing not supported\n",
45			progname);
46			nprocs = 1;
47			newmig = (MIGRATION *)malloc(memlen);
48			#else
49			if (nprocs <= 1) { /* single process? */
50			newmig = (MIGRATION *)malloc(memlen);
51			} else { /* else need to share memory */
52			newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ\|PROT_WRITE,
53			MAP_ANON\|MAP_SHARED, -1, 0);
54			if ((void *)newmig == MAP_FAILED)
55			newmig = NULL;
56			}
57			#endif
58			if (newmig == NULL) {
59			fprintf(stderr, "%s: cannot allocate new migration\n", progname);
60			exit(1);
61			}
62			newmig->rbfv[0] = from_rbf;
63			newmig->rbfv[1] = to_rbf;
64			/* insert in edge lists */
65			newmig->enxt[0] = from_rbf->ejl;
66			from_rbf->ejl = newmig;
67			newmig->enxt[1] = to_rbf->ejl;
68			to_rbf->ejl = newmig;
69			newmig->next = mig_list; /* push onto global list */
70			return(mig_list = newmig);
71			}
72
73			#ifdef _WIN32
74			#define await_children(n) (void)(n)
75			#define run_subprocess() 0
76			#define end_subprocess() (void)0
77			#else
78
79			/* Wait for the specified number of child processes to complete */
80			static void
81			await_children(int n)
82			{
83			int exit_status = 0;
84
85			if (n > nchild)
86			n = nchild;
87			while (n-- > 0) {
88			int status;
89			if (wait(&status) < 0) {
90			fprintf(stderr, "%s: missing child(ren)!\n", progname);
91			nchild = 0;
92			break;
93			}
94			--nchild;
95			if (status) { /* something wrong */
96			if ((status = WEXITSTATUS(status)))
97			exit_status = status;
98			else
99			exit_status += !exit_status;
100			fprintf(stderr, "%s: subprocess died\n", progname);
101			n = nchild; /* wait for the rest */
102			}
103			}
104			if (exit_status)
105			exit(exit_status);
106			}
107
108			/* Start child process if multiprocessing selected */
109			static pid_t
110			run_subprocess(void)
111			{
112			int status;
113			pid_t pid;
114
115			if (nprocs <= 1) /* any children requested? */
116			return(0);
117			await_children(nchild + 1 - nprocs); /* free up child process */
118			if ((pid = fork())) {
119			if (pid < 0) {
120			fprintf(stderr, "%s: cannot fork subprocess\n",
121			progname);
122	greg	2.6	await_children(nchild);
123	greg	2.2	exit(1);
124			}
125			++nchild; /* subprocess started */
126			return(pid);
127			}
128			nchild = -1;
129			return(0); /* put child to work */
130			}
131
132			/* If we are in subprocess, call exit */
133			#define end_subprocess() if (nchild < 0) _exit(0); else
134
135			#endif /* ! _WIN32 */
136
137	greg	2.3	/* Comparison routine needed for sorting price row */
138			static int
139			msrt_cmp(void b, const void p1, const void *p2)
140			{
141			PRICEMAT pm = (PRICEMAT )b;
142			int ri = ((const short *)p1 - pm->sord) / pm->ncols;
143			float c1 = pricerow(pm,ri)[(const short )p1];
144			float c2 = pricerow(pm,ri)[(const short )p2];
145
146			if (c1 > c2) return(1);
147			if (c1 < c2) return(-1);
148			return(0);
149			}
150
151	greg	2.1	/* Compute (and allocate) migration price matrix for optimization */
152	greg	2.3	static void
153			price_routes(PRICEMAT pm, const RBFNODE from_rbf, const RBFNODE *to_rbf)
154	greg	2.1	{
155			FVECT vto = (FVECT )malloc(sizeof(FVECT) * to_rbf->nrbf);
156			int i, j;
157
158	greg	2.3	pm->nrows = from_rbf->nrbf;
159			pm->ncols = to_rbf->nrbf;
160			pm->price = (float )malloc(sizeof(float) pm->nrows*pm->ncols);
161			pm->sord = (short )malloc(sizeof(short) pm->nrows*pm->ncols);
162
163			if ((pm->price == NULL) \| (pm->sord == NULL) \| (vto == NULL)) {
164	greg	2.1	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
165			progname);
166			exit(1);
167			}
168			for (j = to_rbf->nrbf; j--; ) /* save repetitive ops. */
169			ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
170
171			for (i = from_rbf->nrbf; i--; ) {
172			const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
173			FVECT vfrom;
174			ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
175	greg	2.3	for (j = to_rbf->nrbf; j--; ) {
176	greg	2.7	double dprod = DOT(vfrom, vto[j]);
177			pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
178	greg	2.1	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
179	greg	2.3	psortrow(pm,i)[j] = j;
180			}
181			qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
182	greg	2.1	}
183			free(vto);
184			}
185
186	greg	2.3	/* Free price matrix */
187			static void
188			free_routes(PRICEMAT *pm)
189	greg	2.1	{
190	greg	2.3	free(pm->price); pm->price = NULL;
191			free(pm->sord); pm->sord = NULL;
192	greg	2.1	}
193
194			/* Compute minimum (optimistic) cost for moving the given source material */
195			static double
196	greg	2.3	min_cost(double amt2move, const double avail, const PRICEMAT pm, int s)
197	greg	2.1	{
198			double total_cost = 0;
199	greg	2.3	int j;
200	greg	2.1
201			if (amt2move <= FTINY) /* pre-emptive check */
202	greg	2.5	return(.0);
203	greg	2.1	/* move cheapest first */
204	greg	2.3	for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
205			int d = psortrow(pm,s)[j];
206	greg	2.1	double amt = (amt2move < avail[d]) ? amt2move : avail[d];
207
208	greg	2.3	total_cost += amt * pricerow(pm,s)[d];
209	greg	2.1	amt2move -= amt;
210			}
211			return(total_cost);
212			}
213
214			/* Take a step in migration by choosing optimal bucket to transfer */
215			static double
216	greg	2.3	migration_step(MIGRATION mig, double src_rem, double dst_rem, const PRICEMAT pm)
217	greg	2.1	{
218	greg	2.4	const double maxamt = 1./(double)pm->ncols;
219	greg	2.2	const double minamt = maxamt*5e-6;
220	greg	2.5	double *src_cost;
221	greg	2.1	struct {
222			int s, d; /* source and destination */
223			double price; /* price estimate per amount moved */
224			double amt; /* amount we can move */
225			} cur, best;
226			int i;
227	greg	2.5	/* allocate cost array */
228			src_cost = (double )malloc(sizeof(double)pm->nrows);
229			if (src_cost == NULL) {
230			fprintf(stderr, "%s: Out of memory in migration_step()\n",
231			progname);
232			exit(1);
233	greg	2.1	}
234	greg	2.3	for (i = pm->nrows; i--; ) /* starting costs for diff. */
235			src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
236	greg	2.1
237			/* find best source & dest. */
238			best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
239	greg	2.3	for (cur.s = pm->nrows; cur.s--; ) {
240	greg	2.1	double cost_others = 0;
241	greg	2.5
242	greg	2.2	if (src_rem[cur.s] <= minamt)
243	greg	2.1	continue;
244	greg	2.4	/* examine cheapest dest. */
245			for (i = 0; i < pm->ncols; i++)
246	greg	2.5	if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
247	greg	2.4	break;
248			if (i >= pm->ncols)
249	greg	2.5	break;
250	greg	2.4	if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
251	greg	2.5	continue; /* no point checking further */
252	greg	2.1	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
253			src_rem[cur.s] : dst_rem[cur.d];
254			if (cur.amt > maxamt) cur.amt = maxamt;
255			dst_rem[cur.d] -= cur.amt; /* add up differential costs */
256	greg	2.3	for (i = pm->nrows; i--; )
257	greg	2.1	if (i != cur.s)
258	greg	2.3	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
259	greg	2.1	- src_cost[i];
260			dst_rem[cur.d] += cur.amt; /* undo trial move */
261			cur.price += cost_others/cur.amt; /* adjust effective price */
262			if (cur.price < best.price) /* are we better than best? */
263			best = cur;
264			}
265	greg	2.5	free(src_cost); /* finish up */
266
267			if ((best.s < 0) \| (best.d < 0)) /* nothing left to move? */
268	greg	2.1	return(.0);
269	greg	2.5	/* else make the actual move */
270	greg	2.2	mtx_coef(mig,best.s,best.d) += best.amt;
271	greg	2.1	src_rem[best.s] -= best.amt;
272			dst_rem[best.d] -= best.amt;
273			return(best.amt);
274			}
275
276			/* Compute and insert migration along directed edge (may fork child) */
277			static MIGRATION *
278			create_migration(RBFNODE from_rbf, RBFNODE to_rbf)
279			{
280	greg	2.2	const double end_thresh = 5e-6;
281	greg	2.3	PRICEMAT pmtx;
282	greg	2.1	MIGRATION *newmig;
283			double src_rem, dst_rem;
284			double total_rem = 1., move_amt;
285	greg	2.6	int i, j;
286	greg	2.1	/* check if exists already */
287			for (newmig = from_rbf->ejl; newmig != NULL;
288			newmig = nextedge(from_rbf,newmig))
289			if (newmig->rbfv[1] == to_rbf)
290			return(NULL);
291			/* else allocate */
292	greg	2.7	#ifdef DEBUG
293			fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
294			get_theta180(from_rbf->invec),
295			get_phi360(from_rbf->invec));
296			fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
297			get_theta180(to_rbf->invec),
298			get_phi360(to_rbf->invec),
299			from_rbf->nrbf, to_rbf->nrbf);
300			#endif
301	greg	2.1	newmig = new_migration(from_rbf, to_rbf);
302			if (run_subprocess())
303			return(newmig); /* child continues */
304	greg	2.3	price_routes(&pmtx, from_rbf, to_rbf);
305	greg	2.1	src_rem = (double )malloc(sizeof(double)from_rbf->nrbf);
306			dst_rem = (double )malloc(sizeof(double)to_rbf->nrbf);
307			if ((src_rem == NULL) \| (dst_rem == NULL)) {
308			fprintf(stderr, "%s: Out of memory in create_migration()\n",
309			progname);
310			exit(1);
311			}
312			/* starting quantities */
313			memset(newmig->mtx, 0, sizeof(float)from_rbf->nrbfto_rbf->nrbf);
314			for (i = from_rbf->nrbf; i--; )
315			src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
316	greg	2.6	for (j = to_rbf->nrbf; j--; )
317			dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
318
319	greg	2.1	do { /* move a bit at a time */
320	greg	2.3	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
321	greg	2.1	total_rem -= move_amt;
322	greg	2.2	} while ((total_rem > end_thresh) & (move_amt > 0));
323	greg	2.6
324	greg	2.1	for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
325	greg	2.6	double nf = rbf_volume(&from_rbf->rbfa[i]);
326	greg	2.1	if (nf <= FTINY) continue;
327			nf = from_rbf->vtotal / nf;
328			for (j = to_rbf->nrbf; j--; )
329	greg	2.6	mtx_coef(newmig,i,j) = nf; / row now sums to 1.0 */
330	greg	2.1	}
331			end_subprocess(); /* exit here if subprocess */
332	greg	2.3	free_routes(&pmtx); /* free working arrays */
333	greg	2.1	free(src_rem);
334			free(dst_rem);
335			return(newmig);
336			}
337
338			/* Check if prospective vertex would create overlapping triangle */
339			static int
340			overlaps_tri(const RBFNODE bv0, const RBFNODE bv1, const RBFNODE *pv)
341			{
342			const MIGRATION *ej;
343			RBFNODE *vother[2];
344			int im_rev;
345			/* find shared edge in mesh */
346			for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
347			const RBFNODE *tv = opp_rbf(pv,ej);
348			if (tv == bv0) {
349			im_rev = is_rev_tri(ej->rbfv[0]->invec,
350			ej->rbfv[1]->invec, bv1->invec);
351			break;
352			}
353			if (tv == bv1) {
354			im_rev = is_rev_tri(ej->rbfv[0]->invec,
355			ej->rbfv[1]->invec, bv0->invec);
356			break;
357			}
358			}
359			if (!get_triangles(vother, ej)) /* triangle on same side? */
360			return(0);
361			return(vother[im_rev] != NULL);
362			}
363
364			/* Find context hull vertex to complete triangle (oriented call) */
365			static RBFNODE *
366			find_chull_vert(const RBFNODE rbf0, const RBFNODE rbf1)
367			{
368			FVECT vmid, vejn, vp;
369			RBFNODE rbf, rbfbest = NULL;
370			double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
371
372			VSUB(vejn, rbf1->invec, rbf0->invec);
373			VADD(vmid, rbf0->invec, rbf1->invec);
374			if (normalize(vejn) == 0 \|\| normalize(vmid) == 0)
375			return(NULL);
376			/* XXX exhaustive search */
377			/* Find triangle with minimum rotation from perpendicular */
378			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
379			if ((rbf == rbf0) \| (rbf == rbf1))
380			continue;
381			tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
382			if (DOT(vp, vmid) <= FTINY)
383			continue; /* wrong orientation */
384			area2 = .25*DOT(vp,vp);
385			VSUB(vp, rbf->invec, rbf0->invec);
386			dprod = -DOT(vp, vejn);
387			VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
388			dprod = DOT(vp, vmid) / VLEN(vp);
389			if (dprod <= bestdprod + FTINY(1 - 2(area2 < bestarea2)))
390			continue; /* found better already */
391			if (overlaps_tri(rbf0, rbf1, rbf))
392			continue; /* overlaps another triangle */
393			rbfbest = rbf;
394			bestdprod = dprod; /* new one to beat */
395			bestarea2 = area2;
396			}
397			return(rbfbest);
398			}
399
400			/* Create new migration edge and grow mesh recursively around it */
401			static void
402			mesh_from_edge(MIGRATION *edge)
403			{
404			MIGRATION ej0, ej1;
405			RBFNODE *tvert[2];
406
407			if (edge == NULL)
408			return;
409			/* triangle on either side? */
410			get_triangles(tvert, edge);
411			if (tvert[0] == NULL) { /* grow mesh on right */
412			tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
413			if (tvert[0] != NULL) {
414			if (tvert[0]->ord > edge->rbfv[0]->ord)
415			ej0 = create_migration(edge->rbfv[0], tvert[0]);
416			else
417			ej0 = create_migration(tvert[0], edge->rbfv[0]);
418			if (tvert[0]->ord > edge->rbfv[1]->ord)
419			ej1 = create_migration(edge->rbfv[1], tvert[0]);
420			else
421			ej1 = create_migration(tvert[0], edge->rbfv[1]);
422			mesh_from_edge(ej0);
423			mesh_from_edge(ej1);
424			}
425			} else if (tvert[1] == NULL) { /* grow mesh on left */
426			tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
427			if (tvert[1] != NULL) {
428			if (tvert[1]->ord > edge->rbfv[0]->ord)
429			ej0 = create_migration(edge->rbfv[0], tvert[1]);
430			else
431			ej0 = create_migration(tvert[1], edge->rbfv[0]);
432			if (tvert[1]->ord > edge->rbfv[1]->ord)
433			ej1 = create_migration(edge->rbfv[1], tvert[1]);
434			else
435			ej1 = create_migration(tvert[1], edge->rbfv[1]);
436			mesh_from_edge(ej0);
437			mesh_from_edge(ej1);
438			}
439			}
440			}
441
442			/* Build our triangle mesh from recorded RBFs */
443			void
444			build_mesh(void)
445			{
446			double best2 = M_PI*M_PI;
447			RBFNODE *shrt_edj[2];
448			RBFNODE rbf0, rbf1;
449			/* check if isotropic */
450			if (single_plane_incident) {
451			for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
452			if (rbf0->next != NULL)
453			create_migration(rbf0, rbf0->next);
454			await_children(nchild);
455			return;
456			}
457			shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
458			for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
459			for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
460			double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
461			if (dist2 < best2) {
462			shrt_edj[0] = rbf0;
463			shrt_edj[1] = rbf1;
464			best2 = dist2;
465			}
466			}
467			if (shrt_edj[0] == NULL) {
468			fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
469			exit(1);
470			}
471			/* build mesh from this edge */
472			if (shrt_edj[0]->ord < shrt_edj[1]->ord)
473			mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
474			else
475			mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
476			/* complete migrations */
477			await_children(nchild);
478			}