Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Raster Polygon Enumerator
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2008, Frank Warmerdam
9 : * Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_alg_priv.h"
16 :
17 : #include <cstddef>
18 : #include <limits>
19 :
20 : #include "cpl_conv.h"
21 : #include "cpl_error.h"
22 :
23 : /*! @cond Doxygen_Suppress */
24 :
25 : /************************************************************************/
26 : /* GDALRasterPolygonEnumeratorT() */
27 : /************************************************************************/
28 :
29 : template <class DataType, class EqualityTest>
30 126 : GDALRasterPolygonEnumeratorT<
31 : DataType, EqualityTest>::GDALRasterPolygonEnumeratorT(int nConnectednessIn)
32 126 : : nConnectedness(nConnectednessIn)
33 :
34 : {
35 126 : CPLAssert(nConnectedness == 4 || nConnectedness == 8);
36 126 : }
37 :
38 : /************************************************************************/
39 : /* ~GDALRasterPolygonEnumeratorT() */
40 : /************************************************************************/
41 :
42 : template <class DataType, class EqualityTest>
43 126 : GDALRasterPolygonEnumeratorT<DataType,
44 : EqualityTest>::~GDALRasterPolygonEnumeratorT()
45 :
46 : {
47 126 : Clear();
48 126 : }
49 :
50 : /************************************************************************/
51 : /* Clear() */
52 : /************************************************************************/
53 :
54 : template <class DataType, class EqualityTest>
55 136 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::Clear()
56 :
57 : {
58 136 : CPLFree(panPolyIdMap);
59 136 : CPLFree(panPolyValue);
60 :
61 136 : panPolyIdMap = nullptr;
62 136 : panPolyValue = nullptr;
63 :
64 136 : nNextPolygonId = 0;
65 136 : nPolyAlloc = 0;
66 136 : }
67 :
68 : /************************************************************************/
69 : /* MergePolygon() */
70 : /* */
71 : /* Update the polygon map to indicate the merger of two polygons. */
72 : /************************************************************************/
73 :
74 : template <class DataType, class EqualityTest>
75 23636 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::MergePolygon(
76 : int nSrcId, int nDstIdInit)
77 :
78 : {
79 : // Figure out the final dest id.
80 23636 : int nDstIdFinal = nDstIdInit;
81 23667 : while (panPolyIdMap[nDstIdFinal] != nDstIdFinal)
82 31 : nDstIdFinal = panPolyIdMap[nDstIdFinal];
83 :
84 : // Map the whole intermediate chain to it.
85 23636 : int nDstIdCur = nDstIdInit;
86 23667 : while (panPolyIdMap[nDstIdCur] != nDstIdCur)
87 : {
88 31 : int nNextDstId = panPolyIdMap[nDstIdCur];
89 31 : panPolyIdMap[nDstIdCur] = nDstIdFinal;
90 31 : nDstIdCur = nNextDstId;
91 : }
92 :
93 : // And map the whole source chain to it too (can be done in one pass).
94 23640 : while (panPolyIdMap[nSrcId] != nSrcId)
95 : {
96 4 : int nNextSrcId = panPolyIdMap[nSrcId];
97 4 : panPolyIdMap[nSrcId] = nDstIdFinal;
98 4 : nSrcId = nNextSrcId;
99 : }
100 23636 : panPolyIdMap[nSrcId] = nDstIdFinal;
101 23636 : }
102 :
103 : /************************************************************************/
104 : /* NewPolygon() */
105 : /* */
106 : /* Allocate a new polygon id, and reallocate the polygon maps */
107 : /* if needed. */
108 : /************************************************************************/
109 :
110 : template <class DataType, class EqualityTest>
111 41869 : int GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::NewPolygon(
112 : DataType nValue)
113 :
114 : {
115 41869 : if (nNextPolygonId == std::numeric_limits<int>::max())
116 : {
117 0 : CPLError(CE_Failure, CPLE_AppDefined,
118 : "GDALRasterPolygonEnumeratorT::NewPolygon(): maximum number "
119 : "of polygons reached");
120 0 : return -1;
121 : }
122 41869 : if (nNextPolygonId >= nPolyAlloc)
123 : {
124 : int nPolyAllocNew;
125 230 : if (nPolyAlloc < (std::numeric_limits<int>::max() - 20) / 2)
126 230 : nPolyAllocNew = nPolyAlloc * 2 + 20;
127 : else
128 0 : nPolyAllocNew = std::numeric_limits<int>::max();
129 : #if SIZEOF_VOIDP == 4
130 : if (nPolyAllocNew >
131 : static_cast<int>(std::numeric_limits<size_t>::max() /
132 : sizeof(GInt32)) ||
133 : nPolyAllocNew >
134 : static_cast<int>(std::numeric_limits<size_t>::max() /
135 : sizeof(DataType)))
136 : {
137 : CPLError(CE_Failure, CPLE_OutOfMemory,
138 : "GDALRasterPolygonEnumeratorT::NewPolygon(): too many "
139 : "polygons");
140 : return -1;
141 : }
142 : #endif
143 : auto panPolyIdMapNew = static_cast<GInt32 *>(
144 230 : VSI_REALLOC_VERBOSE(panPolyIdMap, nPolyAllocNew * sizeof(GInt32)));
145 230 : auto panPolyValueNew = static_cast<DataType *>(VSI_REALLOC_VERBOSE(
146 : panPolyValue, nPolyAllocNew * sizeof(DataType)));
147 230 : if (panPolyIdMapNew == nullptr || panPolyValueNew == nullptr)
148 : {
149 0 : VSIFree(panPolyIdMapNew);
150 0 : VSIFree(panPolyValueNew);
151 0 : return -1;
152 : }
153 230 : panPolyIdMap = panPolyIdMapNew;
154 230 : panPolyValue = panPolyValueNew;
155 230 : nPolyAlloc = nPolyAllocNew;
156 : }
157 :
158 41869 : const int nPolyId = nNextPolygonId;
159 41869 : panPolyIdMap[nPolyId] = nPolyId;
160 41869 : panPolyValue[nPolyId] = nValue;
161 41869 : nNextPolygonId++;
162 :
163 41869 : return nPolyId;
164 : }
165 :
166 : /************************************************************************/
167 : /* CompleteMerges() */
168 : /* */
169 : /* Make a pass through the maps, ensuring every polygon id */
170 : /* points to the final id it should use, not an intermediate */
171 : /* value. */
172 : /************************************************************************/
173 :
174 : template <class DataType, class EqualityTest>
175 64 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::CompleteMerges()
176 :
177 : {
178 64 : int nFinalPolyCount = 0;
179 :
180 17404 : for (int iPoly = 0; iPoly < nNextPolygonId; iPoly++)
181 : {
182 : // Figure out the final id.
183 17340 : int nId = panPolyIdMap[iPoly];
184 26153 : while (nId != panPolyIdMap[nId])
185 : {
186 8813 : nId = panPolyIdMap[nId];
187 : }
188 :
189 : // Then map the whole intermediate chain to it.
190 17340 : int nIdCur = panPolyIdMap[iPoly];
191 17340 : panPolyIdMap[iPoly] = nId;
192 26153 : while (nIdCur != panPolyIdMap[nIdCur])
193 : {
194 8813 : int nNextId = panPolyIdMap[nIdCur];
195 8813 : panPolyIdMap[nIdCur] = nId;
196 8813 : nIdCur = nNextId;
197 : }
198 :
199 17340 : if (panPolyIdMap[iPoly] == iPoly)
200 7269 : nFinalPolyCount++;
201 : }
202 :
203 64 : CPLDebug("GDALRasterPolygonEnumerator",
204 : "Counted %d polygon fragments forming %d final polygons.",
205 : nNextPolygonId, nFinalPolyCount);
206 64 : }
207 :
208 : /************************************************************************/
209 : /* ProcessLine() */
210 : /* */
211 : /* Assign ids to polygons, one line at a time. */
212 : /************************************************************************/
213 :
214 : template <class DataType, class EqualityTest>
215 3433 : bool GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::ProcessLine(
216 : DataType *panLastLineVal, DataType *panThisLineVal, GInt32 *panLastLineId,
217 : GInt32 *panThisLineId, int nXSize)
218 :
219 : {
220 : EqualityTest eq;
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Special case for the first line. */
224 : /* -------------------------------------------------------------------- */
225 3433 : if (panLastLineVal == nullptr)
226 : {
227 2864 : for (int i = 0; i < nXSize; i++)
228 : {
229 2728 : if (panThisLineVal[i] == GP_NODATA_MARKER)
230 : {
231 1089 : panThisLineId[i] = -1;
232 : }
233 3167 : else if (i == 0 ||
234 1528 : !(eq.operator()(panThisLineVal[i], panThisLineVal[i - 1])))
235 : {
236 690 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
237 690 : if (panThisLineId[i] < 0)
238 0 : return false;
239 : }
240 : else
241 : {
242 949 : panThisLineId[i] = panThisLineId[i - 1];
243 : }
244 : }
245 :
246 136 : return true;
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Process each pixel comparing to the previous pixel, and to */
251 : /* the last line. */
252 : /* -------------------------------------------------------------------- */
253 875876 : for (int i = 0; i < nXSize; i++)
254 : {
255 872579 : if (panThisLineVal[i] == GP_NODATA_MARKER)
256 : {
257 104851 : panThisLineId[i] = -1;
258 : }
259 1533864 : else if (i > 0 &&
260 766141 : eq.operator()(panThisLineVal[i], panThisLineVal[i - 1]))
261 : {
262 690132 : panThisLineId[i] = panThisLineId[i - 1];
263 :
264 1335830 : if (eq.operator()(panLastLineVal[i], panThisLineVal[i]) &&
265 645697 : (panPolyIdMap[panLastLineId[i]] !=
266 645697 : panPolyIdMap[panThisLineId[i]]))
267 : {
268 23607 : MergePolygon(panLastLineId[i], panThisLineId[i]);
269 : }
270 :
271 158 : if (nConnectedness == 8 &&
272 690290 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]) &&
273 17 : (panPolyIdMap[panLastLineId[i - 1]] !=
274 17 : panPolyIdMap[panThisLineId[i]]))
275 : {
276 0 : MergePolygon(panLastLineId[i - 1], panThisLineId[i]);
277 : }
278 :
279 297 : if (nConnectedness == 8 && i < nXSize - 1 &&
280 690429 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
281 34 : (panPolyIdMap[panLastLineId[i + 1]] !=
282 34 : panPolyIdMap[panThisLineId[i]]))
283 : {
284 23 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
285 : }
286 : }
287 77596 : else if (eq.operator()(panLastLineVal[i], panThisLineVal[i]))
288 : {
289 36301 : panThisLineId[i] = panLastLineId[i];
290 : }
291 41937 : else if (i > 0 && nConnectedness == 8 &&
292 642 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]))
293 : {
294 49 : panThisLineId[i] = panLastLineId[i - 1];
295 :
296 43 : if (i < nXSize - 1 &&
297 92 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
298 6 : (panPolyIdMap[panLastLineId[i + 1]] !=
299 6 : panPolyIdMap[panThisLineId[i]]))
300 : {
301 6 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
302 : }
303 : }
304 41856 : else if (i < nXSize - 1 && nConnectedness == 8 &&
305 610 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]))
306 : {
307 67 : panThisLineId[i] = panLastLineId[i + 1];
308 : }
309 : else
310 : {
311 41179 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
312 41179 : if (panThisLineId[i] < 0)
313 0 : return false;
314 : }
315 : }
316 3297 : return true;
317 : }
318 :
319 : template class GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest>;
320 :
321 : template class GDALRasterPolygonEnumeratorT<float, FloatEqualityTest>;
322 :
323 : /*! @endcond */
|
