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 224 : GDALRasterPolygonEnumeratorT<
31 : DataType, EqualityTest>::GDALRasterPolygonEnumeratorT(int nConnectednessIn)
32 224 : : nConnectedness(nConnectednessIn)
33 :
34 : {
35 224 : CPLAssert(nConnectedness == 4 || nConnectedness == 8);
36 224 : }
37 :
38 : /************************************************************************/
39 : /* ~GDALRasterPolygonEnumeratorT() */
40 : /************************************************************************/
41 :
42 : template <class DataType, class EqualityTest>
43 224 : GDALRasterPolygonEnumeratorT<DataType,
44 : EqualityTest>::~GDALRasterPolygonEnumeratorT()
45 :
46 : {
47 224 : Clear();
48 224 : }
49 :
50 : /************************************************************************/
51 : /* Clear() */
52 : /************************************************************************/
53 :
54 : template <class DataType, class EqualityTest>
55 243 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::Clear()
56 :
57 : {
58 243 : CPLFree(panPolyIdMap);
59 243 : CPLFree(panPolyValue);
60 :
61 243 : panPolyIdMap = nullptr;
62 243 : panPolyValue = nullptr;
63 :
64 243 : nNextPolygonId = 0;
65 243 : nPolyAlloc = 0;
66 243 : }
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 23879 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::MergePolygon(
76 : int nSrcId, int nDstIdInit)
77 :
78 : {
79 : // Figure out the final dest id.
80 23879 : int nDstIdFinal = nDstIdInit;
81 23912 : while (panPolyIdMap[nDstIdFinal] != nDstIdFinal)
82 33 : nDstIdFinal = panPolyIdMap[nDstIdFinal];
83 :
84 : // Map the whole intermediate chain to it.
85 23879 : int nDstIdCur = nDstIdInit;
86 23912 : 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 23883 : while (panPolyIdMap[nSrcId] != nSrcId)
95 : {
96 4 : int nNextSrcId = panPolyIdMap[nSrcId];
97 4 : panPolyIdMap[nSrcId] = nDstIdFinal;
98 4 : nSrcId = nNextSrcId;
99 : }
100 23879 : panPolyIdMap[nSrcId] = nDstIdFinal;
101 23879 : }
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 48098 : int GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::NewPolygon(
112 : DataType nValue)
113 :
114 : {
115 48098 : 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 48098 : if (nNextPolygonId >= nPolyAlloc)
123 : {
124 : int nPolyAllocNew;
125 393 : if (nPolyAlloc < (std::numeric_limits<int>::max() - 20) / 2)
126 393 : 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 393 : VSI_REALLOC_VERBOSE(panPolyIdMap, nPolyAllocNew * sizeof(GInt32)));
145 393 : auto panPolyValueNew = static_cast<DataType *>(VSI_REALLOC_VERBOSE(
146 : panPolyValue, nPolyAllocNew * sizeof(DataType)));
147 393 : if (panPolyIdMapNew == nullptr || panPolyValueNew == nullptr)
148 : {
149 0 : VSIFree(panPolyIdMapNew);
150 0 : VSIFree(panPolyValueNew);
151 0 : return -1;
152 : }
153 393 : panPolyIdMap = panPolyIdMapNew;
154 393 : panPolyValue = panPolyValueNew;
155 393 : nPolyAlloc = nPolyAllocNew;
156 : }
157 :
158 48098 : const int nPolyId = nNextPolygonId;
159 48098 : panPolyIdMap[nPolyId] = nPolyId;
160 48098 : panPolyValue[nPolyId] = nValue;
161 48098 : nNextPolygonId++;
162 :
163 48098 : 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 113 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::CompleteMerges()
176 :
177 : {
178 113 : int nFinalPolyCount = 0;
179 :
180 20493 : for (int iPoly = 0; iPoly < nNextPolygonId; iPoly++)
181 : {
182 : // Figure out the final id.
183 20380 : int nId = panPolyIdMap[iPoly];
184 29194 : while (nId != panPolyIdMap[nId])
185 : {
186 8814 : nId = panPolyIdMap[nId];
187 : }
188 :
189 : // Then map the whole intermediate chain to it.
190 20380 : int nIdCur = panPolyIdMap[iPoly];
191 20380 : panPolyIdMap[iPoly] = nId;
192 29194 : while (nIdCur != panPolyIdMap[nIdCur])
193 : {
194 8814 : int nNextId = panPolyIdMap[nIdCur];
195 8814 : panPolyIdMap[nIdCur] = nId;
196 8814 : nIdCur = nNextId;
197 : }
198 :
199 20380 : if (panPolyIdMap[iPoly] == iPoly)
200 10194 : nFinalPolyCount++;
201 : }
202 :
203 113 : CPLDebug("GDALRasterPolygonEnumerator",
204 : "Counted %d polygon fragments forming %d final polygons.",
205 : nNextPolygonId, nFinalPolyCount);
206 113 : }
207 :
208 : /************************************************************************/
209 : /* ProcessLine() */
210 : /* */
211 : /* Assign ids to polygons, one line at a time. */
212 : /************************************************************************/
213 :
214 : template <class DataType, class EqualityTest>
215 4486 : 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 4486 : if (panLastLineVal == nullptr)
226 : {
227 4006 : for (int i = 0; i < nXSize; i++)
228 : {
229 3763 : if (panThisLineVal[i] == GP_NODATA_MARKER)
230 : {
231 1154 : panThisLineId[i] = -1;
232 : }
233 5009 : else if (i == 0 ||
234 2400 : !(eq.operator()(panThisLineVal[i], panThisLineVal[i - 1])))
235 : {
236 1186 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
237 1186 : if (panThisLineId[i] < 0)
238 0 : return false;
239 : }
240 : else
241 : {
242 1423 : panThisLineId[i] = panThisLineId[i - 1];
243 : }
244 : }
245 :
246 243 : return true;
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Process each pixel comparing to the previous pixel, and to */
251 : /* the last line. */
252 : /* -------------------------------------------------------------------- */
253 893044 : for (int i = 0; i < nXSize; i++)
254 : {
255 888801 : if (panThisLineVal[i] == GP_NODATA_MARKER)
256 : {
257 105014 : panThisLineId[i] = -1;
258 : }
259 1565086 : else if (i > 0 &&
260 781296 : eq.operator()(panThisLineVal[i], panThisLineVal[i - 1]))
261 : {
262 698718 : panThisLineId[i] = panThisLineId[i - 1];
263 :
264 1351950 : if (eq.operator()(panLastLineVal[i], panThisLineVal[i]) &&
265 653227 : (panPolyIdMap[panLastLineId[i]] !=
266 653227 : panPolyIdMap[panThisLineId[i]]))
267 : {
268 23817 : MergePolygon(panLastLineId[i], panThisLineId[i]);
269 : }
270 :
271 347 : if (nConnectedness == 8 &&
272 699065 : 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 699362 : 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 85069 : else if (eq.operator()(panLastLineVal[i], panThisLineVal[i]))
288 : {
289 37923 : panThisLineId[i] = panLastLineId[i];
290 : }
291 48401 : 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 48241 : 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 46912 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
312 46912 : if (panThisLineId[i] < 0)
313 0 : return false;
314 : }
315 : }
316 4243 : return true;
317 : }
318 :
319 : template class GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest>;
320 :
321 : template class GDALRasterPolygonEnumeratorT<float, FloatEqualityTest>;
322 :
323 : /*! @endcond */
|
