Coverage Report

Created: 2026-03-30 09:00

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/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
}