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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_port.h"
31 : #include "ogr_api.h"
32 :
33 : #include <cmath>
34 : #include <cstddef>
35 : #include <list>
36 : #include <vector>
37 :
38 : #include "ogr_core.h"
39 : #include "ogr_geometry.h"
40 : #include "cpl_conv.h"
41 : #include "cpl_error.h"
42 :
43 : /************************************************************************/
44 : /* CheckPoints() */
45 : /* */
46 : /* Check if two points are closer than the current best */
47 : /* distance. Update the current best distance if they are. */
48 : /************************************************************************/
49 :
50 28692 : static bool CheckPoints(OGRLineString *poLine1, int iPoint1,
51 : OGRLineString *poLine2, int iPoint2,
52 : double *pdfDistance)
53 :
54 : {
55 28692 : if (pdfDistance == nullptr || *pdfDistance == 0)
56 : {
57 35781 : if (poLine1->getX(iPoint1) == poLine2->getX(iPoint2) &&
58 7291 : poLine1->getY(iPoint1) == poLine2->getY(iPoint2))
59 : {
60 5332 : if (pdfDistance)
61 4954 : *pdfDistance = 0.0;
62 5332 : return true;
63 : }
64 23158 : return false;
65 : }
66 :
67 : const double dfDeltaX =
68 202 : std::abs(poLine1->getX(iPoint1) - poLine2->getX(iPoint2));
69 :
70 202 : if (dfDeltaX > *pdfDistance)
71 30 : return false;
72 :
73 : const double dfDeltaY =
74 172 : std::abs(poLine1->getY(iPoint1) - poLine2->getY(iPoint2));
75 :
76 172 : if (dfDeltaY > *pdfDistance)
77 6 : return false;
78 :
79 166 : const double dfDistance = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
80 :
81 166 : if (dfDistance < *pdfDistance)
82 : {
83 166 : *pdfDistance = dfDistance;
84 166 : return true;
85 : }
86 :
87 0 : return false;
88 : }
89 :
90 : /************************************************************************/
91 : /* AddEdgeToRing() */
92 : /************************************************************************/
93 :
94 2717 : static void AddEdgeToRing(OGRLinearRing *poRing, OGRLineString *poLine,
95 : bool bReverse, double dfTolerance)
96 :
97 : {
98 : /* -------------------------------------------------------------------- */
99 : /* Establish order and range of traverse. */
100 : /* -------------------------------------------------------------------- */
101 2717 : const int nVertToAdd = poLine->getNumPoints();
102 :
103 2717 : int iStart = bReverse ? nVertToAdd - 1 : 0;
104 2717 : const int iEnd = bReverse ? 0 : nVertToAdd - 1;
105 2717 : const int iStep = bReverse ? -1 : 1;
106 :
107 : /* -------------------------------------------------------------------- */
108 : /* Should we skip a repeating vertex? */
109 : /* -------------------------------------------------------------------- */
110 5083 : if (poRing->getNumPoints() > 0 &&
111 2366 : CheckPoints(poRing, poRing->getNumPoints() - 1, poLine, iStart,
112 : &dfTolerance))
113 : {
114 2366 : iStart += iStep;
115 : }
116 :
117 2717 : poRing->addSubLineString(poLine, iStart, iEnd);
118 2717 : }
119 :
120 : /************************************************************************/
121 : /* OGRBuildPolygonFromEdges() */
122 : /************************************************************************/
123 :
124 : /**
125 : * Build a ring from a bunch of arcs.
126 : *
127 : * @param hLines handle to an OGRGeometryCollection (or OGRMultiLineString)
128 : * containing the line string geometries to be built into rings.
129 : * @param bBestEffort not yet implemented???.
130 : * @param bAutoClose indicates if the ring should be close when first and
131 : * last points of the ring are the same.
132 : * @param dfTolerance tolerance into which two arcs are considered
133 : * close enough to be joined.
134 : * @param peErr OGRERR_NONE on success, or OGRERR_FAILURE on failure.
135 : * @return a handle to the new geometry, a polygon.
136 : *
137 : */
138 :
139 343 : OGRGeometryH OGRBuildPolygonFromEdges(OGRGeometryH hLines,
140 : CPL_UNUSED int bBestEffort,
141 : int bAutoClose, double dfTolerance,
142 : OGRErr *peErr)
143 :
144 : {
145 343 : if (hLines == nullptr)
146 : {
147 0 : if (peErr != nullptr)
148 0 : *peErr = OGRERR_NONE;
149 0 : return nullptr;
150 : }
151 :
152 : /* -------------------------------------------------------------------- */
153 : /* Check for the case of a geometrycollection that can be */
154 : /* promoted to MultiLineString. */
155 : /* -------------------------------------------------------------------- */
156 343 : OGRGeometry *poGeom = OGRGeometry::FromHandle(hLines);
157 343 : if (wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection)
158 : {
159 3055 : for (auto &&poMember : poGeom->toGeometryCollection())
160 : {
161 2715 : if (wkbFlatten(poMember->getGeometryType()) != wkbLineString)
162 : {
163 1 : if (peErr != nullptr)
164 1 : *peErr = OGRERR_FAILURE;
165 1 : CPLError(CE_Failure, CPLE_NotSupported,
166 : "The geometry collection contains "
167 : "non-line string geometries");
168 1 : return nullptr;
169 : }
170 : }
171 : }
172 2 : else if (wkbFlatten(poGeom->getGeometryType()) != wkbMultiLineString)
173 : {
174 1 : if (peErr != nullptr)
175 1 : *peErr = OGRERR_FAILURE;
176 1 : CPLError(CE_Failure, CPLE_NotSupported,
177 : "The passed geometry is not an OGRGeometryCollection "
178 : "(or OGRMultiLineString) "
179 : "containing line string geometries");
180 1 : return nullptr;
181 : }
182 :
183 341 : bool bSuccess = true;
184 341 : OGRGeometryCollection *poLines = poGeom->toGeometryCollection();
185 682 : std::vector<OGRLinearRing *> apoRings;
186 :
187 : /* -------------------------------------------------------------------- */
188 : /* Setup array of line markers indicating if they have been */
189 : /* added to a ring yet. */
190 : /* -------------------------------------------------------------------- */
191 341 : const int nEdges = poLines->getNumGeometries();
192 682 : std::list<OGRLineString *> oListEdges;
193 3060 : for (int i = 0; i < nEdges; i++)
194 : {
195 2719 : OGRLineString *poLine = poLines->getGeometryRef(i)->toLineString();
196 2719 : if (poLine->getNumPoints() >= 2)
197 : {
198 2717 : oListEdges.push_back(poLine);
199 : }
200 : }
201 :
202 : /* ==================================================================== */
203 : /* Loop generating rings. */
204 : /* ==================================================================== */
205 692 : while (!oListEdges.empty())
206 : {
207 : /* --------------------------------------------------------------------
208 : */
209 : /* Find the first unconsumed edge. */
210 : /* --------------------------------------------------------------------
211 : */
212 351 : OGRLineString *poLine = oListEdges.front();
213 351 : oListEdges.erase(oListEdges.begin());
214 :
215 : /* --------------------------------------------------------------------
216 : */
217 : /* Start a new ring, copying in the current line directly */
218 : /* --------------------------------------------------------------------
219 : */
220 351 : OGRLinearRing *poRing = new OGRLinearRing();
221 :
222 351 : AddEdgeToRing(poRing, poLine, FALSE, 0);
223 :
224 : /* ====================================================================
225 : */
226 : /* Loop adding edges to this ring until we make a whole pass */
227 : /* within finding anything to add. */
228 : /* ====================================================================
229 : */
230 351 : bool bWorkDone = true;
231 351 : double dfBestDist = dfTolerance;
232 :
233 2717 : while (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
234 2371 : nullptr) &&
235 2717 : !oListEdges.empty() && bWorkDone)
236 : {
237 2366 : bool bReverse = false;
238 :
239 2366 : bWorkDone = false;
240 2366 : dfBestDist = dfTolerance;
241 :
242 : // We consider linking the end to the beginning. If this is
243 : // closer than any other option we will just close the loop.
244 :
245 : // CheckPoints(poRing, 0, poRing, poRing->getNumPoints()-1,
246 : // &dfBestDist);
247 :
248 : // Find unused edge with end point closest to our loose end.
249 2366 : OGRLineString *poBestEdge = nullptr;
250 2366 : std::list<OGRLineString *>::iterator oBestIter;
251 11604 : for (auto oIter = oListEdges.begin(); oIter != oListEdges.end();
252 9238 : ++oIter)
253 : {
254 11592 : poLine = *oIter;
255 :
256 11592 : if (CheckPoints(poLine, 0, poRing, poRing->getNumPoints() - 1,
257 : &dfBestDist))
258 : {
259 2031 : poBestEdge = poLine;
260 2031 : oBestIter = oIter;
261 2031 : bReverse = false;
262 : }
263 11592 : if (CheckPoints(poLine, poLine->getNumPoints() - 1, poRing,
264 11592 : poRing->getNumPoints() - 1, &dfBestDist))
265 : {
266 335 : poBestEdge = poLine;
267 335 : oBestIter = oIter;
268 335 : bReverse = true;
269 : }
270 :
271 : // If we found an exact match, jump now.
272 11592 : if (dfBestDist == 0.0 && poBestEdge != nullptr)
273 2354 : break;
274 : }
275 :
276 : // We found one within tolerance - add it.
277 2366 : if (poBestEdge)
278 : {
279 2366 : AddEdgeToRing(poRing, poBestEdge, bReverse, dfTolerance);
280 :
281 2366 : oListEdges.erase(oBestIter);
282 2366 : bWorkDone = true;
283 : }
284 : }
285 :
286 : /* --------------------------------------------------------------------
287 : */
288 : /* Did we fail to complete the ring? */
289 : /* --------------------------------------------------------------------
290 : */
291 351 : dfBestDist = dfTolerance;
292 :
293 351 : if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
294 : &dfBestDist))
295 : {
296 0 : CPLDebug("OGR",
297 : "Failed to close ring %d.\n"
298 : "End Points are: (%.8f,%.7f) and (%.7f,%.7f)",
299 0 : static_cast<int>(apoRings.size()), poRing->getX(0),
300 0 : poRing->getY(0), poRing->getX(poRing->getNumPoints() - 1),
301 0 : poRing->getY(poRing->getNumPoints() - 1));
302 :
303 0 : bSuccess = false;
304 : }
305 :
306 : /* --------------------------------------------------------------------
307 : */
308 : /* Do we need to auto-close this ring? */
309 : /* --------------------------------------------------------------------
310 : */
311 351 : dfBestDist = dfTolerance;
312 :
313 351 : if (bAutoClose)
314 : {
315 37 : if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
316 : &dfBestDist))
317 : {
318 0 : poRing->addPoint(poRing->getX(0), poRing->getY(0),
319 : poRing->getZ(0));
320 : }
321 37 : else if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
322 : nullptr))
323 : {
324 : // The endpoints are very close but do not exactly match.
325 : // Alter the last one so it is equal to the first, to prevent
326 : // invalid self-intersecting rings.
327 5 : poRing->setPoint(poRing->getNumPoints() - 1, poRing->getX(0),
328 : poRing->getY(0), poRing->getZ(0));
329 : }
330 : }
331 :
332 351 : apoRings.push_back(poRing);
333 : } // Next ring.
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Identify exterior ring - it will be the largest. #3610 */
337 : /* -------------------------------------------------------------------- */
338 341 : double maxarea = 0.0;
339 341 : int maxring = -1;
340 341 : OGREnvelope tenv;
341 :
342 692 : for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn)
343 : {
344 351 : apoRings[rn]->getEnvelope(&tenv);
345 351 : const double tarea = (tenv.MaxX - tenv.MinX) * (tenv.MaxY - tenv.MinY);
346 351 : if (tarea > maxarea)
347 : {
348 342 : maxarea = tarea;
349 342 : maxring = rn;
350 : }
351 : }
352 :
353 341 : OGRPolygon *poPolygon = new OGRPolygon();
354 :
355 341 : if (maxring != -1)
356 : {
357 341 : poPolygon->addRingDirectly(apoRings[maxring]);
358 692 : for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn)
359 : {
360 351 : if (rn == maxring)
361 341 : continue;
362 10 : poPolygon->addRingDirectly(apoRings[rn]);
363 : }
364 : }
365 : else
366 : {
367 0 : for (auto &poRing : apoRings)
368 0 : delete poRing;
369 : }
370 :
371 341 : if (peErr != nullptr)
372 : {
373 341 : *peErr = bSuccess ? OGRERR_NONE : OGRERR_FAILURE;
374 : }
375 :
376 341 : return OGRGeometry::ToHandle(poPolygon);
377 : }
|
