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 220 : GDALRasterPolygonEnumeratorT<
31 : DataType, EqualityTest>::GDALRasterPolygonEnumeratorT(int nConnectednessIn)
32 220 : : nConnectedness(nConnectednessIn)
33 :
34 : {
35 220 : CPLAssert(nConnectedness == 4 || nConnectedness == 8);
36 220 : }
37 :
38 : /************************************************************************/
39 : /* ~GDALRasterPolygonEnumeratorT() */
40 : /************************************************************************/
41 :
42 : template <class DataType, class EqualityTest>
43 220 : GDALRasterPolygonEnumeratorT<DataType,
44 : EqualityTest>::~GDALRasterPolygonEnumeratorT()
45 :
46 : {
47 220 : Clear();
48 220 : }
49 :
50 : /************************************************************************/
51 : /* Clear() */
52 : /************************************************************************/
53 :
54 : template <class DataType, class EqualityTest>
55 239 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::Clear()
56 :
57 : {
58 239 : CPLFree(panPolyIdMap);
59 239 : CPLFree(panPolyValue);
60 :
61 239 : panPolyIdMap = nullptr;
62 239 : panPolyValue = nullptr;
63 :
64 239 : nNextPolygonId = 0;
65 239 : nPolyAlloc = 0;
66 239 : }
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 23859 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::MergePolygon(
76 : int nSrcId, int nDstIdInit)
77 :
78 : {
79 : // Figure out the final dest id.
80 23859 : int nDstIdFinal = nDstIdInit;
81 23892 : while (panPolyIdMap[nDstIdFinal] != nDstIdFinal)
82 33 : nDstIdFinal = panPolyIdMap[nDstIdFinal];
83 :
84 : // Map the whole intermediate chain to it.
85 23859 : int nDstIdCur = nDstIdInit;
86 23892 : while (panPolyIdMap[nDstIdCur] != nDstIdCur)
87 : {
88 33 : int nNextDstId = panPolyIdMap[nDstIdCur];
89 33 : panPolyIdMap[nDstIdCur] = nDstIdFinal;
90 33 : nDstIdCur = nNextDstId;
91 : }
92 :
93 : // And map the whole source chain to it too (can be done in one pass).
94 23863 : while (panPolyIdMap[nSrcId] != nSrcId)
95 : {
96 4 : int nNextSrcId = panPolyIdMap[nSrcId];
97 4 : panPolyIdMap[nSrcId] = nDstIdFinal;
98 4 : nSrcId = nNextSrcId;
99 : }
100 23859 : panPolyIdMap[nSrcId] = nDstIdFinal;
101 23859 : }
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 47514 : int GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::NewPolygon(
112 : DataType nValue)
113 :
114 : {
115 47514 : 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 47514 : if (nNextPolygonId >= nPolyAlloc)
123 : {
124 : int nPolyAllocNew;
125 383 : if (nPolyAlloc < (std::numeric_limits<int>::max() - 20) / 2)
126 383 : 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 383 : VSI_REALLOC_VERBOSE(panPolyIdMap, nPolyAllocNew * sizeof(GInt32)));
145 383 : auto panPolyValueNew = static_cast<DataType *>(VSI_REALLOC_VERBOSE(
146 : panPolyValue, nPolyAllocNew * sizeof(DataType)));
147 383 : if (panPolyIdMapNew == nullptr || panPolyValueNew == nullptr)
148 : {
149 0 : VSIFree(panPolyIdMapNew);
150 0 : VSIFree(panPolyValueNew);
151 0 : return -1;
152 : }
153 383 : panPolyIdMap = panPolyIdMapNew;
154 383 : panPolyValue = panPolyValueNew;
155 383 : nPolyAlloc = nPolyAllocNew;
156 : }
157 :
158 47514 : const int nPolyId = nNextPolygonId;
159 47514 : panPolyIdMap[nPolyId] = nPolyId;
160 47514 : panPolyValue[nPolyId] = nValue;
161 47514 : nNextPolygonId++;
162 :
163 47514 : 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 111 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::CompleteMerges()
176 :
177 : {
178 111 : int nFinalPolyCount = 0;
179 :
180 20199 : for (int iPoly = 0; iPoly < nNextPolygonId; iPoly++)
181 : {
182 : // Figure out the final id.
183 20088 : int nId = panPolyIdMap[iPoly];
184 28902 : while (nId != panPolyIdMap[nId])
185 : {
186 8814 : nId = panPolyIdMap[nId];
187 : }
188 :
189 : // Then map the whole intermediate chain to it.
190 20088 : int nIdCur = panPolyIdMap[iPoly];
191 20088 : panPolyIdMap[iPoly] = nId;
192 28902 : while (nIdCur != panPolyIdMap[nIdCur])
193 : {
194 8814 : int nNextId = panPolyIdMap[nIdCur];
195 8814 : panPolyIdMap[nIdCur] = nId;
196 8814 : nIdCur = nNextId;
197 : }
198 :
199 20088 : if (panPolyIdMap[iPoly] == iPoly)
200 9912 : nFinalPolyCount++;
201 : }
202 :
203 111 : CPLDebug("GDALRasterPolygonEnumerator",
204 : "Counted %d polygon fragments forming %d final polygons.",
205 : nNextPolygonId, nFinalPolyCount);
206 111 : }
207 :
208 : /************************************************************************/
209 : /* ProcessLine() */
210 : /* */
211 : /* Assign ids to polygons, one line at a time. */
212 : /************************************************************************/
213 :
214 : template <class DataType, class EqualityTest>
215 4406 : 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 4406 : if (panLastLineVal == nullptr)
226 : {
227 3922 : for (int i = 0; i < nXSize; i++)
228 : {
229 3683 : if (panThisLineVal[i] == GP_NODATA_MARKER)
230 : {
231 1154 : panThisLineId[i] = -1;
232 : }
233 4853 : else if (i == 0 ||
234 2324 : !(eq.operator()(panThisLineVal[i], panThisLineVal[i - 1])))
235 : {
236 1152 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
237 1152 : if (panThisLineId[i] < 0)
238 0 : return false;
239 : }
240 : else
241 : {
242 1377 : panThisLineId[i] = panThisLineId[i - 1];
243 : }
244 : }
245 :
246 239 : return true;
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Process each pixel comparing to the previous pixel, and to */
251 : /* the last line. */
252 : /* -------------------------------------------------------------------- */
253 891448 : for (int i = 0; i < nXSize; i++)
254 : {
255 887281 : if (panThisLineVal[i] == GP_NODATA_MARKER)
256 : {
257 105014 : panThisLineId[i] = -1;
258 : }
259 1562116 : else if (i > 0 &&
260 779852 : eq.operator()(panThisLineVal[i], panThisLineVal[i - 1]))
261 : {
262 697888 : panThisLineId[i] = panThisLineId[i - 1];
263 :
264 1350370 : if (eq.operator()(panLastLineVal[i], panThisLineVal[i]) &&
265 652483 : (panPolyIdMap[panLastLineId[i]] !=
266 652483 : panPolyIdMap[panThisLineId[i]]))
267 : {
268 23797 : MergePolygon(panLastLineId[i], panThisLineId[i]);
269 : }
270 :
271 347 : if (nConnectedness == 8 &&
272 698235 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]) &&
273 38 : (panPolyIdMap[panLastLineId[i - 1]] !=
274 38 : panPolyIdMap[panThisLineId[i]]))
275 : {
276 0 : MergePolygon(panLastLineId[i - 1], panThisLineId[i]);
277 : }
278 :
279 644 : if (nConnectedness == 8 && i < nXSize - 1 &&
280 698532 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
281 74 : (panPolyIdMap[panLastLineId[i + 1]] !=
282 74 : panPolyIdMap[panThisLineId[i]]))
283 : {
284 50 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
285 : }
286 : }
287 84379 : else if (eq.operator()(panLastLineVal[i], panThisLineVal[i]))
288 : {
289 37783 : panThisLineId[i] = panLastLineId[i];
290 : }
291 47851 : else if (i > 0 && nConnectedness == 8 &&
292 1255 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]))
293 : {
294 98 : panThisLineId[i] = panLastLineId[i - 1];
295 :
296 88 : if (i < nXSize - 1 &&
297 186 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
298 12 : (panPolyIdMap[panLastLineId[i + 1]] !=
299 12 : panPolyIdMap[panThisLineId[i]]))
300 : {
301 12 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
302 : }
303 : }
304 47691 : else if (i < nXSize - 1 && nConnectedness == 8 &&
305 1193 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]))
306 : {
307 136 : panThisLineId[i] = panLastLineId[i + 1];
308 : }
309 : else
310 : {
311 46362 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
312 46362 : if (panThisLineId[i] < 0)
313 0 : return false;
314 : }
315 : }
316 4167 : return true;
317 : }
318 :
319 : template class GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest>;
320 :
321 : template class GDALRasterPolygonEnumeratorT<float, FloatEqualityTest>;
322 :
323 : /*! @endcond */
|
