Coverage Report

Created: 2025-06-09 08:44

/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
}