Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Raster to Polygon Converter
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2008, Frank Warmerdam
9 : * Copyright (c) 2009-2011, Even Rouault <even dot 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.h"
32 :
33 : #include <cstring>
34 :
35 : #include <algorithm>
36 : #include <cassert>
37 : #include <set>
38 : #include <vector>
39 : #include <utility>
40 :
41 : #include "cpl_conv.h"
42 : #include "cpl_error.h"
43 : #include "cpl_progress.h"
44 : #include "cpl_vsi.h"
45 : #include "gdal.h"
46 : #include "gdal_alg_priv.h"
47 :
48 : #define MY_MAX_INT 2147483647
49 :
50 : /*
51 : * General Plan
52 : *
53 : * 1) make a pass with the polygon enumerator to build up the
54 : * polygon map array. Also accumulate polygon size information.
55 : *
56 : * 2) Identify the polygons that need to be merged.
57 : *
58 : * 3) Make a pass with the polygon enumerator. For each "to be merged"
59 : * polygon keep track of its largest neighbour.
60 : *
61 : * 4) Fix up remappings that would go to polygons smaller than the seive
62 : * size. Ensure these in term map to the largest neighbour of the
63 : * "to be sieved" polygons.
64 : *
65 : * 5) Make another pass with the polygon enumerator. This time we remap
66 : * the actual pixel values of all polygons to be merged.
67 : */
68 :
69 : /************************************************************************/
70 : /* GPMaskImageData() */
71 : /* */
72 : /* Mask out image pixels to a special nodata value if the mask */
73 : /* band is zero. */
74 : /************************************************************************/
75 :
76 123 : static CPLErr GPMaskImageData(GDALRasterBandH hMaskBand, GByte *pabyMaskLine,
77 : int iY, int nXSize, std::int64_t *panImageLine)
78 :
79 : {
80 123 : const CPLErr eErr = GDALRasterIO(hMaskBand, GF_Read, 0, iY, nXSize, 1,
81 : pabyMaskLine, nXSize, 1, GDT_Byte, 0, 0);
82 123 : if (eErr == CE_None)
83 : {
84 1722 : for (int i = 0; i < nXSize; i++)
85 : {
86 1599 : if (pabyMaskLine[i] == 0)
87 48 : panImageLine[i] = GP_NODATA_MARKER;
88 : }
89 : }
90 :
91 123 : return eErr;
92 : }
93 :
94 : // TODO: What is "eaches" supposed to be?
95 :
96 : /************************************************************************/
97 : /* CompareNeighbour() */
98 : /* */
99 : /* Compare two neighbouring polygons, and update eaches */
100 : /* "biggest neighbour" if the other is larger than its current */
101 : /* largest neighbour. */
102 : /* */
103 : /* Note that this should end up with each polygon knowing the */
104 : /* id of its largest neighbour. No attempt is made to */
105 : /* restrict things to small polygons that we will be merging, */
106 : /* nor to exclude assigning "biggest neighbours" that are still */
107 : /* smaller than our sieve threshold. */
108 : /************************************************************************/
109 :
110 22263 : static inline void CompareNeighbour(int nPolyId1, int nPolyId2,
111 : int *panPolyIdMap,
112 : std::int64_t * /* panPolyValue */,
113 : const std::vector<int> &anPolySizes,
114 : std::vector<int> &anBigNeighbour)
115 :
116 : {
117 : // Nodata polygon do not need neighbours, and cannot be neighbours
118 : // to valid polygons.
119 22263 : if (nPolyId1 < 0 || nPolyId2 < 0)
120 34 : return;
121 :
122 : // Make sure we are working with the final merged polygon ids.
123 22229 : nPolyId1 = panPolyIdMap[nPolyId1];
124 22229 : nPolyId2 = panPolyIdMap[nPolyId2];
125 :
126 22229 : if (nPolyId1 == nPolyId2)
127 7680 : return;
128 :
129 : // Nodata polygon do not need neighbours, and cannot be neighbours
130 : // to valid polygons.
131 : // Should no longer happen with r28826 optimization.
132 : // if( panPolyValue[nPolyId1] == GP_NODATA_MARKER
133 : // || panPolyValue[nPolyId2] == GP_NODATA_MARKER )
134 : // return;
135 :
136 25407 : if (anBigNeighbour[nPolyId1] == -1 ||
137 10858 : anPolySizes[anBigNeighbour[nPolyId1]] < anPolySizes[nPolyId2])
138 3811 : anBigNeighbour[nPolyId1] = nPolyId2;
139 :
140 29088 : if (anBigNeighbour[nPolyId2] == -1 ||
141 14539 : anPolySizes[anBigNeighbour[nPolyId2]] < anPolySizes[nPolyId1])
142 175 : anBigNeighbour[nPolyId2] = nPolyId1;
143 : }
144 :
145 : /************************************************************************/
146 : /* GDALSieveFilter() */
147 : /************************************************************************/
148 :
149 : /**
150 : * Removes small raster polygons.
151 : *
152 : * The function removes raster polygons smaller than a provided
153 : * threshold size (in pixels) and replaces them with the pixel value
154 : * of the largest neighbour polygon.
155 : *
156 : * Polygon are determined (per GDALRasterPolygonEnumerator) as regions of
157 : * the raster where the pixels all have the same value, and that are contiguous
158 : * (connected).
159 : *
160 : * Pixels determined to be "nodata" per hMaskBand will not be treated as part
161 : * of a polygon regardless of their pixel values. Nodata areas will never be
162 : * changed nor affect polygon sizes.
163 : *
164 : * Polygons smaller than the threshold with no neighbours that are as large
165 : * as the threshold will not be altered. Polygons surrounded by nodata areas
166 : * will therefore not be altered.
167 : *
168 : * The algorithm makes three passes over the input file to enumerate the
169 : * polygons and collect limited information about them. Memory use is
170 : * proportional to the number of polygons (roughly 24 bytes per polygon), but
171 : * is not directly related to the size of the raster. So very large raster
172 : * files can be processed effectively if there aren't too many polygons. But
173 : * extremely noisy rasters with many one pixel polygons will end up being
174 : * expensive (in memory) to process.
175 : *
176 : * @param hSrcBand the source raster band to be processed.
177 : * @param hMaskBand an optional mask band. All pixels in the mask band with a
178 : * value other than zero will be considered suitable for inclusion in polygons.
179 : * @param hDstBand the output raster band. It may be the same as hSrcBand
180 : * to update the source in place.
181 : * @param nSizeThreshold raster polygons with sizes smaller than this will
182 : * be merged into their largest neighbour.
183 : * @param nConnectedness either 4 indicating that diagonal pixels are not
184 : * considered directly adjacent for polygon membership purposes or 8
185 : * indicating they are.
186 : * @param papszOptions algorithm options in name=value list form. None
187 : * currently supported.
188 : * @param pfnProgress callback for reporting algorithm progress matching the
189 : * GDALProgressFunc() semantics. May be NULL.
190 : * @param pProgressArg callback argument passed to pfnProgress.
191 : *
192 : * @return CE_None on success or CE_Failure if an error occurs.
193 : */
194 :
195 10 : CPLErr CPL_STDCALL GDALSieveFilter(GDALRasterBandH hSrcBand,
196 : GDALRasterBandH hMaskBand,
197 : GDALRasterBandH hDstBand, int nSizeThreshold,
198 : int nConnectedness,
199 : CPL_UNUSED char **papszOptions,
200 : GDALProgressFunc pfnProgress,
201 : void *pProgressArg)
202 : {
203 10 : VALIDATE_POINTER1(hSrcBand, "GDALSieveFilter", CE_Failure);
204 10 : VALIDATE_POINTER1(hDstBand, "GDALSieveFilter", CE_Failure);
205 :
206 10 : if (pfnProgress == nullptr)
207 8 : pfnProgress = GDALDummyProgress;
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Allocate working buffers. */
211 : /* -------------------------------------------------------------------- */
212 10 : int nXSize = GDALGetRasterBandXSize(hSrcBand);
213 10 : int nYSize = GDALGetRasterBandYSize(hSrcBand);
214 : auto *panLastLineVal = static_cast<std::int64_t *>(
215 10 : VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
216 : auto *panThisLineVal = static_cast<std::int64_t *>(
217 10 : VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
218 : auto *panLastLineId =
219 10 : static_cast<GInt32 *>(VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
220 : auto *panThisLineId =
221 10 : static_cast<GInt32 *>(VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
222 : auto *panThisLineWriteVal = static_cast<std::int64_t *>(
223 10 : VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
224 : GByte *pabyMaskLine = hMaskBand != nullptr
225 10 : ? static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize))
226 10 : : nullptr;
227 10 : if (panLastLineVal == nullptr || panThisLineVal == nullptr ||
228 10 : panLastLineId == nullptr || panThisLineId == nullptr ||
229 10 : panThisLineWriteVal == nullptr ||
230 4 : (hMaskBand != nullptr && pabyMaskLine == nullptr))
231 : {
232 0 : CPLFree(panThisLineId);
233 0 : CPLFree(panLastLineId);
234 0 : CPLFree(panThisLineVal);
235 0 : CPLFree(panLastLineVal);
236 0 : CPLFree(panThisLineWriteVal);
237 0 : CPLFree(pabyMaskLine);
238 0 : return CE_Failure;
239 : }
240 :
241 : /* -------------------------------------------------------------------- */
242 : /* The first pass over the raster is only used to build up the */
243 : /* polygon id map so we will know in advance what polygons are */
244 : /* what on the second pass. */
245 : /* -------------------------------------------------------------------- */
246 20 : GDALRasterPolygonEnumerator oFirstEnum(nConnectedness);
247 20 : std::vector<int> anPolySizes;
248 :
249 10 : CPLErr eErr = CE_None;
250 189 : for (int iY = 0; eErr == CE_None && iY < nYSize; iY++)
251 : {
252 179 : eErr = GDALRasterIO(hSrcBand, GF_Read, 0, iY, nXSize, 1, panThisLineVal,
253 : nXSize, 1, GDT_Int64, 0, 0);
254 :
255 179 : if (eErr == CE_None && hMaskBand != nullptr)
256 41 : eErr = GPMaskImageData(hMaskBand, pabyMaskLine, iY, nXSize,
257 : panThisLineVal);
258 179 : if (eErr != CE_None)
259 0 : break;
260 :
261 179 : if (iY == 0)
262 20 : eErr = oFirstEnum.ProcessLine(nullptr, panThisLineVal, nullptr,
263 : panThisLineId, nXSize)
264 10 : ? CE_None
265 : : CE_Failure;
266 : else
267 338 : eErr = oFirstEnum.ProcessLine(panLastLineVal, panThisLineVal,
268 : panLastLineId, panThisLineId, nXSize)
269 169 : ? CE_None
270 : : CE_Failure;
271 179 : if (eErr != CE_None)
272 0 : break;
273 :
274 : /* --------------------------------------------------------------------
275 : */
276 : /* Accumulate polygon sizes. */
277 : /* --------------------------------------------------------------------
278 : */
279 179 : if (oFirstEnum.nNextPolygonId > static_cast<int>(anPolySizes.size()))
280 169 : anPolySizes.resize(oFirstEnum.nNextPolygonId);
281 :
282 11438 : for (int iX = 0; iX < nXSize; iX++)
283 : {
284 11259 : const int iPoly = panThisLineId[iX];
285 :
286 11259 : if (iPoly >= 0 && anPolySizes[iPoly] < MY_MAX_INT)
287 11243 : anPolySizes[iPoly] += 1;
288 : }
289 :
290 : /* --------------------------------------------------------------------
291 : */
292 : /* swap this/last lines. */
293 : /* --------------------------------------------------------------------
294 : */
295 179 : std::swap(panLastLineVal, panThisLineVal);
296 179 : std::swap(panLastLineId, panThisLineId);
297 :
298 : /* --------------------------------------------------------------------
299 : */
300 : /* Report progress, and support interrupts. */
301 : /* --------------------------------------------------------------------
302 : */
303 179 : if (!pfnProgress(0.25 * ((iY + 1) / static_cast<double>(nYSize)), "",
304 : pProgressArg))
305 : {
306 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
307 0 : eErr = CE_Failure;
308 : }
309 : }
310 :
311 : /* -------------------------------------------------------------------- */
312 : /* Make a pass through the maps, ensuring every polygon id */
313 : /* points to the final id it should use, not an intermediate */
314 : /* value. */
315 : /* -------------------------------------------------------------------- */
316 10 : if (eErr == CE_None)
317 10 : oFirstEnum.CompleteMerges();
318 :
319 : /* -------------------------------------------------------------------- */
320 : /* Push the sizes of merged polygon fragments into the */
321 : /* merged polygon id's count. */
322 : /* -------------------------------------------------------------------- */
323 10 : assert(oFirstEnum.panPolyIdMap != nullptr); // for Coverity
324 10 : assert(oFirstEnum.panPolyValue != nullptr); // for Coverity
325 7199 : for (int iPoly = 0; iPoly < oFirstEnum.nNextPolygonId; iPoly++)
326 : {
327 7189 : if (oFirstEnum.panPolyIdMap[iPoly] != iPoly)
328 : {
329 3488 : GIntBig nSize = anPolySizes[oFirstEnum.panPolyIdMap[iPoly]];
330 :
331 3488 : nSize += anPolySizes[iPoly];
332 :
333 3488 : if (nSize > MY_MAX_INT)
334 0 : nSize = MY_MAX_INT;
335 :
336 3488 : anPolySizes[oFirstEnum.panPolyIdMap[iPoly]] =
337 : static_cast<int>(nSize);
338 3488 : anPolySizes[iPoly] = 0;
339 : }
340 : }
341 :
342 : /* -------------------------------------------------------------------- */
343 : /* We will use a new enumerator for the second pass primarily */
344 : /* so we can preserve the first pass map. */
345 : /* -------------------------------------------------------------------- */
346 20 : GDALRasterPolygonEnumerator oSecondEnum(nConnectedness);
347 :
348 10 : std::vector<int> anBigNeighbour;
349 10 : anBigNeighbour.resize(anPolySizes.size());
350 :
351 7199 : for (int iPoly = 0; iPoly < static_cast<int>(anPolySizes.size()); iPoly++)
352 7189 : anBigNeighbour[iPoly] = -1;
353 :
354 : /* ==================================================================== */
355 : /* Second pass ... identify the largest neighbour for each */
356 : /* polygon. */
357 : /* ==================================================================== */
358 189 : for (int iY = 0; eErr == CE_None && iY < nYSize; iY++)
359 : {
360 : /* --------------------------------------------------------------------
361 : */
362 : /* Read the image data. */
363 : /* --------------------------------------------------------------------
364 : */
365 179 : eErr = GDALRasterIO(hSrcBand, GF_Read, 0, iY, nXSize, 1, panThisLineVal,
366 : nXSize, 1, GDT_Int64, 0, 0);
367 :
368 179 : if (eErr == CE_None && hMaskBand != nullptr)
369 41 : eErr = GPMaskImageData(hMaskBand, pabyMaskLine, iY, nXSize,
370 : panThisLineVal);
371 :
372 179 : if (eErr != CE_None)
373 0 : continue;
374 :
375 : /* --------------------------------------------------------------------
376 : */
377 : /* Determine what polygon the various pixels belong to (redoing */
378 : /* the same thing done in the first pass above). */
379 : /* --------------------------------------------------------------------
380 : */
381 179 : if (iY == 0)
382 20 : eErr = oSecondEnum.ProcessLine(nullptr, panThisLineVal, nullptr,
383 : panThisLineId, nXSize)
384 10 : ? CE_None
385 : : CE_Failure;
386 : else
387 338 : eErr = oSecondEnum.ProcessLine(panLastLineVal, panThisLineVal,
388 : panLastLineId, panThisLineId, nXSize)
389 169 : ? CE_None
390 : : CE_Failure;
391 :
392 179 : if (eErr != CE_None)
393 0 : continue;
394 :
395 : /* --------------------------------------------------------------------
396 : */
397 : /* Check our neighbours, and update our biggest neighbour map */
398 : /* as appropriate. */
399 : /* --------------------------------------------------------------------
400 : */
401 11438 : for (int iX = 0; iX < nXSize; iX++)
402 : {
403 11259 : if (iY > 0)
404 : {
405 11087 : CompareNeighbour(panThisLineId[iX], panLastLineId[iX],
406 : oFirstEnum.panPolyIdMap,
407 : oFirstEnum.panPolyValue, anPolySizes,
408 : anBigNeighbour);
409 :
410 11087 : if (iX > 0 && nConnectedness == 8)
411 48 : CompareNeighbour(panThisLineId[iX], panLastLineId[iX - 1],
412 : oFirstEnum.panPolyIdMap,
413 : oFirstEnum.panPolyValue, anPolySizes,
414 : anBigNeighbour);
415 :
416 11087 : if (iX < nXSize - 1 && nConnectedness == 8)
417 48 : CompareNeighbour(panThisLineId[iX], panLastLineId[iX + 1],
418 : oFirstEnum.panPolyIdMap,
419 : oFirstEnum.panPolyValue, anPolySizes,
420 : anBigNeighbour);
421 : }
422 :
423 11259 : if (iX > 0)
424 11080 : CompareNeighbour(panThisLineId[iX], panThisLineId[iX - 1],
425 : oFirstEnum.panPolyIdMap,
426 : oFirstEnum.panPolyValue, anPolySizes,
427 : anBigNeighbour);
428 :
429 : // We don't need to compare to next pixel or next line
430 : // since they will be compared to us.
431 : }
432 :
433 : /* --------------------------------------------------------------------
434 : */
435 : /* Swap pixel value, and polygon id lines to be ready for the */
436 : /* next line. */
437 : /* --------------------------------------------------------------------
438 : */
439 179 : std::swap(panLastLineVal, panThisLineVal);
440 179 : std::swap(panLastLineId, panThisLineId);
441 :
442 : /* --------------------------------------------------------------------
443 : */
444 : /* Report progress, and support interrupts. */
445 : /* --------------------------------------------------------------------
446 : */
447 179 : if (!pfnProgress(0.25 + 0.25 * ((iY + 1) / static_cast<double>(nYSize)),
448 : "", pProgressArg))
449 : {
450 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
451 0 : eErr = CE_Failure;
452 : }
453 : }
454 :
455 : /* -------------------------------------------------------------------- */
456 : /* If our biggest neighbour is still smaller than the */
457 : /* threshold, then try tracking to that polygons biggest */
458 : /* neighbour, and so forth. */
459 : /* -------------------------------------------------------------------- */
460 10 : int nFailedMerges = 0;
461 10 : int nIsolatedSmall = 0;
462 10 : int nSieveTargets = 0;
463 :
464 7199 : for (int iPoly = 0; iPoly < static_cast<int>(anPolySizes.size()); iPoly++)
465 : {
466 7189 : if (oFirstEnum.panPolyIdMap[iPoly] != iPoly)
467 3817 : continue;
468 :
469 : // Ignore nodata polygons.
470 3701 : if (oFirstEnum.panPolyValue[iPoly] == GP_NODATA_MARKER)
471 0 : continue;
472 :
473 : // Don't try to merge polygons larger than the threshold.
474 3701 : if (anPolySizes[iPoly] >= nSizeThreshold)
475 : {
476 309 : anBigNeighbour[iPoly] = -1;
477 309 : continue;
478 : }
479 :
480 3392 : nSieveTargets++;
481 :
482 : // if we have no neighbours but we are small, what shall we do?
483 3392 : if (anBigNeighbour[iPoly] == -1)
484 : {
485 0 : nIsolatedSmall++;
486 0 : continue;
487 : }
488 :
489 3392 : std::set<int> oSetVisitedPoly;
490 3392 : oSetVisitedPoly.insert(iPoly);
491 :
492 : // Walk through our neighbours until we find a polygon large enough.
493 3392 : int iFinalId = iPoly;
494 3392 : bool bFoundBigEnoughPoly = false;
495 : while (true)
496 : {
497 3420 : iFinalId = anBigNeighbour[iFinalId];
498 3420 : if (iFinalId < 0)
499 : {
500 19 : break;
501 : }
502 : // If the biggest neighbour is larger than the threshold
503 : // then we are golden.
504 3401 : if (anPolySizes[iFinalId] >= nSizeThreshold)
505 : {
506 3372 : bFoundBigEnoughPoly = true;
507 3372 : break;
508 : }
509 : // Check that we don't cycle on an already visited polygon.
510 29 : if (oSetVisitedPoly.find(iFinalId) != oSetVisitedPoly.end())
511 1 : break;
512 28 : oSetVisitedPoly.insert(iFinalId);
513 : }
514 :
515 3392 : if (!bFoundBigEnoughPoly)
516 : {
517 20 : nFailedMerges++;
518 20 : anBigNeighbour[iPoly] = -1;
519 20 : continue;
520 : }
521 :
522 : // Map the whole intermediate chain to it.
523 3372 : int iPolyCur = iPoly;
524 3380 : while (anBigNeighbour[iPolyCur] != iFinalId)
525 : {
526 8 : int iNextPoly = anBigNeighbour[iPolyCur];
527 8 : anBigNeighbour[iPolyCur] = iFinalId;
528 8 : iPolyCur = iNextPoly;
529 : }
530 : }
531 :
532 10 : CPLDebug("GDALSieveFilter",
533 : "Small Polygons: %d, Isolated: %d, Unmergable: %d", nSieveTargets,
534 : nIsolatedSmall, nFailedMerges);
535 :
536 : /* ==================================================================== */
537 : /* Make a third pass over the image, actually applying the */
538 : /* merges. We reuse the second enumerator but preserve the */
539 : /* "final maps" from the first. */
540 : /* ==================================================================== */
541 10 : oSecondEnum.Clear();
542 :
543 378 : for (int iY = 0; oFirstEnum.panPolyIdMap != nullptr && // for Coverity
544 378 : eErr == CE_None && iY < nYSize;
545 : iY++)
546 : {
547 : /* --------------------------------------------------------------------
548 : */
549 : /* Read the image data. */
550 : /* --------------------------------------------------------------------
551 : */
552 179 : eErr = GDALRasterIO(hSrcBand, GF_Read, 0, iY, nXSize, 1, panThisLineVal,
553 : nXSize, 1, GDT_Int64, 0, 0);
554 :
555 179 : memcpy(panThisLineWriteVal, panThisLineVal,
556 179 : sizeof(panThisLineVal[0]) * nXSize);
557 :
558 179 : if (eErr == CE_None && hMaskBand != nullptr)
559 41 : eErr = GPMaskImageData(hMaskBand, pabyMaskLine, iY, nXSize,
560 : panThisLineVal);
561 :
562 179 : if (eErr != CE_None)
563 0 : continue;
564 :
565 : /* --------------------------------------------------------------------
566 : */
567 : /* Determine what polygon the various pixels belong to (redoing */
568 : /* the same thing done in the first pass above). */
569 : /* --------------------------------------------------------------------
570 : */
571 179 : if (iY == 0)
572 10 : oSecondEnum.ProcessLine(nullptr, panThisLineVal, nullptr,
573 : panThisLineId, nXSize);
574 : else
575 169 : oSecondEnum.ProcessLine(panLastLineVal, panThisLineVal,
576 : panLastLineId, panThisLineId, nXSize);
577 :
578 : /* --------------------------------------------------------------------
579 : */
580 : /* Reprocess the actual pixel values according to the polygon */
581 : /* merging, and write out this line of image data. */
582 : /* --------------------------------------------------------------------
583 : */
584 11438 : for (int iX = 0; iX < nXSize; iX++)
585 : {
586 11259 : int iThisPoly = panThisLineId[iX];
587 11259 : if (iThisPoly >= 0)
588 : {
589 11243 : iThisPoly = oFirstEnum.panPolyIdMap[iThisPoly];
590 :
591 11243 : if (anBigNeighbour[iThisPoly] != -1)
592 : {
593 3377 : panThisLineWriteVal[iX] =
594 3377 : oFirstEnum.panPolyValue[anBigNeighbour[iThisPoly]];
595 : }
596 : }
597 : }
598 :
599 : /* --------------------------------------------------------------------
600 : */
601 : /* Write the update data out. */
602 : /* --------------------------------------------------------------------
603 : */
604 179 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, iY, nXSize, 1,
605 : panThisLineWriteVal, nXSize, 1, GDT_Int64, 0, 0);
606 :
607 : /* --------------------------------------------------------------------
608 : */
609 : /* Swap pixel value, and polygon id lines to be ready for the */
610 : /* next line. */
611 : /* --------------------------------------------------------------------
612 : */
613 179 : std::swap(panLastLineVal, panThisLineVal);
614 179 : std::swap(panLastLineId, panThisLineId);
615 :
616 : /* --------------------------------------------------------------------
617 : */
618 : /* Report progress, and support interrupts. */
619 : /* --------------------------------------------------------------------
620 : */
621 358 : if (eErr == CE_None &&
622 179 : !pfnProgress(0.5 + 0.5 * ((iY + 1) / static_cast<double>(nYSize)),
623 : "", pProgressArg))
624 : {
625 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
626 0 : eErr = CE_Failure;
627 : }
628 : }
629 :
630 : /* -------------------------------------------------------------------- */
631 : /* Cleanup */
632 : /* -------------------------------------------------------------------- */
633 10 : CPLFree(panThisLineId);
634 10 : CPLFree(panLastLineId);
635 10 : CPLFree(panThisLineVal);
636 10 : CPLFree(panLastLineVal);
637 10 : CPLFree(panThisLineWriteVal);
638 10 : CPLFree(pabyMaskLine);
639 :
640 10 : return eErr;
641 : }
|
