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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_port.h"
31 : #include "gdal_alg_priv.h"
32 :
33 : #include <cstddef>
34 : #include <limits>
35 :
36 : #include "cpl_conv.h"
37 : #include "cpl_error.h"
38 :
39 : /*! @cond Doxygen_Suppress */
40 :
41 : /************************************************************************/
42 : /* GDALRasterPolygonEnumeratorT() */
43 : /************************************************************************/
44 :
45 : template <class DataType, class EqualityTest>
46 124 : GDALRasterPolygonEnumeratorT<
47 : DataType, EqualityTest>::GDALRasterPolygonEnumeratorT(int nConnectednessIn)
48 124 : : nConnectedness(nConnectednessIn)
49 :
50 : {
51 124 : CPLAssert(nConnectedness == 4 || nConnectedness == 8);
52 124 : }
53 :
54 : /************************************************************************/
55 : /* ~GDALRasterPolygonEnumeratorT() */
56 : /************************************************************************/
57 :
58 : template <class DataType, class EqualityTest>
59 124 : GDALRasterPolygonEnumeratorT<DataType,
60 : EqualityTest>::~GDALRasterPolygonEnumeratorT()
61 :
62 : {
63 124 : Clear();
64 124 : }
65 :
66 : /************************************************************************/
67 : /* Clear() */
68 : /************************************************************************/
69 :
70 : template <class DataType, class EqualityTest>
71 134 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::Clear()
72 :
73 : {
74 134 : CPLFree(panPolyIdMap);
75 134 : CPLFree(panPolyValue);
76 :
77 134 : panPolyIdMap = nullptr;
78 134 : panPolyValue = nullptr;
79 :
80 134 : nNextPolygonId = 0;
81 134 : nPolyAlloc = 0;
82 134 : }
83 :
84 : /************************************************************************/
85 : /* MergePolygon() */
86 : /* */
87 : /* Update the polygon map to indicate the merger of two polygons. */
88 : /************************************************************************/
89 :
90 : template <class DataType, class EqualityTest>
91 23636 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::MergePolygon(
92 : int nSrcId, int nDstIdInit)
93 :
94 : {
95 : // Figure out the final dest id.
96 23636 : int nDstIdFinal = nDstIdInit;
97 23667 : while (panPolyIdMap[nDstIdFinal] != nDstIdFinal)
98 31 : nDstIdFinal = panPolyIdMap[nDstIdFinal];
99 :
100 : // Map the whole intermediate chain to it.
101 23636 : int nDstIdCur = nDstIdInit;
102 23667 : while (panPolyIdMap[nDstIdCur] != nDstIdCur)
103 : {
104 31 : int nNextDstId = panPolyIdMap[nDstIdCur];
105 31 : panPolyIdMap[nDstIdCur] = nDstIdFinal;
106 31 : nDstIdCur = nNextDstId;
107 : }
108 :
109 : // And map the whole source chain to it too (can be done in one pass).
110 23640 : while (panPolyIdMap[nSrcId] != nSrcId)
111 : {
112 4 : int nNextSrcId = panPolyIdMap[nSrcId];
113 4 : panPolyIdMap[nSrcId] = nDstIdFinal;
114 4 : nSrcId = nNextSrcId;
115 : }
116 23636 : panPolyIdMap[nSrcId] = nDstIdFinal;
117 23636 : }
118 :
119 : /************************************************************************/
120 : /* NewPolygon() */
121 : /* */
122 : /* Allocate a new polygon id, and reallocate the polygon maps */
123 : /* if needed. */
124 : /************************************************************************/
125 :
126 : template <class DataType, class EqualityTest>
127 41869 : int GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::NewPolygon(
128 : DataType nValue)
129 :
130 : {
131 41869 : if (nNextPolygonId == std::numeric_limits<int>::max())
132 : {
133 0 : CPLError(CE_Failure, CPLE_AppDefined,
134 : "GDALRasterPolygonEnumeratorT::NewPolygon(): maximum number "
135 : "of polygons reached");
136 0 : return -1;
137 : }
138 41869 : if (nNextPolygonId >= nPolyAlloc)
139 : {
140 : int nPolyAllocNew;
141 230 : if (nPolyAlloc < (std::numeric_limits<int>::max() - 20) / 2)
142 230 : nPolyAllocNew = nPolyAlloc * 2 + 20;
143 : else
144 0 : nPolyAllocNew = std::numeric_limits<int>::max();
145 : #if SIZEOF_VOIDP == 4
146 : if (nPolyAllocNew >
147 : static_cast<int>(std::numeric_limits<size_t>::max() /
148 : sizeof(GInt32)) ||
149 : nPolyAllocNew >
150 : static_cast<int>(std::numeric_limits<size_t>::max() /
151 : sizeof(DataType)))
152 : {
153 : CPLError(CE_Failure, CPLE_OutOfMemory,
154 : "GDALRasterPolygonEnumeratorT::NewPolygon(): too many "
155 : "polygons");
156 : return -1;
157 : }
158 : #endif
159 : auto panPolyIdMapNew = static_cast<GInt32 *>(
160 230 : VSI_REALLOC_VERBOSE(panPolyIdMap, nPolyAllocNew * sizeof(GInt32)));
161 230 : auto panPolyValueNew = static_cast<DataType *>(VSI_REALLOC_VERBOSE(
162 : panPolyValue, nPolyAllocNew * sizeof(DataType)));
163 230 : if (panPolyIdMapNew == nullptr || panPolyValueNew == nullptr)
164 : {
165 0 : VSIFree(panPolyIdMapNew);
166 0 : VSIFree(panPolyValueNew);
167 0 : return -1;
168 : }
169 230 : panPolyIdMap = panPolyIdMapNew;
170 230 : panPolyValue = panPolyValueNew;
171 230 : nPolyAlloc = nPolyAllocNew;
172 : }
173 :
174 41869 : const int nPolyId = nNextPolygonId;
175 41869 : panPolyIdMap[nPolyId] = nPolyId;
176 41869 : panPolyValue[nPolyId] = nValue;
177 41869 : nNextPolygonId++;
178 :
179 41869 : return nPolyId;
180 : }
181 :
182 : /************************************************************************/
183 : /* CompleteMerges() */
184 : /* */
185 : /* Make a pass through the maps, ensuring every polygon id */
186 : /* points to the final id it should use, not an intermediate */
187 : /* value. */
188 : /************************************************************************/
189 :
190 : template <class DataType, class EqualityTest>
191 62 : void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::CompleteMerges()
192 :
193 : {
194 62 : int nFinalPolyCount = 0;
195 :
196 17402 : for (int iPoly = 0; iPoly < nNextPolygonId; iPoly++)
197 : {
198 : // Figure out the final id.
199 17340 : int nId = panPolyIdMap[iPoly];
200 26153 : while (nId != panPolyIdMap[nId])
201 : {
202 8813 : nId = panPolyIdMap[nId];
203 : }
204 :
205 : // Then map the whole intermediate chain to it.
206 17340 : int nIdCur = panPolyIdMap[iPoly];
207 17340 : panPolyIdMap[iPoly] = nId;
208 26153 : while (nIdCur != panPolyIdMap[nIdCur])
209 : {
210 8813 : int nNextId = panPolyIdMap[nIdCur];
211 8813 : panPolyIdMap[nIdCur] = nId;
212 8813 : nIdCur = nNextId;
213 : }
214 :
215 17340 : if (panPolyIdMap[iPoly] == iPoly)
216 7269 : nFinalPolyCount++;
217 : }
218 :
219 62 : CPLDebug("GDALRasterPolygonEnumerator",
220 : "Counted %d polygon fragments forming %d final polygons.",
221 : nNextPolygonId, nFinalPolyCount);
222 62 : }
223 :
224 : /************************************************************************/
225 : /* ProcessLine() */
226 : /* */
227 : /* Assign ids to polygons, one line at a time. */
228 : /************************************************************************/
229 :
230 : template <class DataType, class EqualityTest>
231 3413 : bool GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::ProcessLine(
232 : DataType *panLastLineVal, DataType *panThisLineVal, GInt32 *panLastLineId,
233 : GInt32 *panThisLineId, int nXSize)
234 :
235 : {
236 : EqualityTest eq;
237 :
238 : /* -------------------------------------------------------------------- */
239 : /* Special case for the first line. */
240 : /* -------------------------------------------------------------------- */
241 3413 : if (panLastLineVal == nullptr)
242 : {
243 2842 : for (int i = 0; i < nXSize; i++)
244 : {
245 2708 : if (panThisLineVal[i] == GP_NODATA_MARKER)
246 : {
247 1069 : panThisLineId[i] = -1;
248 : }
249 3167 : else if (i == 0 ||
250 1528 : !(eq.operator()(panThisLineVal[i], panThisLineVal[i - 1])))
251 : {
252 690 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
253 690 : if (panThisLineId[i] < 0)
254 0 : return false;
255 : }
256 : else
257 : {
258 949 : panThisLineId[i] = panThisLineId[i - 1];
259 : }
260 : }
261 :
262 134 : return true;
263 : }
264 :
265 : /* -------------------------------------------------------------------- */
266 : /* Process each pixel comparing to the previous pixel, and to */
267 : /* the last line. */
268 : /* -------------------------------------------------------------------- */
269 875678 : for (int i = 0; i < nXSize; i++)
270 : {
271 872399 : if (panThisLineVal[i] == GP_NODATA_MARKER)
272 : {
273 104671 : panThisLineId[i] = -1;
274 : }
275 1533864 : else if (i > 0 &&
276 766141 : eq.operator()(panThisLineVal[i], panThisLineVal[i - 1]))
277 : {
278 690132 : panThisLineId[i] = panThisLineId[i - 1];
279 :
280 1335830 : if (eq.operator()(panLastLineVal[i], panThisLineVal[i]) &&
281 645697 : (panPolyIdMap[panLastLineId[i]] !=
282 645697 : panPolyIdMap[panThisLineId[i]]))
283 : {
284 23607 : MergePolygon(panLastLineId[i], panThisLineId[i]);
285 : }
286 :
287 158 : if (nConnectedness == 8 &&
288 690290 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]) &&
289 17 : (panPolyIdMap[panLastLineId[i - 1]] !=
290 17 : panPolyIdMap[panThisLineId[i]]))
291 : {
292 0 : MergePolygon(panLastLineId[i - 1], panThisLineId[i]);
293 : }
294 :
295 297 : if (nConnectedness == 8 && i < nXSize - 1 &&
296 690429 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
297 34 : (panPolyIdMap[panLastLineId[i + 1]] !=
298 34 : panPolyIdMap[panThisLineId[i]]))
299 : {
300 23 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
301 : }
302 : }
303 77596 : else if (eq.operator()(panLastLineVal[i], panThisLineVal[i]))
304 : {
305 36301 : panThisLineId[i] = panLastLineId[i];
306 : }
307 41937 : else if (i > 0 && nConnectedness == 8 &&
308 642 : eq.operator()(panLastLineVal[i - 1], panThisLineVal[i]))
309 : {
310 49 : panThisLineId[i] = panLastLineId[i - 1];
311 :
312 43 : if (i < nXSize - 1 &&
313 92 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]) &&
314 6 : (panPolyIdMap[panLastLineId[i + 1]] !=
315 6 : panPolyIdMap[panThisLineId[i]]))
316 : {
317 6 : MergePolygon(panLastLineId[i + 1], panThisLineId[i]);
318 : }
319 : }
320 41856 : else if (i < nXSize - 1 && nConnectedness == 8 &&
321 610 : eq.operator()(panLastLineVal[i + 1], panThisLineVal[i]))
322 : {
323 67 : panThisLineId[i] = panLastLineId[i + 1];
324 : }
325 : else
326 : {
327 41179 : panThisLineId[i] = NewPolygon(panThisLineVal[i]);
328 41179 : if (panThisLineId[i] < 0)
329 0 : return false;
330 : }
331 : }
332 3279 : return true;
333 : }
334 :
335 : template class GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest>;
336 :
337 : template class GDALRasterPolygonEnumeratorT<float, FloatEqualityTest>;
338 :
339 : /*! @endcond */
|
