1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2010 Sandro Santilli <strk@kbt.io>
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@kbt.io)
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
sharedPathsOp(const Geometry & g1,const Geometry & g2,PathList & sameDirection,PathList & oppositeDirection)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 */
SharedPathsOp(const geom::Geometry & g1,const geom::Geometry & g2)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
checkLinealInput(const geom::Geometry & g)66 SharedPathsOp::checkLinealInput(const geom::Geometry& g)
67 {
68 if(! dynamic_cast<const LineString*>(&g) &&
69 ! dynamic_cast<const MultiLineString*>(&g)) {
70 throw util::IllegalArgumentException("Geometry is not lineal");
71 }
72 }
73
74 /* public */
75 void
getSharedPaths(PathList & forwDir,PathList & backDir)76 SharedPathsOp::getSharedPaths(PathList& forwDir, PathList& backDir)
77 {
78 PathList paths;
79 findLinearIntersections(paths);
80 for(size_t i = 0, n = paths.size(); i < n; ++i) {
81 LineString* path = paths[i];
82 if(isSameDirection(*path)) {
83 forwDir.push_back(path);
84 }
85 else {
86 backDir.push_back(path);
87 }
88 }
89 }
90
91 /* static private */
92 void
clearEdges(PathList & edges)93 SharedPathsOp::clearEdges(PathList& edges)
94 {
95 for(PathList::const_iterator
96 i = edges.begin(), e = edges.end();
97 i != e; ++i) {
98 delete *i;
99 }
100 edges.clear();
101 }
102
103 /* private */
104 void
findLinearIntersections(PathList & to)105 SharedPathsOp::findLinearIntersections(PathList& to)
106 {
107 using geos::operation::overlay::OverlayOp;
108
109 // TODO: optionally use the tolerance,
110 // snapping _g2 over _g1 ?
111
112 std::unique_ptr<Geometry> full(OverlayOp::overlayOp(
113 &_g1, &_g2, OverlayOp::opINTERSECTION));
114
115 // NOTE: intersection of equal lines yelds splitted lines,
116 // should we sew them back ?
117
118 for(size_t i = 0, n = full->getNumGeometries(); i < n; ++i) {
119 const Geometry* sub = full->getGeometryN(i);
120 const LineString* path = dynamic_cast<const LineString*>(sub);
121 if(path && ! path->isEmpty()) {
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
isForward(const geom::LineString & edge,const geom::Geometry & geom)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
172