diff options
Diffstat (limited to 'blossom5-v2.05.src/GEOM/GPMinit.cpp')
| -rw-r--r-- | blossom5-v2.05.src/GEOM/GPMinit.cpp | 200 |
1 files changed, 200 insertions, 0 deletions
diff --git a/blossom5-v2.05.src/GEOM/GPMinit.cpp b/blossom5-v2.05.src/GEOM/GPMinit.cpp new file mode 100644 index 0000000..0e31231 --- /dev/null +++ b/blossom5-v2.05.src/GEOM/GPMinit.cpp @@ -0,0 +1,200 @@ +#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+
+
+
+
+// greedy procedure to make sure that a perfect matching exists:
+// 1. construct a matching among existing edges
+// (greedy procedure: pick a node, check whether there are edges leading
+// to unmatched nodes, if there are pick the edge with the smallest length).
+// 2. take remaining unmatched nodes, construct kd-tree for them,
+// assign an ordering to nodes (last visited time during left-most depth-first search),
+// add edges between consecutive nodes (2*i,2*i+1)
+void GeomPerfectMatching::CompleteInitialMatching()
+{
+ if (options.verbose) printf("adding edges to make sure that a perfect matching exists...");
+ PointId p, q;
+ Edge* e;
+ double len, len_min;
+ int unmatched_num = 0, edge_num0 = edge_num;
+
+ // construct greedy matching
+ for (p=0; p<node_num; p++)
+ {
+ if (nodes[p].is_marked) continue;
+ q = -1;
+ for (e=nodes[p].first[0]; e; e=e->next[0])
+ {
+ if (nodes[e->head[0]].is_marked) continue;
+ len = Dist2(p, e->head[0]);
+ if (q < 0 || len_min > len)
+ {
+ q = e->head[0];
+ len_min = len;
+ }
+ }
+ if (q >= 0)
+ {
+ nodes[p].is_marked = nodes[q].is_marked = 1;
+ }
+ else unmatched_num ++;
+ }
+
+ if (unmatched_num == 0)
+ {
+ for (p=0; p<node_num; p++) nodes[p].is_marked = 0;
+ return;
+ }
+
+ //printf("%d unmatched\n", unmatched_num);
+
+ REAL* unmatched_coords = new REAL[unmatched_num*DIM];
+ int* rev_mapping = new int[unmatched_num];
+ unmatched_num = 0;
+ for (p=0; p<node_num; p++)
+ {
+ if (nodes[p].is_marked) nodes[p].is_marked = 0;
+ else
+ {
+ memcpy(unmatched_coords+unmatched_num*DIM, coords+p*DIM, DIM*sizeof(REAL));
+ rev_mapping[unmatched_num ++] = p;
+ }
+ }
+
+ GPMKDTree* kd_tree = new GPMKDTree(DIM, unmatched_num, unmatched_coords, this);
+ kd_tree->AddPerfectMatching(rev_mapping);
+
+ delete kd_tree;
+ delete [] unmatched_coords;
+ delete [] rev_mapping;
+
+ if (options.verbose) printf("done (%d edges)\n", edge_num-edge_num0);
+}
+
+void GeomPerfectMatching::InitKNN(int K)
+{
+ if (node_num != node_num_max) { printf("InitKNN() cannot be called before all points have been added!\n"); exit(1); }
+ if (options.verbose) printf("adding K nearest neighbors (K=%d)\n", K);
+
+ int dir, k;
+ PointId p;
+ Edge* e;
+
+ if (K > node_num - 1) K = node_num - 1;
+
+ GPMKDTree* kd_tree = new GPMKDTree(DIM, node_num, coords, this);
+ PointId* neighbors = new PointId[K];
+
+ for (p=0; p<node_num; p++)
+ {
+ for (dir=0; dir<2; dir++)
+ for (e=nodes[p].first[dir]; e; e=e->next[dir])
+ {
+ nodes[e->head[dir]].is_marked = 1;
+ }
+
+ kd_tree->ComputeKNN(p, K, neighbors);
+ for (k=0; k<K; k++)
+ {
+ if (nodes[neighbors[k]].is_marked) continue;
+ AddInitialEdge(p, neighbors[k]);
+ nodes[neighbors[k]].is_marked = 1;
+ }
+
+ for (dir=0; dir<2; dir++)
+ for (e=nodes[p].first[dir]; e; e=e->next[dir])
+ {
+ nodes[e->head[dir]].is_marked = 0;
+ }
+ }
+
+ delete kd_tree;
+ delete [] neighbors;
+}
+
+#ifdef DELAUNAY_TRIANGLE
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4311)
+#pragma warning(disable: 4312)
+#endif
+
+extern "C" {
+#define ANSI_DECLARATORS
+#define TRILIBRARY
+#define NO_TIMER
+#define main NO_MAIN_FUNCTION
+#include "../triangle/triangle.c"
+}
+
+void GeomPerfectMatching::InitDelaunay()
+{
+ if (node_num < 16) return;
+ if (options.verbose) printf("adding edges in Delaunay triangulation\n");
+
+ int k;
+
+ struct triangulateio in, out, vorout;
+ in.numberofpoints = node_num;
+ in.numberofpointattributes = 0;
+ in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
+ for (k=0; k<2*node_num; k++) in.pointlist[k] = coords[k];
+ in.pointattributelist = NULL;
+ in.pointmarkerlist = NULL;
+ in.numberofsegments = 0;
+ in.numberofholes = 0;
+ in.numberofregions = 0;
+ in.regionlist = 0;
+
+ out.pointlist = (REAL *) NULL;
+ out.pointattributelist = (REAL *) NULL;
+ out.pointmarkerlist = (int *) NULL;
+ out.trianglelist = (int *) NULL;
+ out.triangleattributelist = (REAL *) NULL;
+ out.neighborlist = (int *) NULL;
+ out.segmentlist = (int *) NULL;
+ out.segmentmarkerlist = (int *) NULL;
+ out.edgelist = (int *) NULL;
+ out.edgemarkerlist = (int *) NULL;
+
+ vorout.pointlist = (REAL *) NULL;
+ vorout.pointattributelist = (REAL *) NULL;
+ vorout.edgelist = (int *) NULL;
+ vorout.normlist = (REAL *) NULL;
+
+ triangulate("pczAevn", &in, &out, &vorout);
+
+ free(in.pointlist);
+ free(out.pointlist);
+ free(out.pointmarkerlist);
+ free(out.trianglelist);
+ free(out.neighborlist);
+ free(out.segmentlist);
+ free(out.segmentmarkerlist);
+ free(out.edgemarkerlist);
+ free(vorout.pointlist);
+ free(vorout.pointattributelist);
+ free(vorout.edgelist);
+ free(vorout.normlist);
+
+ for (k=0; k<out.numberofedges; k++) AddInitialEdge(out.edgelist[2*k], out.edgelist[2*k+1]);
+
+ free(out.edgelist);
+}
+
+
+
+#else
+
+void GeomPerfectMatching::InitDelaunay()
+{
+ printf("You need to download the 'Triangle' software from \n\thttp://www.cs.cmu.edu/~quake/triangle.html ,\nextract it to the directory GeomPerfectMatching and define DELAUNARY_TRIANG in GeomPerfectMatching.h\n");
+ exit(1);
+}
+
+#endif
+
|