src/operation/sharedpaths/SharedPathsOp.cpp

   1 /**********************************************************************
   2  *
   3  * GEOS - Geometry Engine Open Source
   4  * http://geos.osgeo.org
   5  *
   6  * Copyright (C) 2010      Sandro Santilli <strk@keybit.net>
   7  *
   8  * This is free software; you can redistribute and/or modify it under
   9  * the terms of the GNU Lesser General Public Licence as published
  10  * by the Free Software Foundation.
  11  * See the COPYING file for more information.
  12  *
  13  **********************************************************************
  14  *
  15  * Last port: original work
  16  *
  17  * Developed by Sandro Santilli (strk@keybit.net)
  18  * for Faunalia (http://www.faunalia.it)
  19  * with funding from Regione Toscana - Settore SISTEMA INFORMATIVO
  20  * TERRITORIALE ED AMBIENTALE - for the project: "Sviluppo strumenti
  21  * software per il trattamento di dati geografici basati su QuantumGIS
  22  * e Postgis (CIG 0494241492)"
  23  *
  24  **********************************************************************/
  25
  26 #include <geos/operation/sharedpaths/SharedPathsOp.h>
  27 #include <geos/geom/GeometryFactory.h>
  28 #include <geos/geom/Geometry.h>
  29 #include <geos/geom/LineString.h>
  30 #include <geos/geom/MultiLineString.h>
  31 #include <geos/linearref/LinearLocation.h>
  32 #include <geos/linearref/LocationIndexOfPoint.h>
  33 #include <geos/operation/overlay/OverlayOp.h>
  34 #include <geos/util/IllegalArgumentException.h>
  35
  36 using namespace geos::geom;
  37
  38 namespace geos {
  39 namespace operation { // geos.operation
  40 namespace sharedpaths { // geos.operation.sharedpaths
  41
  42 /* public static */
  43 void
  44 SharedPathsOp::sharedPathsOp(const Geometry& g1, const Geometry& g2,
  45     PathList& sameDirection,
  46     PathList& oppositeDirection)
  47 {
  48   SharedPathsOp sp(g1, g2);
  49   sp.getSharedPaths(sameDirection, oppositeDirection);
  50 }
  51
  52 /* public */
  53 SharedPathsOp::SharedPathsOp(
  54     const geom::Geometry& g1, const geom::Geometry& g2)
  55   :
  56   _g1(g1),
  57   _g2(g2),
  58   _gf(*g1.getFactory())
  59 {
  60   checkLinealInput(_g1);
  61   checkLinealInput(_g2);
  62 }
  63
  64 /* private */
  65 void
  66 SharedPathsOp::checkLinealInput(const geom::Geometry& g)
  67 {
  68   if ( ! dynamic_cast<const LineString*>(&g) &&
  69        ! dynamic_cast<const MultiLineString*>(&g) )
  70   {
  71     throw util::IllegalArgumentException("Geometry is not lineal");
  72   }
  73 }
  74
  75 /* public */
  76 void
  77 SharedPathsOp::getSharedPaths(PathList& forwDir, PathList& backDir)
  78 {
  79   PathList paths;
  80   findLinearIntersections(paths);
  81   for (size_t i=0, n=paths.size(); i<n; ++i)
  82   {
  83     LineString* path = paths[i];
  84     if ( isSameDirection(*path) ) forwDir.push_back(path);
  85     else backDir.push_back(path);
  86   }
  87 }
  88
  89 /* static private */
  90 void
  91 SharedPathsOp::clearEdges(PathList& edges)
  92 {
  93   for (PathList::const_iterator
  94         i=edges.begin(), e=edges.end();
  95         i!=e; ++i)
  96   {
  97     delete *i;
  98   }
  99   edges.clear();
 100 }
 101
 102 /* private */
 103 void
 104 SharedPathsOp::findLinearIntersections(PathList& to)
 105 {
 106   using geos::operation::overlay::OverlayOp;
 107
 108   // TODO: optionally use the tolerance,
 109   //       snapping _g2 over _g1 ?
 110
 111   std::auto_ptr<Geometry> full ( OverlayOp::overlayOp(
 112     &_g1, &_g2, OverlayOp::opINTERSECTION) );
 113
 114   // NOTE: intersection of equal lines yelds splitted lines,
 115   //       should we sew them back ?
 116
 117   for (size_t i=0, n=full->getNumGeometries(); i<n; ++i)
 118   {
 119     const Geometry* sub = full->getGeometryN(i);
 120     const LineString* path = dynamic_cast<const LineString*>(sub);
 121     if ( path ) {
 122       // NOTE: we're making a copy here, wouldn't be needed
 123       //       for a simple predicate
 124       to.push_back(_gf.createLineString(*path).release());
 125     }
 126   }
 127 }
 128
 129 /* private */
 130 bool
 131 SharedPathsOp::isForward(const geom::LineString& edge,
 132                        const geom::Geometry& geom)
 133 {
 134   using namespace geos::linearref;
 135
 136   /*
 137      ALGO:
 138       1. find first point of edge on geom (linearref)
 139       2. find second point of edge on geom (linearref)
 140       3. if first < second, we're forward
 141
 142      PRECONDITIONS:
 143       1. edge has at least 2 points
 144       2. edge first two points are not equal
 145       3. geom is simple
 146    */
 147
 148   const Coordinate& pt1 = edge.getCoordinateN(0);
 149   const Coordinate& pt2 = edge.getCoordinateN(1);
 150
 151   /*
 152    * We move the coordinate somewhat closer, to avoid
 153    * vertices of the geometry being checked (geom).
 154    *
 155    * This is mostly only needed when one of the two points
 156    * of the edge is an endpoint of a _closed_ geom.
 157    * We have an unit test for this...
 158    */
 159   Coordinate pt1i = LinearLocation::pointAlongSegmentByFraction(pt1, pt2, 0.1);
 160   Coordinate pt2i = LinearLocation::pointAlongSegmentByFraction(pt1, pt2, 0.9);
 161
 162   LinearLocation l1 = LocationIndexOfPoint::indexOf(&geom, pt1i);
 163   LinearLocation l2 = LocationIndexOfPoint::indexOf(&geom, pt2i);
 164
 165   return l1.compareTo(l2) < 0;
 166 }
 167
 168 } // namespace geos.operation.sharedpaths
 169 } // namespace geos::operation
 170 } // namespace geos
 171