src/simplify/TaggedLineStringSimplifier.cpp

   1 /**********************************************************************
   2  *
   3  * GEOS - Geometry Engine Open Source
   4  * http://geos.osgeo.org
   5  *
   6  * Copyright (C) 2006 Refractions Research Inc.
   7  *
   8  * This is free software; you can redistribute and/or modify it under
   9  * the terms of the GNU Lesser General Licence as published
  10  * by the Free Software Foundation.
  11  * See the COPYING file for more information.
  12  *
  13  **********************************************************************
  14  *
  15  * Last port: simplify/TaggedLineStringSimplifier.java r536 (JTS-1.12+)
  16  *
  17  **********************************************************************/
  18
  19 #include <geos/simplify/TaggedLineStringSimplifier.h>
  20 #include <geos/simplify/LineSegmentIndex.h>
  21 #include <geos/algorithm/LineIntersector.h>
  22 #include <geos/simplify/TaggedLineString.h>
  23 #include <geos/simplify/TaggedLineSegment.h>
  24 #include <geos/geom/CoordinateSequence.h>
  25 #include <geos/geom/LineString.h>
  26 //#include <geos/geom/Geometry.h> // for auto_ptr destructor
  27 //#include <geos/geom/GeometryFactory.h>
  28 //#include <geos/geom/CoordinateSequenceFactory.h>
  29
  30 #include <cassert>
  31 #include <memory>
  32
  33 #ifndef GEOS_DEBUG
  34 #define GEOS_DEBUG 0
  35 #endif
  36
  37 #ifdef GEOS_DEBUG
  38 #include <iostream>
  39 #endif
  40
  41 using namespace geos::geom;
  42 using namespace std;
  43
  44 namespace geos {
  45 namespace simplify { // geos::simplify
  46
  47 /*public*/
  48 TaggedLineStringSimplifier::TaggedLineStringSimplifier(
  49                 LineSegmentIndex* nInputIndex,
  50                 LineSegmentIndex* nOutputIndex)
  51         :
  52         inputIndex(nInputIndex),
  53         outputIndex(nOutputIndex),
  54         li(new algorithm::LineIntersector()),
  55         line(NULL),
  56         linePts(NULL),
  57         distanceTolerance(0.0)
  58 {
  59 }
  60
  61 /*public*/
  62 void
  63 TaggedLineStringSimplifier::simplify(TaggedLineString* nLine)
  64 {
  65         assert(nLine);
  66         line = nLine;
  67
  68         linePts = line->getParentCoordinates();
  69         assert(linePts);
  70
  71 #if GEOS_DEBUG
  72         std::cerr << "TaggedLineStringSimplifier[" << this << "] "
  73              << " TaggedLineString[" << line << "] "
  74              << " has " << linePts->size() << " coords in input"
  75              << std::endl;
  76 #endif
  77
  78         if ( ! linePts->size() ) return;
  79         simplifySection(0, linePts->size() - 1, 0);
  80
  81 }
  82
  83
  84 /*private*/
  85 void
  86 TaggedLineStringSimplifier::simplifySection(std::size_t i,
  87                 std::size_t j, std::size_t depth)
  88 {
  89         depth += 1;
  90
  91 #if GEOS_DEBUG
  92         std::cerr << "TaggedLineStringSimplifier[" << this << "] "
  93                   << " simplifying section " << i << "-" << j
  94                   << std::endl;
  95 #endif
  96
  97         vector<std::size_t> sectionIndex(2);
  98
  99         if((i+1) == j)
 100         {
 101
 102 #if GEOS_DEBUG
 103                 std::cerr << "single segment, no flattening"
 104                           << std::endl;
 105 #endif
 106
 107                 auto_ptr<TaggedLineSegment> newSeg(new
 108                         TaggedLineSegment(*(line->getSegment(i))));
 109
 110                 line->addToResult(newSeg);
 111                 // leave this segment in the input index, for efficiency
 112                 return;
 113         }
 114
 115         bool isValidToSimplify = true;
 116
 117         /**
 118          * Following logic ensures that there is enough points in the
 119          * output line.
 120          * If there is already more points than the minimum, there's
 121          * nothing to check.
 122          * Otherwise, if in the worst case there wouldn't be enough points,
 123          * don't flatten this segment (which avoids the worst case scenario)
 124          */
 125         if (line->getResultSize() < line->getMinimumSize())
 126         {
 127                 std::size_t worstCaseSize = depth + 1;
 128                 if (worstCaseSize < line->getMinimumSize())
 129                         isValidToSimplify = false;
 130         }
 131
 132         double distance;
 133
 134         // pass distance by ref
 135         std::size_t furthestPtIndex = findFurthestPoint(linePts, i, j, distance);
 136
 137 #if GEOS_DEBUG
 138         std::cerr << "furthest point " << furthestPtIndex
 139                   << " at distance " << distance
 140                   << std::endl;
 141 #endif
 142
 143         // flattening must be less than distanceTolerance
 144         if ( distance > distanceTolerance ) isValidToSimplify = false;
 145
 146         // test if flattened section would cause intersection
 147         LineSegment candidateSeg(linePts->getAt(i), linePts->getAt(j));
 148
 149         sectionIndex[0] = i;
 150         sectionIndex[1] = j;
 151
 152         if (hasBadIntersection(line, sectionIndex, candidateSeg))
 153                         isValidToSimplify = false;
 154
 155         if (isValidToSimplify)
 156         {
 157
 158                 auto_ptr<TaggedLineSegment> newSeg = flatten(i, j);
 159
 160 #if GEOS_DEBUG
 161                 std::cerr << "isValidToSimplify, adding seg "
 162                           << newSeg->p0 << ", " << newSeg->p1
 163                           << " to TaggedLineSegment["<<line<<"] result "
 164                           << std::endl;
 165 #endif
 166
 167                 line->addToResult(newSeg);
 168                 return;
 169         }
 170
 171         simplifySection(i, furthestPtIndex, depth);
 172         simplifySection(furthestPtIndex, j, depth);
 173
 174 }
 175
 176
 177 /*private*/
 178 auto_ptr<TaggedLineSegment>
 179 TaggedLineStringSimplifier::flatten(std::size_t start, std::size_t end)
 180 {
 181         // make a new segment for the simplified geometry
 182         const Coordinate& p0 = linePts->getAt(start);
 183         const Coordinate& p1 = linePts->getAt(end);
 184         auto_ptr<TaggedLineSegment> newSeg(new TaggedLineSegment(p0, p1));
 185         // update the indexes
 186         remove(line, start, end);
 187         outputIndex->add(newSeg.get());
 188         return newSeg;
 189 }
 190
 191 /*private*/
 192 bool
 193 TaggedLineStringSimplifier::hasBadIntersection(
 194                 const TaggedLineString* parentLine,
 195                 const vector<std::size_t>& sectionIndex,
 196                 const LineSegment& candidateSeg)
 197 {
 198         if (hasBadOutputIntersection(candidateSeg))
 199                 return true;
 200
 201         if (hasBadInputIntersection(parentLine, sectionIndex, candidateSeg))
 202                 return true;
 203
 204         return false;
 205 }
 206
 207 /*private*/
 208 bool
 209 TaggedLineStringSimplifier::hasBadOutputIntersection(
 210                 const LineSegment& candidateSeg)
 211 {
 212         auto_ptr< vector<LineSegment*> > querySegs =
 213                 outputIndex->query(&candidateSeg);
 214
 215         for (vector<LineSegment*>::iterator
 216                         it = querySegs->begin(), iEnd = querySegs->end();
 217                         it != iEnd;
 218                         ++it)
 219         {
 220                 LineSegment* querySeg = *it;
 221                 assert(querySeg);
 222                 if (hasInteriorIntersection(*querySeg, candidateSeg))
 223                 {
 224                         return true;
 225                 }
 226         }
 227
 228         return false;
 229 }
 230
 231 /*private*/
 232 bool
 233 TaggedLineStringSimplifier::hasInteriorIntersection(
 234                         const LineSegment& seg0,
 235                         const LineSegment& seg1) const
 236 {
 237         li->computeIntersection(seg0.p0, seg0.p1, seg1.p0, seg1.p1);
 238         return li->isInteriorIntersection();
 239 }
 240
 241 /*private*/
 242 bool
 243 TaggedLineStringSimplifier::hasBadInputIntersection(
 244                 const TaggedLineString* parentLine,
 245                 const vector<std::size_t>& sectionIndex,
 246                 const LineSegment& candidateSeg)
 247 {
 248         auto_ptr< vector<LineSegment*> > querySegs =
 249                 inputIndex->query(&candidateSeg);
 250
 251         for (vector<LineSegment*>::iterator
 252                         it = querySegs->begin(), iEnd = querySegs->end();
 253                         it != iEnd;
 254                         ++it)
 255         {
 256                 assert(*it);
 257                 assert(dynamic_cast<TaggedLineSegment*>(*it));
 258                 TaggedLineSegment* querySeg =
 259                         static_cast<TaggedLineSegment*>(*it);
 260
 261                 if (hasInteriorIntersection(*querySeg, candidateSeg))
 262                 {
 263
 264                         if ( isInLineSection(parentLine,
 265                                         sectionIndex, querySeg) )
 266                         {
 267                                 continue;
 268                         }
 269
 270                         return true;
 271                 }
 272         }
 273
 274         return false;
 275 }
 276
 277 /*static private*/
 278 bool
 279 TaggedLineStringSimplifier::isInLineSection(
 280                 const TaggedLineString* line,
 281                 const vector<std::size_t>& sectionIndex,
 282                 const TaggedLineSegment* seg)
 283 {
 284         // not in this line
 285         if (seg->getParent() != line->getParent())
 286                 return false;
 287
 288         std::size_t segIndex = seg->getIndex();
 289         if (segIndex >= sectionIndex[0] && segIndex < sectionIndex[1])
 290                 return true;
 291
 292         return false;
 293 }
 294
 295 /*private*/
 296 void
 297 TaggedLineStringSimplifier::remove(const TaggedLineString* line,
 298                 std::size_t start,
 299                 std::size_t end)
 300 {
 301         assert(end <= line->getSegments().size() );
 302         assert(start < end); // I'm not sure this should always be true
 303
 304         for (std::size_t i = start; i < end; i++)
 305         {
 306                 const TaggedLineSegment* seg = line->getSegment(i);
 307                 inputIndex->remove(seg);
 308         }
 309 }
 310
 311 /*private static*/
 312 std::size_t
 313 TaggedLineStringSimplifier::findFurthestPoint(
 314                 const geom::CoordinateSequence* pts,
 315                 std::size_t i, std::size_t j,
 316                 double& maxDistance)
 317 {
 318         LineSegment seg(pts->getAt(i), pts->getAt(j));
 319 #if GEOS_DEBUG
 320         std::cerr << __FUNCTION__ << "segment " << seg
 321                   << std::endl;
 322 #endif
 323         double maxDist = -1.0;
 324         std::size_t maxIndex = i;
 325         for (std::size_t k = i + 1; k < j; k++)
 326         {
 327                 const Coordinate& midPt = pts->getAt(k);
 328                 double distance = seg.distance(midPt);
 329 #if GEOS_DEBUG
 330                 std::cerr << "dist to " << midPt
 331                           << ": " << distance
 332                           << std::endl;
 333 #endif
 334                 if (distance > maxDist) {
 335 #if GEOS_DEBUG
 336                         std::cerr << "this is max"
 337                                   << std::endl;
 338 #endif
 339                         maxDist = distance;
 340                         maxIndex = k;
 341                 }
 342         }
 343         maxDistance = maxDist;
 344         return maxIndex;
 345
 346 }
 347
 348 } // namespace geos::simplify
 349 } // namespace geos