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