src/operation/overlay/snap/LineStringSnapper.cpp

   1 /**********************************************************************
   2  * $Id$
   3  *
   4  * GEOS - Geometry Engine Open Source
   5  * http://geos.refractions.net
   6  *
   7  * Copyright (C) 2009-2010  Sandro Santilli <strk@keybit.net>
   8  * Copyright (C) 2006 Refractions Research Inc.
   9  *
  10  * This is free software; you can redistribute and/or modify it under
  11  * the terms of the GNU Lesser General Public Licence as published
  12  * by the Free Software Foundation.
  13  * See the COPYING file for more information.
  14  *
  15  ***********************************************************************
  16  *
  17  * Last port: operation/overlay/snap/LineStringSnapper.java r320 (JTS-1.12)
  18  *
  19  **********************************************************************/
  20
  21 #include <geos/operation/overlay/snap/LineStringSnapper.h>
  22 #include <geos/geom/CoordinateSequence.h>
  23 #include <geos/geom/Coordinate.h>
  24 #include <geos/geom/CoordinateList.h>
  25 #include <geos/util/UniqueCoordinateArrayFilter.h>
  26 #include <geos/geom/LineSegment.h>
  27
  28 #include <vector>
  29 #include <memory>
  30
  31 #ifndef GEOS_DEBUG
  32 #define GEOS_DEBUG 0
  33 #endif
  34
  35 #if GEOS_DEBUG
  36 #include <iostream>
  37 using std::cerr;
  38 using std::endl;
  39 #endif
  40
  41 //using namespace std;
  42 using namespace geos::geom;
  43
  44 namespace geos {
  45 namespace operation { // geos.operation
  46 namespace overlay { // geos.operation.overlay
  47 namespace snap { // geos.operation.overlay.snap
  48
  49 /*public*/
  50 std::auto_ptr<Coordinate::Vect>
  51 LineStringSnapper::snapTo(const geom::Coordinate::ConstVect& snapPts)
  52 {
  53         geom::CoordinateList coordList(srcPts);
  54
  55         snapVertices(coordList, snapPts);
  56         snapSegments(coordList, snapPts);
  57
  58         return coordList.toCoordinateArray();
  59 }
  60
  61 /*private*/
  62 void
  63 LineStringSnapper::snapVertices(geom::CoordinateList& srcCoords,
  64                         const geom::Coordinate::ConstVect& snapPts)
  65 {
  66   if ( srcCoords.empty() ) return;
  67
  68 #if GEOS_DEBUG
  69 cerr << "Snapping vertices of: " << srcCoords << endl;
  70 #endif
  71
  72         using geom::CoordinateList;
  73
  74         geom::Coordinate::ConstVect::const_iterator not_found = snapPts.end();
  75
  76         // try snapping vertices
  77         // if src is a ring then don't snap final vertex
  78         CoordinateList::iterator it = srcCoords.begin();
  79         CoordinateList::iterator end = srcCoords.end();
  80         CoordinateList::iterator last = end; --last;
  81   if ( isClosed ) --end;
  82         for ( ; it != end; ++it )
  83         {
  84                 Coordinate& srcPt = *it;
  85
  86 #if GEOS_DEBUG
  87 cerr << "Checking for a snap for source coordinate " << srcPt << endl;
  88 #endif
  89
  90                 geom::Coordinate::ConstVect::const_iterator found = findSnapForVertex(srcPt, snapPts);
  91                 if ( found == not_found )
  92                 {       // no snaps found (or no need to snap)
  93 #if GEOS_DEBUG
  94 cerr << " no snap found" << endl;
  95 #endif
  96                         continue;
  97                 }
  98
  99                 assert(*found);
 100                 const Coordinate& snapPt = *(*found);
 101
 102 #if GEOS_DEBUG
 103 cerr << " found snap point " << snapPt << endl;
 104 #endif
 105
 106                 // update src with snap pt
 107                 *it = snapPt;
 108
 109 #if GEOS_DEBUG
 110 cerr << " source point became " << srcPt << endl;
 111 #endif
 112
 113                 // keep final closing point in synch (rings only)
 114                 if (it == srcCoords.begin() && isClosed)
 115                 {
 116                         *last = snapPt;
 117                 }
 118         }
 119 }
 120
 121 /*private*/
 122 Coordinate::ConstVect::const_iterator
 123 LineStringSnapper::findSnapForVertex(const Coordinate& pt,
 124                         const Coordinate::ConstVect& snapPts)
 125 {
 126
 127         // TODO: use std::find_if
 128         Coordinate::ConstVect::const_iterator end=snapPts.end();
 129         for ( Coordinate::ConstVect::const_iterator
 130                         it=snapPts.begin();
 131                         it != end;
 132                         ++it)
 133         {
 134                 assert(*it);
 135                 const Coordinate& snapPt = *(*it);
 136
 137                 // shouldn't we look for *all* segments to be snapped rather then a single one?
 138                 if ( snapPt.equals2D(pt) )
 139                 {
 140 #if GEOS_DEBUG
 141 cerr << " points are equal, returning not-found " << endl;
 142 #endif
 143                         return end;
 144                         //continue;
 145                 }
 146
 147                 double dist = snapPt.distance(pt);
 148 #if GEOS_DEBUG
 149 cerr << " distance from snap point " << snapPt << ": " << dist << endl;
 150 #endif
 151
 152                 if ( dist < snapTolerance )
 153                 {
 154 #if GEOS_DEBUG
 155 cerr << " snap point within tolerance, returning iterator to it" << endl;
 156 #endif
 157                         return it;
 158                 }
 159         }
 160
 161 #if GEOS_DEBUG
 162 cerr << " no snap point within distance, returning not-found" << endl;
 163 #endif
 164
 165         return end;
 166 }
 167
 168
 169 /*private*/
 170 void
 171 LineStringSnapper::snapSegments(geom::CoordinateList& srcCoords,
 172                         const geom::Coordinate::ConstVect& snapPts)
 173 {
 174
 175   // nothing to do if there are no source coords..
 176   if ( srcCoords.empty() ) return;
 177
 178 #if GEOS_DEBUG
 179 cerr << "Snapping segments of: " << srcCoords << endl;
 180 #endif
 181
 182         for ( Coordinate::ConstVect::const_iterator
 183                         it=snapPts.begin(), end=snapPts.end();
 184                         it != end;
 185                         ++it)
 186         {
 187                 assert(*it);
 188                 const Coordinate& snapPt = *(*it);
 189
 190 #if GEOS_DEBUG
 191 cerr << "Checking for a segment to snap to snapPt " << snapPt << endl;
 192 #endif
 193
 194                 // shouldn't we look for *all* segments to be snapped
 195                 // rather then a single one?
 196                 CoordinateList::iterator too_far = srcCoords.end(); --too_far;
 197                 CoordinateList::iterator segpos =
 198                         findSegmentToSnap(snapPt, srcCoords.begin(), too_far);
 199                 if ( segpos == too_far)
 200                 {
 201 #if GEOS_DEBUG
 202 cerr << " No segment to snap" << endl;
 203 #endif
 204                         continue;
 205                 }
 206 #if GEOS_DEBUG
 207 cerr << " Segment to be snapped found, inserting point" << endl;
 208 #endif
 209                 // insert must happen one-past first point (before next point)
 210                 ++segpos;
 211                 srcCoords.insert(segpos, snapPt);
 212         }
 213
 214 #if GEOS_DEBUG
 215 cerr << " After segment snapping, srcCoors are: " << srcCoords << endl;
 216 #endif
 217
 218 }
 219
 220 /*private*/
 221 /* NOTE: this is called findSegmentIndexToSnap in JTS */
 222 CoordinateList::iterator
 223 LineStringSnapper::findSegmentToSnap(
 224                         const Coordinate& snapPt,
 225                         CoordinateList::iterator from,
 226                         CoordinateList::iterator too_far)
 227 {
 228         LineSegment seg;
 229         double minDist = snapTolerance+1; // make sure the first closer then
 230                                           // snapTolerance is accepted
 231         CoordinateList::iterator match=too_far;
 232
 233         // TODO: use std::find_if
 234         for ( ; from != too_far; ++from)
 235         {
 236                 seg.p0 = *from;
 237                 CoordinateList::iterator to = from;
 238                 ++to;
 239                 seg.p1 = *to;
 240
 241 #if GEOS_DEBUG
 242 cerr << " Checking segment " << seg << endl;
 243 #endif
 244
 245                 /**                                                                              * Check if the snap pt is equal to one of
 246                  * the segment endpoints.
 247                  *
 248                  * If the snap pt is already in the src list,
 249                  * don't snap at all (unless allowSnappingToSourceVertices
 250                  * is set to true)
 251                  */
 252                 if ( seg.p0.equals2D(snapPt) || seg.p1.equals2D(snapPt) )
 253                 {
 254
 255                         if (allowSnappingToSourceVertices) {
 256 #if GEOS_DEBUG
 257 cerr << "   snap point matches a segment endpoint, checking next segment"
 258      << endl;
 259 #endif
 260                                 continue;
 261                         } else {
 262 #if GEOS_DEBUG
 263 cerr << "   snap point matches a segment endpoint, giving up seek" << endl;
 264 #endif
 265                                 return too_far;
 266                         }
 267                 }
 268
 269                 double dist = seg.distance(snapPt);
 270                 if ( dist < snapTolerance ) {
 271       if ( dist < minDist ) {
 272 #if GEOS_DEBUG
 273 cerr << "   snap point distance " << dist << " within tolerance "
 274      << snapTolerance << " and closer than previous candidate " << minDist
 275      << endl;
 276 #endif
 277         match = from;
 278         minDist = dist;
 279       }
 280 #if GEOS_DEBUG
 281       else {
 282 cerr << "   snap point distance " << dist << " within tolerance "
 283      << snapTolerance << " but not closer than previous candidate " << minDist
 284      << endl;
 285       }
 286 #endif
 287     }
 288 #if GEOS_DEBUG
 289     else {
 290 cerr << "   snap point distance " << dist << " bigger than tolerance "
 291      << snapTolerance << endl;
 292     }
 293 #endif
 294         }
 295
 296         return match;
 297 }
 298
 299 } // namespace geos.operation.snap
 300 } // namespace geos.operation.overlay
 301 } // namespace geos.operation
 302 } // namespace geos
 303