rfst.c

生成直角Steiner树的程序包

			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> y > p2 -> y) OR
			    ((p1 -> y EQ p2 -> y) AND
			     ((p1 -> x > p2 -> x) OR
			      ((p1 -> x EQ p2 -> x) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	/* Now do actual sorting.  Exchange first/last and sift down. */
	while (n > 1) {
		/* Largest is at index [0], swap with index [n-1],	*/
		/* thereby putting it into final position.		*/
		--n;
		i = index [0];
		index [0] = index [n];
		index [n] = i;

		/* Now restore the heap by sifting index [0] down. */
		j = 0;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> y > p1 -> y) OR
				    ((p2 -> y EQ p1 -> y) AND
				     ((p2 -> x > p1 -> x) OR
				      ((p2 -> x EQ p1 -> x) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> y > p2 -> y) OR
			    ((p1 -> y EQ p2 -> y) AND
			     ((p1 -> x > p2 -> x) OR
			      ((p1 -> x EQ p2 -> x) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	return (index);
}

/*
 * This routine finds all pairs of points whose coordinates are exactly
 * identical.  This is a degeneracy that can cause successive algorithms
 * much heartburn, so we take care of them immediately by keeping only
 * the earliest (lowest pnum) point, and marking the others as duplicates.
 * We generate a partition of such terminals into subsets -- if terminals
 * A and B reside in the same subset then they have identical coordinates.
 * We refer to such a subset as a Duplicate Terminal Group.
 *
 * For each terminal group we retain ONLY THE FIRST member -- the terminal
 * map bits are turned OFF for all other members of a terminal group,
 * effectively eliminating those terminals from the problem.
 */

	static
	int
generate_duplicate_terminal_groups (

struct pset *		pts,	/* IN - original terminal set */
int *			xorder, /* IN - terminal numbers sorted by X coord */
int ***			grps	/* OUT - duplicate terminal groups */
)
{
int			i;
int			j;
int			n;
int			n_grps;
struct point *		p0;
struct point *		p1;
struct point *		p2;
int *			ip;
int **			real_ptrs;
int *			real_terms;
int **			ptrs;
int *			terms;

	n = pts -> n;

	n_grps = 0;

	ptrs	= NEWA (n + 1, int *);
	terms	= NEWA (n, int);

	ip = &terms [0];
	for (i = 1; i < n; ) {
		p0 = &(pts -> a [xorder [i - 1]]);
		p1 = &(pts -> a [xorder [i]]);

		if ((p0 -> y NE p1 -> y) OR (p0 -> x NE p1 -> x)) {
			/* Not identical. */
			++i;
			continue;
		}

		/* Terminals xorder[i-1] and xorder[i] are identical. */

		for (j = i + 1; j < n; j++) {
			p2 = &(pts -> a [xorder [j]]);
			if (p0 -> y NE p2 -> y) break;
			if (p0 -> x NE p2 -> x) break;
		}
		/* j now points to the first non-equal terminal */

		/* Make a new duplicate terminal group... */
		ptrs [n_grps++] = ip;
		*ip++ = xorder [i - 1];
		while (i < j) {
			*ip++ = xorder [i++];
		}

		/* Skip whole group of coincident points and continue */
	}
	ptrs [n_grps] = ip;

	if (n_grps <= 0) {
		*grps = NULL;
	}
	else {
		/* Transfer to permanent memory of proper size... */
		real_ptrs = NEWA (n_grps + 1, int *);
		real_terms = NEWA (ip - terms, int);
		(void) memcpy ((char *) real_terms,
			       (char *) terms,
			       (ip - terms) * sizeof (int));
		for (i = 0; i <= n_grps; i++) {
			real_ptrs [i] = &real_terms [ptrs [i] - ptrs [0]];
		}
		*grps = real_ptrs;
	}

	free ((char *) terms);
	free ((char *) ptrs);

	return (n_grps);
}

/*
 * This routine removes all but the first terminal in each duplicate
 * terminal group.  We also prepare arrays for mapping forward from
 * old to new terminal numbers, and backward from new terminal numbers
 * to old.
 */

	static
	struct pset *
remove_duplicates (

struct pset *		pts,		/* IN - original point set */
int			ndg,		/* IN - number of duplicate groups */
int **			list,		/* IN - list of duplicate groups */
int **			fwd_map_ptr,	/* OUT - map to renumber old to new */
int **			rev_map_ptr	/* OUT - map to renumber new to old */
)
{
int		i;
int		j;
int		n;
int		kmasks;
int		numdel;
int		new_n;
int *		ip1;
int *		ip2;
int *		fwd;
int *		rev;
bitmap_t *	deleted;
struct pset *	newpts;

	n = pts -> n;

	kmasks = BMAP_ELTS (n);
	deleted = NEWA (kmasks, bitmap_t);
	for (i = 0; i < kmasks; i++) {
		deleted [i] = 0;
	}

	numdel = 0;
	for (i = 0; i < ndg; i++) {
		ip1 = list [i];
		ip2 = list [i + 1];

		/* Retain the first in this group, exclude all	*/
		/* of the others...				*/
		while (++ip1 < ip2) {
			if (BITON (deleted, *ip1)) {
				fatal ("remove_duplicates: Bug 1.");
			}
			++numdel;
			SETBIT (deleted, *ip1);
		}
	}

	new_n = n - numdel;

	fwd = NEWA (n, int);
	rev = NEWA (new_n, int);
	newpts = NEW_PSET (new_n);
	memset (newpts, 0, PSET_SIZE (new_n));
	newpts -> n = new_n;
	j = 0;
	for (i = 0; i < n; i++) {
		if (BITON (deleted, i)) {
			fwd [i] = -1;
		}
		else {
			newpts -> a [j].x		= pts -> a [i].x;
			newpts -> a [j].y		= pts -> a [i].y;
			newpts -> a [j].pnum		= j;
			rev [j] = i;
			fwd [i] = j;
			++j;
		}
	}

	free ((char *) deleted);

	*fwd_map_ptr = fwd;
	*rev_map_ptr = rev;

	return (newpts);
}

/*
 * Compute the RFSTs for the given set of terminals, which are now
 * guaranteed to be unique.
 */

	static
	void
compute_rfsts_for_unique_terminals (

struct rinfo *		rip	/* IN - global RFST info */
)
{
int			i;
int			n;
int			dir;
int			nedges;
struct pset *		pts;
struct edge *		ep;
struct edge *		mst_edges;
struct rlist *		rp;
dist_t			mst_len;
double			ub_shortleg [2];
char			buf1 [32];

	pts = rip -> pts;
	n = pts -> n;

	compute_successors (rip);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute Successors:     %s\n", buf1);
	}

	rip -> empty_rect = init_empty_rectangles (pts, rip -> succ [0]);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Empty Rectangles:       %s\n", buf1);
	}

	compute_ub0 (rip);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute UB0:            %s\n", buf1);
	}

	mst_edges = NEWA (n - 1, struct edge);
	nedges = rect_mst (pts, mst_edges, rip -> empty_rect);
	if (nedges NE n - 1) {
		fatal ("compute_rfsts_for_unique_terminals: Bug 1.");
	}

	mst_len = 0.0;
	ep = mst_edges;
	for (i = 0; i < nedges; ep++, i++) {
		mst_len += ep -> len;
	}
	rip -> mst_length = mst_len;

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute MST:            %s\n", buf1);
	}

	rip -> bsd = compute_bsd (nedges, mst_edges, 0);

	free ((char *) mst_edges);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute BSD:            %s\n", buf1);
	}

	compute_zt (rip);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute ZT:             %s\n", buf1);
	}

	compute_ub1 (rip);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute UB1:            %s\n", buf1);
	}

	/*
	 * Preprocessing finished.  Start growing FSTs!
	 */

	rip -> terms		= NEWA (n, int);
	rip -> longterms	= NEWA (n + 1, int);
	rip -> maxedges		= NEWA (n, dist_t);
	rip -> shortterm	= NEWA (n, int);
	rip -> lrindex		= NEWA (n, int);
	rip -> term_check	= NEWA (n, bool);
	rip -> hash		= NEWA (n, struct rlist *);

	for (i = 0; i < n; i++) {
		rip -> lrindex [i] = 0;
		rip -> term_check [i] = FALSE;
		rip -> hash [i] = NULL;
	}
	rip -> list.forw	= &(rip -> list);
	rip -> list.back	= &(rip -> list);

	rip -> fsts_checked = 0;

	ub_shortleg [0] = INF_DISTANCE;
	ub_shortleg [1] = INF_DISTANCE;
	if (MaxFSTSize EQ 0) MaxFSTSize = n;

	for (dir = 0; dir < 2; dir++) {
		for (i = 0; i < n; i++) {
			rip -> terms [0]	= i;
			rip -> maxedges [i]	= 0.0;
			/* Long leg candidate list is initially empty,	*/
			/* add candidates on demand.			*/
			rip -> longterms [0]	= i;
			rip -> longterms [1]	= -1;
			grow_RFST (rip,
				   1,		/* size */
				   0.0,		/* length */
				   dir,
				   0.0,		/* ub_length */
				   ub_shortleg,
				   0);		/* longindex */

#if DO_STATISTICS
			for (j = 0; longterms [j] >= 0; j++) {
				++long_leg_count;
			}
#endif
		}
	}

	/* Finally add MST-edges */

	ep = &(rip -> bsd -> mst_edges [1]);
	for (i = 1; i < n; i++) {
		rip -> terms [0] = ep -> p1;
		rip -> terms [1] = ep -> p2;
		test_and_save_fst (rip,
				   2,
				   ep -> len,
				   0,
				   TYPE_1);
		++ep;
	}

#if DO_STATISTICS
	fprintf (stderr, "Total FSTs checked: %d\n", rip -> fsts_checked);
#endif

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Grow FSTs:              %s\n", buf1);
	}

	/* Disconnect FSTs from hash table. */

	for (rp = rip -> list.forw;
	     rp NE &(rip -> list);
	     rp = rp -> forw) {
		rp -> next = NULL;
	}

	/* Clean up. */

	free ((char *) (rip -> hash));		rip -> hash = NULL;
	free ((char *) (rip -> term_check));	rip -> term_check = NULL;
	free ((char *) (rip -> lrindex));	rip -> lrindex = NULL;
	free ((char *) (rip -> shortterm));	rip -> shortterm = NULL;
	free ((char *) (rip -> maxedges));	rip -> maxedges = NULL;
	free ((char *) (rip -> longterms));	rip -> longterms = NULL;
	free ((char *) (rip -> terms));		rip -> terms = NULL;

	free ((char *) (rip -> ub1 [0]));
	free ((char *) (rip -> zt [0] [0]));
	free ((char *) (rip -> zt [0]));

	shutdown_bsd (rip -> bsd);

	free ((char *) (rip -> ub0 [0]));
	free ((char *) (rip -> empty_rect));	rip -> empty_rect = NULL;
	free ((char *) (rip -> succ [0]));

	memset (&(rip -> succ), 0, sizeof (rip -> succ));
	memset (&(rip -> ub0),	0, sizeof (rip -> ub0));
	memset (&(rip -> ub1),	0, sizeof (rip -> ub1));
	memset (&(rip -> zt),	0, sizeof (rip -> zt));
}

/*
 * Sort the terminals by both X and Y coordinates, and then create the
 * successor lists.  These permit us to start from a random terminal and
 * traverse all remaining terminals in order by one of the four compass
 * directions.
 */

	static
	void
compute_successors (

struct rinfo *		rip		/* The global RFST info */
)
{
int			i;
int			n;
struct pset *		pts;
int *			index;
int *			succ;

	pts = rip -> pts;
	n = pts -> n;

	/* Allocate permanent successor lists. */
	succ = NEWA (4 * n, int);
	rip -> succ [0] = succ;		succ += n;
	rip -> succ [1] = succ;		succ += n;
	rip -> succ [2] = succ;		succ += n;
	rip -> succ [3] = succ;
	rip -> succ [4] = rip -> succ [0];
	rip -> succ [5] = rip -> succ [1];
	rip -> succ [6] = rip -> succ [2];

	/* Get terminal numbers sorted by X coordinate. */
	index = rip -> x_order;

	/* Set direction 0 (due east) successors. */
	succ = rip -> succ [0];
	for (i = 1; i < n; i++) {
		succ [index [i - 1]] = index [i];
	}
	succ [index [i - 1]] = -1;

	/* Set direction 2 (due west) successors. */
	succ = rip -> succ [2];
	for (i = 1; i < n; i++) {
		succ [index [i]] = index [i - 1];
	}
	succ [index [0]] = -1;

	/* Get terminal numbers sorted by Y coordinate. */
	index = rip -> y_order;

	/* Set direction 1 (due north) successors. */
	succ = rip -> succ [1];
	for (i = 1; i < n; i++) {
		succ [index [i - 1]] = index [i];
	}
	succ [index [i - 1]] = -1;

	/* Set direction 3 (due south) successors. */
	succ = rip -> succ [3];
	for (i = 1; i < n; i++) {
		succ [index [i]] = index [i - 1];
	}
	succ [index [0]] = -1;
}

/*
 * Compute the UB0 bounds.
 */

	static
	void
compute_ub0 (

struct rinfo *		rip		/* IN/OUT - global RFST info */
)
{
int			i;
int			j;
int			n;
int			dir;
struct pset *		pts;
struct point *		p1;
struct point *		p2;
int *			succ;
dist_t *		array;
dstdiroff_t		dir_off;
dstdiroff_t		dirp_off;
double			bound;
double			d1, d2, d3;

	pts = rip -> pts;
	n = pts -> n;

	array = NEWA (4 * n, dist_t);
	rip -> ub0 [0]	= array;	array += n;
	rip -> ub0 [1]	= array;	array += n;
	rip -> ub0 [2]	= array;	array += n;
	rip -> ub0 [3]	= array;
	rip -> ub0 [4]	= rip -> ub0 [0];
	rip -> ub0 [5]	= rip -> ub0 [1];
	rip -> ub0 [6]	= rip -> ub0 [2];

	for (dir = 0; dir < 4; dir++) {
		array =	 rip -> ub0 [dir];
		succ = rip -> succ [dir];
		dir_off = DSTDIR_GETOFFSET (dir);
		dirp_off = DSTDIR_GETOFFSET (dir + 1);

		p1 = &(pts -> a [0]);
		for (i = 0; i < n; i++, p1++) {
			bound = INF_DISTANCE;
			for (j = succ [i]; j >= 0; j = succ [j]) {
				p2 = &(pts -> a [j]);
				d1 = DSTDIR_OFFSET (p1, p2, dir_off);
				if (d1 > bound) break;
				d2 = DSTDIR_OFFSET (p1, p2, dirp_off);
				if (d1 > d2) {
					d3 = d1 + d2;	/* RDIST (p1,p2) */
					if (d3 < bound) {
						bound = d3;
					}
				}
			}
			array [i] = bound;
		}
	}
}

/*
 * Find the short leg candidates for each terminal i and direction.
 */

	static
	void
compute_zt (

struct rinfo *		rip		/* IN/OUT - global RFST info */
)
{
int			i;
int			j;
int			n;
int			dir;
int			dirp;
int			lr;
int			index;
struct pset *		pts;
struct point *		p1;
struct point *		p2;
int *			succ;
int *			ip;
int *			ip1;
int **			zt;
dstdiroff_t		dir_off;
dstdiroff_t		dir1_off;
lrdiroff_t		lrdir_off;
dist_t *		ub0;
dist_t *		ub0p;
dist_t			d1, d2;
dist_t			limit;
dist_t			b;
struct ibuf *		root;
struct ibuf *		ibp;
struct ibuf **		hookp;
int *			bufp;
int			bleft;
int *			offset;
int *			op;

static int		lr_table [] = {0, 1, 1, 0};
static int		dirp_table [] = {3, 2, 3, 2};


	pts = rip -> pts;
	n = pts -> n;

	/* Compute the short leg candidate lists.  We save them in a	*/
	/* dynamic buffer first.  Later we massage the data structure	*/
	/* into something more appropriate for access.			*/

	root = NULL;
	hookp = &root;
	bufp = NULL;
	bleft = 0;
	index = 0;
	offset = NEWA (4 * (n + 1), int);
	op = offset;

	for (dir = 0; dir < 4; dir++) {
		lr	= lr_table [dir];
		dirp	= dirp_table [dir];

		succ = rip -> succ [dir];
		dir_off = DSTDIR_GETOFFSET (dir);
		dir1_off = DSTDIR_GETOFFSET (dir + 1);
		lrdir_off = LRINDEX_GETOFFSET (dir);
		ub0 = rip -> ub0 [dir];
		ub0p = rip -> ub0 [dirp];

		p1 = &(pts -> a [0]);
		for (i = 0; i < n; i++, p1++) {
			*op++ = index;
			limit = ub0 [i];
			for (j = succ [i]; j >= 0; j = succ [j]) {
				p2 = &(pts -> a [j]);
				d1 = DSTDIR_OFFSET (p1, p2, dir_off);
				if (d1 == 0.0) continue;
				if (d1 > limit) break;

				d2 = DSTDIR_OFFSET (p1, p2, dir1_off);
				if (d2 EQ 0.0) break;
				if (LRINDEX_OFFSET (p1, p2, lrdir_off) NE lr) {
					continue;
				}
				if (d2 > ub0p [j]) continue;

				b = bsd (rip -> bsd, i, j);
				if (d1 > b) continue;
				if (d2 > b) continue;
				if (is_empty_rectangle (rip -> empty_rect, i, j)) {
					/* j is a short leg candidate for i */
					/* Store it away! */
					if (bleft <= 0) {
						ibp = (struct ibuf *)
							new (offsetof (struct ibuf,
								       buf [n]));
						ibp -> next = NULL;
						ibp -> count = 0;
						*hookp = ibp;