/src/gdal/ogr/ograssemblepolygon.cpp
Line | Count | Source |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: S-57 Reader |
4 | | * Purpose: Implements polygon assembly from a bunch of arcs. |
5 | | * Author: Frank Warmerdam, warmerdam@pobox.com |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 1999, Frank Warmerdam |
9 | | * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com> |
10 | | * |
11 | | * SPDX-License-Identifier: MIT |
12 | | ****************************************************************************/ |
13 | | |
14 | | #include "cpl_port.h" |
15 | | #include "ogr_api.h" |
16 | | |
17 | | #include <cmath> |
18 | | #include <cstddef> |
19 | | #include <list> |
20 | | #include <vector> |
21 | | |
22 | | #include "ogr_core.h" |
23 | | #include "ogr_geometry.h" |
24 | | #include "cpl_conv.h" |
25 | | #include "cpl_error.h" |
26 | | |
27 | | /************************************************************************/ |
28 | | /* CheckPoints() */ |
29 | | /* */ |
30 | | /* Check if two points are closer than the current best */ |
31 | | /* distance. Update the current best distance if they are. */ |
32 | | /************************************************************************/ |
33 | | |
34 | | static bool CheckPoints(const OGRLineString *poLine1, int iPoint1, |
35 | | const OGRLineString *poLine2, int iPoint2, |
36 | | double *pdfDistance) |
37 | | |
38 | 22.0M | { |
39 | 22.0M | if (pdfDistance == nullptr || *pdfDistance == 0) |
40 | 1.75M | { |
41 | 1.75M | if (poLine1->getX(iPoint1) == poLine2->getX(iPoint2) && |
42 | 940k | poLine1->getY(iPoint1) == poLine2->getY(iPoint2)) |
43 | 721k | { |
44 | 721k | if (pdfDistance) |
45 | 239k | *pdfDistance = 0.0; |
46 | 721k | return true; |
47 | 721k | } |
48 | 1.03M | return false; |
49 | 1.75M | } |
50 | | |
51 | 20.2M | const double dfDeltaX = |
52 | 20.2M | std::abs(poLine1->getX(iPoint1) - poLine2->getX(iPoint2)); |
53 | | |
54 | 20.2M | if (dfDeltaX > *pdfDistance) |
55 | 14.6M | return false; |
56 | | |
57 | 5.67M | const double dfDeltaY = |
58 | 5.67M | std::abs(poLine1->getY(iPoint1) - poLine2->getY(iPoint2)); |
59 | | |
60 | 5.67M | if (dfDeltaY > *pdfDistance) |
61 | 1.00M | return false; |
62 | | |
63 | 4.66M | const double dfDistance = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY); |
64 | | |
65 | 4.66M | if (dfDistance < *pdfDistance) |
66 | 1.00M | { |
67 | 1.00M | *pdfDistance = dfDistance; |
68 | 1.00M | return true; |
69 | 1.00M | } |
70 | | |
71 | 3.66M | return false; |
72 | 4.66M | } |
73 | | |
74 | | /************************************************************************/ |
75 | | /* AddEdgeToRing() */ |
76 | | /************************************************************************/ |
77 | | |
78 | | static void AddEdgeToRing(OGRLinearRing *poRing, const OGRLineString *poLine, |
79 | | bool bReverse, double dfTolerance) |
80 | | |
81 | 592k | { |
82 | | /* -------------------------------------------------------------------- */ |
83 | | /* Establish order and range of traverse. */ |
84 | | /* -------------------------------------------------------------------- */ |
85 | 592k | const int nVertToAdd = poLine->getNumPoints(); |
86 | | |
87 | 592k | int iStart = bReverse ? nVertToAdd - 1 : 0; |
88 | 592k | const int iEnd = bReverse ? 0 : nVertToAdd - 1; |
89 | 592k | const int iStep = bReverse ? -1 : 1; |
90 | | |
91 | | /* -------------------------------------------------------------------- */ |
92 | | /* Should we skip a repeating vertex? */ |
93 | | /* -------------------------------------------------------------------- */ |
94 | 592k | if (poRing->getNumPoints() > 0 && |
95 | 255k | CheckPoints(poRing, poRing->getNumPoints() - 1, poLine, iStart, |
96 | 255k | &dfTolerance)) |
97 | 255k | { |
98 | 255k | iStart += iStep; |
99 | 255k | } |
100 | | |
101 | 592k | poRing->addSubLineString(poLine, iStart, iEnd); |
102 | 592k | } |
103 | | |
104 | | /************************************************************************/ |
105 | | /* OGRBuildPolygonFromEdges() */ |
106 | | /************************************************************************/ |
107 | | |
108 | | /** |
109 | | * Build a ring from a bunch of arcs. |
110 | | * |
111 | | * @param hLines handle to an OGRGeometryCollection (or OGRMultiLineString) |
112 | | * containing the line string geometries to be built into rings. |
113 | | * @param bBestEffort not yet implemented???. |
114 | | * @param bAutoClose indicates if the ring should be close when first and |
115 | | * last points of the ring are the same. |
116 | | * @param dfTolerance tolerance into which two arcs are considered |
117 | | * close enough to be joined. |
118 | | * @param peErr OGRERR_NONE on success, or OGRERR_FAILURE on failure. |
119 | | * @return a handle to the new geometry, a polygon or multipolygon. |
120 | | * |
121 | | */ |
122 | | |
123 | | OGRGeometryH OGRBuildPolygonFromEdges(OGRGeometryH hLines, |
124 | | CPL_UNUSED int bBestEffort, |
125 | | int bAutoClose, double dfTolerance, |
126 | | OGRErr *peErr) |
127 | | |
128 | 249k | { |
129 | 249k | if (hLines == nullptr) |
130 | 0 | { |
131 | 0 | if (peErr != nullptr) |
132 | 0 | *peErr = OGRERR_NONE; |
133 | 0 | return nullptr; |
134 | 0 | } |
135 | | |
136 | | /* -------------------------------------------------------------------- */ |
137 | | /* Check for the case of a geometrycollection that can be */ |
138 | | /* promoted to MultiLineString. */ |
139 | | /* -------------------------------------------------------------------- */ |
140 | 249k | OGRGeometry *poGeom = OGRGeometry::FromHandle(hLines); |
141 | 249k | if (wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection) |
142 | 114k | { |
143 | 114k | for (auto &&poMember : poGeom->toGeometryCollection()) |
144 | 246k | { |
145 | 246k | if (wkbFlatten(poMember->getGeometryType()) != wkbLineString) |
146 | 6.18k | { |
147 | 6.18k | if (peErr != nullptr) |
148 | 6.18k | *peErr = OGRERR_FAILURE; |
149 | 6.18k | CPLError(CE_Failure, CPLE_NotSupported, |
150 | 6.18k | "The geometry collection contains " |
151 | 6.18k | "non-line string geometries"); |
152 | 6.18k | return nullptr; |
153 | 6.18k | } |
154 | 246k | } |
155 | 114k | } |
156 | 135k | else if (wkbFlatten(poGeom->getGeometryType()) != wkbMultiLineString) |
157 | 0 | { |
158 | 0 | if (peErr != nullptr) |
159 | 0 | *peErr = OGRERR_FAILURE; |
160 | 0 | CPLError(CE_Failure, CPLE_NotSupported, |
161 | 0 | "The passed geometry is not an OGRGeometryCollection " |
162 | 0 | "(or OGRMultiLineString) " |
163 | 0 | "containing line string geometries"); |
164 | 0 | return nullptr; |
165 | 0 | } |
166 | | |
167 | 243k | bool bSuccess = true; |
168 | 243k | const OGRGeometryCollection *poLines = poGeom->toGeometryCollection(); |
169 | 243k | std::vector<std::unique_ptr<OGRGeometry>> apoPolys; |
170 | | |
171 | | /* -------------------------------------------------------------------- */ |
172 | | /* Setup array of line markers indicating if they have been */ |
173 | | /* added to a ring yet. */ |
174 | | /* -------------------------------------------------------------------- */ |
175 | 243k | const int nEdges = poLines->getNumGeometries(); |
176 | 243k | std::list<const OGRLineString *> oListEdges; |
177 | 848k | for (int i = 0; i < nEdges; i++) |
178 | 604k | { |
179 | 604k | const OGRLineString *poLine = |
180 | 604k | poLines->getGeometryRef(i)->toLineString(); |
181 | 604k | if (poLine->getNumPoints() >= 2) |
182 | 592k | { |
183 | 592k | oListEdges.push_back(poLine); |
184 | 592k | } |
185 | 604k | } |
186 | | |
187 | | /* ==================================================================== */ |
188 | | /* Loop generating rings. */ |
189 | | /* ==================================================================== */ |
190 | 580k | while (!oListEdges.empty()) |
191 | 336k | { |
192 | | /* -------------------------------------------------------------------- |
193 | | */ |
194 | | /* Find the first unconsumed edge. */ |
195 | | /* -------------------------------------------------------------------- |
196 | | */ |
197 | 336k | const OGRLineString *poLine = oListEdges.front(); |
198 | 336k | oListEdges.erase(oListEdges.begin()); |
199 | | |
200 | | /* -------------------------------------------------------------------- |
201 | | */ |
202 | | /* Start a new ring, copying in the current line directly */ |
203 | | /* -------------------------------------------------------------------- |
204 | | */ |
205 | 336k | auto poRing = std::make_unique<OGRLinearRing>(); |
206 | | |
207 | 336k | AddEdgeToRing(poRing.get(), poLine, FALSE, 0); |
208 | | |
209 | | /* ==================================================================== |
210 | | */ |
211 | | /* Loop adding edges to this ring until we make a whole pass */ |
212 | | /* within finding anything to add. */ |
213 | | /* ==================================================================== |
214 | | */ |
215 | 336k | bool bWorkDone = true; |
216 | 336k | double dfBestDist = dfTolerance; |
217 | | |
218 | 615k | while (!CheckPoints(poRing.get(), 0, poRing.get(), |
219 | 615k | poRing->getNumPoints() - 1, nullptr) && |
220 | 336k | !oListEdges.empty() && bWorkDone) |
221 | 278k | { |
222 | 278k | bool bReverse = false; |
223 | | |
224 | 278k | bWorkDone = false; |
225 | 278k | dfBestDist = dfTolerance; |
226 | | |
227 | | // We consider linking the end to the beginning. If this is |
228 | | // closer than any other option we will just close the loop. |
229 | | |
230 | | // CheckPoints(poRing, 0, poRing, poRing->getNumPoints()-1, |
231 | | // &dfBestDist); |
232 | | |
233 | | // Find unused edge with end point closest to our loose end. |
234 | 278k | const OGRLineString *poBestEdge = nullptr; |
235 | 278k | std::list<const OGRLineString *>::iterator oBestIter; |
236 | 10.2M | for (auto oIter = oListEdges.begin(); oIter != oListEdges.end(); |
237 | 9.93M | ++oIter) |
238 | 10.1M | { |
239 | 10.1M | poLine = *oIter; |
240 | | |
241 | 10.1M | if (CheckPoints(poLine, 0, poRing.get(), |
242 | 10.1M | poRing->getNumPoints() - 1, &dfBestDist)) |
243 | 237k | { |
244 | 237k | poBestEdge = poLine; |
245 | 237k | oBestIter = oIter; |
246 | 237k | bReverse = false; |
247 | 237k | } |
248 | 10.1M | if (CheckPoints(poLine, poLine->getNumPoints() - 1, |
249 | 10.1M | poRing.get(), poRing->getNumPoints() - 1, |
250 | 10.1M | &dfBestDist)) |
251 | 251k | { |
252 | 251k | poBestEdge = poLine; |
253 | 251k | oBestIter = oIter; |
254 | 251k | bReverse = true; |
255 | 251k | } |
256 | | |
257 | | // If we found an exact match, jump now. |
258 | 10.1M | if (dfBestDist == 0.0 && poBestEdge != nullptr) |
259 | 247k | break; |
260 | 10.1M | } |
261 | | |
262 | | // We found one within tolerance - add it. |
263 | 278k | if (poBestEdge) |
264 | 255k | { |
265 | 255k | AddEdgeToRing(poRing.get(), poBestEdge, bReverse, dfTolerance); |
266 | | |
267 | 255k | oListEdges.erase(oBestIter); |
268 | 255k | bWorkDone = true; |
269 | 255k | } |
270 | 278k | } |
271 | | |
272 | | /* -------------------------------------------------------------------- |
273 | | */ |
274 | | /* Did we fail to complete the ring? */ |
275 | | /* -------------------------------------------------------------------- |
276 | | */ |
277 | 336k | dfBestDist = dfTolerance; |
278 | | |
279 | 336k | if (!CheckPoints(poRing.get(), 0, poRing.get(), |
280 | 336k | poRing->getNumPoints() - 1, &dfBestDist)) |
281 | 47.7k | { |
282 | 47.7k | CPLDebug("OGR", |
283 | 47.7k | "Failed to close ring %d.\n" |
284 | 47.7k | "End Points are: (%.8f,%.7f) and (%.7f,%.7f)", |
285 | 47.7k | static_cast<int>(apoPolys.size()), poRing->getX(0), |
286 | 47.7k | poRing->getY(0), poRing->getX(poRing->getNumPoints() - 1), |
287 | 47.7k | poRing->getY(poRing->getNumPoints() - 1)); |
288 | | |
289 | 47.7k | bSuccess = false; |
290 | 47.7k | } |
291 | | |
292 | | /* -------------------------------------------------------------------- |
293 | | */ |
294 | | /* Do we need to auto-close this ring? */ |
295 | | /* -------------------------------------------------------------------- |
296 | | */ |
297 | 336k | dfBestDist = dfTolerance; |
298 | | |
299 | 336k | if (bAutoClose) |
300 | 243k | { |
301 | 243k | if (!CheckPoints(poRing.get(), 0, poRing.get(), |
302 | 243k | poRing->getNumPoints() - 1, &dfBestDist)) |
303 | 30.4k | { |
304 | 30.4k | poRing->addPoint(poRing->getX(0), poRing->getY(0), |
305 | 30.4k | poRing->getZ(0)); |
306 | 30.4k | } |
307 | 213k | else if (!CheckPoints(poRing.get(), 0, poRing.get(), |
308 | 213k | poRing->getNumPoints() - 1, nullptr)) |
309 | 10.1k | { |
310 | | // The endpoints are very close but do not exactly match. |
311 | | // Alter the last one so it is equal to the first, to prevent |
312 | | // invalid self-intersecting rings. |
313 | 10.1k | poRing->setPoint(poRing->getNumPoints() - 1, poRing->getX(0), |
314 | 10.1k | poRing->getY(0), poRing->getZ(0)); |
315 | 10.1k | } |
316 | 243k | } |
317 | | |
318 | 336k | auto poPoly = std::make_unique<OGRPolygon>(); |
319 | 336k | poPoly->addRing(std::move(poRing)); |
320 | 336k | apoPolys.push_back(std::move(poPoly)); |
321 | 336k | } // Next ring. |
322 | | |
323 | 243k | if (peErr != nullptr) |
324 | 108k | { |
325 | 108k | *peErr = bSuccess ? OGRERR_NONE : OGRERR_FAILURE; |
326 | 108k | } |
327 | | |
328 | 243k | return OGRGeometry::ToHandle( |
329 | 243k | OGRGeometryFactory::organizePolygons(apoPolys).release()); |
330 | 249k | } |