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