efst.c

生成直角Steiner树的程序包

			if (sqr_dist(&(eip -> eqp[r].E), &(eqpj -> E))	>= dist_qpi) continue;
			if (sqr_dist(&(eip -> eqp[r].E), &(eqpk -> RP)) >= dist_qpi) continue;

			rotate2(&(eqpj -> E), &(eqpk -> DC),
				-EDIST(&(eip -> eqp[r].E), &(eqpj -> E))/(2.0*eqpk -> DR), 2.0, &CRP);

			a  = sqrt(aa);
			b  = (a + 2.0 * EDIST(&(eqpi -> E), &(eqpk -> RP))) / SQRT3;
			bb = b*b;
			dist_opqr = EDIST(&(eip -> eqp[r].E), &(eqpk -> DC));

			c = dist_opqr*dist_opqr + eqpk -> DR2;
			d = 2.0*((eqpk -> RP.x - eqpk -> DC.x) * (eqpk -> DC.x - eip -> eqp[r].E.x) +
				 (eqpk -> RP.y - eqpk -> DC.y) * (eqpk -> DC.y - eip -> eqp[r].E.y));
			e = 2.0*((eqpk -> RP.y - eqpk -> DC.y) * (eqpk -> DC.x - eip -> eqp[r].E.x) -
				 (eqpk -> RP.x - eqpk -> DC.x) * (eqpk -> DC.y - eip -> eqp[r].E.y));
			f = (aa+bb)/2.0-c; ff = f*f;
			g = e-a*b;	   gg = g*g;
			h = d-(aa-bb)/2.0; hh = h*h;
			if (solve_quadratic(gg + hh, f*g, ff - hh, &sin1, &sin2)) {
				if ((fabs(sin1) < eip -> eps) OR (fabs(sin2) < eip -> eps)) continue;
				rotate(&(eqpk -> RP), &(eqpk -> DC), sin1, (f+g*sin1)/h, &CRPP);
				if (left_turn(&(eqpk -> E), &(eqpk -> RP), &CRPP) AND
				    right_turn(&(eqpk -> E), &CRP, &CRPP))
					CRP = CRPP;
				rotate(&(eqpk -> RP), &(eqpk -> DC), sin2, (f+g*sin2)/h, &CRPP);
				if (left_turn(&(eqpk -> E), &(eqpk -> RP), &CRPP) AND
				    right_turn(&(eqpk -> E), &CRP, &CRPP))
					CRP = CRPP;
			}

			if (NOT test_and_save_RP(eip, eqpi, eqpj, eqpk, &CRP)) return FALSE;

			aa	 = sqr_dist(&(eqpj -> E), &(eqpk -> RP));
			dist_qpi = 0.999 * aa;
		}
	}

	return TRUE;
}


/*
 * Upper bound test
 */

	static
	bool
upper_bound_test (

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 */
)
{
	dist_t	low_arc_length, rlb, llb, lower_bound, upper_bound, dist, distr, distl,
		rsmt, lsmt, dx, dy, l,
		ub_ri = INF_DISTANCE, ub_lj = INF_DISTANCE, ub_rk = INF_DISTANCE, ub_lk = INF_DISTANCE,
		distrj, distli, cosp, sinp, a, aa, b, bb, sin1, sin2;
	struct point P, CRP, CRPP, CLP, CLPP, M;

	/* Lower bound */

	if (EDIST(&(eqpk -> RP), &(eqpk -> LP)) < eip -> eps) return FALSE;
	rlb = EDIST(&(eqpk -> RP), &(eqpk -> E)) * (1.0 - eip -> eps * ((double) eqpk -> S));
	llb = EDIST(&(eqpk -> LP), &(eqpk -> E)) * (1.0 - eip -> eps * ((double) eqpk -> S));
	lower_bound = MIN(rlb, llb);

	M.x = (eqpk -> RP.x + eqpk -> LP.x) / 2.0;
	M.y = (eqpk -> RP.y + eqpk -> LP.y) / 2.0;
	P.x = eqpk -> DC.x + (eqpk -> DR / EDIST(&(eqpk -> DC), &M))* (M.x - eqpk -> DC.x);
	P.y = eqpk -> DC.y + (eqpk -> DR / EDIST(&(eqpk -> DC), &M))* (M.y - eqpk -> DC.y);
	low_arc_length = 2.0*EDIST(&P, &(eqpk -> RP));

	/* 1. upper bound */

	distr = closest_terminal(eip, &(eqpk -> RP), eqpi);
	distl = closest_terminal(eip, &(eqpk -> LP), eqpj);
	upper_bound = (eqpi -> UB + distr) + (eqpj -> UB + distl) + low_arc_length;
	if (upper_bound	 < lower_bound) return FALSE;
	rsmt = upper_bound;
	lsmt = upper_bound;

	/* 2. upper bound */

	dist = distl;
	if (distr < dist) dist = distr;
	upper_bound = eqpi -> UB + eqpj -> UB + dist +
		      EDIST(&(eqpk -> RP), &(eqpk -> LP)) + eqpk -> BS;
	if (upper_bound	 < lower_bound) return FALSE;
	if (rsmt > upper_bound) rsmt = upper_bound;
	if (lsmt > upper_bound) lsmt = upper_bound;

	/* 3. upper bound */

	if (eqpi -> R EQ NULL)
		ub_ri = EDIST(&(eqpk -> RP), &(eqpi -> E));
	else {
		eip -> termlist -> a[0] = eqpk -> RP;
		eip -> termlist -> n	= 1;
		eqpoint_terminals(eip, eqpi, eip -> termlist, TRUE);
		ub_ri = upper_bound_heuristic(eip -> termlist, eip -> bsd);
	}
	upper_bound = ub_ri + (eqpj -> UB + distl) + low_arc_length;
	if (upper_bound < lower_bound) return FALSE;
	if (rsmt > upper_bound) rsmt = upper_bound;
	if (lsmt > upper_bound) lsmt = upper_bound;

	/* 4. upper bound */

	if (eqpj -> R EQ NULL)
		ub_lj = EDIST(&(eqpk -> LP), &(eqpj -> E));
	else {
		eip -> termlist -> a[0] = eqpk -> LP;
		eip -> termlist -> n	= 1;
		eqpoint_terminals(eip, eqpj, eip -> termlist, TRUE);
		ub_lj = upper_bound_heuristic(eip -> termlist, eip -> bsd);
	}
	upper_bound = eqpi -> UB + distr;
	if (ub_ri < upper_bound) upper_bound = ub_ri;
	upper_bound += ub_lj + low_arc_length;
	if (upper_bound < lower_bound) return FALSE;
	if (rsmt > upper_bound) rsmt = upper_bound;
	if (lsmt > upper_bound) lsmt = upper_bound;

	/* 5. upper bound */

	upper_bound = 0.0;
	eip -> termlist -> a[0] = eqpk -> RP;
	eip -> termlist -> a[1] = eqpk -> LP;
	eip -> termlist -> n	= 2;
	if (eqpi -> R EQ NULL)
		upper_bound += distr;
	else
		eqpoint_terminals(eip, eqpi, eip -> termlist, TRUE);
	if (eqpj -> R EQ NULL)
		upper_bound += distl;
	else
		eqpoint_terminals(eip, eqpj, eip -> termlist, TRUE);
	if (eip -> termlist -> n EQ 2)
		upper_bound += eqpk -> BS + low_arc_length/2.0;
		/* was arc_length before (changed 12dec98) */
	else
		upper_bound += upper_bound_heuristic(eip -> termlist, eip -> bsd) +
			       low_arc_length/2.0;
	if (upper_bound < lower_bound) return FALSE;
	if (rsmt > upper_bound) rsmt = upper_bound;
	if (lsmt > upper_bound) lsmt = upper_bound;

	/* 6. upper bound */

	if (eqpi -> R) {
		eip -> termlist -> a[0] = eqpk -> RP;
		eip -> termlist -> n	= 1;
		eqpoint_terminals(eip, eqpk, eip -> termlist, TRUE);
		ub_rk = upper_bound_heuristic(eip -> termlist, eip -> bsd);
		upper_bound = ub_rk + low_arc_length;
		if (upper_bound < lower_bound) return FALSE;
		if (rsmt > upper_bound) rsmt = upper_bound;
	}

	/* 7. upper bound */

	if (eqpj -> R) {
		eip -> termlist -> a[0] = eqpk -> LP;
		eip -> termlist -> n	= 1;
		eqpoint_terminals(eip, eqpk, eip -> termlist, TRUE);
		ub_lk = upper_bound_heuristic(eip -> termlist, eip -> bsd);
		upper_bound = ub_lk + low_arc_length;
		if (upper_bound < lower_bound) return FALSE;
		if (lsmt > upper_bound) lsmt = upper_bound;
	}

	/* 8. upper bound */

	upper_bound = eqpj -> UB + low_arc_length + eqpk -> BS;
	if (eqpi -> R EQ NULL)
		upper_bound += EDIST(&(eqpk -> RP), &(eqpi -> E));
	else
		upper_bound += ub_ri;
	if (upper_bound < lower_bound) return FALSE;
	if (rsmt > upper_bound) rsmt = upper_bound;

	/* 9. upper bound */

	upper_bound = eqpi -> UB + low_arc_length + eqpk -> BS;
	if (eqpj -> R EQ NULL)
		upper_bound += EDIST(&(eqpk -> LP), &(eqpj -> E));
	else
		upper_bound += ub_lj;
	if (upper_bound < lower_bound) return FALSE;
	if (lsmt > upper_bound) lsmt = upper_bound;

	/* Do the final push */

	distrj = closest_terminal(eip, &(eqpk -> RP), eqpj);
	upper_bound = eqpi -> UB + eqpj -> UB + distr + distrj;
	if (rsmt > upper_bound) rsmt = upper_bound;
	distli = closest_terminal(eip, &(eqpk -> LP), eqpi);
	upper_bound = eqpi -> UB + eqpj -> UB + distl + distli;
	if (lsmt > upper_bound) lsmt = upper_bound;

	if (eqpi -> R NE NULL)
		upper_bound = ub_rk;
	else {
		if (eqpj -> R EQ NULL)
			upper_bound = EDIST(&(eqpj -> E), &(eqpk -> RP)) +
				      EDIST(&(eqpi -> E), &(eqpk -> RP));
		else {
			eip -> termlist -> a[0] = eqpk -> RP;
			eip -> termlist -> n	= 1;
			eqpoint_terminals(eip, eqpj, eip -> termlist, TRUE);
			upper_bound = upper_bound_heuristic(eip -> termlist, eip -> bsd) +
					EDIST(&(eqpi -> E), &(eqpk -> RP));
		}
	}
	if (upper_bound < rsmt) rsmt = upper_bound;

	if (rsmt < 0.999 * rlb) {
		CRP = eqpj -> E;
		memset (&CRPP, 0, sizeof (CRPP));
		get_angle_vector(&(eqpi -> E), &(eqpk -> E), &(eqpk -> RP), &dx, &dy);
		l = sqrt(dx*dx + dy*dy);
		cosp = dx/l; sinp = dy/l;
		a = eqpk -> DR*(SQRT3*cosp+sinp);     aa = a*a;
		b = eqpk -> DR*(SQRT3*sinp-cosp+2.0); bb = b*b;

		if (solve_quadratic(aa + bb, rsmt*b, rsmt*rsmt - aa, &sin1, &sin2)) {
			if ((fabs(sin1) < eip -> eps) OR (fabs(sin2) < eip -> eps)) return TRUE;
			rotate(&(eqpk -> RP), &(eqpk -> DC), sin1, (b*sin1+rsmt)/a, &CRPP);
			if (left_turn(&(eqpk -> E), &(eqpk -> RP), &CRPP) AND
			    right_turn(&(eqpk -> E), &CRP, &CRPP))
				CRP = CRPP;
			rotate(&(eqpk -> RP), &(eqpk -> DC), sin2, (b*sin2+rsmt)/a, &CRPP);
			if (left_turn(&(eqpk -> E), &(eqpk -> RP), &CRPP) AND
			    right_turn(&(eqpk -> E), &CRP, &CRPP))
				CRP = CRPP;

			if (NOT test_and_save_RP(eip, eqpi, eqpj, eqpk, &CRP)) return FALSE;
		}
	}

	if (eqpj -> R)
		upper_bound = ub_lk;
	else {
		if (eqpi -> R EQ NULL)
			upper_bound = EDIST(&(eqpi -> E), &(eqpk -> LP)) +
				      EDIST(&(eqpj -> E), &(eqpk -> LP));
		else {
			eip -> termlist -> a[0] = eqpk -> LP;
			eip -> termlist -> n	= 1;
			eqpoint_terminals(eip, eqpi, eip -> termlist, TRUE);
			upper_bound = upper_bound_heuristic(eip -> termlist, eip -> bsd) +
				      EDIST(&(eqpj -> E), &(eqpk -> LP));
		}
	}
	if (upper_bound < lsmt) lsmt = upper_bound;

	if (lsmt < 0.999 * llb) {
		CLP = eqpi -> E;
		memset (&CLPP, 0, sizeof (CLPP));
		get_angle_vector(&(eqpk -> LP), &(eqpk -> E), &(eqpj -> E), &dx, &dy);
		l = sqrt(dx*dx + dy*dy);
		cosp = dx/l; sinp = dy/l;
		a = eqpk -> DR*(SQRT3*cosp+sinp);     aa=a*a;
		b = eqpk -> DR*(SQRT3*sinp-cosp+2.0); bb=b*b;

		if (solve_quadratic(aa + bb, lsmt*b, lsmt*lsmt - aa, &sin1, &sin2)) {
			if ((fabs(sin1) < eip -> eps) OR (fabs(sin2) < eip -> eps)) return TRUE;
			rotate(&(eqpk -> LP), &(eqpk -> DC), -sin1, (b*sin1+lsmt)/a, &CLPP);
			if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
			    left_turn(&(eqpk -> E), &CLP, &CLPP))
				CLP = CLPP;
			rotate(&(eqpk -> LP), &(eqpk -> DC), -sin2, (b*sin2+lsmt)/a, &CLPP);
			if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
			    left_turn(&(eqpk -> E), &CLP, &CLPP))
				CLP = CLPP;

			if (NOT test_and_save_LP(eip, eqpi, eqpj, eqpk, &CLP)) return FALSE;
		}
	}

	return TRUE;
}


/*
 * Wedge property test
 */

	static
	bool
wedge_test (

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 */
)
{
	int r, t;
	int right_counter = 0;
	int middle_counter = 0;
	int left_counter = 0;
	int top = highest_terminal(eqpk);
	bool flag;
	dist_t dist, dist1, dist2, bsdi, bsdj;
	struct point SP;
	struct eqp_t * eqpt;
	struct eqp_t * other_eqp;

	other_eqp = eqpj;
	if (NOT disjoint(eip, other_eqp)) other_eqp = eqpi;
	set_member_arr(eip, other_eqp, TRUE);

#ifdef HAVE_GMP
	if (Multiple_Precision > 0) {
		update_eqpoint_and_displacement (eip, eqpk);
	}
#endif

	for (t = 0; t < eip -> pts -> n; t++) {

		/* Terminal should not be part of this eq-point */
		if (eip -> MEMB[t]) continue;

		eqpt = &(eip -> eqp[t]);

		/* Is terminal completely outside feasible area? */
		if ((left_turn(&(eqpi -> E), &(eqpk -> LP), &(eqpt -> E))) OR
		    (right_turn(&(eqpj -> E), &(eqpk -> RP), &(eqpt ->E)))) continue;

		if (left_turn(&(eqpk -> E), &(eqpk -> LP), &(eqpt -> E))) {
			left_counter++; continue;
		}

		if (right_turn(&(eqpk -> E), &(eqpk -> RP), &(eqpt -> E))) {
			right_counter++; continue;
		}

		if (sqr_dist(&(eqpk -> DC), &(eqpt -> E)) > eqpk -> DR2) {
			middle_counter++;
			if (t <= top) continue; /* avoid duplicates */

			project_point(&(eqpk -> E), &(eqpk -> DC),
				      &(eqpt -> E), &SP);
			dist = EDIST(&(eqpt -> E), &SP) *
			       (1.0 - eip -> eps * ((double) eqpk -> S));

			/* Is the last edge too long? */
			if (dist >= getBSD(eip, eqpt, eqpk)) continue;

			flag = TRUE;
			for (r = 0; r < eip -> pts -> n; r++) {
				if (r EQ t) continue;
				if ((EDIST(&SP, &(eip -> eqp[r].E))	     < dist) AND
				    (EDIST(&(eqpt -> E), &(eip -> eqp[r].E)) < dist)) {
					flag = FALSE; break;
				}
			}
			if (NOT flag) continue;

			eip -> termlist -> a[0] = eqpt -> E;
			eip -> termlist -> n	= 1;
			eqpoint_terminals(eip, eqpk, eip -> termlist, TRUE);

			if (eqpk -> S > 2) {
				dist1 = upper_bound_heuristic(eip -> termlist, eip -> bsd);
				dist2 = eq_point_dist(eip, eqpt, eqpk) *
					(1.0 - eip -> eps * ((double) eqpk -> S));
				if (dist1 < dist2) {
					flag = FALSE;
				}
				else {
					bsdi = getBSD(eip, eqpi, eqpt);
					bsdj = getBSD(eip, eqpj, eqpt);
					if ((bsdi + eqpk -> UB < dist2) OR
					    (bsdj + eqpk -> UB < dist2) OR
					    (bsdi + bsdj + eqpi -> UB + eqpj -> UB < dist2))
						flag = FALSE;
				}
			}
			if (NOT flag) continue;

			/* This FST survived all tests. Save it! */
			test_and_save_fst(eip, eqpt, eqpk);
		}
	}

	set_member_arr(eip, other_eqp, FALSE);
	flag = FALSE;
	if (middle_counter >= 1) {
		flag = TRUE;
		if ((middle_counter EQ 1) AND (left_counter + right_counter EQ 0)) return FALSE;
	}
	else {
		if ((left_counter >= 1) AND (right_counter >= 1)) flag = TRUE;
	}
	return(flag);
}

/*
 * Compute the EFSTs for the given set of terminals, which are now
 * guaranteed to be unique.
 */

	static
	void
compute_efsts_for_unique_terminals (

struct einfo *		eip	/* IN - global EFST info */
)
{
int			n;
int			nedges;
struct pset *		pts;
struct edge *		ep;
struct edge *		mst_edges;
dist_t			mst_len;
char			buf1 [32];
eterm_t			*new_Zp;
int			i, j, k, si, l, size, iter, sz, starti, endi;
dist_t			upper_bound;
struct eqp_t		*eqpi, *eqpj, *eqpk, *eqpt, *eqp_old;
struct eqp_t		**eqp_list, **eqpp, **eqppp;
struct elist		*rp;

	pts = eip -> pts;
	n = pts -> n;

	/* Compute minimum spanning tree */

	mst_edges = NEWA (n - 1, struct edge);
	nedges = euclidean_mst (pts, mst_edges);
	if (nedges NE n - 1) {
		fatal ("compute_efsts_for_unique_terminals: Bug 1.");
	}

	mst_len = 0.0;
	eip -> max_mst_edge = 0.0;
	ep = mst_edges;
	for (i = 0; i < nedges; ep++, i++) {
		mst_len += ep -> len;
		/* Get longest MST edge */
		if (eip -> max_mst_edge < ep -> len) eip -> max_mst_edge = ep -> len;
	}
	eip -> mst_length = mst_len;

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

	/* Set relative epsilon for floating point comparisons.
	   We use the constant EpsilonFactor to indicate
	   the maximum relative error that is expected.
	*/

	eip -> eps   = ((double) EpsilonFactor) * DBL_EPSILON;

	/* Compute bottleneck Steiner distances */

	eip -> 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);
	}

	/* Set up hashing data structure */

	eip -> term_check	= NEWA (n, bool);
	eip -> hash		= NEWA (n, struct elist *);

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

	/* Compute the mean of all terminals.  We translate the terminals */
	/* so that the mean is at the origin.  The coordinates of the */
	/* translated instance have smaller magnitude, which permits us */
	/* to compute eq-point coordinates with higher precision. */

	eip -> mean.x = 0.0;
	eip -> mean.y = 0.0;
	for (k = 0; k < n; k++) {
		eip -> mean.x += pts -> a[k].x;
		eip -> mean.y += pts -> a[k].y;
	}
	eip -> mean.x = floor( eip -> mean.x / (double) n)