rfst.c

生成直角Steiner树的程序包

						hookp = &(ibp -> next);
						bufp = &(ibp -> buf [0]);
						bleft = n;
					}
					*bufp++ = j;
					++(ibp -> count);
					--bleft;
					++index;
				}
			}
		}
		*op++ = index;
	}

	/* Copy list data into nice contiguous memory. */

	ip = NEWA (index, int);
	ip1 = ip;
	for (ibp = root; ibp NE NULL; ibp = ibp -> next) {
		bufp = &(ibp -> buf [0]);
		for (i = 0; i < ibp -> count; i++) {
			*ip1++ = *bufp++;
		}
	}

	if (ip1 NE ip + index) {
		fatal ("compute_zt: Bug 1.");
	}

	/* Allocate arrays of pointers into the lists. */

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

	op = offset;
	for (dir = 0; dir < 4; dir++) {
		zt = rip -> zt [dir];
		for (i = 0; i <= n; i++, p1++) {
			zt [i] = &(ip [*op++]);
		}
	}
	if (rip -> zt [0] [0] + index NE rip -> zt [3] [n]) {
		fatal ("compute_zt: Bug 2.");
	}

	while (root NE NULL) {
		ibp = root;
		root = ibp -> next;
		free ((char *) ibp);
	}

	free ((char *) offset);
}

/*
 * Compute the upper bounds UB1.
 */

	static
	void
compute_ub1 (

struct rinfo *		rip		/* IN/OUT - global RFST info */
)
{
int			i;
int			j;
int			n;
int			dir;
int			lr;
int			last;
struct pset *		pts;
struct point *		p1;
struct point *		p2;
struct point *		p3;
int *			succ;
int *			ip1;
int *			ip2;
int **			zt;
dstdiroff_t		dir_off;
dstdiroff_t		dirp_off;
lrdiroff_t		lrdir_off;
dist_t *		array;
dist_t *		ub1;
dist_t			d1, d2, d3, d4;
dist_t			bound;
struct point		s;

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

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

	array = NEWA (4 * n, dist_t);

	rip -> ub1 [0]	= array;	array += n;
	rip -> ub1 [1]	= array;	array += n;
	rip -> ub1 [2]	= array;	array += n;
	rip -> ub1 [3]	= array;
	rip -> ub1 [4]	= rip -> ub1 [0];
	rip -> ub1 [5]	= rip -> ub1 [1];
	rip -> ub1 [6]	= rip -> ub1 [2];

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

		ub1 = rip -> ub1 [dir];
		zt = rip -> zt [dir];

		succ = rip -> succ [dir];
		dir_off = DSTDIR_GETOFFSET (dir);
		dirp_off = DSTDIR_GETOFFSET (dir + 1);
		lrdir_off = LRINDEX_GETOFFSET (dir);

		p1 = &(pts -> a [0]);
		for (i = 0; i < n; i++, p1++) {
			ip1 = zt [i];
			ip2 = zt [i + 1];
			if (ip1 >= ip2) {
				/* No short leg candidates.  UB1 is zero! */
				ub1 [i] = 0.0;
				continue;
			}

			/* Get last short leg candidate in this direction. */
			last = ip2 [-1];

			bound = INF_DISTANCE;
			p3 = &(pts -> a [last]);
			/* Get Steiner point for this last terminal. */

			SPOINT (&s, p1, p3, dir);

			d3 = DSTDIR_OFFSET (p1, p3, dirp_off);
			for (j = succ [last]; j >= 0; j = succ [j]) {
				p2 = &(pts -> a [j]);
				d1 = DSTDIR_OFFSET (&s, p2, dir_off);
				if (d1 > bound) break;

				d2 = DSTDIR_OFFSET (&s, p2, dirp_off);
				if ((LRINDEX_OFFSET (p1, p2, lrdir_off) EQ lr)
				    AND
				    (d3 > d2)) {
					bound = d1;
					break;
				}
				if (d1 > d2) {
					d4 = d1 + d2; /* RDIST(&s, p2) */
					if (d4 < bound) {
						bound = d4;
					}
				}
			}
			ub1 [i] = DSTDIR_OFFSET (p1, p3, dir_off) + bound;
		}
	}
}

/*
 * Recursively grow a rectilinear FST.
 */

	static
	void
grow_RFST (

struct rinfo *	rip,		/* IN/OUT - The global RFST info */
int		size,		/* IN - number of terms in partial tree */
dist_t		length,		/* IN - length of partial tree */
int		dir,		/* IN - growth direction from root */
dist_t		ub_length,	/* IN - upper bound on partial tree length */
dist_t *	ub_shortleg,	/* IN - left/right short leg upper bounds */
int		longindex	/* IN - current long leg terminal */
)
{
int			i, j, k, l, j2;
int			lr;
int			r;
int			lastlr;
int			dirp;
int			lp;
bool			pass_bsd;
bool			need_to_restore;
struct pset *		pts;
struct point *		p;
struct point *		q;
struct point *		root;
struct point *		last;
int *			terms;
int *			succ;
int *			ip1;
int *			ip2;
dist_t			dstdir, dstdirp;
dist_t			d1, d2, d3;
dist_t			ub1;
dist_t			b;
dist_t			max_backbone;
dist_t			min_bsd;
dist_t			new_ub_length;
dist_t			new_ub_shortleg [2];
dist_t			lastdstdirp;
dist_t			long_leg_max;

#define LONGTERMS	rip -> longterms

	terms		= rip -> terms;
	r		= terms [0];
	l		= terms [size - 1];
	lastlr		= rip -> lrindex [l];
	succ		= rip -> succ [dir];

	pts		= rip -> pts;

	root = &(pts -> a [r]);
	last = &(pts -> a [l]);

	max_backbone = INF_DISTANCE;

	lastdstdirp = DSTDIRP (root, last, dir);

	/* We set this flag whenever we change one of the "maxedges" values. */
	need_to_restore = FALSE;

	for (;;) {
		i = LONGTERMS [++longindex];
		if (i < -1) {
			/* No more candidates, and we have already	*/
			/* determined that no more can be found.	*/
			break;
		}
		if (i EQ -1) {
			/* Dynamically add next candidate to longterms. */
			for (i = succ [LONGTERMS [longindex - 1]];
			     i >= 0;
			     i = succ [i]) {
				p = &(pts -> a [i]);
				dstdirp = DSTDIRP (root, p, dir);
				if (dstdirp EQ 0.0) {
					lr = 2;
					rip -> lrindex [i] = 2;
					rip -> shortterm [i] = 0;
					LONGTERMS [longindex] = i;
					LONGTERMS [longindex + 1] = -2;
					break;
				}

				lr = LRINDEX (root, p, dir);
				if (lr EQ 0) {
					dirp = (dir - 1) & 0x03;
				}
				else {
					dirp = (dir + 1) & 0x03;
				}

				/* Find short leg candidate (if it exists) */
				j = -1;
				ip1 = rip -> zt [dirp] [i];
				ip2 = rip -> zt [dirp] [i + 1];
				while (ip1 < ip2) {
					k = *--ip2;
					if (LRINDEX (root, &(pts -> a [k]), dir)
					    EQ lr) {
						j = k;
						break;
					}
				}
				rip -> shortterm [i] = j;

				/* Check upper bounds */
				ub1 = 0.0;
				if (j >= 0) {
					ub1 = rip -> ub1 [dirp] [i];
				}
				d1 = rip -> ub0 [dirp] [i];
				if (d1 < ub1) {
					d1 = ub1;
				}
				if (dstdirp > d1) continue;

				rip -> lrindex [i] = lr;
				LONGTERMS [longindex]	  = i;
				LONGTERMS [longindex + 1] = -1;
				break;
			}
			if (i < 0) {
				/* This long leg has no further candidates! */
				/* Mark it so we never try to do find any   */
				/* more candidates! */
				LONGTERMS [longindex] = -2;
				break;
			}
		}
		else {
			/* Next long leg candidate available in longterms! */
			p = &(pts -> a [i]);
			lr = rip -> lrindex [i];
			dstdirp = DSTDIRP (root, p, dir);
		}

		dstdir	 = DSTDIR (last, p, dir);

		/* Check if consecutive terminals share Steiner point. */
		if ((size >= 3) AND (dstdir EQ 0.0)) continue;

		/* Update maximum backbone length using empty diamond */
		/* property. */
		if (dstdirp < dstdir) {
			d1 = dstdir + dstdirp;
			if (d1 < max_backbone) {
				max_backbone = d1;
			}
		}

		/* Update maximum backbone length using empty corner	*/
		/* rectangle property.					*/
		if ((size >= 2) AND (lr EQ lastlr) AND (dstdirp < lastdstdirp)) {
			if (dstdir < max_backbone) {
				max_backbone = dstdir;
			}
		}

		/* Check length of new backbone segment */
		if (dstdir > max_backbone) break;

		if (lr EQ 2) {
			/* Terminal is ON the backbone!	 Save	*/
			/* immediately as type (i) and break.	*/
			if (size >= 2) {
				terms [size] = i;
				test_and_save_fst (rip,
						   size + 1,
						   length + dstdir + dstdirp,
						   dir,
						   TYPE_1);
			}
			break;
		}

		/* Terminal on the wrong side of the long leg? */
		if ((size >= 2) AND (lr EQ lastlr)) continue;

		/* Empty rectangle with last terminal? */
		if (NOT is_empty_rectangle (rip -> empty_rect, l, i)) continue;

		/* Check if new backbone segment is longer than any BSD. */
		pass_bsd = TRUE;
		min_bsd	 = INF_DISTANCE;
		for (j = 0; j < size; j++) {
			k = terms [j];
			d1 = rip -> maxedges [k];
			if (dstdir > d1) {
				d1 = dstdir;
				rip -> maxedges [k] = d1;
				need_to_restore = TRUE;
			}
			b = bsd (rip -> bsd, i, k);
			if (d1 > b) {
				pass_bsd = FALSE;
				break;
			}
			if (b < min_bsd) {
				min_bsd = b;
			}
		}
		if (NOT pass_bsd) continue;
		new_ub_length = ub_length + min_bsd;

		/* dirp = (lr EQ 0) ? dir+3 : dir+1; */
		dirp = (dir + (lr + lr - 1)) & 3;

		/* Try to make a type (ii) FST */

		j = rip -> shortterm [i];
		if (j < 0) goto try_typeI;	/* no short leg candidate! */

		/* Is backbone rectangle empty? */
		if (NOT is_empty_rectangle (rip -> empty_rect, r, j)) goto try_typeI;
		/* Check UB1. */
		if (dstdirp > rip -> ub1 [dirp] [i]) goto try_typeI;

		/* Check BSD for each terminal in current tree */
		q = &(pts -> a [j]);
		for (j2 = 0; j2 < size; j2++) {
			k = terms [j2];
			d1 = rip -> maxedges [k];
			d2 = DSTDIRP (root, q, dir);
			d3 = DSTDIRP (p, q, dir);
			if (d2 > d1) {
				d1 = d2;
			}
			/* d1 = MAX (maxedges [k], DSTDIRP (root, q, dir)) */
			if (d1 > d3) {
				d3 = d1;
			}
			/* d3 = MAX (d1, DSTDIRP (p, q, dir)) */
			if (d3 > bsd (rip -> bsd, i, k)) goto try_typeI;
			d3 = DSTDIR (p, q, dir);
			if (d1 > d3) {
				d3 = d1;
			}
			/* d3 = MAX (d1, DSTDIR (p, q, dir)) */
			if (d3 > bsd (rip -> bsd, j, k)) goto try_typeI;
		}

#if UB_SHORTLEG
		/* Check short leg upper bound */
		if (DSTDIRP (root, q, dir) > ub_shortleg [lr]) goto try_growing;
#endif

		/* Perform FST specific tests and save type (ii) tree */
		terms [size]		= i;
		terms [size + 1]	= j;
		test_and_save_fst (rip,
				   size + 2,
				   length + DSTDIR (last, q, dir)
					  + DSTDIRP (root, p, dir),
				   dir,
				   TYPE_2);

		/* Try to make a type (i) FST */
try_typeI:

		/* Check UB0 */
		if (dstdirp > rip -> ub0 [dirp] [i]) continue;

		/* Check BSD to each terminal in previous tree */
		pass_bsd = TRUE;
		for (j2 = 0; j2 < size; j2++) {
			k = terms [j2];
			d1 = rip -> maxedges [k];
			if (dstdirp > d1) {
				d1 = dstdirp;
			}
			if (d1 > bsd (rip -> bsd, k, i)) {
				pass_bsd = FALSE;
				break;
			}
		}
		if (NOT pass_bsd) continue;

		/* Check if BSD upper bound is shorter than partial tree */
		/* (including connecting Steiner point).		 */
		if (length + dstdir > new_ub_length) continue;

		/* Check if BSD upper bound is shorter than partial tree */
		if (length + dstdir + dstdirp > new_ub_length) goto try_growing;

		/* Do not make 2-terminal FSTs (MST edges). */
		if (size <= 1) goto try_growing;

		/* Is backbone rectangle empty? */
		if (NOT is_empty_rectangle (rip -> empty_rect, r, i)) goto try_growing;

#if UB_SHORTLEG
		/* Check short leg upper bound */
		if (dstdirp > ub_shortleg [lr]) goto try_growing;
#endif

		/* Perform FST specific tests and save type (i) tree */
		terms [size] = i;
		new_ub_length = test_and_save_fst (rip,
						   size + 1,
						   length + dstdir + dstdirp,
						   dir,
						   TYPE_1);

		/* Try to grow the current tree. */

try_growing:
		/* Should we generate larger FSTs? */
		if (size >= MaxFSTSize) continue;

		/* Upper bound (A). */
		d1 = ub_shortleg [lr];
		if (dstdirp < d1) {
			d1 = dstdirp;
		}
		new_ub_shortleg [lr] = d1;

		d1 = ub_shortleg [1-lr];
		if (size >= 2) {
			/* Upper bound (B). */
			d2 = rip -> ub0 [dirp] [i];
			if (min_bsd < d2) {
				d2 = min_bsd;
			}
			d2 -= dstdirp;
			if (d2 < d1) {
				d1 = d2;
			}

			/* Upper bound (C). */
			if (dstdir < d1) {
				d1 = dstdir;
			}

			if (size >= 3) {
				/* Upper bound (D). */
				lp = terms [size - 2];
				d2 = DSTDIRP (root, &(pts -> a [lp]), dir);
				if (dstdirp < d2) {
					d2 = dstdirp;
				}
				d2 = DSTDIR (&(pts -> a [lp]), p, dir) - d2;
				if (d2 < d1) {
					d1 = d2;
				}
			}
		}
		new_ub_shortleg [1-lr] = d1;

		terms [size] = i;
		rip -> maxedges [i] = dstdirp;
		grow_RFST (rip,
			   size + 1,
			   length + dstdir + dstdirp,
			   dir,
			   new_ub_length,
			   new_ub_shortleg,
			   longindex);
	}

	if (need_to_restore) {
		/* Restore caller's maxedges values.  We do this by	*/
		/* recomputing them from scratch in a scan back toward	*/
		/* the root.						*/
		long_leg_max = 0.0;
		for (i = size - 1; i > 0; i--) {
			k = terms [i];
			p = &(pts -> a [k]);
			d1 = DSTDIRP (root, p, dir);
			if (long_leg_max > d1) {
				d1 = long_leg_max;
			}
			rip -> maxedges [k] = d1;
			j = terms [i - 1];
			q = &(pts -> a [j]);
			d1 = DSTDIR (q, p, dir);
			if (d1 > long_leg_max) {
				long_leg_max = d1;
			}
		}
		rip -> maxedges [r] = long_leg_max;
	}
}

/*
 * This routine performs all of the FST specific screening tests.
 * If all are passed, the FST is saved.
 */

	static
	dist_t
test_and_save_fst (

struct rinfo *	rip,		/* IN/OUT - The global RFST info */
int		size,		/* IN - number of terminals in RFST */
dist_t		length,		/* IN - length of RFST */
int		dir,		/* IN - backbone growth direction from root */
int		type		/* IN - type of tree */
)
{
int			i, j, k;
int			z12, z13, z23;
int			nstein;
int			nedges;
int			last;
int *			terms;
struct pset *		pts;
struct point *		p1;
struct point *		p2;
struct point *		p3;
struct point *		p4;
struct rlist *		rp;
struct rlist **		hookp;
struct rlist *		rp1;
struct rlist *		rp2;
struct pset *		new_terms;
struct pset *		new_steiners;
dist_t			b;
dist_t			d1;
struct full_set *	fsp;
struct edge *		edges;
struct point		s1;
struct point		s2;

	++(rip -> fsts_checked);

	terms = rip -> terms;

	/* Is this FST too large? */
	if (size > MaxFSTSize) return (length);

	if (size > 2) {
		b = bmst_terms_length (terms, size, rip -> bsd);
		if (length >= b) return (b);
	}

	/* Simple duplicate tests for size 3 and 4 */
	if (dir EQ 1) {
		if (size EQ 3) return (length);
		if ((size EQ 4) AND (type NE TYPE_1)) return (length);
	}

	pts = rip -> pts;

	if (size > 4) {
		/* No two terminals on the long leg should share a	*/
		/* Steiner point.					*/

		for (i = 1; i < size; i++) {
			d1 = DSTDIR (&(pts -> a [terms [i-1]]),
				     &(pts -> a [terms [i]]),
				     dir);
			if (d1 EQ 0.0) return (length);
		}
	}
	else if (size EQ 4) {
		p1 = &(pts -> a [terms [0]]);
		p2 = &(pts -> a [terms [1]]);
		if (DSTDIR (p1, p2, dir) EQ 0.0) return (length);
		p3 = &(pts -> a [terms [2]]);
		p4 = &(pts -> a [terms [3]]);
		if (DSTDIR (p3, p4, dir) EQ 0.0) return (length);

		if (DSTDIR (p2, p3, dir) EQ 0.0) {
			if (DSTDIRP (p1, p4, dir) NE 0.0) {
				return (length);
			}
			type = CROSS;
		}
	}
	else if (size EQ 3) {
		/* Make sure that 3-terminal FST is not degenerate */
		p1 = &(pts -> a [terms [0]]);
		p2 = &(pts -> a [terms [1]]);
		p3 = &(pts -> a [terms [2]]);
		z12 = (DSTDIR  (p1, p2, dir) EQ 0.0) + (DSTDIRP (p1, p2, dir) EQ 0.0);
		z13 = (DSTDIR  (p1, p3, dir) EQ 0.0) + (DSTDIRP (p1, p3, dir) EQ 0.0);
		z23 = (DSTDIR  (p2, p3, dir) EQ 0.0) + (DSTDIRP (p2, p3, dir) EQ 0.0);
		if (z12 + z13 > 1) return (length);
		if (z12 + z23 > 1) return (length);
		if (z13 + z23 > 1) return (length);
	}

#if EMPTY_DIAMOND_PROPERTY

	/* Check empty diamond property for transformed FST */
	i = 0;
	last = size - 1;
	if (type EQ TYPE_2) {
		last = size - 2;
		if ((size & 1) EQ 0) {
			i = 1;		/* type (ii), even */
		}
	}
	else if ((size & 1) NE 0) {
		i = 1;			/* type (i), odd */
	}
	p1 = &(pts -> a [terms [0]]);
	p2 = &(pts -> a [terms [size-1]]);
	while (i < last) {
		/* Check skew diamond... */
		SPOINT (&s1, p1, &(pts -> a [terms [i]]), dir);
		SPOINT (&s2, p2, &(pts -> a [terms [i+1]]), dir);
		if (NOT diamond_empty (rip, &s1, &s2, terms [i], dir)) {
			return (length);
		}
		++i;
		if (i >= last) break;

		/* Check diamond for flat segment... */
		SPOINT (&s1, p2, &(pts -> a [terms [i]]), dir);
		SPOINT (&s2, p2, &(pts -> a [terms [i+1]]), dir);
		if (NOT diamond_empty (rip, &s1, &s2, terms [i], dir)) {
			return (length);
		}
		++i;
	}