efst.c

生成直角Steiner树的程序包

 * Merge two disjoint ordered lists of terminal numbers.  The result
 * is of course also ordered.
 */

	static
	eterm_t *
merge_terminal_lists (

struct einfo * eip,  /* IN/OUT - global EFST info */
struct eqp_t * eqpi, /* IN - first eq-point */
struct eqp_t * eqpj, /* IN - second eq-point */
struct eqp_t * eqpk  /* IN - new eq-point */
)
{
	eterm_t *p1, *endp1, *p2, *endp2, *Zp, *eqpZ_old;
	int t1, t2;
	struct eqp_t *eqpt;

	/* Set new eq-point terminal list pointer */
	eqpk -> Z = eip -> eqpZ_curr;

	if (eip -> eqpZ_curr + eqpi -> S + eqpj -> S
		>  eip -> eqpZ + eip -> eqpZ_size)  {

		/* Terminal list space exhausted - double array */

		if (Print_Detailed_Timings)
			fprintf(stderr, "- doubling terminal list array\n");
		eqpZ_old = eip -> eqpZ;
		eip -> eqpZ = NEWA ( eip -> eqpZ_size * 2, eterm_t );
		memcpy ( eip -> eqpZ, eqpZ_old, eip -> eqpZ_size * sizeof(eterm_t) );
		eip -> eqpZ_size = 2 * eip -> eqpZ_size;
		eip -> eqpZ_curr = UPDATE_PTR( eip -> eqpZ_curr, eqpZ_old, eip -> eqpZ );

		/* Update pointers from eq-point array */
		for (eqpt = eip -> eqp; eqpt <= eqpk; eqpt++)
			eqpt -> Z = UPDATE_PTR( eqpt -> Z, eqpZ_old, eip -> eqpZ );
		free( eqpZ_old );
	}

	p1 = eqpi -> Z;
	p2 = eqpj -> Z;
	endp1 = p1 + eqpi -> S;
	endp2 = p2 + eqpj -> S;

	/* We know each list contains at least one item! */
	t1 = *p1++;
	t2 = *p2++;

	/* I tried it lots of different ways and discovered that
	   these goto's actually produce the MOST readable form! */
	Zp = eip -> eqpZ_curr;
	for (;;) {
		if (t1 < t2) {
			*Zp++ = t1;
			if (p1 >= endp1) goto finish2;
			t1 = *p1++;
		}
		else {
			/* Disjoint implies t2 < t1 right here. */
			*Zp++ = t2;
			if (p2 >= endp2) goto finish1;
			t2 = *p2++;
		}
	}

finish1:
	for (;;) {
		*Zp++ = t1;
		if (p1 >= endp1) break;
		t1 = *p1++;
	}
	return (Zp);

finish2:
	for (;;) {
		*Zp++ = t2;
		if (p2 >= endp2) break;
		t2 = *p2++;
	}
	return (Zp);
}

/*
 * Test if terminals involved in eq-point j are disjoint from those
 * already flagged.
 */

	static
	bool
disjoint (

struct einfo * eip, /* IN - global EFST info */
struct eqp_t * eqp  /* IN - eq-point */
)
{
	eterm_t *p = eqp -> Z;
	eterm_t *endp = p + eqp -> S;

	while (p < endp) {
		int t = *p++;
		if (eip -> MEMB [t]) return FALSE;
	}
	return TRUE;
}

/*
 * Return terminal with highest index involved in eq-point k
 */

	static
	int
highest_terminal (

struct eqp_t * eqp  /* IN - eq-point */
)
{
	return ( (eqp -> Z) [eqp -> S  -  1]);
}


/*
 * Return the closest terminal to Q involved in the construction of
 * equilateral point number k
 */

	static
	dist_t
closest_terminal (

struct einfo * eip,  /* IN - global EFST info */
struct point * Q,    /* IN - query point */
struct eqp_t * eqpk  /* IN - eq-point */
)
{
	eterm_t *p, *endp;
	int t;
	dist_t min_dist2, dist2;
	int min_t;

	p = eqpk -> Z;
	endp = p + eqpk -> S;

	min_t = *p++;
	min_dist2 = sqr_dist(Q, &(eip -> eqp[min_t].E));
	while (p < endp) {
		t = *p++;
		dist2 = sqr_dist(Q, &(eip -> eqp[t].E));
		if (dist2 < min_dist2) {
			min_dist2 = dist2;
			min_t = t;
		}
	}

	return (sqrt (min_dist2));
}

/*
 * Get bottleneck Steiner distance between equilateral points i and j
 */

	static
	double
getBSD (

struct einfo * eip,  /* IN - global EFST info */
struct eqp_t * eqpi, /* IN - first eq-point */
struct eqp_t * eqpj  /* IN - second eq-point */
)
{
	dist_t rbsd, lbsd;
	if (eqpi -> R EQ NULL) {
		if (eqpj -> R EQ NULL) return bsd(eip -> bsd, eqpi -> E.pnum, eqpj -> E.pnum);
		rbsd = getBSD(eip, eqpi, eqpj -> R);
		lbsd = getBSD(eip, eqpi, eqpj -> L);
	}
	else {
		rbsd = getBSD(eip, eqpi -> R, eqpj);
		lbsd = getBSD(eip, eqpi -> L, eqpj);
	}
	if (rbsd < lbsd) return rbsd;
	return lbsd;
}


/*
 * Computation of eq-point displacement vector
 */

	static
	void
eq_point_disp_vector (

struct einfo * eip,   /* IN - global EFST info */
struct eqp_t * eqpi,  /* IN - first eq-point */
struct eqp_t * eqpj,  /* IN - second eq-point */
struct eqp_t * eqpk   /* IN - new eq-point */
)
{
pnum_t		ti, tj;
double		dx, dy, rdx, rdy;

	/* First compute vector between old eq-points terminal references */
	ti = eqpi -> DV.pnum;
	tj = eqpj -> DV.pnum;

	dx = eip -> eqp [tj].E.x - eip -> eqp [ti].E.x;
	dy = eip -> eqp [tj].E.y - eip -> eqp [ti].E.y;

	/* Add displacements */
	dx -= eqpi -> DV.x;
	dy -= eqpi -> DV.y;
	dx += eqpj -> DV.x;
	dy += eqpj -> DV.y;

	/* Rotate 60 degrees and save */
	rotate60(dx, dy, &rdx, &rdy);
	eqpk -> DV.x	= rdx + eqpi -> DV.x;
	eqpk -> DV.y	= rdy + eqpi -> DV.y;
	eqpk -> DV.pnum = eqpi -> DV.pnum;
}

/*
 * Computation of distance between eq-points.
 * We do this in a numerically careful way by computing displacements
 * and not absolute coordinates.
 */

	static
	dist_t
eq_point_dist (

struct einfo * eip,   /* IN - global EFST info */
struct eqp_t * eqpi,  /* IN - first eq-point */
struct eqp_t * eqpj   /* IN - second eq-point */
)
{
pnum_t		ti, tj;
double		dx, dy, rdx, rdy;

	/* First compute vector between eq-points terminal references */
	ti = eqpi -> DV.pnum;
	tj = eqpj -> DV.pnum;

	dx = eip -> eqp [tj].E.x - eip -> eqp [ti].E.x;
	dy = eip -> eqp [tj].E.y - eip -> eqp [ti].E.y;

	/* Add displacements */
	dx -= eqpi -> DV.x;
	dy -= eqpi -> DV.y;
	dx += eqpj -> DV.x;
	dy += eqpj -> DV.y;

	/* Finally return length of vector */
	return hypot (dx, dy);
}

/*
 * List of terminals involved in the construction of eq-point
 */

	static
	void
eqpoint_terminals (

struct einfo * eip,  /* IN - global EFST info */
struct eqp_t * eqpk, /* IN - eq-point */
struct pset  * pts,  /* OUT - list of point */
bool   append	     /* IN - append to list? */
)
{
	eterm_t *p = eqpk -> Z;
	eterm_t *endp = p + eqpk -> S;
	struct point *pt;

	if (NOT append) pts -> n = 0;

	pt = pts -> a + pts -> n;
	while (p < endp)
		*(pt++) = eip -> eqp[ *p++ ].E;
	pts -> n += eqpk -> S;
}

/*
 * Initialize data structure for eq-point rectangles.
 * We use a hashing like data structure; the rectangle enclosing all
 * the terminals is divided into squares with side length of the
 * longest MST edge. For each such square the list of eq-point
 * rectangles overlapping with that square are stored in an array.

 * The idea of using eq-point rectangles was proposed by Ernst Althaus,
 * Max-Planck-Institut fur Informatik, Germany.
 * He used a 2D-search trees to store the rectangles; we apply
 * a simpler alternative here.
 */

	static
	void
initialize_eqp_rectangles (

struct einfo * eip /* IN/OUT - global EFST info */
)
{
	int i;
	struct point * p;

	/* Allocate square data structure */
	eip -> eqp_squares = NEWA(eip -> pts -> n, struct eqp_square_t *);

	/* Find range of terminal coordinates */
	eip -> minx =  INF_DISTANCE;
	eip -> maxx = -INF_DISTANCE;
	eip -> miny =  INF_DISTANCE;
	eip -> maxy = -INF_DISTANCE;
	for (i = 0; i < eip -> pts -> n; i++) {
		p = &(eip -> pts -> a[i]);
		eip -> minx = MIN( eip -> minx, p -> x - eip -> mean.x);
		eip -> maxx = MAX( eip -> maxx, p -> x - eip -> mean.x);
		eip -> miny = MIN( eip -> miny, p -> y - eip -> mean.y);
		eip -> maxy = MAX( eip -> maxy, p -> y - eip -> mean.y);
		eip -> eqp_squares[i] = NULL;
	}

	/* Set up basic paramaters */
	eip -> eqp_square_size = eip -> max_mst_edge;
	eip -> srangex = floor( (eip -> maxx - eip -> minx) / eip -> eqp_square_size) + 1;
	eip -> srangey = floor( (eip -> maxy - eip -> miny) / eip -> eqp_square_size) + 1;
}


/*
 * Free memory used by eq-point rectangle data structure.
 */

	static
	void
destroy_eqp_rectangles (

struct einfo * eip /* IN/OUT - global EFST info */
)
{
	int i;

	for (i = 0; i < eip -> pts -> n; i++)
	if (eip -> eqp_squares[i] NE NULL) {
		free( eip -> eqp_squares[i][0].eqp ); /* eq-point lists */
		free( eip -> eqp_squares[i]	   ); /* pointers to eq-point lists */
	}

	free( eip -> eqp_squares );
}

/*
 * Compute boundary of rectangle in which the Steiner point joining
 * the current eq-point to another eq-point can be placed
 */

	static
	void
compute_eqp_rectangle (

struct einfo * eip,  /* IN - global EFST info */
struct eqp_t * eqpk, /* IN - eq-point */
dist_t * minx,	     /* OUT - minimum x-coordinate of boundary */
dist_t * maxx,	     /* OUT - maximum x-coordinate of boundary */
dist_t * miny,	     /* OUT - minimum y-coordinate of boundary */
dist_t * maxy	     /* OUT - maximum x-coordinate of boundary */
)
{
	struct point p;
	dist_t alpha;

	if (eqpk -> S EQ 1) {

		/* This is a terminal */
		*minx = eqpk -> E.x - eip -> max_mst_edge;
		*maxx = eqpk -> E.x + eip -> max_mst_edge;
		*miny = eqpk -> E.y - eip -> max_mst_edge;
		*maxy = eqpk -> E.y + eip -> max_mst_edge;
	}
	else {

		/* This is a "normal" eq-point */
		*minx = eqpk -> LP.x; *maxx = eqpk -> LP.x;
		*miny = eqpk -> LP.y; *maxy = eqpk -> LP.y;
		UPDATE_RECTANGLE_BOUNDS(eqpk -> RP);

		/* Projection of eqpk -> LP (from eqpk -> E) to a point
		   at distance max_mst_edge away. */
		alpha = 1.0 + eip -> max_mst_edge / EDIST(&(eqpk -> E), &(eqpk -> LP));
		p.x   = eqpk -> E.x + alpha * (eqpk -> LP.x - eqpk -> E.x);
		p.y   = eqpk -> E.y + alpha * (eqpk -> LP.y - eqpk -> E.y);
		UPDATE_RECTANGLE_BOUNDS(p);

		/* Projection of eqpk -> RP (from eqpk -> E) to a point
		   at distance max_mst_edge away. */
		alpha = 1.0 + eip -> max_mst_edge / EDIST(&(eqpk -> E), &(eqpk -> RP));
		p.x   = eqpk -> E.x + alpha * (eqpk -> RP.x - eqpk -> E.x);
		p.y   = eqpk -> E.y + alpha * (eqpk -> RP.y - eqpk -> E.y);
		UPDATE_RECTANGLE_BOUNDS(p);

		/* Now check all intersections with axis parallel lines
		   through  eqpk -> DC */
		p.x = eqpk -> DC.x + eqpk -> DR + eip -> max_mst_edge;
		p.y = eqpk -> DC.y;
		if (right_turn(&(eqpk -> E), &(eqpk -> LP), &p) AND
		    left_turn(&(eqpk -> E), &(eqpk -> RP), &p))
			UPDATE_RECTANGLE_BOUNDS(p);

		p.x = eqpk -> DC.x - eqpk -> DR - eip -> max_mst_edge;
		p.y = eqpk -> DC.y;
		if (right_turn(&(eqpk -> E), &(eqpk -> LP), &p) AND
		    left_turn(&(eqpk -> E), &(eqpk -> RP), &p))
			UPDATE_RECTANGLE_BOUNDS(p);

		p.x = eqpk -> DC.x;
		p.y = eqpk -> DC.y + eqpk -> DR + eip -> max_mst_edge;
		if (right_turn(&(eqpk -> E), &(eqpk -> LP), &p) AND
		    left_turn(&(eqpk -> E), &(eqpk -> RP), &p))
			UPDATE_RECTANGLE_BOUNDS(p);

		p.x = eqpk -> DC.x;
		p.y = eqpk -> DC.y - eqpk -> DR - eip -> max_mst_edge;
		if (right_turn(&(eqpk -> E), &(eqpk -> LP), &p) AND
		    left_turn(&(eqpk -> E), &(eqpk -> RP), &p))
			UPDATE_RECTANGLE_BOUNDS(p);
	}

	*minx -= (fabs(*minx) * eip -> eps * ((double) eqpk -> S));
	*miny -= (fabs(*miny) * eip -> eps * ((double) eqpk -> S));
	*maxx += (fabs(*maxx) * eip -> eps * ((double) eqpk -> S));
	*maxy += (fabs(*maxy) * eip -> eps * ((double) eqpk -> S));
}

/*
 * Save all eq-point rectangles for a given eq-point size
 */

	static
	void
save_eqp_rectangles (

struct einfo * eip,  /* IN/OUT - global EFST info */
int first_eqp,	     /* IN - index of first eq-point to be saved */
int last_eqp	     /* IN - index of last eq-point to be saved */
)
{
	int i, j, k, si, size;
	struct eqp_t * eqpk;
	dist_t minx, maxx, miny, maxy;
	int total_squares;
	int * curr_count;
	struct eqp_square_t* sqr;
	struct eqp_t ** eqp_alloc;

	if (first_eqp > last_eqp) return; /* no eq-points to save */

	/* Go through all equilateral points.
	 * Find min/max square indices in each dimension and count the total
	 * number of squares that are overlapped by these eq-points.
	 */

	size = eip -> eqp[first_eqp].S;
	total_squares = 0;
	for (k = first_eqp; k <= last_eqp; k++) {
		eqpk = &(eip -> eqp[k]);

		/* Compute geometric range of rectangle for this eq-point */
		compute_eqp_rectangle(eip, eqpk, &minx, &maxx, &miny, &maxy);

		/* Compute square index range for this eq-point */

		eqpk -> SMINX = MAX(0,		      floor( (minx - eip -> minx) /
						       eip -> eqp_square_size));
		eqpk -> SMAXX = MIN(eip -> srangex-1, floor( (maxx - eip -> minx) /
						       eip -> eqp_square_size));
		eqpk -> SMINY = MAX(0,		      floor( (miny - eip -> miny) /
						       eip -> eqp_square_size));
		eqpk -> SMAXY = MIN(eip -> srangey-1, floor( (maxy - eip -> miny) /
						       eip -> eqp_square_size));

		total_squares += (eqpk -> SMAXX - eqpk -> SMINX + 1) *
				 (eqpk -> SMAXY - eqpk -> SMINY + 1);
	}

	/* Allocate space for these rectangles */
	eip -> eqp_squares[size] = NEWA(eip -> srangex * eip -> srangey,
						struct eqp_square_t);
	sqr = eip -> eqp_squares[size];

	/* Initialize eq-point counts */
	curr_count = NEWA(eip -> srangex * eip -> srangey, int);
	for (i = 0; i < eip -> srangex; i++)
		for (j = 0; j < eip -> srangey; j++) {
			si = i * eip -> srangey + j;
			sqr[si].n = 0; curr_count[si] = 0;