Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Interpolate in nodata areas.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2008, Frank Warmerdam
9 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2015, Sean Gillies <sean@mapbox.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ***************************************************************************/
14 :
15 : #include "cpl_port.h"
16 : #include "gdal_alg.h"
17 :
18 : #include <cmath>
19 : #include <cstring>
20 :
21 : #include <algorithm>
22 : #include <string>
23 : #include <utility>
24 :
25 : #include "cpl_conv.h"
26 : #include "cpl_error.h"
27 : #include "cpl_progress.h"
28 : #include "cpl_string.h"
29 : #include "cpl_vsi.h"
30 : #include "gdal.h"
31 : #include "gdal_priv.h"
32 :
33 : /************************************************************************/
34 : /* GDALFilterLine() */
35 : /* */
36 : /* Apply 3x3 filtering one one scanline with masking for which */
37 : /* pixels are to be interpolated (ThisFMask) and which window */
38 : /* pixels are valid to include in the interpolation (TMask). */
39 : /************************************************************************/
40 :
41 2280 : static void GDALFilterLine(const float *pafLastLine, const float *pafThisLine,
42 : const float *pafNextLine, float *pafOutLine,
43 : const GByte *pabyLastTMask,
44 : const GByte *pabyThisTMask,
45 : const GByte *pabyNextTMask,
46 : const GByte *pabyThisFMask, int nXSize)
47 :
48 : {
49 631866 : for (int iX = 0; iX < nXSize; iX++)
50 : {
51 629586 : if (!pabyThisFMask[iX])
52 : {
53 579740 : pafOutLine[iX] = pafThisLine[iX];
54 579740 : continue;
55 : }
56 :
57 49846 : CPLAssert(pabyThisTMask[iX]);
58 :
59 49846 : double dfValSum = 0.0;
60 49846 : double dfWeightSum = 0.0;
61 :
62 : // Previous line.
63 49846 : if (pafLastLine != nullptr)
64 : {
65 49846 : if (iX > 0 && pabyLastTMask[iX - 1])
66 : {
67 47442 : dfValSum += pafLastLine[iX - 1];
68 47442 : dfWeightSum += 1.0;
69 : }
70 49846 : if (pabyLastTMask[iX])
71 : {
72 48642 : dfValSum += pafLastLine[iX];
73 48642 : dfWeightSum += 1.0;
74 : }
75 49846 : if (iX < nXSize - 1 && pabyLastTMask[iX + 1])
76 : {
77 47282 : dfValSum += pafLastLine[iX + 1];
78 47282 : dfWeightSum += 1.0;
79 : }
80 : }
81 :
82 : // Current Line.
83 49846 : if (iX > 0 && pabyThisTMask[iX - 1])
84 : {
85 47846 : dfValSum += pafThisLine[iX - 1];
86 47846 : dfWeightSum += 1.0;
87 : }
88 49846 : if (pabyThisTMask[iX])
89 : {
90 49846 : dfValSum += pafThisLine[iX];
91 49846 : dfWeightSum += 1.0;
92 : }
93 49846 : if (iX < nXSize - 1 && pabyThisTMask[iX + 1])
94 : {
95 47644 : dfValSum += pafThisLine[iX + 1];
96 47644 : dfWeightSum += 1.0;
97 : }
98 :
99 : // Next line.
100 49846 : if (pafNextLine != nullptr)
101 : {
102 49846 : if (iX > 0 && pabyNextTMask[iX - 1])
103 : {
104 41958 : dfValSum += pafNextLine[iX - 1];
105 41958 : dfWeightSum += 1.0;
106 : }
107 49846 : if (pabyNextTMask[iX])
108 : {
109 43084 : dfValSum += pafNextLine[iX];
110 43084 : dfWeightSum += 1.0;
111 : }
112 49846 : if (iX < nXSize - 1 && pabyNextTMask[iX + 1])
113 : {
114 41720 : dfValSum += pafNextLine[iX + 1];
115 41720 : dfWeightSum += 1.0;
116 : }
117 : }
118 :
119 49846 : pafOutLine[iX] = static_cast<float>(dfValSum / dfWeightSum);
120 : }
121 2280 : }
122 :
123 : /************************************************************************/
124 : /* GDALMultiFilter() */
125 : /* */
126 : /* Apply multiple iterations of a 3x3 smoothing filter over a */
127 : /* band with masking controlling what pixels should be */
128 : /* filtered (FiltMaskBand non zero) and which pixels can be */
129 : /* considered valid contributors to the filter */
130 : /* (TargetMaskBand non zero). */
131 : /* */
132 : /* This implementation attempts to apply many iterations in */
133 : /* one IO pass by managing the filtering over a rolling buffer */
134 : /* of nIterations+2 scanlines. While possibly clever this */
135 : /* makes the algorithm implementation largely */
136 : /* incomprehensible. */
137 : /************************************************************************/
138 :
139 92 : static CPLErr GDALMultiFilter(GDALRasterBandH hTargetBand,
140 : GDALRasterBandH hTargetMaskBand,
141 : GDALRasterBandH hFiltMaskBand, int nIterations,
142 : GDALProgressFunc pfnProgress, void *pProgressArg)
143 :
144 : {
145 92 : const int nXSize = GDALGetRasterBandXSize(hTargetBand);
146 92 : const int nYSize = GDALGetRasterBandYSize(hTargetBand);
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Report starting progress value. */
150 : /* -------------------------------------------------------------------- */
151 92 : if (!pfnProgress(0.0, "Smoothing Filter...", pProgressArg))
152 : {
153 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
154 0 : return CE_Failure;
155 : }
156 :
157 : /* -------------------------------------------------------------------- */
158 : /* Allocate rotating buffers. */
159 : /* -------------------------------------------------------------------- */
160 92 : const int nBufLines = nIterations + 2;
161 :
162 : GByte *pabyTMaskBuf =
163 92 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nBufLines));
164 : GByte *pabyFMaskBuf =
165 92 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nBufLines));
166 :
167 : float *paf3PassLineBuf = static_cast<float *>(
168 92 : VSI_MALLOC3_VERBOSE(nXSize, nBufLines, 3 * sizeof(float)));
169 92 : if (pabyTMaskBuf == nullptr || pabyFMaskBuf == nullptr ||
170 : paf3PassLineBuf == nullptr)
171 : {
172 0 : CPLFree(pabyTMaskBuf);
173 0 : CPLFree(pabyFMaskBuf);
174 0 : CPLFree(paf3PassLineBuf);
175 :
176 0 : return CE_Failure;
177 : }
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Process rotating buffers. */
181 : /* -------------------------------------------------------------------- */
182 :
183 92 : CPLErr eErr = CE_None;
184 :
185 92 : int iPassCounter = 0;
186 :
187 2612 : for (int nNewLine = 0; // Line being loaded (zero based scanline).
188 2612 : eErr == CE_None && nNewLine < nYSize + nIterations; nNewLine++)
189 : {
190 : /* --------------------------------------------------------------------
191 : */
192 : /* Rotate pass buffers. */
193 : /* --------------------------------------------------------------------
194 : */
195 2520 : iPassCounter = (iPassCounter + 1) % 3;
196 :
197 2520 : float *const pafSLastPass =
198 2520 : paf3PassLineBuf + ((iPassCounter + 0) % 3) * nXSize * nBufLines;
199 2520 : float *const pafLastPass =
200 2520 : paf3PassLineBuf + ((iPassCounter + 1) % 3) * nXSize * nBufLines;
201 2520 : float *const pafThisPass =
202 2520 : paf3PassLineBuf + ((iPassCounter + 2) % 3) * nXSize * nBufLines;
203 :
204 : /* --------------------------------------------------------------------
205 : */
206 : /* Where does the new line go in the rotating buffer? */
207 : /* --------------------------------------------------------------------
208 : */
209 2520 : const int iBufOffset = nNewLine % nBufLines;
210 :
211 : /* --------------------------------------------------------------------
212 : */
213 : /* Read the new data line if it is't off the bottom of the */
214 : /* image. */
215 : /* --------------------------------------------------------------------
216 : */
217 2520 : if (nNewLine < nYSize)
218 : {
219 4820 : eErr = GDALRasterIO(hTargetMaskBand, GF_Read, 0, nNewLine, nXSize,
220 2410 : 1, pabyTMaskBuf + nXSize * iBufOffset, nXSize,
221 : 1, GDT_Byte, 0, 0);
222 :
223 2410 : if (eErr != CE_None)
224 0 : break;
225 :
226 4820 : eErr = GDALRasterIO(hFiltMaskBand, GF_Read, 0, nNewLine, nXSize, 1,
227 2410 : pabyFMaskBuf + nXSize * iBufOffset, nXSize, 1,
228 : GDT_Byte, 0, 0);
229 :
230 2410 : if (eErr != CE_None)
231 0 : break;
232 :
233 4820 : eErr = GDALRasterIO(hTargetBand, GF_Read, 0, nNewLine, nXSize, 1,
234 2410 : pafThisPass + nXSize * iBufOffset, nXSize, 1,
235 : GDT_Float32, 0, 0);
236 :
237 2410 : if (eErr != CE_None)
238 0 : break;
239 : }
240 :
241 : /* --------------------------------------------------------------------
242 : */
243 : /* Loop over the loaded data, applying the filter to all loaded */
244 : /* lines with neighbours. */
245 : /* --------------------------------------------------------------------
246 : */
247 5310 : for (int iFLine = nNewLine - 1;
248 5310 : eErr == CE_None && iFLine >= nNewLine - nIterations; iFLine--)
249 : {
250 2790 : const int iLastOffset = (iFLine - 1) % nBufLines;
251 2790 : const int iThisOffset = (iFLine) % nBufLines;
252 2790 : const int iNextOffset = (iFLine + 1) % nBufLines;
253 :
254 : // Default to preserving the old value.
255 2790 : if (iFLine >= 0)
256 2590 : memcpy(pafThisPass + iThisOffset * nXSize,
257 2590 : pafLastPass + iThisOffset * nXSize,
258 2590 : sizeof(float) * nXSize);
259 :
260 : // TODO: Enable first and last line.
261 : // Skip the first and last line.
262 2790 : if (iFLine < 1 || iFLine >= nYSize - 1)
263 : {
264 510 : continue;
265 : }
266 :
267 2280 : GDALFilterLine(pafSLastPass + iLastOffset * nXSize,
268 2280 : pafLastPass + iThisOffset * nXSize,
269 2280 : pafThisPass + iNextOffset * nXSize,
270 2280 : pafThisPass + iThisOffset * nXSize,
271 2280 : pabyTMaskBuf + iLastOffset * nXSize,
272 2280 : pabyTMaskBuf + iThisOffset * nXSize,
273 2280 : pabyTMaskBuf + iNextOffset * nXSize,
274 2280 : pabyFMaskBuf + iThisOffset * nXSize, nXSize);
275 : }
276 :
277 : /* --------------------------------------------------------------------
278 : */
279 : /* Write out the top data line that will be rolling out of our */
280 : /* buffer. */
281 : /* --------------------------------------------------------------------
282 : */
283 2520 : const int iLineToSave = nNewLine - nIterations;
284 :
285 2520 : if (iLineToSave >= 0 && eErr == CE_None)
286 : {
287 2410 : const int iBufOffset2 = iLineToSave % nBufLines;
288 :
289 2410 : eErr = GDALRasterIO(hTargetBand, GF_Write, 0, iLineToSave, nXSize,
290 2410 : 1, pafThisPass + nXSize * iBufOffset2, nXSize,
291 : 1, GDT_Float32, 0, 0);
292 : }
293 :
294 : /* --------------------------------------------------------------------
295 : */
296 : /* Report progress. */
297 : /* --------------------------------------------------------------------
298 : */
299 5040 : if (eErr == CE_None &&
300 2520 : !pfnProgress((nNewLine + 1) /
301 2520 : static_cast<double>(nYSize + nIterations),
302 : "Smoothing Filter...", pProgressArg))
303 : {
304 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
305 0 : eErr = CE_Failure;
306 : }
307 : }
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Cleanup */
311 : /* -------------------------------------------------------------------- */
312 92 : CPLFree(pabyTMaskBuf);
313 92 : CPLFree(pabyFMaskBuf);
314 92 : CPLFree(paf3PassLineBuf);
315 :
316 92 : return eErr;
317 : }
318 :
319 : /************************************************************************/
320 : /* QUAD_CHECK() */
321 : /* */
322 : /* macro for checking whether a point is nearer than the */
323 : /* existing closest point. */
324 : /************************************************************************/
325 :
326 4370290 : inline void QUAD_CHECK(double &dfQuadDist, float &fQuadValue, int target_x,
327 : GUInt32 target_y, int origin_x, int origin_y,
328 : float fTargetValue, GUInt32 nNoDataVal)
329 : {
330 4370290 : if (target_y != nNoDataVal)
331 : {
332 468800 : const double dfDx =
333 468800 : static_cast<double>(target_x) - static_cast<double>(origin_x);
334 468800 : const double dfDy =
335 468800 : static_cast<double>(target_y) - static_cast<double>(origin_y);
336 468800 : double dfDistSq = dfDx * dfDx + dfDy * dfDy;
337 :
338 468800 : if (dfDistSq < dfQuadDist * dfQuadDist)
339 : {
340 170775 : CPLAssert(dfDistSq > 0.0);
341 170775 : dfQuadDist = sqrt(dfDistSq);
342 170775 : fQuadValue = fTargetValue;
343 : }
344 : }
345 4370290 : }
346 :
347 : /************************************************************************/
348 : /* GDALFillNodata() */
349 : /************************************************************************/
350 :
351 : /**
352 : * Fill selected raster regions by interpolation from the edges.
353 : *
354 : * This algorithm will interpolate values for all designated
355 : * nodata pixels (marked by zeros in hMaskBand). For each pixel
356 : * a four direction conic search is done to find values to interpolate
357 : * from (using inverse distance weighting by default). Once all values are
358 : * interpolated, zero or more smoothing iterations (3x3 average
359 : * filters on interpolated pixels) are applied to smooth out
360 : * artifacts.
361 : *
362 : * This algorithm is generally suitable for interpolating missing
363 : * regions of fairly continuously varying rasters (such as elevation
364 : * models for instance). It is also suitable for filling small holes
365 : * and cracks in more irregularly varying images (like airphotos). It
366 : * is generally not so great for interpolating a raster from sparse
367 : * point data - see the algorithms defined in gdal_grid.h for that case.
368 : *
369 : * @param hTargetBand the raster band to be modified in place.
370 : * @param hMaskBand a mask band indicating pixels to be interpolated
371 : * (zero valued). If hMaskBand is set to NULL, this method will internally use
372 : * the mask band returned by GDALGetMaskBand(hTargetBand).
373 : * @param dfMaxSearchDist the maximum number of pixels to search in all
374 : * directions to find values to interpolate from.
375 : * @param bDeprecatedOption unused argument, should be zero.
376 : * @param nSmoothingIterations the number of 3x3 smoothing filter passes to
377 : * run (0 or more).
378 : * @param papszOptions additional name=value options in a string list.
379 : * <ul>
380 : * <li>TEMP_FILE_DRIVER=gdal_driver_name. For example MEM.</li>
381 : * <li>NODATA=value (starting with GDAL 2.4).
382 : * Source pixels at that value will be ignored by the interpolator. Warning:
383 : * currently this will not be honored by smoothing passes.</li>
384 : * <li>INTERPOLATION=INV_DIST/NEAREST (GDAL >= 3.9). By default, pixels are
385 : * interpolated using an inverse distance weighting (INV_DIST). It is also
386 : * possible to choose a nearest neighbour (NEAREST) strategy.</li>
387 : * </ul>
388 : * @param pfnProgress the progress function to report completion.
389 : * @param pProgressArg callback data for progress function.
390 : *
391 : * @return CE_None on success or CE_Failure if something goes wrong.
392 : */
393 :
394 117 : CPLErr CPL_STDCALL GDALFillNodata(GDALRasterBandH hTargetBand,
395 : GDALRasterBandH hMaskBand,
396 : double dfMaxSearchDist,
397 : CPL_UNUSED int bDeprecatedOption,
398 : int nSmoothingIterations, char **papszOptions,
399 : GDALProgressFunc pfnProgress,
400 : void *pProgressArg)
401 :
402 : {
403 117 : VALIDATE_POINTER1(hTargetBand, "GDALFillNodata", CE_Failure);
404 :
405 117 : const int nXSize = GDALGetRasterBandXSize(hTargetBand);
406 117 : const int nYSize = GDALGetRasterBandYSize(hTargetBand);
407 :
408 117 : if (dfMaxSearchDist == 0.0)
409 0 : dfMaxSearchDist = std::max(nXSize, nYSize) + 1;
410 :
411 117 : const int nMaxSearchDist = static_cast<int>(floor(dfMaxSearchDist));
412 :
413 : const char *pszInterpolation =
414 117 : CSLFetchNameValueDef(papszOptions, "INTERPOLATION", "INV_DIST");
415 117 : const bool bNearest = EQUAL(pszInterpolation, "NEAREST");
416 117 : if (!EQUAL(pszInterpolation, "INV_DIST") &&
417 6 : !EQUAL(pszInterpolation, "NEAREST"))
418 : {
419 0 : CPLError(CE_Failure, CPLE_NotSupported,
420 : "Unsupported interpolation method: %s", pszInterpolation);
421 0 : return CE_Failure;
422 : }
423 :
424 : // Special "x" pixel values identifying pixels as special.
425 117 : GDALDataType eType = GDT_UInt16;
426 117 : GUInt32 nNoDataVal = 65535;
427 :
428 117 : if (nXSize > 65533 || nYSize > 65533)
429 : {
430 0 : eType = GDT_UInt32;
431 0 : nNoDataVal = 4000002;
432 : }
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Determine format driver for temp work files. */
436 : /* -------------------------------------------------------------------- */
437 : CPLString osTmpFileDriver =
438 234 : CSLFetchNameValueDef(papszOptions, "TEMP_FILE_DRIVER", "GTiff");
439 117 : GDALDriverH hDriver = GDALGetDriverByName(osTmpFileDriver.c_str());
440 :
441 117 : if (hDriver == nullptr)
442 : {
443 0 : CPLError(CE_Failure, CPLE_AppDefined,
444 : "TEMP_FILE_DRIVER=%s driver is not registered",
445 : osTmpFileDriver.c_str());
446 0 : return CE_Failure;
447 : }
448 :
449 117 : if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr)
450 : {
451 0 : CPLError(CE_Failure, CPLE_AppDefined,
452 : "TEMP_FILE_DRIVER=%s driver is incapable of creating "
453 : "temp work files",
454 : osTmpFileDriver.c_str());
455 0 : return CE_Failure;
456 : }
457 :
458 234 : CPLStringList aosWorkFileOptions;
459 117 : if (osTmpFileDriver == "GTiff")
460 : {
461 117 : aosWorkFileOptions.SetNameValue("COMPRESS", "LZW");
462 117 : aosWorkFileOptions.SetNameValue("BIGTIFF", "IF_SAFER");
463 : }
464 :
465 234 : const CPLString osTmpFile = CPLGenerateTempFilename("");
466 :
467 117 : std::unique_ptr<GDALDataset> poTmpMaskDS;
468 117 : if (hMaskBand == nullptr)
469 : {
470 112 : hMaskBand = GDALGetMaskBand(hTargetBand);
471 : }
472 7 : else if (nSmoothingIterations > 0 &&
473 2 : hMaskBand != GDALGetMaskBand(hTargetBand))
474 : {
475 : // If doing smoothing operations and the user provided its own
476 : // mask band, we must make a copy of it to be able to update it
477 : // when we fill pixels during the initial pass.
478 1 : const CPLString osMaskTmpFile = osTmpFile + "fill_mask_work.tif";
479 1 : poTmpMaskDS.reset(GDALDataset::FromHandle(
480 : GDALCreate(hDriver, osMaskTmpFile, nXSize, nYSize, 1, GDT_Byte,
481 1 : aosWorkFileOptions.List())));
482 1 : if (poTmpMaskDS == nullptr)
483 : {
484 0 : CPLError(CE_Failure, CPLE_AppDefined,
485 : "Could not create poTmpMaskDS work file. Check driver "
486 : "capabilities.");
487 0 : return CE_Failure;
488 : }
489 1 : poTmpMaskDS->MarkSuppressOnClose();
490 : auto hTmpMaskBand =
491 1 : GDALRasterBand::ToHandle(poTmpMaskDS->GetRasterBand(1));
492 1 : if (GDALRasterBandCopyWholeRaster(hMaskBand, hTmpMaskBand, nullptr,
493 1 : nullptr, nullptr) != CE_None)
494 : {
495 0 : return CE_Failure;
496 : }
497 1 : hMaskBand = hTmpMaskBand;
498 : }
499 :
500 : // If there are smoothing iterations, reserve 10% of the progress for them.
501 117 : const double dfProgressRatio = nSmoothingIterations > 0 ? 0.9 : 1.0;
502 :
503 117 : const char *pszNoData = CSLFetchNameValue(papszOptions, "NODATA");
504 117 : bool bHasNoData = false;
505 117 : float fNoData = 0.0f;
506 117 : if (pszNoData)
507 : {
508 2 : bHasNoData = true;
509 2 : fNoData = static_cast<float>(CPLAtof(pszNoData));
510 : }
511 :
512 : /* -------------------------------------------------------------------- */
513 : /* Initialize progress counter. */
514 : /* -------------------------------------------------------------------- */
515 117 : if (pfnProgress == nullptr)
516 113 : pfnProgress = GDALDummyProgress;
517 :
518 117 : if (!pfnProgress(0.0, "Filling...", pProgressArg))
519 : {
520 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
521 0 : return CE_Failure;
522 : }
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Create a work file to hold the Y "last value" indices. */
526 : /* -------------------------------------------------------------------- */
527 234 : const CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";
528 :
529 : auto poYDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
530 : GDALCreate(hDriver, osYTmpFile, nXSize, nYSize, 1, eType,
531 234 : aosWorkFileOptions.List())));
532 :
533 117 : if (poYDS == nullptr)
534 : {
535 0 : CPLError(
536 : CE_Failure, CPLE_AppDefined,
537 : "Could not create Y index work file. Check driver capabilities.");
538 0 : return CE_Failure;
539 : }
540 117 : poYDS->MarkSuppressOnClose();
541 :
542 : GDALRasterBandH hYBand =
543 117 : GDALRasterBand::FromHandle(poYDS->GetRasterBand(1));
544 :
545 : /* -------------------------------------------------------------------- */
546 : /* Create a work file to hold the pixel value associated with */
547 : /* the "last xy value" pixel. */
548 : /* -------------------------------------------------------------------- */
549 234 : const CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";
550 :
551 : auto poValDS =
552 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALCreate(
553 : hDriver, osValTmpFile, nXSize, nYSize, 1,
554 234 : GDALGetRasterDataType(hTargetBand), aosWorkFileOptions.List())));
555 :
556 117 : if (poValDS == nullptr)
557 : {
558 0 : CPLError(
559 : CE_Failure, CPLE_AppDefined,
560 : "Could not create XY value work file. Check driver capabilities.");
561 0 : return CE_Failure;
562 : }
563 117 : poValDS->MarkSuppressOnClose();
564 :
565 : GDALRasterBandH hValBand =
566 117 : GDALRasterBand::FromHandle(poValDS->GetRasterBand(1));
567 :
568 : /* -------------------------------------------------------------------- */
569 : /* Create a mask file to make it clear what pixels can be filtered */
570 : /* on the filtering pass. */
571 : /* -------------------------------------------------------------------- */
572 234 : const CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";
573 :
574 : auto poFiltMaskDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
575 : GDALCreate(hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1, GDT_Byte,
576 234 : aosWorkFileOptions.List())));
577 :
578 117 : if (poFiltMaskDS == nullptr)
579 : {
580 0 : CPLError(CE_Failure, CPLE_AppDefined,
581 : "Could not create mask work file. Check driver capabilities.");
582 0 : return CE_Failure;
583 : }
584 117 : poFiltMaskDS->MarkSuppressOnClose();
585 :
586 : GDALRasterBandH hFiltMaskBand =
587 117 : GDALRasterBand::FromHandle(poFiltMaskDS->GetRasterBand(1));
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Allocate buffers for last scanline and this scanline. */
591 : /* -------------------------------------------------------------------- */
592 :
593 : GUInt32 *panLastY =
594 117 : static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));
595 : GUInt32 *panThisY =
596 117 : static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));
597 : GUInt32 *panTopDownY =
598 117 : static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));
599 : float *pafLastValue =
600 117 : static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));
601 : float *pafThisValue =
602 117 : static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));
603 : float *pafTopDownValue =
604 117 : static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));
605 : float *pafScanline =
606 117 : static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));
607 117 : GByte *pabyMask = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nXSize, 1));
608 117 : GByte *pabyFiltMask = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nXSize, 1));
609 :
610 117 : CPLErr eErr = CE_None;
611 :
612 117 : if (panLastY == nullptr || panThisY == nullptr || panTopDownY == nullptr ||
613 117 : pafLastValue == nullptr || pafThisValue == nullptr ||
614 117 : pafTopDownValue == nullptr || pafScanline == nullptr ||
615 117 : pabyMask == nullptr || pabyFiltMask == nullptr)
616 : {
617 0 : eErr = CE_Failure;
618 0 : goto end;
619 : }
620 :
621 11318 : for (int iX = 0; iX < nXSize; iX++)
622 : {
623 11201 : panLastY[iX] = nNoDataVal;
624 : }
625 :
626 : /* ==================================================================== */
627 : /* Make first pass from top to bottom collecting the "last */
628 : /* known value" for each column and writing it out to the work */
629 : /* files. */
630 : /* ==================================================================== */
631 :
632 2624 : for (int iY = 0; iY < nYSize && eErr == CE_None; iY++)
633 : {
634 : /* --------------------------------------------------------------------
635 : */
636 : /* Read data and mask for this line. */
637 : /* --------------------------------------------------------------------
638 : */
639 2507 : eErr = GDALRasterIO(hMaskBand, GF_Read, 0, iY, nXSize, 1, pabyMask,
640 : nXSize, 1, GDT_Byte, 0, 0);
641 :
642 2507 : if (eErr != CE_None)
643 0 : break;
644 :
645 2507 : eErr = GDALRasterIO(hTargetBand, GF_Read, 0, iY, nXSize, 1, pafScanline,
646 : nXSize, 1, GDT_Float32, 0, 0);
647 :
648 2507 : if (eErr != CE_None)
649 0 : break;
650 :
651 : /* --------------------------------------------------------------------
652 : */
653 : /* Figure out the most recent pixel for each column. */
654 : /* --------------------------------------------------------------------
655 : */
656 :
657 654985 : for (int iX = 0; iX < nXSize; iX++)
658 : {
659 652478 : if (pabyMask[iX])
660 : {
661 341578 : pafThisValue[iX] = pafScanline[iX];
662 341578 : panThisY[iX] = iY;
663 : }
664 310900 : else if (iY <= dfMaxSearchDist + panLastY[iX])
665 : {
666 303766 : pafThisValue[iX] = pafLastValue[iX];
667 303766 : panThisY[iX] = panLastY[iX];
668 : }
669 : else
670 : {
671 7134 : panThisY[iX] = nNoDataVal;
672 : }
673 : }
674 :
675 : /* --------------------------------------------------------------------
676 : */
677 : /* Write out best index/value to working files. */
678 : /* --------------------------------------------------------------------
679 : */
680 2507 : eErr = GDALRasterIO(hYBand, GF_Write, 0, iY, nXSize, 1, panThisY,
681 : nXSize, 1, GDT_UInt32, 0, 0);
682 2507 : if (eErr != CE_None)
683 0 : break;
684 :
685 2507 : eErr = GDALRasterIO(hValBand, GF_Write, 0, iY, nXSize, 1, pafThisValue,
686 : nXSize, 1, GDT_Float32, 0, 0);
687 2507 : if (eErr != CE_None)
688 0 : break;
689 :
690 : /* --------------------------------------------------------------------
691 : */
692 : /* Flip this/last buffers. */
693 : /* --------------------------------------------------------------------
694 : */
695 2507 : std::swap(pafThisValue, pafLastValue);
696 2507 : std::swap(panThisY, panLastY);
697 :
698 : /* --------------------------------------------------------------------
699 : */
700 : /* report progress. */
701 : /* --------------------------------------------------------------------
702 : */
703 2507 : if (!pfnProgress(dfProgressRatio *
704 2507 : (0.5 * (iY + 1) / static_cast<double>(nYSize)),
705 : "Filling...", pProgressArg))
706 : {
707 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
708 0 : eErr = CE_Failure;
709 : }
710 : }
711 :
712 11318 : for (int iX = 0; iX < nXSize; iX++)
713 : {
714 11201 : panLastY[iX] = nNoDataVal;
715 : }
716 :
717 : /* ==================================================================== */
718 : /* Now we will do collect similar this/last information from */
719 : /* bottom to top and use it in combination with the top to */
720 : /* bottom search info to interpolate. */
721 : /* ==================================================================== */
722 2624 : for (int iY = nYSize - 1; iY >= 0 && eErr == CE_None; iY--)
723 : {
724 2507 : eErr = GDALRasterIO(hMaskBand, GF_Read, 0, iY, nXSize, 1, pabyMask,
725 : nXSize, 1, GDT_Byte, 0, 0);
726 :
727 2507 : if (eErr != CE_None)
728 0 : break;
729 :
730 2507 : eErr = GDALRasterIO(hTargetBand, GF_Read, 0, iY, nXSize, 1, pafScanline,
731 : nXSize, 1, GDT_Float32, 0, 0);
732 :
733 2507 : if (eErr != CE_None)
734 0 : break;
735 :
736 : /* --------------------------------------------------------------------
737 : */
738 : /* Figure out the most recent pixel for each column. */
739 : /* --------------------------------------------------------------------
740 : */
741 :
742 654985 : for (int iX = 0; iX < nXSize; iX++)
743 : {
744 652478 : if (pabyMask[iX])
745 : {
746 341578 : pafThisValue[iX] = pafScanline[iX];
747 341578 : panThisY[iX] = iY;
748 : }
749 310900 : else if (panLastY[iX] - iY <= dfMaxSearchDist)
750 : {
751 25254 : pafThisValue[iX] = pafLastValue[iX];
752 25254 : panThisY[iX] = panLastY[iX];
753 : }
754 : else
755 : {
756 285646 : panThisY[iX] = nNoDataVal;
757 : }
758 : }
759 :
760 : /* --------------------------------------------------------------------
761 : */
762 : /* Load the last y and corresponding value from the top down pass.
763 : */
764 : /* --------------------------------------------------------------------
765 : */
766 2507 : eErr = GDALRasterIO(hYBand, GF_Read, 0, iY, nXSize, 1, panTopDownY,
767 : nXSize, 1, GDT_UInt32, 0, 0);
768 :
769 2507 : if (eErr != CE_None)
770 0 : break;
771 :
772 2507 : eErr = GDALRasterIO(hValBand, GF_Read, 0, iY, nXSize, 1,
773 : pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0);
774 :
775 2507 : if (eErr != CE_None)
776 0 : break;
777 :
778 : /* --------------------------------------------------------------------
779 : */
780 : /* Attempt to interpolate any pixels that are nodata. */
781 : /* --------------------------------------------------------------------
782 : */
783 2507 : memset(pabyFiltMask, 0, nXSize);
784 654985 : for (int iX = 0; iX < nXSize; iX++)
785 : {
786 652478 : int nThisMaxSearchDist = nMaxSearchDist;
787 :
788 : // If this was a valid target - no change.
789 652478 : if (pabyMask[iX])
790 341578 : continue;
791 :
792 : enum Quadrants
793 : {
794 : QUAD_TOP_LEFT = 0,
795 : QUAD_BOTTOM_LEFT = 1,
796 : QUAD_TOP_RIGHT = 2,
797 : QUAD_BOTTOM_RIGHT = 3,
798 : };
799 :
800 310900 : constexpr int QUAD_COUNT = 4;
801 310900 : double adfQuadDist[QUAD_COUNT] = {};
802 310900 : float afQuadValue[QUAD_COUNT] = {};
803 :
804 1554500 : for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)
805 : {
806 1243600 : adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;
807 1243600 : afQuadValue[iQuad] = 0.0;
808 : }
809 :
810 : // Step left and right by one pixel searching for the closest
811 : // target value for each quadrant.
812 1558920 : for (int iStep = 0; iStep <= nThisMaxSearchDist; iStep++)
813 : {
814 1248020 : const int iLeftX = std::max(0, iX - iStep);
815 1248020 : const int iRightX = std::min(nXSize - 1, iX + iStep);
816 :
817 : // Top left includes current line.
818 1248020 : QUAD_CHECK(adfQuadDist[QUAD_TOP_LEFT],
819 : afQuadValue[QUAD_TOP_LEFT], iLeftX,
820 1248020 : panTopDownY[iLeftX], iX, iY, pafTopDownValue[iLeftX],
821 : nNoDataVal);
822 :
823 : // Bottom left.
824 1248020 : QUAD_CHECK(adfQuadDist[QUAD_BOTTOM_LEFT],
825 : afQuadValue[QUAD_BOTTOM_LEFT], iLeftX,
826 1248020 : panLastY[iLeftX], iX, iY, pafLastValue[iLeftX],
827 : nNoDataVal);
828 :
829 : // Top right and bottom right do no include center pixel.
830 1248020 : if (iStep == 0)
831 310900 : continue;
832 :
833 : // Top right includes current line.
834 937122 : QUAD_CHECK(adfQuadDist[QUAD_TOP_RIGHT],
835 : afQuadValue[QUAD_TOP_RIGHT], iRightX,
836 937122 : panTopDownY[iRightX], iX, iY,
837 937122 : pafTopDownValue[iRightX], nNoDataVal);
838 :
839 : // Bottom right.
840 937122 : QUAD_CHECK(adfQuadDist[QUAD_BOTTOM_RIGHT],
841 : afQuadValue[QUAD_BOTTOM_RIGHT], iRightX,
842 937122 : panLastY[iRightX], iX, iY, pafLastValue[iRightX],
843 : nNoDataVal);
844 :
845 : // Every four steps, recompute maximum distance.
846 937122 : if ((iStep & 0x3) == 0)
847 1171 : nThisMaxSearchDist = static_cast<int>(floor(
848 : std::max(std::max(adfQuadDist[0], adfQuadDist[1]),
849 1171 : std::max(adfQuadDist[2], adfQuadDist[3]))));
850 : }
851 :
852 310900 : bool bHasSrcValues = false;
853 310900 : if (bNearest)
854 : {
855 12 : double dfNearestDist = dfMaxSearchDist + 1;
856 12 : float fNearestValue = 0.0f;
857 :
858 60 : for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)
859 : {
860 48 : if (adfQuadDist[iQuad] < dfNearestDist)
861 : {
862 11 : bHasSrcValues = true;
863 11 : if (!bHasNoData || afQuadValue[iQuad] != fNoData)
864 : {
865 10 : fNearestValue = afQuadValue[iQuad];
866 10 : dfNearestDist = adfQuadDist[iQuad];
867 : }
868 : }
869 : }
870 :
871 12 : if (bHasSrcValues)
872 : {
873 10 : pabyFiltMask[iX] = 255;
874 10 : if (dfNearestDist <= dfMaxSearchDist)
875 : {
876 7 : pabyMask[iX] = 255;
877 7 : pafScanline[iX] = fNearestValue;
878 : }
879 : else
880 3 : pafScanline[iX] = fNoData;
881 : }
882 : }
883 : else
884 : {
885 310888 : double dfWeightSum = 0.0;
886 310888 : double dfValueSum = 0.0;
887 :
888 1554440 : for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)
889 : {
890 1243550 : if (adfQuadDist[iQuad] <= dfMaxSearchDist)
891 : {
892 132924 : bHasSrcValues = true;
893 132924 : if (!bHasNoData || afQuadValue[iQuad] != fNoData)
894 : {
895 132923 : const double dfWeight = 1.0 / adfQuadDist[iQuad];
896 132923 : dfWeightSum += dfWeight;
897 132923 : dfValueSum += afQuadValue[iQuad] * dfWeight;
898 : }
899 : }
900 : }
901 :
902 310888 : if (bHasSrcValues)
903 : {
904 56598 : pabyFiltMask[iX] = 255;
905 56598 : if (dfWeightSum > 0.0)
906 : {
907 56598 : pabyMask[iX] = 255;
908 56598 : pafScanline[iX] =
909 56598 : static_cast<float>(dfValueSum / dfWeightSum);
910 : }
911 : else
912 0 : pafScanline[iX] = fNoData;
913 : }
914 : }
915 : }
916 :
917 : /* --------------------------------------------------------------------
918 : */
919 : /* Write out the updated data and mask information. */
920 : /* --------------------------------------------------------------------
921 : */
922 2507 : eErr = GDALRasterIO(hTargetBand, GF_Write, 0, iY, nXSize, 1,
923 : pafScanline, nXSize, 1, GDT_Float32, 0, 0);
924 :
925 2507 : if (eErr != CE_None)
926 0 : break;
927 :
928 2507 : if (poTmpMaskDS != nullptr)
929 : {
930 : // Update (copy of) mask band when it has been provided by the
931 : // user
932 5 : eErr = GDALRasterIO(hMaskBand, GF_Write, 0, iY, nXSize, 1, pabyMask,
933 : nXSize, 1, GDT_Byte, 0, 0);
934 :
935 5 : if (eErr != CE_None)
936 0 : break;
937 : }
938 :
939 2507 : eErr = GDALRasterIO(hFiltMaskBand, GF_Write, 0, iY, nXSize, 1,
940 : pabyFiltMask, nXSize, 1, GDT_Byte, 0, 0);
941 :
942 2507 : if (eErr != CE_None)
943 0 : break;
944 :
945 : /* --------------------------------------------------------------------
946 : */
947 : /* Flip this/last buffers. */
948 : /* --------------------------------------------------------------------
949 : */
950 2507 : std::swap(pafThisValue, pafLastValue);
951 2507 : std::swap(panThisY, panLastY);
952 :
953 : /* --------------------------------------------------------------------
954 : */
955 : /* report progress. */
956 : /* --------------------------------------------------------------------
957 : */
958 2507 : if (!pfnProgress(
959 : dfProgressRatio *
960 2507 : (0.5 + 0.5 * (nYSize - iY) / static_cast<double>(nYSize)),
961 : "Filling...", pProgressArg))
962 : {
963 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
964 0 : eErr = CE_Failure;
965 : }
966 : }
967 :
968 : /* ==================================================================== */
969 : /* Now we will do iterative average filters over the */
970 : /* interpolated values to smooth things out and make linear */
971 : /* artifacts less obvious. */
972 : /* ==================================================================== */
973 117 : if (eErr == CE_None && nSmoothingIterations > 0)
974 : {
975 92 : if (poTmpMaskDS == nullptr)
976 : {
977 : // Force masks to be to flushed and recomputed when the user
978 : // didn't pass a user-provided hMaskBand, and we assigned it
979 : // to be the mask band of hTargetBand.
980 91 : GDALFlushRasterCache(hMaskBand);
981 : }
982 :
983 92 : void *pScaledProgress = GDALCreateScaledProgress(
984 : dfProgressRatio, 1.0, pfnProgress, pProgressArg);
985 :
986 92 : eErr = GDALMultiFilter(hTargetBand, hMaskBand, hFiltMaskBand,
987 : nSmoothingIterations, GDALScaledProgress,
988 : pScaledProgress);
989 :
990 92 : GDALDestroyScaledProgress(pScaledProgress);
991 : }
992 :
993 : /* -------------------------------------------------------------------- */
994 : /* Close and clean up temporary files. Free working buffers */
995 : /* -------------------------------------------------------------------- */
996 25 : end:
997 117 : CPLFree(panLastY);
998 117 : CPLFree(panThisY);
999 117 : CPLFree(panTopDownY);
1000 117 : CPLFree(pafLastValue);
1001 117 : CPLFree(pafThisValue);
1002 117 : CPLFree(pafTopDownValue);
1003 117 : CPLFree(pafScanline);
1004 117 : CPLFree(pabyMask);
1005 117 : CPLFree(pabyFiltMask);
1006 :
1007 117 : return eErr;
1008 : }
|
