Permalink
Browse files

initial commit

  • Loading branch information...
0 parents commit a51f7facf12431faee83361d1223fe1e4e634b41 Aaron Straup Cope committed Nov 25, 2009
Showing with 814 additions and 0 deletions.
  1. +19 −0 Makefile
  2. +67 −0 README
  3. +235 −0 clustr.cpp
  4. +188 −0 clustr.h
  5. +45 −0 component.cpp
  6. +41 −0 polygon.cpp
  7. +154 −0 shapefile.cpp
  8. +65 −0 shapefile.h
@@ -0,0 +1,19 @@
+# DEBUG := -ggdb
+CFLAGS := $(shell gdal-config --cflags) $(DEBUG) -O3 -funroll-loops
+LIBS := $(shell gdal-config --libs) -lCGAL
+objects = clustr.o component.o shapefile.o polygon.o
+
+clustr: $(objects)
+ g++ $(CFLAGS) -o $@ $+ $(LIBS)
+
+$(objects): %.o: %.cpp
+ g++ $(CFLAGS) -c $+
+
+shape_test: stuff/ogr_test.cpp shapefile.cpp
+ g++ -I. $(CFLAGS) -o $@ $+ $(LIBS)
+
+all: clustr
+
+clean:
+ -rm -f clustr
+ -rm -f *.o *.shp *.shx *.dbf *.prj *.gch *.rpo
@@ -0,0 +1,67 @@
+clustr - construct polygons from tagged points
+-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+written by Schuyler Erle <schuyler@nocat.net>
+(c) 2008 Yahoo! Inc.
+
+Overview
+--------
+Clustr takes a text file containing longitude/latitude points, tagged with
+a bit of text, and attempts to generate minimal polygons that "cover" those
+points, using an "alpha" parameter to determine the notion of "coverage".
+The polygons are written to an ESRI Shapefile, suitable for use in GIS
+software.
+
+Building clustr
+---------------
+You will need to have the GDAL library (http://www.gdal.org/) installed,
+with OGR support enabled, for writing the Shapefiles. You will also need
+the CGAL library (http://www.cgal.org/) installed.
+
+Depending on your system configuration, you may need to edit the Makefile
+to point to the location of the GDAL and CGAL libraries and C/C++ headers.
+
+To build clustr, simply run 'make'. The 'clustr' executable is
+self-contained (aside from shared libraries) and does not rely on any
+supporting files. The code has been tested with GCC 3.4 and 4.2.
+
+Using clustr
+------------
+
+$ clustr [-a <n>] [-p] [-v] <input> <output>
+
+ -h, -? this help message
+ -v be verbose (default: off)
+ -a <n> set alpha value (default: use optimal value)
+ -p output points to shapefile, instead of polygons
+
+If <input> is missing or given as "-", stdin is used.
+
+ * The input file should be formatted as: <tag> <lon> <lat>\n
+
+ * The input file *must* be sorted. Piping it through sort(1) will
+ suffice.
+
+ * Tags must not contain spaces.
+
+If <output> is missing, output is written to "clustr.shp". The output is
+always an ESRI Shapefile.
+
+The -p option simply converts the input file to a Shapefile containing the
+input points. This can be useful for comparing the polygon output of clustr
+to its input dataset in a GIS browser.
+
+The -a option sets the value of "alpha", which in turn determines the
+convexity of the output polygons. Larger values make more convex polygons,
+smaller values make less convex polygons. Smaller values will also filter
+out outlying points, which may be what you want.
+
+Omitting -a or setting it to zero will cause clustr to use the CGAL
+library's idea of "optimal alpha", which, if there are significant outliers
+in the input file, may not actually be what you want. In general, we
+recommend starting with a value of 0.001.
+
+Bugs
+----
+At present, clustr does not handle polygon holes correctly.
+
+- Fin -
@@ -0,0 +1,235 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <list>
+#include <stdexcept>
+#include <cmath>
+
+#include <unistd.h>
+#include "clustr.h"
+#include "shapefile.h"
+
+using namespace std;
+using namespace Clustr;
+
+void write_points_to_shapefile (Shapefile &shape, vector<Point> &pts,
+ string text) {
+ vector<Point>::iterator begin = pts.begin(),
+ end = pts.end();
+ for (vector<Point>::iterator it = begin; it != end; it++) {
+ Feature pt_feat(shape);
+ Geometry pt_geom(wkbPoint);
+ pt_geom.push_back(*it);
+ pt_feat.set(pt_geom);
+ pt_feat.set("tag", text.c_str());
+ shape.add_feature(pt_feat);
+ }
+}
+
+void write_polygon_to_shapefile (Shapefile &shape, Polygon &poly, string text, long int count) {
+ Feature feat(shape);
+ Geometry geom(wkbPolygon);
+ geom.insert_rings(poly.begin(), poly.end());
+ feat.set(geom);
+ feat.set("tag", text.c_str());
+ feat.set("count", count);
+ feat.set("area", poly.area());
+ feat.set("perimeter", poly.perimeter());
+ feat.set("density", ((float)count)/poly.area());
+ shape.add_feature(feat);
+}
+
+bool extract_alpha_component (Mesh::component &c, Polygon &poly, bool verbose)
+{
+ Mesh::vertex_circulator v(c);
+ Vertex_handle v0(*v);
+ Ring ring;
+
+ do {
+ ring.push_back(v->point());
+ } while (*(++v) != v0);
+
+ if (ring.size() < 4) {
+ if (verbose) cerr << "Discarding degenerate shape." << endl;
+ return false;
+ }
+ if (verbose)
+ cerr << "- " << ring.size() << " vertices, "
+ << "area: " << fabs(ring.area())
+ << ", perimeter: " << ring.perimeter() << endl;
+
+ poly.push_back(ring);
+ return true;
+}
+
+void construct_alpha_shape (Config &config, Shapefile &shape,
+ vector<Point> &pts, string text)
+{
+ coord_type alpha = config.alpha;
+ Mesh mesh(pts.begin(), pts.end(), alpha);
+
+ if (alpha == 0) {
+ if (config.verbose) cerr << "Computing optimal alpha... ";
+ alpha = *mesh.find_optimal_alpha(1);
+ if (config.verbose) cerr << alpha << endl;
+ mesh.set_alpha(alpha);
+ }
+
+ if (config.verbose)
+ cerr << mesh.number_of_solid_components() <<
+ " component(s) found for alpha value " << alpha << "." << endl;
+
+ for (Mesh::component c = mesh.components_begin();
+ c != mesh.components_end(); c++) {
+ Polygon poly;
+ if (extract_alpha_component(c, poly, config.verbose)) {
+ if (config.verbose)
+ cerr << "Writing polygon for tag '" << text << "'." << endl;
+ write_polygon_to_shapefile(shape, poly, text, pts.size());
+ }
+ }
+}
+
+void construct_output (Config &config, Shapefile &shape,
+ vector<Point> &pts, string tag) {
+ cerr << "Got " << pts.size() << " points for tag '"
+ << tag << "'." << endl;
+ if (config.points_only) {
+ write_points_to_shapefile(shape, pts, tag);
+ } else {
+ if (pts.size() >= 3) {
+ construct_alpha_shape(config, shape, pts, tag);
+ } else {
+ if (config.verbose)
+ cerr << "Need more than 3 points for tag '"
+ << tag << "'" << endl;
+ }
+ }
+}
+
+void display_usage (void) {
+ cerr << endl
+ << "clustr " << VERSION << " - construct polygons from tagged points"
+ << endl
+ << "written by " << AUTHOR << endl
+ << COPYRIGHT << endl
+ << endl
+ << "Usage: clustr [-a <n>] [-p] [-v] <input> <output>" << endl
+ << " -h, -? this help message" << endl
+ << " -v be verbose (default: off)" << endl
+ << " -a <n> set alpha value (default: use \"optimal\" value)"
+ << endl
+ << " -p output points to shapefile, instead of polygons"
+ << endl
+ << endl
+ << "If <input> is missing or given as \"-\", stdin is used."
+ << endl
+ << "If <output> is missing, output is written to clustr.shp."
+ << endl
+ << "Input file should be formatted as: <tag> <lon> <lat>\\n" << endl
+ << "Tags must not contain spaces." << endl
+ << endl;
+}
+
+bool parse_options (Config &config, int argc, char **argv) {
+ istringstream iss;
+ string optstring("a:pvh?");
+
+ if (argc <= 1) {
+ display_usage();
+ return false;
+ }
+
+ int opt = getopt( argc, argv, optstring.c_str() );
+ while( opt != -1 ) {
+ switch( opt ) {
+ case 'a':
+ iss.str((string) optarg);
+ iss >> config.alpha;
+ break;
+ case 'p':
+ config.points_only = true;
+ break;
+ case 'v':
+ config.verbose = true;
+ break;
+ case 'h': /* fall-through is intentional */
+ case '?':
+ display_usage();
+ return false;
+ }
+ opt = getopt( argc, argv, optstring.c_str() );
+ }
+ if (optind < argc) {
+ config.in_file = argv[optind++];
+ }
+ if (optind < argc) {
+ config.out_file = argv[optind];
+ }
+ return true;
+}
+
+int main(int argc, char **argv) {
+ vector<Point> points;
+ istream *input;
+ Config config;
+ Shapefile *shape;
+
+ if (!parse_options(config, argc, argv))
+ return 1;
+
+ if (config.in_file == "-") {
+ input = &cin;
+ } else {
+ input = new ifstream(config.in_file.c_str());
+ }
+
+ if (config.points_only) {
+ if (config.verbose) {
+ cerr << "Only converting points to Shapefile,"
+ << " no clustering performed." << endl;
+ }
+ shape = new Shapefile(config.out_file, wkbPoint);
+ shape->add_field("tag", OFTString, 64);
+
+ } else {
+ shape = new Shapefile(config.out_file, wkbPolygon);
+ shape->add_field("tag", OFTString, 64);
+ shape->add_field("count", OFTInteger, 8);
+ shape->add_field("area", OFTReal, 10, 2);
+ shape->add_field("perimeter", OFTReal, 10, 2);
+ shape->add_field("density", OFTReal, 10, 2);
+ }
+
+ string line, tag, previous = "";
+ vector<Point> pts;
+ Point pt;
+
+ if (config.verbose)
+ cerr << "Reading points from input." << endl;
+
+ while (!getline(*input,line).eof()) {
+ istringstream is(line);
+ if (!(is >> tag >> pt)) {
+ cerr << "Bad input line: " << line;
+ continue;
+ }
+ if (previous == "") {
+ previous = tag;
+ }
+ if (tag != previous) {
+ construct_output(config, *shape, pts, tag);
+ pts.clear();
+ }
+ pts.push_back(pt);
+ }
+ construct_output(config, *shape, pts, tag);
+
+ // if you want to do memory leak testing, include <ogr_api.h> and
+ // uncomment the following line:
+ // OGRCleanupAll();
+ delete shape;
+ return 0;
+}
+
Oops, something went wrong.

0 comments on commit a51f7fa

Please sign in to comment.