/src/gdal/ogr/ograssemblepolygon.cpp
Line | Count | Source (jump to first uncovered line) |
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(OGRLineString *poLine1, int iPoint1, |
35 | | OGRLineString *poLine2, int iPoint2, |
36 | | double *pdfDistance) |
37 | | |
38 | 0 | { |
39 | 0 | if (pdfDistance == nullptr || *pdfDistance == 0) |
40 | 0 | { |
41 | 0 | if (poLine1->getX(iPoint1) == poLine2->getX(iPoint2) && |
42 | 0 | poLine1->getY(iPoint1) == poLine2->getY(iPoint2)) |
43 | 0 | { |
44 | 0 | if (pdfDistance) |
45 | 0 | *pdfDistance = 0.0; |
46 | 0 | return true; |
47 | 0 | } |
48 | 0 | return false; |
49 | 0 | } |
50 | | |
51 | 0 | const double dfDeltaX = |
52 | 0 | std::abs(poLine1->getX(iPoint1) - poLine2->getX(iPoint2)); |
53 | |
|
54 | 0 | if (dfDeltaX > *pdfDistance) |
55 | 0 | return false; |
56 | | |
57 | 0 | const double dfDeltaY = |
58 | 0 | std::abs(poLine1->getY(iPoint1) - poLine2->getY(iPoint2)); |
59 | |
|
60 | 0 | if (dfDeltaY > *pdfDistance) |
61 | 0 | return false; |
62 | | |
63 | 0 | const double dfDistance = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY); |
64 | |
|
65 | 0 | if (dfDistance < *pdfDistance) |
66 | 0 | { |
67 | 0 | *pdfDistance = dfDistance; |
68 | 0 | return true; |
69 | 0 | } |
70 | | |
71 | 0 | return false; |
72 | 0 | } |
73 | | |
74 | | /************************************************************************/ |
75 | | /* AddEdgeToRing() */ |
76 | | /************************************************************************/ |
77 | | |
78 | | static void AddEdgeToRing(OGRLinearRing *poRing, OGRLineString *poLine, |
79 | | bool bReverse, double dfTolerance) |
80 | | |
81 | 0 | { |
82 | | /* -------------------------------------------------------------------- */ |
83 | | /* Establish order and range of traverse. */ |
84 | | /* -------------------------------------------------------------------- */ |
85 | 0 | const int nVertToAdd = poLine->getNumPoints(); |
86 | |
|
87 | 0 | int iStart = bReverse ? nVertToAdd - 1 : 0; |
88 | 0 | const int iEnd = bReverse ? 0 : nVertToAdd - 1; |
89 | 0 | const int iStep = bReverse ? -1 : 1; |
90 | | |
91 | | /* -------------------------------------------------------------------- */ |
92 | | /* Should we skip a repeating vertex? */ |
93 | | /* -------------------------------------------------------------------- */ |
94 | 0 | if (poRing->getNumPoints() > 0 && |
95 | 0 | CheckPoints(poRing, poRing->getNumPoints() - 1, poLine, iStart, |
96 | 0 | &dfTolerance)) |
97 | 0 | { |
98 | 0 | iStart += iStep; |
99 | 0 | } |
100 | |
|
101 | 0 | poRing->addSubLineString(poLine, iStart, iEnd); |
102 | 0 | } |
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. |
120 | | * |
121 | | */ |
122 | | |
123 | | OGRGeometryH OGRBuildPolygonFromEdges(OGRGeometryH hLines, |
124 | | CPL_UNUSED int bBestEffort, |
125 | | int bAutoClose, double dfTolerance, |
126 | | OGRErr *peErr) |
127 | | |
128 | 0 | { |
129 | 0 | 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 | 0 | OGRGeometry *poGeom = OGRGeometry::FromHandle(hLines); |
141 | 0 | if (wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection) |
142 | 0 | { |
143 | 0 | for (auto &&poMember : poGeom->toGeometryCollection()) |
144 | 0 | { |
145 | 0 | if (wkbFlatten(poMember->getGeometryType()) != wkbLineString) |
146 | 0 | { |
147 | 0 | if (peErr != nullptr) |
148 | 0 | *peErr = OGRERR_FAILURE; |
149 | 0 | CPLError(CE_Failure, CPLE_NotSupported, |
150 | 0 | "The geometry collection contains " |
151 | 0 | "non-line string geometries"); |
152 | 0 | return nullptr; |
153 | 0 | } |
154 | 0 | } |
155 | 0 | } |
156 | 0 | 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 | 0 | bool bSuccess = true; |
168 | 0 | OGRGeometryCollection *poLines = poGeom->toGeometryCollection(); |
169 | 0 | std::vector<OGRLinearRing *> apoRings; |
170 | | |
171 | | /* -------------------------------------------------------------------- */ |
172 | | /* Setup array of line markers indicating if they have been */ |
173 | | /* added to a ring yet. */ |
174 | | /* -------------------------------------------------------------------- */ |
175 | 0 | const int nEdges = poLines->getNumGeometries(); |
176 | 0 | std::list<OGRLineString *> oListEdges; |
177 | 0 | for (int i = 0; i < nEdges; i++) |
178 | 0 | { |
179 | 0 | OGRLineString *poLine = poLines->getGeometryRef(i)->toLineString(); |
180 | 0 | if (poLine->getNumPoints() >= 2) |
181 | 0 | { |
182 | 0 | oListEdges.push_back(poLine); |
183 | 0 | } |
184 | 0 | } |
185 | | |
186 | | /* ==================================================================== */ |
187 | | /* Loop generating rings. */ |
188 | | /* ==================================================================== */ |
189 | 0 | while (!oListEdges.empty()) |
190 | 0 | { |
191 | | /* -------------------------------------------------------------------- |
192 | | */ |
193 | | /* Find the first unconsumed edge. */ |
194 | | /* -------------------------------------------------------------------- |
195 | | */ |
196 | 0 | OGRLineString *poLine = oListEdges.front(); |
197 | 0 | oListEdges.erase(oListEdges.begin()); |
198 | | |
199 | | /* -------------------------------------------------------------------- |
200 | | */ |
201 | | /* Start a new ring, copying in the current line directly */ |
202 | | /* -------------------------------------------------------------------- |
203 | | */ |
204 | 0 | OGRLinearRing *poRing = new OGRLinearRing(); |
205 | |
|
206 | 0 | AddEdgeToRing(poRing, poLine, FALSE, 0); |
207 | | |
208 | | /* ==================================================================== |
209 | | */ |
210 | | /* Loop adding edges to this ring until we make a whole pass */ |
211 | | /* within finding anything to add. */ |
212 | | /* ==================================================================== |
213 | | */ |
214 | 0 | bool bWorkDone = true; |
215 | 0 | double dfBestDist = dfTolerance; |
216 | |
|
217 | 0 | while (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1, |
218 | 0 | nullptr) && |
219 | 0 | !oListEdges.empty() && bWorkDone) |
220 | 0 | { |
221 | 0 | bool bReverse = false; |
222 | |
|
223 | 0 | bWorkDone = false; |
224 | 0 | dfBestDist = dfTolerance; |
225 | | |
226 | | // We consider linking the end to the beginning. If this is |
227 | | // closer than any other option we will just close the loop. |
228 | | |
229 | | // CheckPoints(poRing, 0, poRing, poRing->getNumPoints()-1, |
230 | | // &dfBestDist); |
231 | | |
232 | | // Find unused edge with end point closest to our loose end. |
233 | 0 | OGRLineString *poBestEdge = nullptr; |
234 | 0 | std::list<OGRLineString *>::iterator oBestIter; |
235 | 0 | for (auto oIter = oListEdges.begin(); oIter != oListEdges.end(); |
236 | 0 | ++oIter) |
237 | 0 | { |
238 | 0 | poLine = *oIter; |
239 | |
|
240 | 0 | if (CheckPoints(poLine, 0, poRing, poRing->getNumPoints() - 1, |
241 | 0 | &dfBestDist)) |
242 | 0 | { |
243 | 0 | poBestEdge = poLine; |
244 | 0 | oBestIter = oIter; |
245 | 0 | bReverse = false; |
246 | 0 | } |
247 | 0 | if (CheckPoints(poLine, poLine->getNumPoints() - 1, poRing, |
248 | 0 | poRing->getNumPoints() - 1, &dfBestDist)) |
249 | 0 | { |
250 | 0 | poBestEdge = poLine; |
251 | 0 | oBestIter = oIter; |
252 | 0 | bReverse = true; |
253 | 0 | } |
254 | | |
255 | | // If we found an exact match, jump now. |
256 | 0 | if (dfBestDist == 0.0 && poBestEdge != nullptr) |
257 | 0 | break; |
258 | 0 | } |
259 | | |
260 | | // We found one within tolerance - add it. |
261 | 0 | if (poBestEdge) |
262 | 0 | { |
263 | 0 | AddEdgeToRing(poRing, poBestEdge, bReverse, dfTolerance); |
264 | |
|
265 | 0 | oListEdges.erase(oBestIter); |
266 | 0 | bWorkDone = true; |
267 | 0 | } |
268 | 0 | } |
269 | | |
270 | | /* -------------------------------------------------------------------- |
271 | | */ |
272 | | /* Did we fail to complete the ring? */ |
273 | | /* -------------------------------------------------------------------- |
274 | | */ |
275 | 0 | dfBestDist = dfTolerance; |
276 | |
|
277 | 0 | if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1, |
278 | 0 | &dfBestDist)) |
279 | 0 | { |
280 | 0 | CPLDebug("OGR", |
281 | 0 | "Failed to close ring %d.\n" |
282 | 0 | "End Points are: (%.8f,%.7f) and (%.7f,%.7f)", |
283 | 0 | static_cast<int>(apoRings.size()), poRing->getX(0), |
284 | 0 | poRing->getY(0), poRing->getX(poRing->getNumPoints() - 1), |
285 | 0 | poRing->getY(poRing->getNumPoints() - 1)); |
286 | |
|
287 | 0 | bSuccess = false; |
288 | 0 | } |
289 | | |
290 | | /* -------------------------------------------------------------------- |
291 | | */ |
292 | | /* Do we need to auto-close this ring? */ |
293 | | /* -------------------------------------------------------------------- |
294 | | */ |
295 | 0 | dfBestDist = dfTolerance; |
296 | |
|
297 | 0 | if (bAutoClose) |
298 | 0 | { |
299 | 0 | if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1, |
300 | 0 | &dfBestDist)) |
301 | 0 | { |
302 | 0 | poRing->addPoint(poRing->getX(0), poRing->getY(0), |
303 | 0 | poRing->getZ(0)); |
304 | 0 | } |
305 | 0 | else if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1, |
306 | 0 | nullptr)) |
307 | 0 | { |
308 | | // The endpoints are very close but do not exactly match. |
309 | | // Alter the last one so it is equal to the first, to prevent |
310 | | // invalid self-intersecting rings. |
311 | 0 | poRing->setPoint(poRing->getNumPoints() - 1, poRing->getX(0), |
312 | 0 | poRing->getY(0), poRing->getZ(0)); |
313 | 0 | } |
314 | 0 | } |
315 | |
|
316 | 0 | apoRings.push_back(poRing); |
317 | 0 | } // Next ring. |
318 | | |
319 | | /* -------------------------------------------------------------------- */ |
320 | | /* Identify exterior ring - it will be the largest. #3610 */ |
321 | | /* -------------------------------------------------------------------- */ |
322 | 0 | double maxarea = 0.0; |
323 | 0 | int maxring = -1; |
324 | 0 | OGREnvelope tenv; |
325 | |
|
326 | 0 | for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn) |
327 | 0 | { |
328 | 0 | apoRings[rn]->getEnvelope(&tenv); |
329 | 0 | const double tarea = (tenv.MaxX - tenv.MinX) * (tenv.MaxY - tenv.MinY); |
330 | 0 | if (tarea > maxarea) |
331 | 0 | { |
332 | 0 | maxarea = tarea; |
333 | 0 | maxring = rn; |
334 | 0 | } |
335 | 0 | } |
336 | |
|
337 | 0 | OGRPolygon *poPolygon = new OGRPolygon(); |
338 | |
|
339 | 0 | if (maxring != -1) |
340 | 0 | { |
341 | 0 | poPolygon->addRingDirectly(apoRings[maxring]); |
342 | 0 | for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn) |
343 | 0 | { |
344 | 0 | if (rn == maxring) |
345 | 0 | continue; |
346 | 0 | poPolygon->addRingDirectly(apoRings[rn]); |
347 | 0 | } |
348 | 0 | } |
349 | 0 | else |
350 | 0 | { |
351 | 0 | for (auto &poRing : apoRings) |
352 | 0 | delete poRing; |
353 | 0 | } |
354 | |
|
355 | 0 | if (peErr != nullptr) |
356 | 0 | { |
357 | 0 | *peErr = bSuccess ? OGRERR_NONE : OGRERR_FAILURE; |
358 | 0 | } |
359 | |
|
360 | 0 | return OGRGeometry::ToHandle(poPolygon); |
361 | 0 | } |