greedy.c

生成直角Steiner树的程序包

	}

	in.numberofpointattributes	= 0;
	in.pointattributelist		= NULL;
	in.pointmarkerlist		= NULL;
	in.numberofsegments		= 0;
	in.numberofholes		= 0;
	in.numberofregions		= 0;
	in.regionlist			= NULL;
	out.pointlist			= NULL;
	out.trianglelist		= NULL;
	out.neighborlist		= NULL;

	triangulate ("znNBQ", &in, &out, &vorout);

	if (out.numberoftriangles EQ 0) {

		/* There are no triangles (all input points co-linear).
		   Compute SMT as straight line segment between points */
		minp.x = INF_DISTANCE; maxp.x = -INF_DISTANCE;
		minp.y = INF_DISTANCE; maxp.y = -INF_DISTANCE;
		for (i = 0; i < pts -> n; i++) {
			minp.x = MIN(minp.x, pts -> a [i].x);
			maxp.x = MAX(maxp.x, pts -> a [i].x);
			minp.y = MIN(minp.y, pts -> a [i].y);
			maxp.y = MAX(maxp.y, pts -> a [i].y);
		}

		free (in.pointlist);
		free (out.pointlist);
		free (out.trianglelist);
		free (out.neighborlist);
		return (EDIST(&minp, &maxp));
	}

	/* Construct edge information */
	greedy_edges	= NEWA (out.numberofedges, struct greedy_edge);
	triangleedges	= NEWA (3*out.numberoftriangles, int);
	for (i = 0; i < 3*out.numberoftriangles; i++) {
		triangleedges [i] = -1;
	}

	ei = 0;
	for (i = 0; i < out.numberoftriangles; i++) {
		for (j = 0; j < 3; j++) {
			p1 = out.trianglelist [3*i + j];
			p2 = out.trianglelist [3*i + ((j+1) % 3)];
			if (triangleedges [3*i + j] EQ -1) {
				/* only add once */
				greedy_edges [ei].len	 = EDIST (&(pts -> a [p1]),
								  &(pts -> a [p2]));
				/* Get BSD info if it is there */
				if ((bsdp NE NULL) AND
				    (pts -> a [p1].pnum >= 0) AND
				    (pts -> a [p2].pnum >= 0)) {
					greedy_edges [ei].bsd = bsd (bsdp,
								     pts -> a [p1].pnum,
								     pts -> a [p2].pnum);
				}
				else {
					greedy_edges [ei].bsd = greedy_edges [ei].len;
				}

				greedy_edges [ei].p1	 = p1;
				greedy_edges [ei].p2	 = p2;
				greedy_edges [ei].mst	 = FALSE;
				triangleedges [3*i + j]	 = ei;

				/* Now go through neighbouring triangles */
				/* and add information about the new edge */
				for (ni = 0; ni < 3; ni++) {
					nt = out.neighborlist [3*i + ni];
					if (nt EQ -1) continue;
					for (nj = 0; nj < 3; nj++) {
						np1 = out.trianglelist [3*nt + nj];
						np2 = out.trianglelist [3*nt + ((nj+1) % 3)];
						if ((np1 EQ p2) AND
						    (np2 EQ p1)) {
							 /* found reverse edge */
							triangleedges [3*nt + nj] = ei;
						}
					}
				}
				++ei;
			}
		}
	}

	/* Sort edges */

	sortededges = NEWA (out.numberofedges, int);
	for (i = 0; i < out.numberofedges; i++) {
		sortededges [i] = i;
	}
	qsort (sortededges, out.numberofedges, sizeof(int), comp_edges);

	/* Use Kruskal to find MST */

	dsuf_create (&mstsets, pts -> n);
	for (i = 0; i < pts -> n; i++) {
		dsuf_makeset (&mstsets, i);
	}
	total_mst_l = 0.0;
	for (i = 0; i < out.numberofedges; i++) {
		/* go through edges in sorted order */
		ei = sortededges [i];
		root [1] = dsuf_find (&mstsets, greedy_edges [ei].p1);
		root [2] = dsuf_find (&mstsets, greedy_edges [ei].p2);

		if (root [1] NE root [2]) {
			dsuf_unite (&mstsets, root [1], root [2]);
			total_mst_l += greedy_edges [ei].len;
			greedy_edges [ei].mst = TRUE; /* remember this edge */
		}
	}

	/* Generate FSTs with 3 and 4 terminals */

	greedy_fsts = NEWA (4*out.numberoftriangles, struct greedy_fst);
	fst_count = 0;
	for (i = 0; i < out.numberoftriangles; i++) {

		/* Count the number of MST edges and find their length sum */
		mst_count1 = 0;
		mst_l1	   = 0.0;
		for (j = 0; j < 3; j++) {
			ei  = triangleedges [3*i + j];
			l[j] = greedy_edges [ei].len;
			if (greedy_edges [ei].mst) {
				mst_count1++;
				mst_l1 += l[j];
			}
		}

		for (j = 0; j < 3; j++) {
			terms [j] = out.trianglelist [3*i + j];
		}
		terms [3] = -1;

		/* Do the three terminals induced a connected dubgraph of the MST? */
		if (mst_count1 EQ 2) {
			/* Yes, so we know the length of the corresponding MST */
			mst_l = mst_l1;
		}
		else {
			/* No, they do not. Compute MST for these three terminals */
			mst_l = MIN(l[0]+l[1], MIN(l[0]+l[2], l[1]+l[2]));
		}

		/* Add 3-terminal FST */
		add_fst (pts, terms, 3, greedy_fsts, &fst_count, mst_l, (mst_count1 EQ 2));

		/* Now go through neighbouring triangles and add */
		/* valid FSTs */
		for (ni = 0; ni < 3; ni++) {

			nt = out.neighborlist [3*i + ni]; /* get triangle index */
			if (nt <= i) continue;		  /* only consider once... */

			/* Find neighbouring edge and outgoing MST */
			/* edge (note that there can be at most one) */
			neighbour_edge	   = -1;
			neighbour_edge_idx = -1;
			mst_count2	   =  0;

			for (nj = 0; nj < 3; nj++) {
				ei = triangleedges [3*nt + nj];
				if (greedy_edges [ei].mst) {
					mst_count2++;
				}

				/* Is this the neighouring edge? */
				for (j = 0; j < 3; j++) {
					if (triangleedges [3*i + j] EQ ei) {
						neighbour_edge	   = ei;
						neighbour_edge_idx = nj;
					}
				}
			}


			/* Idenfify fourth terminal */

			p4 = out.trianglelist [3*nt + ((neighbour_edge_idx + 2) % 3)];

			/* Make list of points on border of triangles */

			j = -1;
			jj = -1;
			do {
				terms [++jj] = out.trianglelist [3*i + (++j)];
			} while (triangleedges [3*i + j] NE neighbour_edge);
			terms [++jj] = p4;
			while (j < 2) {
				terms [++jj] = out.trianglelist [3*i + (++j)];
			}

			/* Make sure that the two triangles make up a */
			/* convex region */

			convex_region = TRUE;
			for (j = 0; j < 4; j++) {
				if (right_turn (&(pts -> a [terms [j]]),
						&(pts -> a [terms [(j+1) % 4]]),
						&(pts -> a [terms [(j+2) % 4]]))) {
					convex_region = FALSE;
					break;
				}
			}
			if (NOT convex_region) continue; /* do not try to construct FST */

			/* Compute MST for terminals */
			mst_l = mst_length(pts, terms, 4);

			/* Add 4-terminal FST */
			add_fst (pts,
				 terms,
				 4,
				 greedy_fsts,
				 &fst_count,
				 mst_l,
				 (mst_count1 + mst_count2 >= 3));
		}
	}

	/* Sort FSTs */

	sortedfsts = NEWA (fst_count, int);
	for (i = 0; i < fst_count; i++) {
		sortedfsts [i] = i;
	}
	qsort (sortedfsts, fst_count, sizeof (int), comp_fsts);

	/* Use greedy concatenation algorithm:
	   1. Use FSTs that span connected subgraphs of MST
	   2. Insert not yet tried FSTs into existing tree
	*/

	new_chosen = NEWA (fst_count, bool);
	best_smt_l = INF_DISTANCE;
	for (k = -1; k < fst_count; k++) { /* k=-1 means use MST FSTs */

		/* Check that FST is not already in current tree */

		if ((k >= 0) AND greedy_fsts [ sortedfsts[k] ].chosen) continue;

		/* Build union-find structure */
		dsuf_create (&fstsets, pts -> n);
		for (i = 0; i < pts -> n; i++) {
			dsuf_makeset (&fstsets, i);
		}

		for (i = 0; i < fst_count; i++)
			new_chosen[i] = greedy_fsts[i].chosen;

		if (k >= 0) {
			/* Insert FST number k */
			fk = sortedfsts [k];
			term_count = greedy_fsts [fk].count;

			for (j = 1; j < term_count; j++) {
				dsuf_unite (&fstsets,
					    dsuf_find (&fstsets, greedy_fsts [fk].p [0]),
					    dsuf_find (&fstsets, greedy_fsts [fk].p [j]));
			}
			total_smt_l = greedy_fsts [fk].len;
			new_chosen[fk] = TRUE;
		}
		else
			total_smt_l = 0.0;

		/* Go through chosen FSTs (in sorted order) */
		for (i = 0; i < fst_count; i++) {
			fi = sortedfsts [i];
			if (NOT greedy_fsts [fi].chosen) continue;
			term_count = greedy_fsts [fi].count;

			/* check that no two terminals are in the same block */
			in_same_block = FALSE;
			for (j = 0; j < term_count; j++) {
				root [j] = dsuf_find (&fstsets, greedy_fsts [fi].p [j]);
			}
			for (j = 0; j < term_count-1; j++) {
				for (jj = j+1; jj < term_count; jj++) {
					if (root [j] EQ root [jj]) {
						in_same_block = TRUE;
						break;
					}
				}
			}
			if (in_same_block) {
				new_chosen [fi] = FALSE;
				continue;
			}
			for (j = 1; j < term_count; j++) {
				dsuf_unite (&fstsets,
					    dsuf_find (&fstsets, greedy_fsts [fi].p [0]),
					    dsuf_find (&fstsets, greedy_fsts [fi].p [j]));
			}
			total_smt_l += greedy_fsts [fi].len;
		}

		/* Go through edges in sorted order */
		for (i = 0; i < out.numberofedges; i++) {
			ei = sortededges [i];
			root [1] = dsuf_find (&fstsets, greedy_edges [ei].p1);
			root [2] = dsuf_find (&fstsets, greedy_edges [ei].p2);
			if (root [1] NE root [2]) {
				dsuf_unite (&fstsets, root [1], root [2]);
				/* Use BSD length here */
				total_smt_l += greedy_edges [ei].bsd;
			}
		}
		dsuf_destroy (&fstsets);

		if (total_smt_l < best_smt_l) {
			best_smt_l = total_smt_l;
			for (i = 0; i < fst_count; i++)
				greedy_fsts [i].chosen = new_chosen [i];
		}

	}

	/* Free all allocated arrays, including those allocated by Triangle */

	dsuf_destroy (&mstsets);

	free (pts);
	free (in.pointlist);
	free (out.pointlist);
	free (out.trianglelist);
	free (out.neighborlist);
	free (triangleedges);
	free (greedy_edges);
	free (greedy_fsts);
	free (sortededges);
	free (sortedfsts);
	free (new_chosen);

	return (best_smt_l);
}