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