efst.c

生成直角Steiner树的程序包

		}

	/* Find counts for each square */
	for (k = first_eqp; k <= last_eqp; k++) {
		eqpk = &(eip -> eqp[k]);

		for (i = eqpk -> SMINX; i <= eqpk -> SMAXX; i++)
			for (j = eqpk -> SMINY; j <= eqpk -> SMAXY; j++)
				sqr[ i * eip -> srangey + j ].n++;
	}

	/* Set up pointers */
	eqp_alloc = NEWA(total_squares, struct eqp_t *);
	for (i = 0; i < eip -> srangex; i++)
		for (j = 0; j < eip -> srangey; j++) {
			si = i * eip -> srangey + j;
			sqr[si].eqp = eqp_alloc; eqp_alloc += sqr[si].n;
		}

	/* Fill arrays */
	for (k = first_eqp; k <= last_eqp; k++) {
		eqpk = &(eip -> eqp[k]);

		for (i = eqpk -> SMINX; i <= eqpk -> SMAXX; i++)
			for (j = eqpk -> SMINY; j <= eqpk -> SMAXY; j++) {
				si = i * eip -> srangey + j;
				sqr[si].eqp[ curr_count[si]++ ] = eqpk;
			}
	}

	free(curr_count);
}


/*
 * Generate compatible eq-points, i.e., eq-points that are
 * close enough to be joined to a given eq-point.
 */

	static
	void
generate_compatible_eqp (

struct einfo * eip,	  /* IN - global EFST info */
int size,		  /* IN - prespecified size of eq-points */
struct eqp_t * eqpi,	  /* IN - given eq-point */
struct eqp_t ** eqp_list  /* OUT - list of compatible eq-points */
)
{
	struct eqp_square_t * sqr;
	int i, j, l, si;
	struct eqp_t ** eqpp;
	struct eqp_t *	eqpj;

	sqr  = eip -> eqp_squares[size];
	eqpp = eqp_list;

	if (sqr) {
		for (i = eqpi -> SMINX; i <= eqpi -> SMAXX; i++)
			for (j = eqpi -> SMINY; j <= eqpi -> SMAXY; j++) {
				si = i * eip -> srangey + j;
				for (l = 0; l < sqr[si].n; l++) {
					eqpj = sqr[si].eqp[l];
					if (NOT (eqpj -> CHOSEN)) {
						*(eqpp++) = eqpj;
						eqpj -> CHOSEN = TRUE;
					}
				}
			}
	}
	*eqpp = NULL;
}



/*
 * Basic projection test (case I)
 */

	static
	bool
projection_test_case_I (

struct einfo * eip,   /* IN - global EFST info */
struct eqp_t * eqpi,  /* IN - first eq-point */
struct eqp_t * eqpj   /* IN - second eq-point */
)
{
	if (eqpi -> R) {
		get_angle_vector(&(eqpi -> RP), &(eqpi -> E), &(eqpj -> E),
				 &(eip -> dxi), &(eip -> dyi));
		if (angle_geq_120(eip -> dxi, eip -> dyi))
			return FALSE;
	}

	if (eqpj -> R) {
		get_angle_vector(&(eqpi -> E), &(eqpj -> E), &(eqpj -> LP),
				 &(eip -> dxj), &(eip -> dyj));
		if (angle_geq_120(eip -> dxj, eip -> dyj))
			return FALSE;
	}

	return TRUE;
}

/*
 * Projection test (cases II - VI)
 */

	static
	bool
projection_test_cases_II_VI (

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 */
)
{
	struct point CLP, CRP;

	if (eqpi -> R) {
		if (angle_geq_60(eip -> dxi, eip -> dyi)) {
			if (sqr_dist(&(eqpk -> DC), &(eqpi -> LP)) >= eqpk -> DR2)
				return FALSE;				     /* case II	 */
			project_point(&(eqpi -> E),  &(eqpk -> DC),
				      &(eqpi -> LP), &(eqpk -> LP));	     /* case III */
			project_point_perp(&(eqpi -> E),  &(eqpk -> DC),
					   &(eqpi -> DC), &(eqpk -> RP));
		}
		else {
			if (sqr_dist(&(eqpi -> DC), &(eqpk -> LP)) <= eqpi -> DR2)
				return FALSE;				     /* case IV	 */
			if (sqr_dist(&(eqpk -> DC), &(eqpi -> RP)) >= eqpk -> DR2) {
				if (sqr_dist(&(eqpk -> DC), &(eqpi -> LP)) >= eqpk -> DR2)
					return FALSE;			     /* case V	 */
				project_point_perp(&(eqpi -> E),  &(eqpk -> DC),
						   &(eqpi -> DC), &(eqpk -> RP));
			}
			else
				project_point(&(eqpi -> E),  &(eqpk -> DC),
					      &(eqpi -> RP), &(eqpk -> RP)); /* case VI	 */

			if (right_turn(&(eqpi -> E), &(eqpk -> LP), &(eqpi -> LP)))
				project_point(&(eqpi -> E),  &(eqpk -> DC),
					      &(eqpi -> LP), &(eqpk -> LP));
		}
	}


	if (eqpj -> R) {
		memset (&CLP, 0, sizeof (CLP));
		memset (&CRP, 0, sizeof (CRP));
		if (angle_geq_60(eip -> dxj, eip -> dyj)) {
			if (sqr_dist(&(eqpk -> DC), &(eqpj -> RP)) >= eqpk -> DR2)
				return FALSE;				     /* case II	 */
			project_point(&(eqpj -> E), &(eqpk -> DC),
				      &(eqpj -> RP), &CRP);		     /* case III */
			if (right_turn(&(eqpk -> E), &CRP, &(eqpk -> LP)))
				return FALSE;
			if (right_turn(&(eqpk -> E), &CRP, &(eqpk -> RP)))
				eqpk -> RP = CRP;
			project_point_perp(&(eqpj -> E),  &(eqpk -> DC),
					   &(eqpj -> DC), &CLP);
			if (right_turn(&(eqpk -> E), &(eqpk -> RP), &CLP))
				return FALSE;
			if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLP))
				eqpk -> LP = CLP;
		}
		else {
			if (sqr_dist(&(eqpj -> DC), &(eqpk -> RP)) <= eqpj -> DR2)
				return FALSE;				     /* case IV */
			if (sqr_dist(&(eqpk -> DC), &(eqpj -> LP)) >= eqpk -> DR2)  {
				if (sqr_dist(&(eqpk -> DC), &(eqpj -> RP)) >= eqpk -> DR2)
					return FALSE;			     /* case V	*/
				project_point_perp(&(eqpj -> E),  &(eqpk -> DC),
						   &(eqpj -> DC), &CLP);
			}
			else
				project_point(&(eqpj -> E),  &(eqpk -> DC),
					      &(eqpj -> LP), &CLP);	     /* case VI */

			if (right_turn(&(eqpk -> E), &(eqpk -> RP), &CLP))
				return FALSE;
			if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLP))
				eqpk -> LP = CLP;
			if (right_turn(&(eqpj -> E), &(eqpj -> RP), &(eqpk -> RP))) {
				project_point(&(eqpj ->E),   &(eqpk -> DC),
					      &(eqpj -> RP), &CRP);
				if (right_turn(&(eqpk -> E), &CRP, &(eqpk -> LP)))
					return FALSE;
				if (right_turn(&(eqpk -> E), &CRP, &(eqpk -> RP)))
					eqpk -> RP = CRP;
			}
		}
	}

	return TRUE;
}


/*
 * Bottleneck property test
 */

	static
	bool
bsd_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 */
)
{
	double bsbs, a, b, aa, sin1, sin2;
	struct point CLP, CRP, CLP1, CRP1, CLP2, CRP2, ai;

	eqpk -> BS = getBSD(eip, eqpi, eqpj);
	eqpk -> UB = eqpi -> UB + eqpj -> UB + eqpk -> BS;

	memset (&CLP, 0, sizeof (CLP));

	if (eqpi -> R EQ NULL) {
		rotate2(&(eqpi -> E), &(eqpk -> DC),
			eqpk -> BS/(2.0*eqpk -> DR), 2.0, &CLP);

		if (NOT test_and_save_LP(eip, eqpi, eqpj, eqpk, &CLP)) return FALSE;
	}
	else {
		bsbs = eqpk -> BS * eqpk -> BS;
		project_point(&(eqpi -> E),  &(eqpi -> DC),
			      &(eqpk -> LP), &ai);
		aa   = sqr_dist(&ai, &(eqpk -> LP));
		if (bsbs < 0.999 * aa) {
			a    = sqrt(aa);
			b    = (a + 2.0 * (  EDIST(&(eqpj -> E), &(eqpk -> LP))
				+ EDIST(&(eqpi -> R -> E), &ai))) / SQRT3;

			if (solve_quadratic(aa + b*b, eqpk -> BS * b, bsbs - aa, &sin1, &sin2)) {
				memset (&CLP1, 0, sizeof (CLP1));
				memset (&CLP2, 0, sizeof (CLP2));
				rotate2(&(eqpk -> LP), &(eqpk -> DC),
					-sin1, (b * sin1 + eqpk -> BS)/a, &CLP1);
				rotate2(&(eqpk -> LP), &(eqpk -> DC),
					-sin2, (b * sin2 + eqpk -> BS)/a, &CLP2);
				if (right_turn(&(eqpk -> E), &CLP2, &CLP1))
					CLP1 = CLP2;

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

	memset (&CRP, 0, sizeof (CRP));

	if (eqpj -> R EQ NULL) {
		rotate2(&(eqpj -> E), &(eqpk -> DC),
			-eqpk -> BS/(2.0*eqpk -> DR), 2.0, &CRP);

		if (NOT test_and_save_RP(eip, eqpi, eqpj, eqpk, &CRP)) return FALSE;
	}
	else {
		bsbs = eqpk -> BS * eqpk -> BS;
		project_point(&(eqpj -> E),  &(eqpj -> DC),
			      &(eqpk -> RP), &ai);
		aa   = sqr_dist(&ai, &(eqpk -> RP));
		if (bsbs < 0.999 * aa) {
			a    = sqrt(aa);
			b    = (a + 2.0 * (  EDIST(&(eqpi -> E), &(eqpk -> RP))
				+ EDIST(&(eqpj -> L -> E), &ai))) / SQRT3;

			if (solve_quadratic(aa + b*b, eqpk -> BS * b, bsbs - aa, &sin1, &sin2)) {
				memset (&CRP1, 0, sizeof (CRP1));
				memset (&CRP2, 0, sizeof (CRP2));
				rotate2(&(eqpk -> RP), &(eqpk -> DC),
					sin1, (b * sin1 + eqpk -> BS)/a, &CRP1);
				rotate2(&(eqpk -> RP), &(eqpk -> DC),
					sin2, (b * sin2 + eqpk -> BS)/a, &CRP2);
				if (right_turn(&(eqpk -> E), &CRP1, &CRP2))
					CRP1 = CRP2;

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

	return TRUE;
}


/*
 * Lune property test
 */

	static
	bool
lune_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;
	struct point CJ, AI, CLP, CLPP, CRP, CRPP;
	dist_t a, aa, b, bb, c, d, e, f, ff, g, gg, h, hh;
	dist_t dist_qicj, dist_oacr, dist_opqr, dist_pqi, dist_piai, dist_qpi;
	dist_t sin1, sin2;

	memset (&CLPP, 0, sizeof (CLPP));
	memset (&CRPP, 0, sizeof (CRPP));

	if (eqpi -> R) {
		project_point(&(eqpi -> E),  &(eqpi -> DC),
			      &(eqpk -> LP), &CJ);

		aa	  = sqr_dist(&CJ, &(eqpk -> LP));
		dist_qicj = 0.999 * aa; /* only look at terminals which really are inside lune */

		for (r = 0; r < eip -> pts -> n; r++) {
			if (sqr_dist(&(eip -> eqp[r].E), &CJ)		>= dist_qicj) continue;
			if (sqr_dist(&(eip -> eqp[r].E), &(eqpk -> LP)) >= dist_qicj) continue;

			CLP = eqpi -> E;
			a   = sqrt(aa);
			b   = (a + 2.0 * (  EDIST(&(eqpj -> E), &(eqpk -> LP))
				 + EDIST(&(eqpi -> R -> E), &CJ))) / SQRT3;
			bb  = b*b;
			dist_oacr = EDIST(&(eip -> eqp[r].E), &(eqpi -> DC));
			c = dist_oacr*dist_oacr + eqpi -> DR2;
			d = 2.0*((CJ.x - eqpi -> DC.x) * (eqpi -> DC.x - eip -> eqp[r].E.x) +
				 (CJ.y - eqpi -> DC.y) * (eqpi -> DC.y - eip -> eqp[r].E.y));
			e = 2.0*((CJ.y - eqpi -> DC.y) * (eqpi -> DC.x - eip -> eqp[r].E.x) -
				 (CJ.x - eqpi -> DC.x) * (eqpi -> 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 -> LP), &(eqpk -> DC), -sin1, (f+g*sin1)/h, &CLPP);
				if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
				    left_turn (&(eqpk -> E), &CLP, &CLPP))
					CLP = CLPP;
				rotate(&(eqpk -> LP), &(eqpk -> DC), -sin2, (f+g*sin2)/h, &CLPP);
				if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
				    left_turn (&(eqpk -> E), &CLP, &CLPP))
					CLP = CLPP;
			}

			dist_opqr = EDIST(&(eip -> eqp[r].E), &(eqpk -> DC));
			c = dist_opqr*dist_opqr + eqpk -> DR2;
			d = 2.0*((eqpk -> LP.x - eqpk -> DC.x) * (eqpk -> DC.x - eip -> eqp[r].E.x) +
				 (eqpk -> LP.y - eqpk -> DC.y) * (eqpk -> DC.y - eip -> eqp[r].E.y));
			e = 2.0*((eqpk -> LP.y - eqpk -> DC.y) * (eqpk -> DC.x - eip -> eqp[r].E.x) -
				 (eqpk -> LP.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 -> LP), &(eqpk -> DC), -sin1, (f+g*sin1)/h, &CLPP);
				if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
				    left_turn (&(eqpk -> E), &CLP, &CLPP))
					CLP = CLPP;
				rotate(&(eqpk -> LP), &(eqpk -> DC), -sin2, (f+g*sin2)/h, &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;

			project_point(&(eqpi -> E),  &(eqpi -> DC),
				      &(eqpk -> LP), &CJ);
			aa	  = sqr_dist(&CJ, &(eqpk -> LP));
			dist_qicj = 0.999 * aa;
		}
	}
	else {
		aa	 = sqr_dist(&(eqpi -> E), &(eqpk -> LP));
		dist_pqi = 0.999 * aa;

		for (r = 0; r < eip -> pts -> n; r++) {
			if (r EQ eqpi -> E.pnum) continue;
			if (sqr_dist(&(eip -> eqp[r].E), &(eqpi -> E))	>= dist_pqi) continue;
			if (sqr_dist(&(eip -> eqp[r].E), &(eqpk -> LP)) >= dist_pqi) continue;

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

			a  = sqrt(aa);
			b  = (a + 2.0 * EDIST(&(eqpj -> E), &(eqpk -> LP))) / SQRT3;
			bb = b*b;
			dist_opqr = EDIST(&(eip -> eqp[r].E), &(eqpk -> DC));
			c = dist_opqr*dist_opqr + eqpk -> DR2;
			d = 2.0*((eqpk -> LP.x - eqpk -> DC.x) * (eqpk -> DC.x - eip -> eqp[r].E.x) +
				 (eqpk -> LP.y - eqpk -> DC.y) * (eqpk -> DC.y - eip -> eqp[r].E.y));
			e = 2.0*((eqpk -> LP.y - eqpk -> DC.y) * (eqpk -> DC.x - eip -> eqp[r].E.x) -
				 (eqpk -> LP.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 -> LP), &(eqpk -> DC), -sin1, (f+g*sin1)/h, &CLPP);
				if (right_turn(&(eqpk -> E), &(eqpk -> LP), &CLPP) AND
				    left_turn (&(eqpk -> E), &CLP, &CLPP))
					CLP = CLPP;
				rotate(&(eqpk -> LP), &(eqpk -> DC), -sin2, (f+g*sin2)/h, &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;

			aa	 = sqr_dist(&(eqpi -> E), &(eqpk -> LP));
			dist_pqi = 0.999 * aa;
		}
	}


	if (eqpj -> R) {
		project_point(&(eqpj -> E),  &(eqpj -> DC),
			      &(eqpk -> RP), &AI);

		aa	  = sqr_dist(&AI, &(eqpk -> RP));
		dist_piai = 0.999 * aa;

		for (r = 0; r < eip -> pts -> n; r++) {
			if (sqr_dist(&(eip -> eqp[r].E), &AI)		>= dist_piai) continue;
			if (sqr_dist(&(eip -> eqp[r].E), &(eqpk -> RP)) >= dist_piai) continue;

			CRP = eqpj -> E;
			a   = sqrt(aa);
			b   = (a + 2.0 * (  EDIST(&(eqpi -> E), &(eqpk -> RP))
				 + EDIST(&(eqpj -> L -> E), &AI))) / SQRT3;
			bb  = b*b;
			dist_oacr = EDIST(&(eip -> eqp[r].E), &(eqpj -> DC));
			c = dist_oacr*dist_oacr + eqpj -> DR2;
			d = 2.0*((AI.x - eqpj -> DC.x) * (eqpj -> DC.x - eip -> eqp[r].E.x) +
				 (AI.y - eqpj -> DC.y) * (eqpj -> DC.y - eip -> eqp[r].E.y));
			e = 2.0*((AI.y - eqpj -> DC.y) * (eqpj -> DC.x - eip -> eqp[r].E.x) -
				 (AI.x - eqpj -> DC.x) * (eqpj -> 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;
			}

			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;

			project_point(&(eqpj -> E),  &(eqpj -> DC),
				      &(eqpk -> RP), &AI);
			aa	  = sqr_dist(&AI, &(eqpk -> RP));
			dist_piai = 0.999 * aa;
		}
	}
	else {
		aa	 = sqr_dist(&(eqpj -> E), &(eqpk -> RP));
		dist_qpi = 0.999 * aa;

		for (r = 0; r < eip -> pts -> n; r++) {
			if (r EQ eqpj -> E.pnum) continue;