[prepair] 02/06: New upstream version 0.7+20160519-417f1bc
Bas Couwenberg
sebastic at debian.org
Sat Jul 1 16:06:09 UTC 2017
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository prepair.
commit f90374434bd617a92a4ec94ff8ff303336a9d93f
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Sat Jul 1 17:43:01 2017 +0200
New upstream version 0.7+20160519-417f1bc
---
CMakeLists.txt | 14 +++--
PolygonRepair.cpp | 179 +++++++++++-------------------------------------------
PolygonRepair.h | 7 +--
README.md | 45 +++++++++-----
TriangleInfo.cpp | 29 ++++-----
TriangleInfo.h | 9 +--
definitions.h | 39 ++++++------
prepair.cpp | 53 ++++++++++------
8 files changed, 154 insertions(+), 221 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 730e8d6..b1fff2a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,4 +1,4 @@
-# Copyright (c) 2009-2013,
+# Copyright (c) 2009-2014,
# Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
# Hugo Ledoux h.ledoux at tudelft.nl
# Martijn Meijers b.m.meijers at tudelft.nl
@@ -47,13 +47,15 @@ endif()
include( ${CGAL_USE_FILE} )
# Boost
-find_package( Boost REQUIRED )
+find_package( Boost REQUIRED COMPONENTS program_options)
if ( NOT Boost_FOUND )
message(SEND_ERROR "prepair requires the Boost library")
return()
endif()
+INCLUDE_DIRECTORIES( ${Boost_INCLUDE_DIR} )
+
# GDAL
find_package( GDAL )
@@ -63,16 +65,20 @@ endif()
include_directories( ${GDAL_INCLUDE_DIR} )
+# Hardening buildflags
+set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} $ENV{CXXFLAGS}")
+set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} $ENV{LDFLAGS}")
+
# Creating entries for target: prepair
# ############################
if ( AS_LIBRARY )
add_library( prepair SHARED TriangleInfo.cpp PolygonRepair.cpp )
else()
- add_executable( prepair TriangleInfo.cpp PolygonRepair.cpp prepair.cpp )
+ add_executable( prepair TriangleInfo.cpp PolygonRepair.cpp prepair.cpp )
endif()
add_to_cached_list( CGAL_EXECUTABLE_TARGETS prepair )
# Link to CGAL and third-party libraries
-target_link_libraries(prepair ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY})
\ No newline at end of file
+target_link_libraries(prepair ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY} ${Boost_LIBRARIES})
\ No newline at end of file
diff --git a/PolygonRepair.cpp b/PolygonRepair.cpp
index 4aa948a..0f84165 100644
--- a/PolygonRepair.cpp
+++ b/PolygonRepair.cpp
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -106,146 +106,6 @@ void PolygonRepair::removeSmallPolygons(OGRMultiPolygon *outPolygons, double min
}
}
-OGRMultiPolygon *PolygonRepair::isr(OGRGeometry *geometry, double tolerance) {
-
- // Convert to multipolygon
- OGRMultiPolygon *multipolygon = new OGRMultiPolygon();
- switch (geometry->getGeometryType()) {
- case wkbLineString: {
- OGRLineString *oldRing = static_cast<OGRLineString *>(geometry);
- OGRPolygon *polygon = new OGRPolygon();
- OGRLinearRing *ring = new OGRLinearRing();
- OGRPoint lastPoint, point;
- for (int currentPoint = 0; currentPoint < oldRing->getNumPoints(); ++currentPoint) {
- oldRing->getPoint(currentPoint, &point);
- if (currentPoint > 0) {
- if (!point.Equals(&lastPoint)) ring->addPoint(&point);
- } lastPoint = point;
- } polygon->addRingDirectly(ring);
- multipolygon->addGeometryDirectly(polygon);
- break;
- }
-
- case wkbPolygon: {
- OGRPolygon *oldPolygon = static_cast<OGRPolygon *>(geometry);
- OGRPolygon *polygon = new OGRPolygon();
- OGRLinearRing *oldRing = oldPolygon->getExteriorRing();
- OGRLinearRing *ring = new OGRLinearRing();
- OGRPoint lastPoint, point;
- for (int currentPoint = 0; currentPoint < oldRing->getNumPoints(); ++currentPoint) {
- oldRing->getPoint(currentPoint, &point);
- if (currentPoint > 0) {
- if (!point.Equals(&lastPoint)) ring->addPoint(&point);
- } lastPoint = point;
- } polygon->addRingDirectly(ring);
- for (int currentRing = 0; currentRing < oldPolygon->getNumInteriorRings(); ++currentRing) {
- oldRing = oldPolygon->getInteriorRing(currentRing);
- ring = new OGRLinearRing();
- for (int currentPoint = 0; currentPoint < oldRing->getNumPoints(); ++currentPoint) {
- oldRing->getPoint(currentPoint, &point);
- if (currentPoint > 0) {
- if (!point.Equals(&lastPoint)) ring->addPoint(&point);
- } lastPoint = point;
- } polygon->addRingDirectly(ring);
- } multipolygon->addGeometryDirectly(polygon);
- break;
- }
-
- case wkbMultiPolygon: {
- OGRMultiPolygon *oldMultipolygon = static_cast<OGRMultiPolygon *>(geometry);
- OGRPolygon *oldPolygon, *polygon;
- OGRLinearRing *oldRing, *ring;
- OGRPoint lastPoint, point;
- for (int currentPolygon = 0; currentPolygon < oldMultipolygon->getNumGeometries(); ++currentPolygon) {
- oldPolygon = static_cast<OGRPolygon *>(oldMultipolygon->getGeometryRef(currentPolygon));
- polygon = new OGRPolygon();
- oldRing = oldPolygon->getExteriorRing();
- ring = new OGRLinearRing();
- for (int currentPoint = 0; currentPoint < oldRing->getNumPoints(); ++currentPoint) {
- oldRing->getPoint(currentPoint, &point);
- if (currentPoint > 0) {
- if (!point.Equals(&lastPoint)) ring->addPoint(&point);
- } lastPoint = point;
- } polygon->addRingDirectly(ring);
- for (int currentRing = 0; currentRing < oldPolygon->getNumInteriorRings(); ++currentRing) {
- oldRing = oldPolygon->getInteriorRing(currentRing);
- ring = new OGRLinearRing();
- for (int currentPoint = 0; currentPoint < oldRing->getNumPoints(); ++currentPoint) {
- oldRing->getPoint(currentPoint, &point);
- if (currentPoint > 0) {
- if (!point.Equals(&lastPoint)) ring->addPoint(&point);
- } lastPoint = point;
- } polygon->addRingDirectly(ring);
- } multipolygon->addGeometryDirectly(polygon);
- } break;
- }
-
- default:
- std::cerr << "PolygonRepair::isr: Cannot understand input." << std::endl;
- return new OGRMultiPolygon();
- break;
- } multipolygon->closeRings();
-
- // Put all line segments in a list (in order)
- // TODO: Only do point location once for segment.
- std::list<Segment> inputSegments;
- for (int currentPolygon = 0; currentPolygon < multipolygon->getNumGeometries(); ++currentPolygon) {
- OGRPolygon *polygon = static_cast<OGRPolygon *>(multipolygon->getGeometryRef(currentPolygon));
- for (int currentPoint = 1; currentPoint < polygon->getExteriorRing()->getNumPoints(); ++currentPoint) {
- Segment thisSegment(Point(polygon->getExteriorRing()->getX(currentPoint-1), polygon->getExteriorRing()->getY(currentPoint-1)),
- Point(polygon->getExteriorRing()->getX(currentPoint), polygon->getExteriorRing()->getY(currentPoint)));
- if (thisSegment.is_degenerate()) {
- std::cout << "This should not happen!!!" << std::endl;
- return new OGRMultiPolygon();
- } inputSegments.push_back(thisSegment);
- } for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
- for (int currentPoint = 1; currentPoint < polygon->getInteriorRing(currentRing)->getNumPoints(); ++currentPoint) {
- Segment thisSegment(Point(polygon->getExteriorRing()->getX(currentPoint-1), polygon->getExteriorRing()->getY(currentPoint-1)),
- Point(polygon->getExteriorRing()->getX(currentPoint), polygon->getExteriorRing()->getY(currentPoint)));
- if (thisSegment.is_degenerate()) {
- std::cout << "This should not happen!!!" << std::endl;
- return new OGRMultiPolygon();
- } std::cout << thisSegment << std::endl;
- inputSegments.push_back(thisSegment);
- }
- }
- }
-
- // Snap
- std::list<std::list<Point> > outputPolylines;
- CGAL::snap_rounding_2<SRT, std::list<Segment>::const_iterator, std::list<std::list<Point> > >(inputSegments.begin(), inputSegments.end(), outputPolylines, tolerance, true, false, 2);
-
- // Put it back into an OGRMultiPolygon
- OGRMultiPolygon *snappedMultipolygon = new OGRMultiPolygon();
- std::list<std::list<Point> >::iterator currentPolyline = outputPolylines.begin();
- for (int currentPolygon = 0; currentPolygon < multipolygon->getNumGeometries(); ++currentPolygon) {
- OGRPolygon *polygon = static_cast<OGRPolygon *>(multipolygon->getGeometryRef(currentPolygon));
- OGRPolygon *snappedPolygon = new OGRPolygon();
- OGRLinearRing *snappedOuterRing = new OGRLinearRing();
- for (int currentPoint = 1; currentPoint < polygon->getExteriorRing()->getNumPoints(); ++currentPoint) {
- std::list<Point>::iterator point = currentPolyline->begin();
- ++point;
- while (point != currentPolyline->end()) {
- snappedOuterRing->addPoint(CGAL::to_double(point->x()), CGAL::to_double(point->y()));
- ++point;
- } ++currentPolyline;
- } snappedPolygon->addRingDirectly(snappedOuterRing);
- for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
- OGRLinearRing *snappedInnerRing = new OGRLinearRing();
- for (int currentPoint = 1; currentPoint < polygon->getInteriorRing(currentRing)->getNumPoints(); ++currentPoint) {
- std::list<Point>::iterator point = currentPolyline->begin();
- ++point;
- while (point != currentPolyline->end()) {
- snappedOuterRing->addPoint(CGAL::to_double(point->x()), CGAL::to_double(point->y()));
- ++point;
- } ++currentPolyline;
- } snappedPolygon->addRingDirectly(snappedInnerRing);
- } snappedMultipolygon->OGRGeometryCollection::addGeometryDirectly(snappedPolygon);
- }
-
- return snappedMultipolygon;
-}
-
double PolygonRepair::computeRobustness(OGRGeometry *geometry) {
double smallestdist = 1e99;
@@ -305,22 +165,44 @@ double PolygonRepair::computeRobustness(OGRGeometry *geometry) {
bool PolygonRepair::saveToShp(OGRGeometry* geometry, const char *fileName) {
const char *driverName = "ESRI Shapefile";
+#if GDAL_VERSION_MAJOR < 2
OGRRegisterAll();
OGRSFDriver *driver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(driverName);
+#else
+ GDALAllRegister();
+ GDALDriver *driver = GetGDALDriverManager()->GetDriverByName(driverName);
+#endif
if (driver == NULL) {
std::cout << "\tError: OGR Shapefile driver not found." << std::endl;
return false;
}
+
+#if GDAL_VERSION_MAJOR < 2
OGRDataSource *dataSource = driver->Open(fileName, false);
+#else
+ GDALDataset *dataSource = (GDALDataset*) GDALOpenEx(fileName, GDAL_OF_READONLY, NULL, NULL, NULL);
+#endif
if (dataSource != NULL) {
std::cout << "\tOverwriting file..." << std::endl;
+#if GDAL_VERSION_MAJOR < 2
if (driver->DeleteDataSource(dataSource->GetName())!= OGRERR_NONE) {
+#else
+ if (driver->Delete(fileName)!= CE_None) {
+#endif
std::cout << "\tError: Couldn't erase file with same name." << std::endl;
return false;
+#if GDAL_VERSION_MAJOR < 2
} OGRDataSource::DestroyDataSource(dataSource);
+#else
+ } GDALClose(dataSource);
+#endif
}
std::cout << "\tCreating " << fileName << std::endl;
+#if GDAL_VERSION_MAJOR < 2
dataSource = driver->CreateDataSource(fileName, NULL);
+#else
+ dataSource = driver->Create(fileName,0,0,0,GDT_Unknown,NULL);
+#endif
if (dataSource == NULL) {
std::cout << "\tError: Could not create file." << std::endl;
return false;
@@ -337,7 +219,11 @@ bool PolygonRepair::saveToShp(OGRGeometry* geometry, const char *fileName) {
std::cout << "\tError: Could not create feature." << std::endl;
}
OGRFeature::DestroyFeature(feature);
+#if GDAL_VERSION_MAJOR < 2
OGRDataSource::DestroyDataSource(dataSource);
+#else
+ GDALClose(dataSource);
+#endif
return true;
}
@@ -358,6 +244,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
vb = triangulation.insert(Point(ring->getX(currentPoint),
ring->getY(currentPoint)),
triangulation.incident_faces(va));
+ if (va == vb) continue;
if (removeOverlappingConstraints && triangulation.is_edge(va, vb, faceOfEdge, indexOfEdge)) {
if (triangulation.is_constrained(std::pair<Triangulation::Face_handle, int>(faceOfEdge, indexOfEdge))) {
triangulation.insert_constraint(va, vb); // trick to remove a partially overlapping constraint
@@ -371,12 +258,14 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
OGRPolygon *polygon = static_cast<OGRPolygon *>(geometry);
// Outer
+ polygon->getExteriorRing()->closeRings();
vb = triangulation.insert(Point(polygon->getExteriorRing()->getX(0), polygon->getExteriorRing()->getY(0)));
for (int currentPoint = 1; currentPoint < polygon->getExteriorRing()->getNumPoints(); ++currentPoint) {
va = vb;
vb = triangulation.insert(Point(polygon->getExteriorRing()->getX(currentPoint),
polygon->getExteriorRing()->getY(currentPoint)),
triangulation.incident_faces(va));
+ if (va == vb) continue;
if (removeOverlappingConstraints && triangulation.is_edge(va, vb, faceOfEdge, indexOfEdge)) {
if (triangulation.is_constrained(std::pair<Triangulation::Face_handle, int>(faceOfEdge, indexOfEdge))) {
triangulation.insert_constraint(va, vb);
@@ -394,6 +283,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
// Inner
for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
+ polygon->getInteriorRing(currentRing)->closeRings();
vb = triangulation.insert(Point(polygon->getInteriorRing(currentRing)->getX(0),
polygon->getInteriorRing(currentRing)->getY(0)));
for (int currentPoint = 1; currentPoint < polygon->getInteriorRing(currentRing)->getNumPoints(); ++currentPoint) {
@@ -401,6 +291,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
vb = triangulation.insert(Point(polygon->getInteriorRing(currentRing)->getX(currentPoint),
polygon->getInteriorRing(currentRing)->getY(currentPoint)),
triangulation.incident_faces(va));
+ if (va == vb) continue;
if (removeOverlappingConstraints && triangulation.is_edge(va, vb, faceOfEdge, indexOfEdge)) {
if (triangulation.is_constrained(std::pair<Triangulation::Face_handle, int>(faceOfEdge, indexOfEdge))) {
triangulation.insert_constraint(va, vb);
@@ -425,6 +316,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
OGRPolygon *polygon = static_cast<OGRPolygon *>(multipolygon->getGeometryRef(currentPolygon));
// Outer
+ polygon->getExteriorRing()->closeRings();
vb = triangulation.insert(Point(polygon->getExteriorRing()->getX(0),
polygon->getExteriorRing()->getY(0)));
for (int currentPoint = 1; currentPoint < polygon->getExteriorRing()->getNumPoints(); ++currentPoint) {
@@ -432,6 +324,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
vb = triangulation.insert(Point(polygon->getExteriorRing()->getX(currentPoint),
polygon->getExteriorRing()->getY(currentPoint)),
triangulation.incident_faces(va));
+ if (va == vb) continue;
if (removeOverlappingConstraints && triangulation.is_edge(va, vb, faceOfEdge, indexOfEdge)) {
if (triangulation.is_constrained(std::pair<Triangulation::Face_handle, int>(faceOfEdge, indexOfEdge))) {
triangulation.insert_constraint(va, vb);
@@ -442,6 +335,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
// Inner
for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
+ polygon->getInteriorRing(currentRing)->closeRings();
vb = triangulation.insert(Point(polygon->getInteriorRing(currentRing)->getX(0),
polygon->getInteriorRing(currentRing)->getY(0)));
for (int currentPoint = 1; currentPoint < polygon->getInteriorRing(currentRing)->getNumPoints(); ++currentPoint) {
@@ -449,6 +343,7 @@ void PolygonRepair::insertConstraints(Triangulation &triangulation, OGRGeometry*
vb = triangulation.insert(Point(polygon->getInteriorRing(currentRing)->getX(currentPoint),
polygon->getInteriorRing(currentRing)->getY(currentPoint)),
triangulation.incident_faces(va));
+ if (va == vb) continue;
if (removeOverlappingConstraints && triangulation.is_edge(va, vb, faceOfEdge, indexOfEdge)) {
if (triangulation.is_constrained(std::pair<Triangulation::Face_handle, int>(faceOfEdge, indexOfEdge))) {
triangulation.insert_constraint(va, vb);
diff --git a/PolygonRepair.h b/PolygonRepair.h
index 0f867c5..137251a 100644
--- a/PolygonRepair.h
+++ b/PolygonRepair.h
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -28,7 +28,6 @@ public:
OGRMultiPolygon *repairOddEven(OGRGeometry *geometry, bool timeResults = false);
OGRMultiPolygon *repairPointSet(OGRGeometry *geometry, bool timeResults = false);
void removeSmallPolygons(OGRMultiPolygon *outPolygons, double minArea);
- OGRMultiPolygon *isr(OGRGeometry *geometry, double tolerance);
double computeRobustness(OGRGeometry *geometry = NULL);
bool saveToShp(OGRGeometry* geometry, const char *fileName);
@@ -47,4 +46,4 @@ private:
void printChain(std::list<Triangulation::Vertex_handle> &chain);
};
-#endif
\ No newline at end of file
+#endif
diff --git a/README.md b/README.md
index 7805c7f..a674cf9 100644
--- a/README.md
+++ b/README.md
@@ -1,41 +1,55 @@
## What is prepair?
-prepair -- pronounce 'pee-repair' as in 'polygon repair' -- permits us to easily repair "broken" GIS polygons, and that according to the international standards ISO 19107 (Geographic information -- Spatial schema). In brief, given a polygon, it *automatically* repairs it and returns back a valid polygon (actually a MultiPolygon if the input represents more than one polygon). Note that this project is only concerned with single polygons, if you're interested in validating how different po [...]
+prepair — pronounce 'pee-repair' as in 'polygon repair' — permits us to easily repair "broken" GIS polygons according to the international standard [ISO19107](http://www.iso.org/iso/catalogue_detail.htm?csnumber=26012) (Geographic information — Spatial schema). Given one input polygon, it *automatically* repairs it and returns back a valid polygon (actually a MultiPolygon since the input can represent more than one polygon — think of a 'bowtie' for instance).
Automated repair methods can be considered as interpreting ambiguous or ill-defined polygons and giving a coherent and clearly defined output. Examples of errors are: polygon has a dangling edge; polygon is not closed; polygon self-intersects; an inner ring of the polygon is located outside the outer ring; etc.
-prepair performs more or less the same as the new PostGIS 2.0's function [ST_MakeValid()](http://postgis.org/documentation/manual-svn/ST_MakeValid.html), but is faster, scales better to massive polygons, and predicting its behaviour is simple (so one can guess how her polygons will be repaired).
-We have implemented two repair paradigms: an extension of the [odd-even algorithm](https://en.wikipedia.org/wiki/Even-odd_rule) to handle GIS polygons containing inner rings and degeneracies; and one where we follow a *set difference* rule for the rings.
+prepair performs more or less the same as the PostGIS 2.0's function [ST_MakeValid()](http://postgis.org/documentation/manual-svn/ST_MakeValid.html), but is faster, scales better to massive polygons, and predicting its behaviour is simple (so one can guess how her polygons will be repaired).
+We have implemented two repair paradigms:
+
+ 1. an extension of the [odd-even algorithm](https://en.wikipedia.org/wiki/Even-odd_rule) to handle GIS polygons containing inner rings and degeneracies;
+ 2. setdiff: one where we follow a *point set difference* rule for the rings (outer - inner).
prepair is based on a constrained triangulation ([CGAL](http://www.cgal.org) is used) and [OGR](http://www.gdal.org/ogr/) is used to read/write WKT.
It is available under a dual license scheme: [GPLv3](http://www.gnu.org/copyleft/gpl.html) and commercial. If you are interested in a commercial license, please contact [Ken Arroyo Ohori](mailto:g.a.k.arroyoohori at tudelft.nl).
+Note that prepair is only concerned with single polygons, and if you're interested in validating how different polygons interact with each other (to be precise: to check if they form a planar partition) have a look at our other project [pprepair](https://github.com/tudelft3d/pprepair).
+
+
## Details
Details of how we automatically repair broken polygons, and what results you can expect, are available in this scientific article:
-
-> Ledoux, H., Arroyo Ohori, K., and Meijers, M. (2014). A triangulation-based approach to automatically repair GIS polygons. *Computers & Geosciences* 66:121-–131. [ [DOI] ](http://dx.doi.org/10.1016/j.cageo.2014.01.009) [ [PDF] ](http://homepage.tudelft.nl/23t4p/pdfs/14cgeo_prepair.pdf)
+> Ledoux, H., Arroyo Ohori, K., and Meijers, M. (2014). A triangulation-based approach to automatically repair GIS polygons. *Computers & Geosciences* 66:121–131. [ [DOI] ](http://dx.doi.org/10.1016/j.cageo.2014.01.009) [ [PDF] ](http://3dgeoinfo.bk.tudelft.nl/hledoux/pdfs/14_cgeo_prepair.pdf)
If you use prepair for a scientific project, please cite this article.
-## How to compile?
-prepair is provided as source code, which is very easy to compile it on Linux and Mac. We plan on offering binaries (including for Windows) in the future.
+## How to get it?
+
+prepair is provided as source code or as 64-bit binaries for [Windows](https://github.com/tudelft-gist/prepair/releases/download/v0.7/prepair_win64.zip) and [Mac](https://github.com/tudelft-gist/prepair/releases/download/v0.7/prepair_mac.zip). The Mac binary requires Kyngchaos' [GDAL 1.11 Complete Framework](http://www.kyngchaos.com/software/frameworks#gdal_complete).
-To compile it, you first need to install the following three (free) libraries:
+prepair is also very easy to compile on Mac and Linux using the included CMake file. It should also work on other Unix-like systems. To compile prepair, you need to have a recent version of the following three (free) libraries:
1. [CGAL](http://www.cgal.org)
2. [OGR](http://www.gdal.org/ogr/)
3. [CMake](http://www.cmake.org)
-Afterwards run:
+Under Mac, if you use Kyngchaos' GDAL Complete Framework, which is used by QGIS, you already have OGR installed. If you need them, a good way to install CGAL and OGR is to use [Homebrew](http://brew.sh):
+
+ $ brew install gdal
+ $ brew install cgal
+
+Once all the dependencies are met, just generate the makefile for your system and compile:
$ cmake .
$ make
-## It's a command-line program only
+
+## How to run it?
+
+You can run prepair from the command-line or through our [QGIS plug-in](https://github.com/tudelft-gist/prepair-qgis), which you can get from the official QGIS repository.
A [WKT](http://en.wikipedia.org/wiki/Well-known_text) or an OGR dataset (shapefile, geojson or GML for instance) is read as input, and a WKT or a shapefile (a MultiPolygon) is given as output:
@@ -48,12 +62,14 @@ A [WKT](http://en.wikipedia.org/wiki/Well-known_text) or an OGR dataset (shapefi
$ ./prepair --shpOut --ogr data/CLC2006_180927.geojson
Creating out.shp
-[Snap rounding](http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Snap_rounding_2/Chapter_main.html) of the input segments can be performed:
- $ /prepair --isr 2 --wkt "POLYGON((0 0, 10 0, 15 5, 10 0, 10 10, 0 10, 0 0))"
+[Snap rounding](http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Snap_rounding_2/Chapter_main.html) of the input segments can be performed with the --isr option:
+
+ $ ./prepair --isr 2 --wkt "POLYGON((0 0, 10 0, 15 5, 10 0, 10 10, 0 10, 0 0))"
MULTIPOLYGON (((11 1,11 11,1 11,1 1,11 1)))
-It's possible to remove small (sliver) polygons in the output by giving the smallest area allowed:
+
+It's also possible to remove small (sliver) polygons in the output by giving the smallest area allowed with the --minarea option:
$ ./prepair --wkt 'POLYGON((0 0, 10 0, 10 11, 11 10, 0 10))'
MULTIPOLYGON (((10 0,10 10,0 10,0 0,10 0)),((11 10,10 11,10 10,11 10)))
@@ -62,9 +78,10 @@ It's possible to remove small (sliver) polygons in the output by giving the smal
Removing polygons smaller than 1 unit^2.
MULTIPOLYGON (((10 0,10 10,0 10,0 0,10 0)))
+
## Examples of invalid input you can try
-The folder 'data' contains examples of relatively big invalid polygons. These are from the [Corine Land Cover 2006 dataset.](http://sia.eionet.europa.eu/CLC2006)
+The folder 'data' contains examples of relatively big invalid polygons. These are from the [Corine Land Cover 2006 dataset](http://sia.eionet.europa.eu/CLC2006).
A 'bowtie' polygon:
diff --git a/TriangleInfo.cpp b/TriangleInfo.cpp
index 68cd6d4..32d9511 100644
--- a/TriangleInfo.cpp
+++ b/TriangleInfo.cpp
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -22,45 +22,42 @@
#include "TriangleInfo.h"
TriangleInfo::TriangleInfo() {
- info = 0x00;
+
}
void TriangleInfo::clear() {
- info = 0x00;
+ info.reset();
}
bool TriangleInfo::beenTagged() {
- return (info & 0x01) == 0x01;
+ return info[0];
}
void TriangleInfo::beenTagged(bool tagged) {
- if (tagged) info |= 0x01;
- else info &= 0xfe;
+ info[0] = tagged;
}
bool TriangleInfo::isInInterior() {
- return (info & 0x02) == 0x02;
+ return info[1];
}
void TriangleInfo::isInInterior(bool inInterior) {
- if (inInterior) info |= 0x03;
- else info = (info & 0xfd) | 0x01;
+ info[0] = true;
+ info[1] = inInterior;
}
bool TriangleInfo::isOnBorder() {
- return (info & 0x04) == 0x04;
+ return info[2];
}
void TriangleInfo::isOnBorder(bool onBorder) {
- if (onBorder) info |= 0x04;
- else info &= 0xfb;
+ info[2] = onBorder;
}
bool TriangleInfo::beenReconstructed() {
- return (info & 0x08) == 0x08;
+ return info[3];
}
void TriangleInfo::beenReconstructed(bool reconstructed) {
- if (reconstructed) info |= 0x08;
- else info &= 0xf7;
+ info[3] = reconstructed;
}
\ No newline at end of file
diff --git a/TriangleInfo.h b/TriangleInfo.h
index df2045e..6fe8393 100644
--- a/TriangleInfo.h
+++ b/TriangleInfo.h
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -22,6 +22,8 @@
#ifndef TRIANGLEINFO_H
#define TRIANGLEINFO_H
+#include <bitset>
+
// TODO: Maybe create a marker manager in the future
class TriangleInfo {
public:
@@ -38,8 +40,7 @@ public:
void beenReconstructed(bool reconstructed);
private:
- unsigned char info; // Bitset: unused unused unused unused
- // not reconstructed/reconstructed normal/border exterior/interior untagged/tagged
+ std::bitset<4> info;
};
#endif
diff --git a/definitions.h b/definitions.h
index e2cb447..dddbf88 100644
--- a/definitions.h
+++ b/definitions.h
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -22,26 +22,20 @@
// Compile-time options
// * if the code crashes, try to compile with EXACT_CONSTRUCTIONS so that
// robust arithmetic is used
-// * switch to the setdiff paradigm by defining SET_DIFF
#define EXACT_CONSTRUCTIONS
+//#define COORDS_3D
#ifndef DEFINITIONS_H
#define DEFINITIONS_H
#include "TriangleInfo.h"
+// OGR
+#include <ogrsf_frmts.h>
+
// STL
-#include <iostream>
-#include <stack>
-#include <set>
-#include <list>
#include <fstream>
-#include <math.h>
-#include <time.h>
-
-// OGR
-#include <gdal/ogrsf_frmts.h>
// CGAL
#ifdef EXACT_CONSTRUCTIONS
@@ -49,18 +43,26 @@
#else
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#endif
+#ifdef COORDS_3D
+#include <CGAL/Projection_traits_xy_3.h>
+#endif
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_hierarchy_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
-#include <CGAL/Snap_rounding_traits_2.h>
-#include <CGAL/Snap_rounding_2.h>
// Kernel
#ifdef EXACT_CONSTRUCTIONS
-typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+typedef CGAL::Exact_predicates_exact_constructions_kernel TK;
+typedef CGAL::Exact_intersections_tag ET;
+#else
+typedef CGAL::Exact_predicates_inexact_constructions_kernel TK;
+typedef CGAL::Exact_predicates_tag ET;
+#endif
+#ifdef COORDS_3D
+typedef CGAL::Projection_traits_xy_3<TK> K;
#else
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef TK K;
#endif
typedef CGAL::Triangulation_vertex_base_2<K> TVB;
@@ -68,12 +70,9 @@ typedef CGAL::Triangulation_hierarchy_vertex_base_2<TVB> VB;
typedef CGAL::Constrained_triangulation_face_base_2<K> FB;
typedef CGAL::Triangulation_face_base_with_info_2<TriangleInfo, K, FB> FBWI;
typedef CGAL::Triangulation_data_structure_2<VB, FBWI> TDS;
-typedef CGAL::Exact_predicates_tag PT;
-typedef CGAL::Exact_intersections_tag IT;
-typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, PT> CDT;
+typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, ET> CDT;
typedef CGAL::Triangulation_hierarchy_2<CDT> CDTH;
typedef CGAL::Constrained_triangulation_plus_2<CDTH> Triangulation;
-typedef CGAL::Snap_rounding_traits_2<K> SRT;
typedef Triangulation::Point Point;
typedef K::Segment_2 Segment;
diff --git a/prepair.cpp b/prepair.cpp
index 5528a5c..75b292f 100644
--- a/prepair.cpp
+++ b/prepair.cpp
@@ -1,11 +1,11 @@
/*
- Copyright (c) 2009-2013,
+ Copyright (c) 2009-2014,
Ken Arroyo Ohori g.a.k.arroyoohori at tudelft.nl
Hugo Ledoux h.ledoux at tudelft.nl
Martijn Meijers b.m.meijers at tudelft.nl
All rights reserved.
- // This file is part of prepair: you can redistribute it and/or modify
+ This file is part of prepair: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
@@ -81,12 +81,32 @@ int main (int argc, const char * argv[]) {
//-- reading from WKT passed directly
else if (strcmp(argv[argNum], "--wkt") == 0) {
unsigned int bufferSize = 100000000;
- char *inputWKT = (char *)malloc(bufferSize*sizeof(char *));
+ char *inputWKT = (char *)malloc(bufferSize*sizeof(char));
strcpy(inputWKT, argv[argNum+1]);
++argNum;
- OGRGeometryFactory::createFromWkt(&inputWKT, NULL, &geometry);
- if (geometry == NULL) {
- std::cout << "Error: WKT is not valid" << std::endl;
+ OGRErr err = OGRGeometryFactory::createFromWkt(&inputWKT, NULL, &geometry);
+ if (err != OGRERR_NONE) {
+ switch (err) {
+ case OGRERR_UNSUPPORTED_GEOMETRY_TYPE:
+ std::cerr << "Error: geometry must be Polygon or MultiPolygon" << std::endl;
+ break;
+ case OGRERR_NOT_ENOUGH_DATA:
+ case OGRERR_CORRUPT_DATA:
+ std::cerr << "Error: corrupted input" << std::endl;
+ break;
+ default:
+ std::cerr << "Error: corrupted input" << std::endl;
+ break;
+ }
+ return 1;
+ }
+ if (geometry->IsEmpty() == 1) {
+ std::cerr << "Error: empty geometry" << std::endl;
+ return 1;
+ }
+ if ( (geometry->getGeometryType() != wkbPolygon) &&
+ (geometry->getGeometryType() != wkbMultiPolygon) ) {
+ std::cerr << "Error: geometry must be Polygon or MultiPolygon" << std::endl;
return 1;
}
}
@@ -94,7 +114,7 @@ int main (int argc, const char * argv[]) {
//-- reading from WKT stored in first line of a text file
else if (strcmp(argv[argNum], "-f") == 0) {
unsigned int bufferSize = 100000000;
- char *inputWKT = (char *)malloc(bufferSize*sizeof(char *));
+ char *inputWKT = (char *)malloc(bufferSize*sizeof(char));
if (argNum + 1 <= argc - 1 && argv[argNum+1][0] != '-') {
std::ifstream infile(argv[argNum+1], std::ifstream::in);
infile.getline(inputWKT, bufferSize);
@@ -112,8 +132,13 @@ int main (int argc, const char * argv[]) {
//-- reading from a ogr dataset (most supported: shp, geojson, gml, etc)
else if (strcmp(argv[argNum], "--ogr") == 0) {
+#if GDAL_VERSION_MAJOR < 2
OGRRegisterAll();
OGRDataSource *dataSource = OGRSFDriverRegistrar::Open(argv[argNum+1], false);
+#else
+ GDALAllRegister();
+ GDALDataset *dataSource = (GDALDataset*) GDALOpenEx(argv[argNum+1], GDAL_OF_READONLY, NULL, NULL, NULL);
+#endif
++argNum;
if (dataSource == NULL) {
std::cerr << "Error: Could not open file." << std::endl;
@@ -151,19 +176,12 @@ int main (int argc, const char * argv[]) {
if (computeRobustness == true)
std::cout << "Robustness of input polygon: " << sqrt(prepair.computeRobustness(geometry)) <<std::endl;
- OGRGeometry *snappedGeometry;
- if (isrTolerance != 0) {
- snappedGeometry = prepair.isr(geometry, isrTolerance);
- } else {
- snappedGeometry = geometry;
- }
-
OGRMultiPolygon *outPolygons;
if (pointSet) {
- outPolygons = prepair.repairPointSet(snappedGeometry, startTime);
+ outPolygons = prepair.repairPointSet(geometry, startTime);
} else {
- outPolygons = prepair.repairOddEven(snappedGeometry, startTime);
+ outPolygons = prepair.repairOddEven(geometry, startTime);
}
if (minArea > 0) {
@@ -196,7 +214,7 @@ int main (int argc, const char * argv[]) {
void usage() {
std::cout << "=== prepair Help ===\n" << std::endl;
- std::cout << "Usage: prepair --wkt 'POLYGON(0 0, 1 0, 1 1, 0 0)'" << std::endl;
+ std::cout << "Usage: prepair --wkt 'POLYGON((0 0, 0 10, 10 0, 10 10, 0 0))'" << std::endl;
std::cout << "OR" << std::endl;
std::cout << "Usage: prepair -f infile.txt (infile.txt must contain one WKT on the 1st line)" << std::endl;
std::cout << "OR" << std::endl;
@@ -208,6 +226,7 @@ void usage() {
std::cout << "--minarea AREA Only output polygons larger than AREA" << std::endl;
std::cout << "--isr GRIDSIZE Snap round the input before repairing" << std::endl;
std::cout << "--time Benchmark the different stages of the process" << std::endl;
+ std::cout << "--shpOut Output a shapefile (out.shp) instead of a WKT" << std::endl;
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/prepair.git
More information about the Pkg-grass-devel
mailing list