Line data Source code
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 30186 : static bool CheckPoints(const OGRLineString *poLine1, int iPoint1,
35 : const OGRLineString *poLine2, int iPoint2,
36 : double *pdfDistance)
37 :
38 : {
39 30186 : if (pdfDistance == nullptr || *pdfDistance == 0)
40 : {
41 37633 : if (poLine1->getX(iPoint1) == poLine2->getX(iPoint2) &&
42 7691 : poLine1->getY(iPoint1) == poLine2->getY(iPoint2))
43 : {
44 5634 : if (pdfDistance)
45 5198 : *pdfDistance = 0.0;
46 5634 : return true;
47 : }
48 24308 : return false;
49 : }
50 :
51 : const double dfDeltaX =
52 244 : std::abs(poLine1->getX(iPoint1) - poLine2->getX(iPoint2));
53 :
54 244 : if (dfDeltaX > *pdfDistance)
55 30 : return false;
56 :
57 : const double dfDeltaY =
58 214 : std::abs(poLine1->getY(iPoint1) - poLine2->getY(iPoint2));
59 :
60 214 : if (dfDeltaY > *pdfDistance)
61 6 : return false;
62 :
63 208 : const double dfDistance = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
64 :
65 208 : if (dfDistance < *pdfDistance)
66 : {
67 208 : *pdfDistance = dfDistance;
68 208 : return true;
69 : }
70 :
71 0 : return false;
72 : }
73 :
74 : /************************************************************************/
75 : /* AddEdgeToRing() */
76 : /************************************************************************/
77 :
78 2868 : static void AddEdgeToRing(OGRLinearRing *poRing, const OGRLineString *poLine,
79 : bool bReverse, double dfTolerance)
80 :
81 : {
82 : /* -------------------------------------------------------------------- */
83 : /* Establish order and range of traverse. */
84 : /* -------------------------------------------------------------------- */
85 2868 : const int nVertToAdd = poLine->getNumPoints();
86 :
87 2868 : int iStart = bReverse ? nVertToAdd - 1 : 0;
88 2868 : const int iEnd = bReverse ? 0 : nVertToAdd - 1;
89 2868 : const int iStep = bReverse ? -1 : 1;
90 :
91 : /* -------------------------------------------------------------------- */
92 : /* Should we skip a repeating vertex? */
93 : /* -------------------------------------------------------------------- */
94 5348 : if (poRing->getNumPoints() > 0 &&
95 2480 : CheckPoints(poRing, poRing->getNumPoints() - 1, poLine, iStart,
96 : &dfTolerance))
97 : {
98 2480 : iStart += iStep;
99 : }
100 :
101 2868 : poRing->addSubLineString(poLine, iStart, iEnd);
102 2868 : }
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 375 : OGRGeometryH OGRBuildPolygonFromEdges(OGRGeometryH hLines,
124 : CPL_UNUSED int bBestEffort,
125 : int bAutoClose, double dfTolerance,
126 : OGRErr *peErr)
127 :
128 : {
129 375 : if (hLines == nullptr)
130 : {
131 0 : if (peErr != nullptr)
132 0 : *peErr = OGRERR_NONE;
133 0 : return nullptr;
134 : }
135 :
136 : /* -------------------------------------------------------------------- */
137 : /* Check for the case of a geometrycollection that can be */
138 : /* promoted to MultiLineString. */
139 : /* -------------------------------------------------------------------- */
140 375 : OGRGeometry *poGeom = OGRGeometry::FromHandle(hLines);
141 375 : if (wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection)
142 : {
143 3233 : for (auto &&poMember : poGeom->toGeometryCollection())
144 : {
145 2862 : if (wkbFlatten(poMember->getGeometryType()) != wkbLineString)
146 : {
147 1 : if (peErr != nullptr)
148 1 : *peErr = OGRERR_FAILURE;
149 1 : CPLError(CE_Failure, CPLE_NotSupported,
150 : "The geometry collection contains "
151 : "non-line string geometries");
152 1 : return nullptr;
153 : }
154 : }
155 : }
156 3 : else if (wkbFlatten(poGeom->getGeometryType()) != wkbMultiLineString)
157 : {
158 1 : if (peErr != nullptr)
159 1 : *peErr = OGRERR_FAILURE;
160 1 : CPLError(CE_Failure, CPLE_NotSupported,
161 : "The passed geometry is not an OGRGeometryCollection "
162 : "(or OGRMultiLineString) "
163 : "containing line string geometries");
164 1 : return nullptr;
165 : }
166 :
167 373 : bool bSuccess = true;
168 373 : const OGRGeometryCollection *poLines = poGeom->toGeometryCollection();
169 746 : 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 373 : const int nEdges = poLines->getNumGeometries();
176 373 : std::list<const OGRLineString *> oListEdges;
177 3243 : for (int i = 0; i < nEdges; i++)
178 : {
179 : const OGRLineString *poLine =
180 2870 : poLines->getGeometryRef(i)->toLineString();
181 2870 : if (poLine->getNumPoints() >= 2)
182 : {
183 2868 : oListEdges.push_back(poLine);
184 : }
185 : }
186 :
187 : /* ==================================================================== */
188 : /* Loop generating rings. */
189 : /* ==================================================================== */
190 761 : while (!oListEdges.empty())
191 : {
192 : /* --------------------------------------------------------------------
193 : */
194 : /* Find the first unconsumed edge. */
195 : /* --------------------------------------------------------------------
196 : */
197 388 : const OGRLineString *poLine = oListEdges.front();
198 388 : oListEdges.erase(oListEdges.begin());
199 :
200 : /* --------------------------------------------------------------------
201 : */
202 : /* Start a new ring, copying in the current line directly */
203 : /* --------------------------------------------------------------------
204 : */
205 776 : auto poRing = std::make_unique<OGRLinearRing>();
206 :
207 388 : 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 388 : bool bWorkDone = true;
216 388 : double dfBestDist = dfTolerance;
217 :
218 5348 : while (!CheckPoints(poRing.get(), 0, poRing.get(),
219 2868 : poRing->getNumPoints() - 1, nullptr) &&
220 2868 : !oListEdges.empty() && bWorkDone)
221 : {
222 2480 : bool bReverse = false;
223 :
224 2480 : bWorkDone = false;
225 2480 : 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 2480 : const OGRLineString *poBestEdge = nullptr;
235 2480 : std::list<const OGRLineString *>::iterator oBestIter;
236 12179 : for (auto oIter = oListEdges.begin(); oIter != oListEdges.end();
237 9699 : ++oIter)
238 : {
239 12167 : poLine = *oIter;
240 :
241 12167 : if (CheckPoints(poLine, 0, poRing.get(),
242 12167 : poRing->getNumPoints() - 1, &dfBestDist))
243 : {
244 2130 : poBestEdge = poLine;
245 2130 : oBestIter = oIter;
246 2130 : bReverse = false;
247 : }
248 12167 : if (CheckPoints(poLine, poLine->getNumPoints() - 1,
249 12167 : poRing.get(), poRing->getNumPoints() - 1,
250 : &dfBestDist))
251 : {
252 350 : poBestEdge = poLine;
253 350 : oBestIter = oIter;
254 350 : bReverse = true;
255 : }
256 :
257 : // If we found an exact match, jump now.
258 12167 : if (dfBestDist == 0.0 && poBestEdge != nullptr)
259 2468 : break;
260 : }
261 :
262 : // We found one within tolerance - add it.
263 2480 : if (poBestEdge)
264 : {
265 2480 : AddEdgeToRing(poRing.get(), poBestEdge, bReverse, dfTolerance);
266 :
267 2480 : oListEdges.erase(oBestIter);
268 2480 : bWorkDone = true;
269 : }
270 : }
271 :
272 : /* --------------------------------------------------------------------
273 : */
274 : /* Did we fail to complete the ring? */
275 : /* --------------------------------------------------------------------
276 : */
277 388 : dfBestDist = dfTolerance;
278 :
279 388 : if (!CheckPoints(poRing.get(), 0, poRing.get(),
280 388 : poRing->getNumPoints() - 1, &dfBestDist))
281 : {
282 0 : CPLDebug("OGR",
283 : "Failed to close ring %d.\n"
284 : "End Points are: (%.8f,%.7f) and (%.7f,%.7f)",
285 0 : static_cast<int>(apoPolys.size()), poRing->getX(0),
286 0 : poRing->getY(0), poRing->getX(poRing->getNumPoints() - 1),
287 0 : poRing->getY(poRing->getNumPoints() - 1));
288 :
289 0 : bSuccess = false;
290 : }
291 :
292 : /* --------------------------------------------------------------------
293 : */
294 : /* Do we need to auto-close this ring? */
295 : /* --------------------------------------------------------------------
296 : */
297 388 : dfBestDist = dfTolerance;
298 :
299 388 : if (bAutoClose)
300 : {
301 58 : if (!CheckPoints(poRing.get(), 0, poRing.get(),
302 58 : poRing->getNumPoints() - 1, &dfBestDist))
303 : {
304 0 : poRing->addPoint(poRing->getX(0), poRing->getY(0),
305 0 : poRing->getZ(0));
306 : }
307 58 : else if (!CheckPoints(poRing.get(), 0, poRing.get(),
308 58 : poRing->getNumPoints() - 1, nullptr))
309 : {
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 15 : poRing->setPoint(poRing->getNumPoints() - 1, poRing->getX(0),
314 10 : poRing->getY(0), poRing->getZ(0));
315 : }
316 : }
317 :
318 388 : auto poPoly = std::make_unique<OGRPolygon>();
319 388 : poPoly->addRing(std::move(poRing));
320 388 : apoPolys.push_back(std::move(poPoly));
321 : } // Next ring.
322 :
323 373 : if (peErr != nullptr)
324 : {
325 372 : *peErr = bSuccess ? OGRERR_NONE : OGRERR_FAILURE;
326 : }
327 :
328 373 : return OGRGeometry::ToHandle(
329 746 : OGRGeometryFactory::organizePolygons(apoPolys).release());
330 : }
|
