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 <[email protected]> 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 [email protected] # Hugo Ledoux [email protected] # Martijn Meijers [email protected] @@ -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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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:[email protected]). +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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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 [email protected] Hugo Ledoux [email protected] Martijn Meijers [email protected] 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 _______________________________________________ Pkg-grass-devel mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
