greedy.c

生成直角Steiner树的程序包

/***********************************************************************

	File:	greedy.c
	Rev:	a-1
	Date:	01/31/2000

	Copyright (c) 2000, 2001 by Martin Zachariasen

************************************************************************

	Implementation of greedy O(n^2) heuristic for
	computing Euclidean Steiner trees
	(Algorithmica 25, 418-437, 1999)

************************************************************************

	Modification Log:

	a-1:	01/31/2000	martinz
		: Created.  Derived from SLL heuristic

************************************************************************/

#include "bsd.h"
#include "dsuf.h"
#include "efuncs.h"
#include "steiner.h"

#define ANSI_DECLARATORS
#define REAL coord_t
#include "triangle.h"


/*
 * Global Routines
 */

dist_t		greedy_heuristic (struct pset *, struct bsd *);


/*
 * Local Types
 */

struct greedy_edge {
	dist_t		len;
	dist_t		bsd;
	pnum_t		p1;
	pnum_t		p2;
	bool		mst;
};

struct greedy_fst {
	dist_t		len;	/* Length of FST */
	dist_t		ratio;	/* Ratio of length to length of MST spanning FST-terms */
	int		count;	/* Number of terminals */
	pnum_t		p [4];	/* List of terminals */
	bool		chosen; /* Is this FST chosen by the algorithm? */
};


/*
 * Local Routines
 */

static void		add_fst (struct pset *		pts,
				 int *			terms,
				 int			term_count,
				 struct greedy_fst *	greedy_fsts,
				 int *			fst_count,
				 dist_t			mst_l,
				 bool			mst_connected);
static int		comp_edges (const void *, const void *);
static int		comp_fsts  (const void *, const void *);
static dist_t		fst_length (struct pset *, int *, int);
static dist_t		mst_length (struct pset *, int *, int);


/*
 * Local Variables
 */

static struct greedy_edge *	greedy_edges;
static struct greedy_fst  *	greedy_fsts;

/*
 * Compute the length of a shortest full Steiner tree
 * for a set of terminals with up to 4 terminals.
 *  Special fast version which only computes length.
 * Returns 0.0 if no Steiner tree exists.
 * Assume terminals are ordered as they appear on Steiner polygon.
 */

	static
	dist_t
fst_length (

struct pset *		pts,		/* IN - terminal list */
int *			terms,		/* IN - indices of terminals to consider */
int			term_count	/* IN - number of terminals */
)
{
int			i;
struct point *		a;
struct point *		b;
struct point *		c;
struct point *		d;
struct point		e, ctr, e_ad, e_cb;
dist_t			l;
dist_t			min_length;

	if (term_count EQ 2) {
		return (EDIST (&(pts -> a [terms [0]]),
			       &(pts -> a [terms [1]])));
	}

	if (term_count EQ 3) {
		a = &(pts -> a [terms [0]]);
		b = &(pts -> a [terms [1]]);
		c = &(pts -> a [terms [2]]);

		eq_point (a, b, &e);
		eq_circle_center (a, b, &e, &ctr);

		if (right_turn (&e, a, c) AND
		    left_turn (&e, b, c) AND
		    (sqr_dist (&ctr, c) > sqr_dist (&ctr, a))) {
			return EDIST (&e, c);
		}
	}

	if (term_count EQ 4) {
		min_length = INF_DISTANCE;
		for (i = 0; i <= 1; i++) {
			/* Using Lemma 5.2 p. 64 in Hwang, Richards, Winter */
			a = &(pts -> a [terms [i]]);
			d = &(pts -> a [terms [i+1]]);
			c = &(pts -> a [terms [i+2]]);
			b = &(pts -> a [terms [(i+3) % 4]]);

			/* Find intersetion point between ac and bd.
			   It is the same in both iterations */

			if ((i EQ 0) AND
			    NOT segment_intersection (a, c, b, d, &ctr)) break;

			eq_point (a, d, &e_ad);
			eq_point (c, b, &e_cb);

			if (NOT wedge120 (d, a, &e_cb) AND
			    NOT wedge120 (&e_cb, d, a) AND
			    NOT wedge120 (&e_ad, b, c) AND
			    NOT wedge120 (b, c, &e_ad) AND
			    NOT wedge120 (a, &ctr, d)) {
				l = EDIST (&e_ad, &e_cb);
				if (l < min_length) {
					min_length = l;
				}
			}
		}
		if (min_length < INF_DISTANCE) return (min_length);
	}

	return (0.0);
}

/*
 * Compute the length of a MST for a set of terminals.
 * Simply call general procedure (which might not be
 * the most effective to thing to do)
 */

	static
	dist_t
mst_length (

struct pset *		pts,		/* IN - terminal list */
int *			terms,		/* IN - indices of terminals to consider */
int			term_count	/* IN - number of terminals */
)
{
int			t;
struct pset *		tpts;
dist_t			mst_l;

	tpts = NEW_PSET (term_count);
	tpts -> n = term_count;
	for (t = 0; t < term_count; t++) {
		tpts -> a[t] = pts -> a[terms[t]];
	}
	mst_l  = euclidean_mst_length(tpts);
	free (tpts);
	return mst_l;
}

/*
 * For sorting edges by length (should be changed!!)
 */

	static
	int
comp_edges (

const void *		p1,
const void *		p2
)
{
dist_t			l1;
dist_t			l2;

	l1 = greedy_edges [ *(int*) p1].len;
	l2 = greedy_edges [ *(int*) p2].len;

	if (l1 < l2) return (-1);
	if (l1 > l2) return (1);
	return (0);
}


/*
 * For sorting FSTs by ratio (should be changed!!)
 */

	static
	int
comp_fsts (

const void *		p1,
const void *		p2
)
{
dist_t			l1;
dist_t			l2;

	l1 = greedy_fsts [ *(int*) p1].ratio;
	l2 = greedy_fsts [ *(int*) p2].ratio;

	if (l1 < l2) return (-1);
	if (l1 > l2) return (1);
	return (0);
}

/*
 * Add generated FST to list of FSTs
 */

	static
	void
add_fst (

struct pset *	pts,		/* IN - terminal list */
int *		terms,		/* IN - indices of terminals in FST */
int		term_count,	/* IN - number of terminals */
struct greedy_fst * greedy_fsts,/* IN/OUT - list of FSTs */
int *		fst_count,	/* IN/OUT - current FST count */
dist_t		mst_l,		/* Length of MST spanning FST-terms */
bool		mst_connected	/* Is induced subgraph of MST connected? */
)
{
int		j;
dist_t		smt_l;
dist_t		ratio;

	smt_l = fst_length (pts, terms, term_count);

	if (smt_l > 0) {
		ratio = smt_l / mst_l;
		if (ratio < 1.0) {
			for (j = 0; j < 4; j++) {
				greedy_fsts [*fst_count].p [j] = terms [j];
			}
			greedy_fsts [*fst_count].len	= smt_l;
			greedy_fsts [*fst_count].ratio	= ratio;
			greedy_fsts [*fst_count].count	= term_count;
			greedy_fsts [*fst_count].chosen = mst_connected;
			++(*fst_count);
		}
	}
}

/*
 * Greedy heuristic by Zachariasen and Winter
 */

	dist_t
greedy_heuristic (

struct pset *	ptss,	/* IN - terminals list for which
				heuristic tree should be constructed */
struct bsd *	bsdp	/* IN - BSD data structure */
)
{
int			i, j, k, jj, ei, fi, fk;
int			ni, nj, nt, np1, np2, p1, p2, p4;
int			fst_count, mst_count1, mst_count2;
int			neighbour_edge, neighbour_edge_idx;
int			term_count;
int			root [4];
int			terms [4];
int *			triangleedges;
int *			sortededges;
int *			sortedfsts;
bool *			new_chosen;
bool			convex_region;
bool			in_same_block;
dist_t			mst_l, mst_l1, total_mst_l;
dist_t			total_smt_l, best_smt_l, l[5], mx, my;
struct point		minp, maxp;
struct dsuf		mstsets, fstsets;
struct triangulateio	in, out, vorout;
struct pset *		pts;

	/* Special cases */

	if (ptss -> n <= 1) return (0.0);
	if (ptss -> n EQ 2) return (EDIST (&(ptss -> a [0]),
					   &(ptss -> a [1])));

	/* Compute mean of all terminals and translate in order */
	/* to improve the precision of computed eq-points.	*/
	mx = 0.0;
	my = 0.0;
	for (i = 0; i < ptss -> n; i++) {
		mx += ptss -> a[i].x;
		my += ptss -> a[i].y;
	}
	mx = floor(mx / ((dist_t) ptss -> n));
	my = floor(my / ((dist_t) ptss -> n));

	/* Transfer to triangulate data structure and call triangle */

	in.numberofpoints = ptss -> n;
	in.pointlist = NEWA (2*in.numberofpoints, REAL);

	pts = NEW_PSET(ptss -> n);
	pts -> n = ptss -> n;
	for (i = 0; i < pts -> n; i++) {
		pts -> a[i].x	 = ptss -> a[i].x - mx;
		pts -> a[i].y	 = ptss -> a[i].y - my;
		pts -> a[i].pnum = ptss -> a[i].pnum;
		in.pointlist [2*i  ] = pts -> a [i].x;
		in.pointlist [2*i+1] = pts -> a [i].y;