1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2016 Daniel Baston
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: precision/MinimumClearance.java (f6187ee2 JTS-1.14)
16 *
17 **********************************************************************/
18
19 #include <geos/algorithm/Distance.h>
20 #include <geos/precision/MinimumClearance.h>
21 #include <geos/index/strtree/STRtree.h>
22 #include <geos/geom/GeometryFactory.h>
23 #include <geos/geom/CoordinateArraySequenceFactory.h>
24 #include <geos/operation/distance/FacetSequenceTreeBuilder.h>
25 #include <geos/geom/LineSegment.h>
26
27 using namespace geos::geom;
28 using namespace geos::operation::distance;
29 using namespace geos::index::strtree;
30
31 namespace geos {
32 namespace precision {
33
MinimumClearance(const Geometry * g)34 MinimumClearance::MinimumClearance(const Geometry* g) : inputGeom(g) {}
35
36 double
getDistance()37 MinimumClearance::getDistance()
38 {
39 compute();
40 return minClearance;
41 }
42
43 std::unique_ptr<LineString>
getLine()44 MinimumClearance::getLine()
45 {
46 compute();
47
48 // return empty line string if no min pts were found
49 if(minClearance == std::numeric_limits<double>::infinity()) {
50 return inputGeom->getFactory()->createLineString();
51 }
52
53 return inputGeom->getFactory()->createLineString(minClearancePts->clone());
54 }
55
56 void
compute()57 MinimumClearance::compute()
58 {
59 class MinClearanceDistance : public ItemDistance {
60 private:
61 double minDist;
62 std::vector<Coordinate> minPts;
63
64 void
65 updatePts(const Coordinate& p, const Coordinate& seg0, const Coordinate& seg1)
66 {
67 LineSegment seg(seg0, seg1);
68
69 minPts[0] = p;
70 seg.closestPoint(p, minPts[1]);
71 }
72
73 public:
74 MinClearanceDistance() :
75 minDist(std::numeric_limits<double>::infinity()),
76 minPts(std::vector<Coordinate>(2))
77 {}
78
79 const std::vector<Coordinate>*
80 getCoordinates()
81 {
82 return &minPts;
83 }
84
85 double
86 distance(const ItemBoundable* b1, const ItemBoundable* b2) override
87 {
88 FacetSequence* fs1 = static_cast<FacetSequence*>(b1->getItem());
89 FacetSequence* fs2 = static_cast<FacetSequence*>(b2->getItem());
90
91 minDist = std::numeric_limits<double>::infinity();
92
93 return distance(fs1, fs2);
94 }
95
96 double
97 distance(const FacetSequence* fs1, const FacetSequence* fs2)
98 {
99 // Compute MinClearance distance metric
100
101 vertexDistance(fs1, fs2);
102 if(fs1->size() == 1 && fs2->size() == 1) {
103 return minDist;
104 }
105 if(minDist <= 0.0) {
106 return minDist;
107 }
108
109 segmentDistance(fs1, fs2);
110 if(minDist <= 0.0) {
111 return minDist;
112 }
113
114 segmentDistance(fs2, fs1);
115 return minDist;
116 }
117
118 double
119 vertexDistance(const FacetSequence* fs1, const FacetSequence* fs2)
120 {
121 for(size_t i1 = 0; i1 < fs1->size(); i1++) {
122 for(size_t i2 = 0; i2 < fs2->size(); i2++) {
123 const Coordinate* p1 = fs1->getCoordinate(i1);
124 const Coordinate* p2 = fs2->getCoordinate(i2);
125 if(!p1->equals2D(*p2)) {
126 double d = p1->distance(*p2);
127 if(d < minDist) {
128 minDist = d;
129 minPts[0] = *p1;
130 minPts[1] = *p2;
131 if(d == 0.0) {
132 return d;
133 }
134 }
135 }
136 }
137 }
138 return minDist;
139 }
140
141 double
142 segmentDistance(const FacetSequence* fs1, const FacetSequence* fs2)
143 {
144 for(size_t i1 = 0; i1 < fs1->size(); i1++) {
145 for(size_t i2 = 1; i2 < fs2->size(); i2++) {
146 const Coordinate* p = fs1->getCoordinate(i1);
147
148 const Coordinate* seg0 = fs2->getCoordinate(i2 - 1);
149 const Coordinate* seg1 = fs2->getCoordinate(i2);
150
151 if(!(p->equals2D(*seg0) || p->equals2D(*seg1))) {
152 double d = geos::algorithm::Distance::pointToSegment(*p, *seg0, *seg1);
153 if(d < minDist) {
154 minDist = d;
155 updatePts(*p, *seg0, *seg1);
156 if(d == 0.0) {
157 return d;
158 }
159 }
160 }
161 }
162 }
163 return minDist;
164 }
165 };
166
167 // already computed
168 if(minClearancePts.get() != nullptr) {
169 return;
170 }
171
172 // initialize to "No Distance Exists" state
173 minClearancePts = std::unique_ptr<CoordinateSequence>(inputGeom->getFactory()->getCoordinateSequenceFactory()->create(2,
174 2));
175 minClearance = std::numeric_limits<double>::infinity();
176
177 // handle empty geometries
178 if(inputGeom->isEmpty()) {
179 return;
180 }
181
182 auto tree = FacetSequenceTreeBuilder::build(inputGeom);
183 MinClearanceDistance mcd;
184 std::pair<const void*, const void*> nearest = tree->nearestNeighbour(&mcd);
185
186 minClearance = mcd.distance(
187 static_cast<const FacetSequence*>(nearest.first),
188 static_cast<const FacetSequence*>(nearest.second));
189
190 const std::vector<Coordinate>* minClearancePtsVec = mcd.getCoordinates();
191 minClearancePts->setAt((*minClearancePtsVec)[0], 0);
192 minClearancePts->setAt((*minClearancePtsVec)[1], 1);
193 }
194
195
196 }
197 }
198