1 : /******************************************************************************
2 : * $Id: rasterfill.cpp 22395 2011-05-17 21:27:03Z warmerdam $
3 : *
4 : * Project: GDAL
5 : * Purpose: Interpolate in nodata areas.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2008, Frank Warmerdam
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 "gdal_alg.h"
31 : #include "cpl_conv.h"
32 : #include "cpl_string.h"
33 :
34 : CPL_CVSID("$Id: rasterfill.cpp 22395 2011-05-17 21:27:03Z warmerdam $");
35 :
36 : /************************************************************************/
37 : /* GDALFilterLine() */
38 : /* */
39 : /* Apply 3x3 filtering one one scanline with masking for which */
40 : /* pixels are to be interpolated (ThisFMask) and which window */
41 : /* pixels are valid to include in the interpolation (TMask). */
42 : /************************************************************************/
43 :
44 : static void
45 0 : GDALFilterLine( float *pafLastLine, float *pafThisLine, float *pafNextLine,
46 : float *pafOutLine,
47 : GByte *pabyLastTMask, GByte *pabyThisTMask, GByte*pabyNextTMask,
48 : GByte *pabyThisFMask, int nXSize )
49 :
50 : {
51 : int iX;
52 :
53 0 : for( iX = 0; iX < nXSize; iX++ )
54 : {
55 0 : if( !pabyThisFMask[iX] )
56 : {
57 0 : pafOutLine[iX] = pafThisLine[iX];
58 0 : continue;
59 : }
60 :
61 0 : CPLAssert( pabyThisTMask[iX] );
62 :
63 0 : double dfValSum = 0.0;
64 0 : double dfWeightSum = 0.0;
65 :
66 : // Previous line
67 0 : if( pafLastLine != NULL )
68 : {
69 0 : if( iX > 0 && pabyLastTMask[iX-1] )
70 : {
71 0 : dfValSum += pafLastLine[iX-1];
72 0 : dfWeightSum += 1.0;
73 : }
74 0 : if( pabyLastTMask[iX] )
75 : {
76 0 : dfValSum += pafLastLine[iX];
77 0 : dfWeightSum += 1.0;
78 : }
79 0 : if( iX < nXSize-1 && pabyLastTMask[iX+1] )
80 : {
81 0 : dfValSum += pafLastLine[iX+1];
82 0 : dfWeightSum += 1.0;
83 : }
84 : }
85 :
86 : // Current Line
87 0 : if( iX > 0 && pabyThisTMask[iX-1] )
88 : {
89 0 : dfValSum += pafThisLine[iX-1];
90 0 : dfWeightSum += 1.0;
91 : }
92 0 : if( pabyThisTMask[iX] )
93 : {
94 0 : dfValSum += pafThisLine[iX];
95 0 : dfWeightSum += 1.0;
96 : }
97 0 : if( iX < nXSize-1 && pabyThisTMask[iX+1] )
98 : {
99 0 : dfValSum += pafThisLine[iX+1];
100 0 : dfWeightSum += 1.0;
101 : }
102 :
103 : // Next line
104 0 : if( pafNextLine != NULL )
105 : {
106 0 : if( iX > 0 && pabyNextTMask[iX-1] )
107 : {
108 0 : dfValSum += pafNextLine[iX-1];
109 0 : dfWeightSum += 1.0;
110 : }
111 0 : if( pabyNextTMask[iX] )
112 : {
113 0 : dfValSum += pafNextLine[iX];
114 0 : dfWeightSum += 1.0;
115 : }
116 0 : if( iX < nXSize-1 && pabyNextTMask[iX+1] )
117 : {
118 0 : dfValSum += pafNextLine[iX+1];
119 0 : dfWeightSum += 1.0;
120 : }
121 : }
122 :
123 0 : pafOutLine[iX] = (float) (dfValSum / dfWeightSum);
124 : }
125 0 : }
126 :
127 : /************************************************************************/
128 : /* GDALMultiFilter() */
129 : /* */
130 : /* Apply multiple iterations of a 3x3 smoothing filter over a */
131 : /* band with masking controlling what pixels should be */
132 : /* filtered (FiltMaskBand non zero) and which pixels can be */
133 : /* considered valid contributors to the filter */
134 : /* (TargetMaskBand non zero). */
135 : /* */
136 : /* This implementation attempts to apply many iterations in */
137 : /* one IO pass by managing the filtering over a rolling buffer */
138 : /* of nIternations+2 scanlines. While possibly clever this */
139 : /* makes the algorithm implementation largely */
140 : /* incomprehensible. */
141 : /************************************************************************/
142 :
143 : static CPLErr
144 0 : GDALMultiFilter( GDALRasterBandH hTargetBand,
145 : GDALRasterBandH hTargetMaskBand,
146 : GDALRasterBandH hFiltMaskBand,
147 : int nIterations,
148 : GDALProgressFunc pfnProgress,
149 : void * pProgressArg )
150 :
151 : {
152 : float *paf3PassLineBuf;
153 : GByte *pabyTMaskBuf;
154 : GByte *pabyFMaskBuf;
155 : float *pafThisPass, *pafLastPass, *pafSLastPass;
156 :
157 0 : int nBufLines = nIterations + 2;
158 0 : int iPassCounter = 0;
159 : int nNewLine; // the line being loaded this time (zero based scanline)
160 0 : int nXSize = GDALGetRasterBandXSize( hTargetBand );
161 0 : int nYSize = GDALGetRasterBandYSize( hTargetBand );
162 0 : CPLErr eErr = CE_None;
163 :
164 : /* -------------------------------------------------------------------- */
165 : /* Report starting progress value. */
166 : /* -------------------------------------------------------------------- */
167 0 : 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 0 : pabyTMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
177 0 : pabyFMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
178 :
179 0 : paf3PassLineBuf = (float *) VSIMalloc3(nXSize, nBufLines, 3 * sizeof(float));
180 0 : if (pabyTMaskBuf == NULL || pabyFMaskBuf == NULL || paf3PassLineBuf == NULL)
181 : {
182 : CPLError(CE_Failure, CPLE_OutOfMemory,
183 0 : "Could not allocate enough memory for temporary buffers");
184 0 : eErr = CE_Failure;
185 0 : goto end;
186 : }
187 :
188 : /* -------------------------------------------------------------------- */
189 : /* Process rotating buffers. */
190 : /* -------------------------------------------------------------------- */
191 0 : for( nNewLine = 0;
192 : eErr == CE_None && nNewLine < nYSize+nIterations;
193 : nNewLine++ )
194 : {
195 : /* -------------------------------------------------------------------- */
196 : /* Rotate pass buffers. */
197 : /* -------------------------------------------------------------------- */
198 0 : iPassCounter = (iPassCounter + 1) % 3;
199 :
200 : pafSLastPass = paf3PassLineBuf
201 0 : + ((iPassCounter+0)%3) * nXSize*nBufLines;
202 : pafLastPass = paf3PassLineBuf
203 0 : + ((iPassCounter+1)%3) * nXSize*nBufLines;
204 : pafThisPass = paf3PassLineBuf
205 0 : + ((iPassCounter+2)%3) * nXSize*nBufLines;
206 :
207 : /* -------------------------------------------------------------------- */
208 : /* Where does the new line go in the rotating buffer? */
209 : /* -------------------------------------------------------------------- */
210 0 : int iBufOffset = nNewLine % nBufLines;
211 :
212 : /* -------------------------------------------------------------------- */
213 : /* Read the new data line if it is't off the bottom of the */
214 : /* image. */
215 : /* -------------------------------------------------------------------- */
216 0 : if( nNewLine < nYSize )
217 : {
218 : eErr =
219 : GDALRasterIO( hTargetMaskBand, GF_Read,
220 : 0, nNewLine, nXSize, 1,
221 : pabyTMaskBuf + nXSize * iBufOffset, nXSize, 1,
222 0 : GDT_Byte, 0, 0 );
223 :
224 0 : if( eErr != CE_None )
225 0 : break;
226 :
227 : eErr =
228 : GDALRasterIO( hFiltMaskBand, GF_Read,
229 : 0, nNewLine, nXSize, 1,
230 : pabyFMaskBuf + nXSize * iBufOffset, nXSize, 1,
231 0 : GDT_Byte, 0, 0 );
232 :
233 0 : if( eErr != CE_None )
234 0 : break;
235 :
236 : eErr =
237 : GDALRasterIO( hTargetBand, GF_Read,
238 : 0, nNewLine, nXSize, 1,
239 : pafThisPass + nXSize * iBufOffset, nXSize, 1,
240 0 : GDT_Float32, 0, 0 );
241 :
242 0 : if( eErr != CE_None )
243 0 : break;
244 : }
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* Loop over the loaded data, applying the filter to all loaded */
248 : /* lines with neighbours. */
249 : /* -------------------------------------------------------------------- */
250 : int iFLine;
251 :
252 0 : for( iFLine = nNewLine-1;
253 : eErr == CE_None && iFLine >= nNewLine-nIterations;
254 : iFLine-- )
255 : {
256 : int iLastOffset, iThisOffset, iNextOffset;
257 :
258 0 : iLastOffset = (iFLine-1) % nBufLines;
259 0 : iThisOffset = (iFLine ) % nBufLines;
260 0 : iNextOffset = (iFLine+1) % nBufLines;
261 :
262 : // default to preserving the old value.
263 0 : if( iFLine >= 0 )
264 : memcpy( pafThisPass + iThisOffset * nXSize,
265 : pafLastPass + iThisOffset * nXSize,
266 0 : sizeof(float) * nXSize );
267 :
268 : // currently this skips the first and last line. Eventually
269 : // we will enable these too. TODO
270 0 : if( iFLine < 1 || iFLine >= nYSize-1 )
271 : {
272 0 : continue;
273 : }
274 :
275 : GDALFilterLine(
276 : pafSLastPass + iLastOffset * nXSize,
277 : pafLastPass + iThisOffset * nXSize,
278 : pafThisPass + iNextOffset * nXSize,
279 : pafThisPass + iThisOffset * nXSize,
280 : pabyTMaskBuf + iLastOffset * nXSize,
281 : pabyTMaskBuf + iThisOffset * nXSize,
282 : pabyTMaskBuf + iNextOffset * nXSize,
283 : pabyFMaskBuf + iThisOffset * nXSize,
284 0 : nXSize );
285 : }
286 :
287 : /* -------------------------------------------------------------------- */
288 : /* Write out the top data line that will be rolling out of our */
289 : /* buffer. */
290 : /* -------------------------------------------------------------------- */
291 0 : int iLineToSave = nNewLine - nIterations;
292 :
293 0 : if( iLineToSave >= 0 && eErr == CE_None )
294 : {
295 0 : iBufOffset = iLineToSave % nBufLines;
296 :
297 : eErr =
298 : GDALRasterIO( hTargetBand, GF_Write,
299 : 0, iLineToSave, nXSize, 1,
300 : pafThisPass + nXSize * iBufOffset, nXSize, 1,
301 0 : GDT_Float32, 0, 0 );
302 : }
303 :
304 : /* -------------------------------------------------------------------- */
305 : /* Report progress. */
306 : /* -------------------------------------------------------------------- */
307 0 : if( eErr == CE_None
308 : && !pfnProgress( (nNewLine+1) / (double) (nYSize+nIterations),
309 : "Smoothing Filter...", pProgressArg ) )
310 : {
311 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
312 0 : eErr = CE_Failure;
313 : }
314 : }
315 :
316 : /* -------------------------------------------------------------------- */
317 : /* Cleanup */
318 : /* -------------------------------------------------------------------- */
319 : end:
320 0 : CPLFree( pabyTMaskBuf );
321 0 : CPLFree( pabyFMaskBuf );
322 0 : CPLFree( paf3PassLineBuf );
323 :
324 0 : return eErr;
325 : }
326 :
327 : /************************************************************************/
328 : /* QUAD_CHECK() */
329 : /* */
330 : /* macro for checking whether a point is nearer than the */
331 : /* existing closest point. */
332 : /************************************************************************/
333 : #define QUAD_CHECK(quad_dist, quad_value, \
334 : target_x, target_y, origin_x, origin_y, target_value ) \
335 : \
336 : if( quad_value != nNoDataVal ) \
337 : { \
338 : double dfDistSq = ((target_x-origin_x) * (target_x-origin_x)) \
339 : + ((target_y-origin_y) * (target_y-origin_y)); \
340 : \
341 : if( dfDistSq < quad_dist*quad_dist ) \
342 : { \
343 : CPLAssert( dfDistSq > 0.0 ); \
344 : quad_dist = sqrt(dfDistSq); \
345 : quad_value = target_value; \
346 : } \
347 : }
348 :
349 : /************************************************************************/
350 : /* GDALFillNodata() */
351 : /************************************************************************/
352 :
353 : /**
354 : * Fill selected raster regions by interpolation from the edges.
355 : *
356 : * This algorithm will interpolate values for all designated
357 : * nodata pixels (marked by zeros in hMaskBand). For each pixel
358 : * a four direction conic search is done to find values to interpolate
359 : * from (using inverse distance weighting). Once all values are
360 : * interpolated, zero or more smoothing iterations (3x3 average
361 : * filters on interpolated pixels) are applied to smooth out
362 : * artifacts.
363 : *
364 : * This algorithm is generally suitable for interpolating missing
365 : * regions of fairly continuously varying rasters (such as elevation
366 : * models for instance). It is also suitable for filling small holes
367 : * and cracks in more irregularly varying images (like airphotos). It
368 : * is generally not so great for interpolating a raster from sparse
369 : * point data - see the algorithms defined in gdal_grid.h for that case.
370 : *
371 : * @param hTargetBand the raster band to be modified in place.
372 : * @param hMaskBand a mask band indicating pixels to be interpolated (zero valued
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 (none
379 : * supported at this time - just pass NULL).
380 : * @param pfnProgress the progress function to report completion.
381 : * @param pProgressArg callback data for progress function.
382 : *
383 : * @return CE_None on success or CE_Failure if something goes wrong.
384 : */
385 :
386 : CPLErr CPL_STDCALL
387 2 : GDALFillNodata( GDALRasterBandH hTargetBand,
388 : GDALRasterBandH hMaskBand,
389 : double dfMaxSearchDist,
390 : int bDeprecatedOption,
391 : int nSmoothingIterations,
392 : char **papszOptions,
393 : GDALProgressFunc pfnProgress,
394 : void * pProgressArg )
395 :
396 : {
397 2 : VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );
398 :
399 2 : int nXSize = GDALGetRasterBandXSize( hTargetBand );
400 2 : int nYSize = GDALGetRasterBandYSize( hTargetBand );
401 2 : CPLErr eErr = CE_None;
402 :
403 : // Special "x" pixel values identifying pixels as special.
404 : GUInt32 nNoDataVal;
405 : GDALDataType eType;
406 :
407 2 : if( dfMaxSearchDist == 0.0 )
408 0 : dfMaxSearchDist = MAX(nXSize,nYSize) + 1;
409 :
410 2 : int nMaxSearchDist = (int) floor(dfMaxSearchDist);
411 :
412 2 : if( nXSize > 65533 || nYSize > 65533 )
413 : {
414 0 : eType = GDT_UInt32;
415 0 : nNoDataVal = 4000002;
416 : }
417 : else
418 : {
419 2 : eType = GDT_UInt16;
420 2 : nNoDataVal = 65535;
421 : }
422 :
423 2 : if( hMaskBand == NULL )
424 0 : hMaskBand = GDALGetMaskBand( hTargetBand );
425 :
426 : /* If there are smoothing iterations, reserve 10% of the progress for them */
427 2 : double dfProgressRatio = (nSmoothingIterations > 0) ? 0.9 : 1.0;
428 :
429 : /* -------------------------------------------------------------------- */
430 : /* Initialize progress counter. */
431 : /* -------------------------------------------------------------------- */
432 2 : if( pfnProgress == NULL )
433 0 : pfnProgress = GDALDummyProgress;
434 :
435 2 : if( !pfnProgress( 0.0, "Filling...", pProgressArg ) )
436 : {
437 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
438 0 : return CE_Failure;
439 : }
440 :
441 : /* -------------------------------------------------------------------- */
442 : /* Create a work file to hold the Y "last value" indices. */
443 : /* -------------------------------------------------------------------- */
444 2 : GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
445 2 : if (hDriver == NULL)
446 : {
447 : CPLError(CE_Failure, CPLE_AppDefined,
448 0 : "GDALFillNodata needs GTiff driver");
449 0 : return CE_Failure;
450 : }
451 :
452 : GDALDatasetH hYDS;
453 : GDALRasterBandH hYBand;
454 : static const char *apszOptions[] = { "COMPRESS=LZW", "BIGTIFF=IF_SAFER",
455 : NULL };
456 2 : CPLString osTmpFile = CPLGenerateTempFilename("");
457 2 : CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";
458 :
459 : hYDS = GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,
460 2 : eType, (char **) apszOptions );
461 :
462 2 : if( hYDS == NULL )
463 0 : return CE_Failure;
464 :
465 2 : hYBand = GDALGetRasterBand( hYDS, 1 );
466 :
467 : /* -------------------------------------------------------------------- */
468 : /* Create a work file to hold the pixel value associated with */
469 : /* the "last xy value" pixel. */
470 : /* -------------------------------------------------------------------- */
471 : GDALDatasetH hValDS;
472 : GDALRasterBandH hValBand;
473 2 : CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";
474 :
475 : hValDS = GDALCreate( hDriver, osValTmpFile, nXSize, nYSize, 1,
476 : GDALGetRasterDataType( hTargetBand ),
477 2 : (char **) apszOptions );
478 :
479 2 : if( hValDS == NULL )
480 0 : return CE_Failure;
481 :
482 2 : hValBand = GDALGetRasterBand( hValDS, 1 );
483 :
484 : /* -------------------------------------------------------------------- */
485 : /* Create a mask file to make it clear what pixels can be filtered */
486 : /* on the filtering pass. */
487 : /* -------------------------------------------------------------------- */
488 : GDALDatasetH hFiltMaskDS;
489 : GDALRasterBandH hFiltMaskBand;
490 2 : CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";
491 :
492 : hFiltMaskDS =
493 : GDALCreate( hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1,
494 2 : GDT_Byte, (char **) apszOptions );
495 :
496 2 : if( hFiltMaskDS == NULL )
497 0 : return CE_Failure;
498 :
499 2 : hFiltMaskBand = GDALGetRasterBand( hFiltMaskDS, 1 );
500 :
501 : /* -------------------------------------------------------------------- */
502 : /* Allocate buffers for last scanline and this scanline. */
503 : /* -------------------------------------------------------------------- */
504 : GUInt32 *panLastY, *panThisY, *panTopDownY;
505 : float *pafLastValue, *pafThisValue, *pafScanline, *pafTopDownValue;
506 : GByte *pabyMask, *pabyFiltMask;
507 : int iX;
508 : int iY;
509 :
510 2 : panLastY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
511 2 : panThisY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
512 2 : panTopDownY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
513 2 : pafLastValue = (float *) VSICalloc(nXSize,sizeof(float));
514 2 : pafThisValue = (float *) VSICalloc(nXSize,sizeof(float));
515 2 : pafTopDownValue = (float *) VSICalloc(nXSize,sizeof(float));
516 2 : pafScanline = (float *) VSICalloc(nXSize,sizeof(float));
517 2 : pabyMask = (GByte *) VSICalloc(nXSize,1);
518 2 : pabyFiltMask = (GByte *) VSICalloc(nXSize,1);
519 2 : if (panLastY == NULL || panThisY == NULL || panTopDownY == NULL ||
520 : pafLastValue == NULL || pafThisValue == NULL || pafTopDownValue == NULL ||
521 : pafScanline == NULL || pabyMask == NULL || pabyFiltMask == NULL)
522 : {
523 : CPLError(CE_Failure, CPLE_OutOfMemory,
524 0 : "Could not allocate enough memory for temporary buffers");
525 :
526 0 : eErr = CE_Failure;
527 0 : goto end;
528 : }
529 :
530 42 : for( iX = 0; iX < nXSize; iX++ )
531 : {
532 40 : panLastY[iX] = nNoDataVal;
533 : }
534 :
535 : /* ==================================================================== */
536 : /* Make first pass from top to bottom collecting the "last */
537 : /* known value" for each column and writing it out to the work */
538 : /* files. */
539 : /* ==================================================================== */
540 :
541 42 : for( iY = 0; iY < nYSize && eErr == CE_None; iY++ )
542 : {
543 : /* -------------------------------------------------------------------- */
544 : /* Read data and mask for this line. */
545 : /* -------------------------------------------------------------------- */
546 : eErr =
547 : GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
548 40 : pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
549 :
550 40 : if( eErr != CE_None )
551 0 : break;
552 :
553 : eErr =
554 : GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
555 40 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
556 :
557 40 : if( eErr != CE_None )
558 0 : break;
559 :
560 : /* -------------------------------------------------------------------- */
561 : /* Figure out the most recent pixel for each column. */
562 : /* -------------------------------------------------------------------- */
563 :
564 840 : for( iX = 0; iX < nXSize; iX++ )
565 : {
566 800 : if( pabyMask[iX] )
567 : {
568 800 : pafThisValue[iX] = pafScanline[iX];
569 800 : panThisY[iX] = iY;
570 : }
571 0 : else if( iY - panLastY[iX] <= dfMaxSearchDist )
572 : {
573 0 : pafThisValue[iX] = pafLastValue[iX];
574 0 : panThisY[iX] = panLastY[iX];
575 : }
576 : else
577 : {
578 0 : panThisY[iX] = nNoDataVal;
579 : }
580 : }
581 :
582 : /* -------------------------------------------------------------------- */
583 : /* Write out best index/value to working files. */
584 : /* -------------------------------------------------------------------- */
585 : eErr = GDALRasterIO( hYBand, GF_Write, 0, iY, nXSize, 1,
586 40 : panThisY, nXSize, 1, GDT_UInt32, 0, 0 );
587 40 : if( eErr != CE_None )
588 0 : break;
589 :
590 : eErr = GDALRasterIO( hValBand, GF_Write, 0, iY, nXSize, 1,
591 40 : pafThisValue, nXSize, 1, GDT_Float32, 0, 0 );
592 40 : if( eErr != CE_None )
593 0 : break;
594 :
595 : /* -------------------------------------------------------------------- */
596 : /* Flip this/last buffers. */
597 : /* -------------------------------------------------------------------- */
598 : {
599 40 : float *pafTmp = pafThisValue;
600 40 : pafThisValue = pafLastValue;
601 40 : pafLastValue = pafTmp;
602 :
603 40 : GUInt32 *panTmp = panThisY;
604 40 : panThisY = panLastY;
605 40 : panLastY = panTmp;
606 : }
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* report progress. */
610 : /* -------------------------------------------------------------------- */
611 40 : if( eErr == CE_None
612 : && !pfnProgress( dfProgressRatio * (0.5*(iY+1) / (double)nYSize),
613 : "Filling...", pProgressArg ) )
614 : {
615 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
616 0 : eErr = CE_Failure;
617 : }
618 : }
619 :
620 : /* ==================================================================== */
621 : /* Now we will do collect similar this/last information from */
622 : /* bottom to top and use it in combination with the top to */
623 : /* bottom search info to interpolate. */
624 : /* ==================================================================== */
625 42 : for( iY = nYSize-1; iY >= 0 && eErr == CE_None; iY-- )
626 : {
627 : eErr =
628 : GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
629 40 : pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
630 :
631 40 : if( eErr != CE_None )
632 0 : break;
633 :
634 : eErr =
635 : GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
636 40 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
637 :
638 40 : if( eErr != CE_None )
639 0 : break;
640 :
641 : /* -------------------------------------------------------------------- */
642 : /* Figure out the most recent pixel for each column. */
643 : /* -------------------------------------------------------------------- */
644 :
645 840 : for( iX = 0; iX < nXSize; iX++ )
646 : {
647 800 : if( pabyMask[iX] )
648 : {
649 800 : pafThisValue[iX] = pafScanline[iX];
650 800 : panThisY[iX] = iY;
651 : }
652 0 : else if( panLastY[iX] - iY <= dfMaxSearchDist )
653 : {
654 0 : pafThisValue[iX] = pafLastValue[iX];
655 0 : panThisY[iX] = panLastY[iX];
656 : }
657 : else
658 : {
659 0 : panThisY[iX] = nNoDataVal;
660 : }
661 : }
662 :
663 : /* -------------------------------------------------------------------- */
664 : /* Load the last y and corresponding value from the top down pass. */
665 : /* -------------------------------------------------------------------- */
666 : eErr =
667 : GDALRasterIO( hYBand, GF_Read, 0, iY, nXSize, 1,
668 40 : panTopDownY, nXSize, 1, GDT_UInt32, 0, 0 );
669 :
670 40 : if( eErr != CE_None )
671 0 : break;
672 :
673 : eErr =
674 : GDALRasterIO( hValBand, GF_Read, 0, iY, nXSize, 1,
675 40 : pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0 );
676 :
677 40 : if( eErr != CE_None )
678 0 : break;
679 :
680 : /* -------------------------------------------------------------------- */
681 : /* Attempt to interpolate any pixels that are nodata. */
682 : /* -------------------------------------------------------------------- */
683 40 : memset( pabyFiltMask, 0, nXSize );
684 840 : for( iX = 0; iX < nXSize; iX++ )
685 : {
686 : int iStep, iQuad;
687 800 : int nThisMaxSearchDist = nMaxSearchDist;
688 :
689 : // If this was a valid target - no change.
690 800 : if( pabyMask[iX] )
691 800 : continue;
692 :
693 : // Quadrants 0:topleft, 1:bottomleft, 2:topright, 3:bottomright
694 : double adfQuadDist[4];
695 : double adfQuadValue[4];
696 :
697 0 : for( iQuad = 0; iQuad < 4; iQuad++ )
698 0 : adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;
699 :
700 : // Step left and right by one pixel searching for the closest
701 : // target value for each quadrant.
702 0 : for( iStep = 0; iStep < nThisMaxSearchDist; iStep++ )
703 : {
704 0 : int iLeftX = MAX(0,iX - iStep);
705 0 : int iRightX = MIN(nXSize-1,iX + iStep);
706 :
707 : // top left includes current line
708 0 : QUAD_CHECK(adfQuadDist[0],adfQuadValue[0],
709 : iLeftX, panTopDownY[iLeftX], iX, iY,
710 : pafTopDownValue[iLeftX] );
711 :
712 : // bottom left
713 0 : QUAD_CHECK(adfQuadDist[1],adfQuadValue[1],
714 : iLeftX, panLastY[iLeftX], iX, iY,
715 : pafLastValue[iLeftX] );
716 :
717 : // top right and bottom right do no include center pixel.
718 0 : if( iStep == 0 )
719 0 : continue;
720 :
721 : // top right includes current line
722 0 : QUAD_CHECK(adfQuadDist[2],adfQuadValue[2],
723 : iRightX, panTopDownY[iRightX], iX, iY,
724 : pafTopDownValue[iRightX] );
725 :
726 : // bottom right
727 0 : QUAD_CHECK(adfQuadDist[3],adfQuadValue[3],
728 : iRightX, panLastY[iRightX], iX, iY,
729 : pafLastValue[iRightX] );
730 :
731 : // every four steps, recompute maximum distance.
732 0 : if( (iStep & 0x3) == 0 )
733 : nThisMaxSearchDist = (int) floor(
734 0 : MAX(MAX(adfQuadDist[0],adfQuadDist[1]),
735 0 : MAX(adfQuadDist[2],adfQuadDist[3])) );
736 : }
737 :
738 0 : double dfWeightSum = 0.0;
739 0 : double dfValueSum = 0.0;
740 :
741 0 : for( iQuad = 0; iQuad < 4; iQuad++ )
742 : {
743 0 : if( adfQuadDist[iQuad] <= dfMaxSearchDist )
744 : {
745 0 : double dfWeight = 1.0 / adfQuadDist[iQuad];
746 :
747 0 : dfWeightSum += dfWeight;
748 0 : dfValueSum += adfQuadValue[iQuad] * dfWeight;
749 : }
750 : }
751 :
752 0 : if( dfWeightSum > 0.0 )
753 : {
754 0 : pabyMask[iX] = 255;
755 0 : pabyFiltMask[iX] = 255;
756 0 : pafScanline[iX] = (float) (dfValueSum / dfWeightSum);
757 : }
758 :
759 : }
760 :
761 : /* -------------------------------------------------------------------- */
762 : /* Write out the updated data and mask information. */
763 : /* -------------------------------------------------------------------- */
764 : eErr =
765 : GDALRasterIO( hTargetBand, GF_Write, 0, iY, nXSize, 1,
766 40 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
767 :
768 40 : if( eErr != CE_None )
769 0 : break;
770 :
771 : eErr =
772 : GDALRasterIO( hFiltMaskBand, GF_Write, 0, iY, nXSize, 1,
773 40 : pabyFiltMask, nXSize, 1, GDT_Byte, 0, 0 );
774 :
775 40 : if( eErr != CE_None )
776 0 : break;
777 :
778 : /* -------------------------------------------------------------------- */
779 : /* Flip this/last buffers. */
780 : /* -------------------------------------------------------------------- */
781 : {
782 40 : float *pafTmp = pafThisValue;
783 40 : pafThisValue = pafLastValue;
784 40 : pafLastValue = pafTmp;
785 :
786 40 : GUInt32 *panTmp = panThisY;
787 40 : panThisY = panLastY;
788 40 : panLastY = panTmp;
789 : }
790 :
791 : /* -------------------------------------------------------------------- */
792 : /* report progress. */
793 : /* -------------------------------------------------------------------- */
794 40 : if( eErr == CE_None
795 : && !pfnProgress( dfProgressRatio*(0.5+0.5*(nYSize-iY) / (double)nYSize),
796 : "Filling...", pProgressArg ) )
797 : {
798 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
799 0 : eErr = CE_Failure;
800 : }
801 : }
802 :
803 : /* ==================================================================== */
804 : /* Now we will do iterative average filters over the */
805 : /* interpolated values to smooth things out and make linear */
806 : /* artifacts less obvious. */
807 : /* ==================================================================== */
808 2 : if( eErr == CE_None && nSmoothingIterations > 0 )
809 : {
810 : // force masks to be to flushed and recomputed.
811 0 : GDALFlushRasterCache( hMaskBand );
812 :
813 : void *pScaledProgress;
814 : pScaledProgress =
815 0 : GDALCreateScaledProgress( dfProgressRatio, 1.0, pfnProgress, NULL );
816 :
817 : eErr = GDALMultiFilter( hTargetBand, hMaskBand, hFiltMaskBand,
818 : nSmoothingIterations,
819 0 : GDALScaledProgress, pScaledProgress );
820 :
821 0 : GDALDestroyScaledProgress( pScaledProgress );
822 : }
823 :
824 : /* -------------------------------------------------------------------- */
825 : /* Close and clean up temporary files. Free working buffers */
826 : /* -------------------------------------------------------------------- */
827 : end:
828 2 : CPLFree(panLastY);
829 2 : CPLFree(panThisY);
830 2 : CPLFree(panTopDownY);
831 2 : CPLFree(pafLastValue);
832 2 : CPLFree(pafThisValue);
833 2 : CPLFree(pafTopDownValue);
834 2 : CPLFree(pafScanline);
835 2 : CPLFree(pabyMask);
836 2 : CPLFree(pabyFiltMask);
837 :
838 2 : GDALClose( hYDS );
839 2 : GDALClose( hValDS );
840 2 : GDALClose( hFiltMaskDS );
841 :
842 2 : GDALDeleteDataset( hDriver, osYTmpFile );
843 2 : GDALDeleteDataset( hDriver, osValTmpFile );
844 2 : GDALDeleteDataset( hDriver, osFiltMaskTmpFile );
845 :
846 2 : return eErr;
847 : }
|
