[geos-commits] r3800 - in trunk: include/geos/operation/overlay/snap src/operation/overlay/snap tests/unit/capi
svn_geos at osgeo.org
svn_geos at osgeo.org
Thu May 30 13:16:22 PDT 2013
Author: strk
Date: 2013-05-30 13:16:22 -0700 (Thu, 30 May 2013)
New Revision: 3800
Modified:
trunk/include/geos/operation/overlay/snap/LineStringSnapper.h
trunk/src/operation/overlay/snap/LineStringSnapper.cpp
trunk/tests/unit/capi/GEOSSnapTest.cpp
Log:
Improve snap algorithm reducing likelyhood of invalid output
- Snap input vertices to closest snap point (#629)
- Do not snap segments to external points (#501)
- Never snap multiple vertices to the same snap point
Work funded by Tuscany Region - SITA. Contract "Support to the use
of GFOSS (Geographic Free and Open Source Software) Desktop tools"
(CIG Z3B06FA6D7).
Modified: trunk/include/geos/operation/overlay/snap/LineStringSnapper.h
===================================================================
--- trunk/include/geos/operation/overlay/snap/LineStringSnapper.h 2013-05-23 17:21:38 UTC (rev 3799)
+++ trunk/include/geos/operation/overlay/snap/LineStringSnapper.h 2013-05-30 20:16:22 UTC (rev 3800)
@@ -146,6 +146,11 @@
geom::CoordinateList::iterator from,
geom::CoordinateList::iterator too_far);
+ geom::CoordinateList::iterator findVertexToSnap(
+ const geom::Coordinate& snapPt,
+ geom::CoordinateList::iterator from,
+ geom::CoordinateList::iterator too_far);
+
// Declare type as noncopyable
LineStringSnapper(const LineStringSnapper& other);
LineStringSnapper& operator=(const LineStringSnapper& rhs);
Modified: trunk/src/operation/overlay/snap/LineStringSnapper.cpp
===================================================================
--- trunk/src/operation/overlay/snap/LineStringSnapper.cpp 2013-05-23 17:21:38 UTC (rev 3799)
+++ trunk/src/operation/overlay/snap/LineStringSnapper.cpp 2013-05-30 20:16:22 UTC (rev 3800)
@@ -15,6 +15,9 @@
*
* Last port: operation/overlay/snap/LineStringSnapper.java r320 (JTS-1.12)
*
+ * NOTE: algorithm changed to improve output quality by reducing
+ * probability of self-intersections
+ *
**********************************************************************/
#include <geos/operation/overlay/snap/LineStringSnapper.h>
@@ -33,6 +36,7 @@
#if GEOS_DEBUG
#include <iostream>
+#include <iomanip>
using std::cerr;
using std::endl;
#endif
@@ -58,6 +62,53 @@
}
/*private*/
+CoordinateList::iterator
+LineStringSnapper::findVertexToSnap(
+ const Coordinate& snapPt,
+ CoordinateList::iterator from,
+ CoordinateList::iterator too_far)
+{
+ double minDist = snapTolerance; // make sure the first closer then
+ // snapTolerance is accepted
+ CoordinateList::iterator match=too_far;
+
+ for ( ; from != too_far; ++from)
+ {
+ Coordinate& c0 = *from;
+
+#if GEOS_DEBUG
+cerr << " Checking vertex " << c0 << endl;
+#endif
+
+ double dist = c0.distance(snapPt);
+ if ( dist >= minDist) {
+#if GEOS_DEBUG
+cerr << " snap point distance " << dist
+ << " not smaller than tolerance "
+ << snapTolerance << " or previous closest "
+ << minDist << endl;
+#endif
+ continue;
+ }
+
+#if GEOS_DEBUG
+ cerr << " snap point distance " << dist << " within tolerance "
+ << snapTolerance << " and closer than previous candidate "
+ << minDist << endl;
+#endif
+
+ if ( dist == 0.0 ) return from; // can't find any closer
+
+ match = from;
+ minDist = dist;
+
+ }
+
+ return match;
+}
+
+/*private*/
+#if 0 // changed to start looping from snap points
void
LineStringSnapper::snapVertices(geom::CoordinateList& srcCoords,
const geom::Coordinate::ConstVect& snapPts)
@@ -65,7 +116,7 @@
if ( srcCoords.empty() ) return;
#if GEOS_DEBUG
-cerr << "Snapping vertices of: " << srcCoords << endl;
+cerr << "Snapping vertices of: " << std::setprecision(15) << srcCoords << endl;
#endif
using geom::CoordinateList;
@@ -97,6 +148,18 @@
assert(*found);
const Coordinate& snapPt = *(*found);
+
+#if 1
+ // Check if closest vertex for found snap is ourselves ?
+ CoordinateList::iterator vertpos = findVertexToSnap(snapPt, srcCoords.begin(), srcCoords.end());
+ if ( vertpos != it ) {
+#if GEOS_DEBUG
+cerr << " Snap found but it has a vertex point closer than this one" << endl;
+ continue;
+#endif
+ }
+#endif
+
#if GEOS_DEBUG
cerr << " found snap point " << snapPt << endl;
@@ -116,15 +179,88 @@
}
}
}
+#endif
+#if 1
/*private*/
+void
+LineStringSnapper::snapVertices(geom::CoordinateList& srcCoords,
+ const geom::Coordinate::ConstVect& snapPts)
+{
+ // nothing to do if there are no source coords..
+ if ( srcCoords.empty() ) return;
+
+#if GEOS_DEBUG
+cerr << "Snapping vertices of: " << srcCoords << endl;
+#endif
+
+ for ( Coordinate::ConstVect::const_iterator
+ it=snapPts.begin(), end=snapPts.end();
+ it != end;
+ ++it)
+ {
+ assert(*it);
+ const Coordinate& snapPt = *(*it);
+
+#if GEOS_DEBUG
+cerr << "Checking for a vertex to snap to snapPt " << snapPt << endl;
+#endif
+
+ CoordinateList::iterator too_far = srcCoords.end();
+ if ( isClosed ) --too_far;
+ CoordinateList::iterator vertpos =
+ findVertexToSnap(snapPt, srcCoords.begin(), too_far);
+ if ( vertpos == too_far)
+ {
+#if GEOS_DEBUG
+cerr << " No vertex to snap" << endl;
+#endif
+ continue;
+ }
+
+#if 0
+ // Check if closest snap for found vertex is ourselves ?
+ Coordinate::ConstVect::const_iterator snappos =
+ findSnapForVertex(*vertpos, snapPts);
+ if ( snappos != it ) {
+#if GEOS_DEBUG
+cerr << " Vertex found but it has a snap point closer than this one" << endl;
+ continue;
+#endif
+ }
+#endif
+
+#if GEOS_DEBUG
+cerr << " Vertex to be snapped found, snapping" << endl;
+#endif
+ *vertpos = snapPt;
+
+ // keep final closing point in synch (rings only)
+ if (vertpos == srcCoords.begin() && isClosed)
+ {
+ vertpos = srcCoords.end(); --vertpos;
+ *vertpos = snapPt;
+ }
+
+ }
+
+#if GEOS_DEBUG
+cerr << " After vertex snapping, srcCoors are: " << srcCoords << endl;
+#endif
+
+}
+#endif
+
+/*private*/
Coordinate::ConstVect::const_iterator
LineStringSnapper::findSnapForVertex(const Coordinate& pt,
const Coordinate::ConstVect& snapPts)
{
+ Coordinate::ConstVect::const_iterator end = snapPts.end();
+ Coordinate::ConstVect::const_iterator candidate = end;
+ double minDist;
// TODO: use std::find_if
- Coordinate::ConstVect::const_iterator end=snapPts.end();
for ( Coordinate::ConstVect::const_iterator
it=snapPts.begin();
it != end;
@@ -133,14 +269,12 @@
assert(*it);
const Coordinate& snapPt = *(*it);
- // shouldn't we look for *all* segments to be snapped rather then a single one?
if ( snapPt.equals2D(pt) )
{
#if GEOS_DEBUG
cerr << " points are equal, returning not-found " << endl;
#endif
return end;
- //continue;
}
double dist = snapPt.distance(pt);
@@ -150,18 +284,20 @@
if ( dist < snapTolerance )
{
-#if GEOS_DEBUG
-cerr << " snap point within tolerance, returning iterator to it" << endl;
-#endif
- return it;
+ if ( candidate == end || dist < minDist ) {
+ minDist = dist;
+ candidate = it;
+ }
}
}
#if GEOS_DEBUG
+ if ( candidate == end ) {
cerr << " no snap point within distance, returning not-found" << endl;
+ }
#endif
- return end;
+ return candidate;
}
@@ -190,8 +326,6 @@
cerr << "Checking for a segment to snap to snapPt " << snapPt << endl;
#endif
- // shouldn't we look for *all* segments to be snapped
- // rather then a single one?
CoordinateList::iterator too_far = srcCoords.end(); --too_far;
CoordinateList::iterator segpos =
findSegmentToSnap(snapPt, srcCoords.begin(), too_far);
@@ -202,12 +336,113 @@
#endif
continue;
}
+
+ /* Check if the snap point falls outside of the segment */
+ // If the snap point is outside, this means that an endpoint
+ // was not snap where it should have been
+ // so what we should do is re-snap the endpoint to this
+ // snapPt and then snap the closest between this and
+ // previous (for pf < 0.0) or next (for pf > 1.0) segment
+ // to the old endpoint.
+ // --strk May 2013
+ //
+ // TODO: simplify this code, make more readable
+ //
+ CoordinateList::iterator to = segpos; ++to;
+ LineSegment seg(*segpos, *to);
+ double pf = seg.projectionFactor(snapPt);
+ if ( pf > 1.0 ) {
#if GEOS_DEBUG
-cerr << " Segment to be snapped found, inserting point" << endl;
+ cerr << " Segment to be snapped is closer on his end point" << endl;
#endif
- // insert must happen one-past first point (before next point)
- ++segpos;
- srcCoords.insert(segpos, snapPt);
+ Coordinate newSnapPt = seg.p1;
+ *to = seg.p1 = snapPt;
+ // now snap from-to (segpos) or to-next (segpos++) to newSnapPt
+ if ( to == too_far ) {
+ if ( isClosed ) {
+#if GEOS_DEBUG
+ cerr << " His end point is the last one, but is closed " << endl;
+#endif
+ *(srcCoords.begin()) = snapPt; // sync to start point
+ to = srcCoords.begin();
+ } else {
+#if GEOS_DEBUG
+ cerr << " His end point is the last one, inserting " << newSnapPt << " before it" << endl;
+#endif
+ srcCoords.insert(to, newSnapPt);
+ continue;
+ }
+ }
+ ++to;
+ LineSegment nextSeg(seg.p1, *to);
+ if ( nextSeg.distance(newSnapPt) < seg.distance(newSnapPt) ) {
+#if GEOS_DEBUG
+ cerr << " Next segment closer, inserting " << newSnapPt << " into " << nextSeg << endl;
+#endif
+ // insert into next segment
+ srcCoords.insert(to, newSnapPt);
+ } else {
+#if GEOS_DEBUG
+ cerr << " This segment closer, inserting " << newSnapPt << " into " << seg << endl;
+#endif
+ // insert must happen one-past first point (before next point)
+ ++segpos;
+ srcCoords.insert(segpos, newSnapPt);
+ }
+ }
+ else if ( pf < 0.0 ) {
+#if GEOS_DEBUG
+ cerr << " Segment to be snapped is closer on his start point" << endl;
+#endif
+ Coordinate newSnapPt = seg.p0;
+ *segpos = seg.p0 = snapPt;
+ // now snap prev-from (--segpos) or from-to (segpos) to newSnapPt
+ if ( segpos == srcCoords.begin() ) {
+ if ( isClosed ) {
+#if GEOS_DEBUG
+ cerr << " His start point is the first one, but is closed " << endl;
+#endif
+ segpos = srcCoords.end(); --segpos;
+ *segpos = snapPt; // sync to end point
+ } else {
+#if GEOS_DEBUG
+ cerr << " His start point is the first one, inserting " << newSnapPt << " before it" << endl;
+#endif
+ ++segpos;
+ srcCoords.insert(segpos, newSnapPt);
+ continue;
+ }
+ }
+
+#if GEOS_DEBUG
+cerr << " Before seg-snapping, srcCoors are: " << srcCoords << endl;
+#endif
+
+ --segpos;
+ LineSegment prevSeg(*segpos, seg.p0);
+ if ( prevSeg.distance(newSnapPt) < seg.distance(newSnapPt) ) {
+#if GEOS_DEBUG
+ cerr << " Prev segment closer, inserting " << newSnapPt << " into " << prevSeg << endl;
+#endif
+ // insert into prev segment
+ srcCoords.insert(segpos, newSnapPt);
+ } else {
+#if GEOS_DEBUG
+ cerr << " This segment closer, inserting " << newSnapPt << " into " << seg << endl;
+#endif
+ // insert must happen one-past first point (before next point)
+ srcCoords.insert(to, newSnapPt);
+ }
+ }
+ else {
+ assert(pf != 0.0);
+#if GEOS_DEBUG
+cerr << " Segment to be snapped found, projection factor is " << pf << ", inserting point" << endl;
+#endif
+ // insert must happen one-past first point (before next point)
+ ++segpos;
+ srcCoords.insert(segpos, snapPt);
+ }
}
#if GEOS_DEBUG
@@ -225,8 +460,8 @@
CoordinateList::iterator too_far)
{
LineSegment seg;
- double minDist = snapTolerance+1; // make sure the first closer then
- // snapTolerance is accepted
+ double minDist = snapTolerance; // make sure the first closer then
+ // snapTolerance is accepted
CoordinateList::iterator match=too_far;
// TODO: use std::find_if
@@ -267,30 +502,27 @@
}
double dist = seg.distance(snapPt);
- if ( dist < snapTolerance ) {
- if ( dist < minDist ) {
+ if ( dist >= minDist) {
#if GEOS_DEBUG
-cerr << " snap point distance " << dist << " within tolerance "
- << snapTolerance << " and closer than previous candidate " << minDist
- << endl;
+cerr << " snap point distance " << dist
+ << " not smaller than tolerance "
+ << snapTolerance << " or previous closest "
+ << minDist << endl;
#endif
- match = from;
- minDist = dist;
- }
-#if GEOS_DEBUG
- else {
-cerr << " snap point distance " << dist << " within tolerance "
- << snapTolerance << " but not closer than previous candidate " << minDist
- << endl;
- }
-#endif
+ continue;
}
+
#if GEOS_DEBUG
- else {
-cerr << " snap point distance " << dist << " bigger than tolerance "
- << snapTolerance << endl;
- }
+ cerr << " snap point distance " << dist << " within tolerance "
+ << snapTolerance << " and closer than previous candidate "
+ << minDist << endl;
#endif
+
+ if ( dist == 0.0 ) return from; // can't find any closer
+
+ match = from;
+ minDist = dist;
+
}
return match;
Modified: trunk/tests/unit/capi/GEOSSnapTest.cpp
===================================================================
--- trunk/tests/unit/capi/GEOSSnapTest.cpp 2013-05-23 17:21:38 UTC (rev 3799)
+++ trunk/tests/unit/capi/GEOSSnapTest.cpp 2013-05-30 20:16:22 UTC (rev 3800)
@@ -133,5 +133,88 @@
ensure_equals(out, "LINESTRING (0 0, 9 0)");
}
+ /// See http://trac.osgeo.org/geos/ticket/501
+ template<>
+ template<>
+ void object::test<5>()
+ {
+ geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0)");
+ geom2_ = GEOSGeomFromWKT("LINESTRING(0 0, 9 0, 10 0, 11 0)");
+ geom3_ = GEOSSnap(geom1_, geom2_, 2);
+
+ char* wkt_c = GEOSWKTWriter_write(w_, geom3_);
+ std::string out(wkt_c);
+ free(wkt_c);
+
+ //ensure_equals(out, "LINESTRING (0 0, 9 0, 10 0)");
+ ensure_equals(out, "LINESTRING (0 0, 9 0, 10 0, 11 0)");
+ }
+
+ /// Test snapping of equidistant segments to outlyers snap point
+ template<>
+ template<>
+ void object::test<6>()
+ {
+ geom1_ = GEOSGeomFromWKT("LINESTRING(0 3,4 1,0 1)");
+ geom2_ = GEOSGeomFromWKT("MULTIPOINT(5 0,4 1)");
+ geom3_ = GEOSSnap(geom1_, geom2_, 2);
+
+ char* wkt_c = GEOSWKTWriter_write(w_, geom3_);
+ std::string out(wkt_c);
+ free(wkt_c);
+
+ ensure_equals(out, "LINESTRING (0 3, 4 1, 5 0, 0 1)");
+ //ensure_equals(out, "LINESTRING (0 3, 4 1, 0 1)");
+ }
+
+ /// Test snapping of equidistant segments to outlyers snap point
+ /// Same as the above but with the snap points order reversed
+ template<>
+ template<>
+ void object::test<7>()
+ {
+ geom1_ = GEOSGeomFromWKT("LINESTRING(0 3,4 1,0 1)");
+ geom2_ = GEOSGeomFromWKT("MULTIPOINT(4 1,5 0)");
+ geom3_ = GEOSSnap(geom1_, geom2_, 2);
+
+ char* wkt_c = GEOSWKTWriter_write(w_, geom3_);
+ std::string out(wkt_c);
+ free(wkt_c);
+
+ ensure_equals(out, "LINESTRING (0 3, 4 1, 5 0, 0 1)");
+ //ensure_equals(out, "LINESTRING (0 3, 4 1, 0 1)");
+ }
+
+ /// Test snapping of closed ring to outlyers snap point
+ template<>
+ template<>
+ void object::test<8>()
+ {
+ geom1_ = GEOSGeomFromWKT("LINESTRING(0 0,10 0,10 10,0 10,0 0)");
+ geom2_ = GEOSGeomFromWKT("MULTIPOINT(0 0,-1 0)");
+ geom3_ = GEOSSnap(geom1_, geom2_, 3);
+
+ char* wkt_c = GEOSWKTWriter_write(w_, geom3_);
+ std::string out(wkt_c);
+ free(wkt_c);
+
+ ensure_equals(out, "LINESTRING (-1 0, 0 0, 10 0, 10 10, 0 10, -1 0)");
+ }
+
+ template<>
+ template<>
+ void object::test<9>()
+ {
+ geom1_ = GEOSGeomFromWKT("LINESTRING(0 2,5 2,9 2,5 0)");
+ geom2_ = GEOSGeomFromWKT("POINT(5 0)");
+ geom3_ = GEOSSnap(geom1_, geom2_, 3);
+
+ char* wkt_c = GEOSWKTWriter_write(w_, geom3_);
+ std::string out(wkt_c);
+ free(wkt_c);
+
+ ensure_equals(out, "LINESTRING (0 2, 5 2, 9 2, 5 0)");
+ }
+
} // namespace tut
More information about the geos-commits
mailing list