1 : /******************************************************************************
2 : * $Id: rasterfill.cpp 19220 2010-03-27 23:03:21Z rouault $
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 19220 2010-03-27 23:03:21Z rouault $");
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 : GDALFilterLine( float *pafLastLine, float *pafThisLine, float *pafNextLine,
46 : float *pafOutLine,
47 : GByte *pabyLastTMask, GByte *pabyThisTMask, GByte*pabyNextTMask,
48 0 : 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] = 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 : GDALMultiFilter( GDALRasterBandH hTargetBand,
145 : GDALRasterBandH hTargetMaskBand,
146 : GDALRasterBandH hFiltMaskBand,
147 : int nIterations,
148 : GDALProgressFunc pfnProgress,
149 0 : 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 0 : 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 : GDALFillNodata( GDALRasterBandH hTargetBand,
388 : GDALRasterBandH hMaskBand,
389 : double dfMaxSearchDist,
390 : int bDeprecatedOption,
391 : int nSmoothingIterations,
392 : char **papszOptions,
393 : GDALProgressFunc pfnProgress,
394 1 : void * pProgressArg )
395 :
396 : {
397 1 : VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );
398 :
399 1 : int nXSize = GDALGetRasterBandXSize( hTargetBand );
400 1 : int nYSize = GDALGetRasterBandYSize( hTargetBand );
401 1 : CPLErr eErr = CE_None;
402 :
403 : // Special "x" pixel values identifying pixels as special.
404 : GUInt32 nNoDataVal;
405 : GDALDataType eType;
406 :
407 1 : if( dfMaxSearchDist == 0.0 )
408 0 : dfMaxSearchDist = MAX(nXSize,nYSize) + 1;
409 :
410 1 : int nMaxSearchDist = (int) floor(dfMaxSearchDist);
411 :
412 1 : if( nXSize > 65533 || nYSize > 65533 )
413 : {
414 0 : eType = GDT_UInt32;
415 0 : nNoDataVal = 4000002;
416 : }
417 : else
418 : {
419 1 : eType = GDT_UInt16;
420 1 : nNoDataVal = 65535;
421 : }
422 :
423 1 : if( hMaskBand == NULL )
424 0 : hMaskBand = GDALGetMaskBand( hTargetBand );
425 :
426 : /* If there are smoothing iterations, reserve 10% of the progress for them */
427 1 : double dfProgressRatio = (nSmoothingIterations > 0) ? 0.9 : 1.0;
428 :
429 : /* -------------------------------------------------------------------- */
430 : /* Initialize progress counter. */
431 : /* -------------------------------------------------------------------- */
432 1 : if( pfnProgress == NULL )
433 0 : pfnProgress = GDALDummyProgress;
434 :
435 1 : 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 1 : GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
445 1 : 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", NULL };
455 1 : CPLString osTmpFile = CPLGenerateTempFilename("");
456 1 : CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";
457 :
458 : hYDS = GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,
459 1 : eType, (char **) apszOptions );
460 :
461 1 : if( hYDS == NULL )
462 1 : return CE_Failure;
463 :
464 1 : hYBand = GDALGetRasterBand( hYDS, 1 );
465 :
466 : /* -------------------------------------------------------------------- */
467 : /* Create a work file to hold the pixel value associated with */
468 : /* the "last xy value" pixel. */
469 : /* -------------------------------------------------------------------- */
470 : GDALDatasetH hValDS;
471 : GDALRasterBandH hValBand;
472 1 : CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";
473 :
474 : hValDS = GDALCreate( hDriver, osValTmpFile, nXSize, nYSize, 1,
475 : GDALGetRasterDataType( hTargetBand ),
476 1 : (char **) apszOptions );
477 :
478 1 : if( hValDS == NULL )
479 0 : return CE_Failure;
480 :
481 1 : hValBand = GDALGetRasterBand( hValDS, 1 );
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Create a mask file to make it clear what pixels can be filtered */
485 : /* on the filtering pass. */
486 : /* -------------------------------------------------------------------- */
487 : GDALDatasetH hFiltMaskDS;
488 : GDALRasterBandH hFiltMaskBand;
489 1 : CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";
490 :
491 : hFiltMaskDS =
492 : GDALCreate( hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1,
493 1 : GDT_Byte, (char **) apszOptions );
494 :
495 1 : if( hFiltMaskDS == NULL )
496 0 : return CE_Failure;
497 :
498 1 : hFiltMaskBand = GDALGetRasterBand( hFiltMaskDS, 1 );
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Allocate buffers for last scanline and this scanline. */
502 : /* -------------------------------------------------------------------- */
503 : GUInt32 *panLastY, *panThisY, *panTopDownY;
504 : float *pafLastValue, *pafThisValue, *pafScanline, *pafTopDownValue;
505 : GByte *pabyMask, *pabyFiltMask;
506 : int iX;
507 : int iY;
508 :
509 1 : panLastY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
510 1 : panThisY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
511 1 : panTopDownY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
512 1 : pafLastValue = (float *) VSICalloc(nXSize,sizeof(float));
513 1 : pafThisValue = (float *) VSICalloc(nXSize,sizeof(float));
514 1 : pafTopDownValue = (float *) VSICalloc(nXSize,sizeof(float));
515 1 : pafScanline = (float *) VSICalloc(nXSize,sizeof(float));
516 1 : pabyMask = (GByte *) VSICalloc(nXSize,1);
517 1 : pabyFiltMask = (GByte *) VSICalloc(nXSize,1);
518 1 : if (panLastY == NULL || panThisY == NULL || panTopDownY == NULL ||
519 : pafLastValue == NULL || pafThisValue == NULL || pafTopDownValue == NULL ||
520 : pafScanline == NULL || pabyMask == NULL || pabyFiltMask == NULL)
521 : {
522 : CPLError(CE_Failure, CPLE_OutOfMemory,
523 0 : "Could not allocate enough memory for temporary buffers");
524 :
525 0 : eErr = CE_Failure;
526 0 : goto end;
527 : }
528 :
529 21 : for( iX = 0; iX < nXSize; iX++ )
530 : {
531 20 : panLastY[iX] = nNoDataVal;
532 : }
533 :
534 : /* ==================================================================== */
535 : /* Make first pass from top to bottom collecting the "last */
536 : /* known value" for each column and writing it out to the work */
537 : /* files. */
538 : /* ==================================================================== */
539 :
540 21 : for( iY = 0; iY < nYSize && eErr == CE_None; iY++ )
541 : {
542 : /* -------------------------------------------------------------------- */
543 : /* Read data and mask for this line. */
544 : /* -------------------------------------------------------------------- */
545 : eErr =
546 : GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
547 20 : pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
548 :
549 20 : if( eErr != CE_None )
550 0 : break;
551 :
552 : eErr =
553 : GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
554 20 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
555 :
556 20 : if( eErr != CE_None )
557 0 : break;
558 :
559 : /* -------------------------------------------------------------------- */
560 : /* Figure out the most recent pixel for each column. */
561 : /* -------------------------------------------------------------------- */
562 :
563 420 : for( iX = 0; iX < nXSize; iX++ )
564 : {
565 400 : if( pabyMask[iX] )
566 : {
567 400 : pafThisValue[iX] = pafScanline[iX];
568 400 : panThisY[iX] = iY;
569 : }
570 0 : else if( iY - panLastY[iX] <= dfMaxSearchDist )
571 : {
572 0 : pafThisValue[iX] = pafLastValue[iX];
573 0 : panThisY[iX] = panLastY[iX];
574 : }
575 : else
576 : {
577 0 : panThisY[iX] = nNoDataVal;
578 : }
579 : }
580 :
581 : /* -------------------------------------------------------------------- */
582 : /* Write out best index/value to working files. */
583 : /* -------------------------------------------------------------------- */
584 : eErr = GDALRasterIO( hYBand, GF_Write, 0, iY, nXSize, 1,
585 20 : panThisY, nXSize, 1, GDT_UInt32, 0, 0 );
586 20 : if( eErr != CE_None )
587 0 : break;
588 :
589 : eErr = GDALRasterIO( hValBand, GF_Write, 0, iY, nXSize, 1,
590 20 : pafThisValue, nXSize, 1, GDT_Float32, 0, 0 );
591 20 : if( eErr != CE_None )
592 0 : break;
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Flip this/last buffers. */
596 : /* -------------------------------------------------------------------- */
597 : {
598 20 : float *pafTmp = pafThisValue;
599 20 : pafThisValue = pafLastValue;
600 20 : pafLastValue = pafTmp;
601 :
602 20 : GUInt32 *panTmp = panThisY;
603 20 : panThisY = panLastY;
604 20 : panLastY = panTmp;
605 : }
606 :
607 : /* -------------------------------------------------------------------- */
608 : /* report progress. */
609 : /* -------------------------------------------------------------------- */
610 20 : if( eErr == CE_None
611 : && !pfnProgress( dfProgressRatio * (0.5*(iY+1) / (double)nYSize),
612 : "Filling...", pProgressArg ) )
613 : {
614 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
615 0 : eErr = CE_Failure;
616 : }
617 : }
618 :
619 : /* ==================================================================== */
620 : /* Now we will do collect similar this/last information from */
621 : /* bottom to top and use it in combination with the top to */
622 : /* bottom search info to interpolate. */
623 : /* ==================================================================== */
624 21 : for( iY = nYSize-1; iY >= 0 && eErr == CE_None; iY-- )
625 : {
626 : eErr =
627 : GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
628 20 : pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
629 :
630 20 : if( eErr != CE_None )
631 0 : break;
632 :
633 : eErr =
634 : GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
635 20 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
636 :
637 20 : if( eErr != CE_None )
638 0 : break;
639 :
640 : /* -------------------------------------------------------------------- */
641 : /* Figure out the most recent pixel for each column. */
642 : /* -------------------------------------------------------------------- */
643 :
644 420 : for( iX = 0; iX < nXSize; iX++ )
645 : {
646 400 : if( pabyMask[iX] )
647 : {
648 400 : pafThisValue[iX] = pafScanline[iX];
649 400 : panThisY[iX] = iY;
650 : }
651 0 : else if( panLastY[iX] - iY <= dfMaxSearchDist )
652 : {
653 0 : pafThisValue[iX] = pafLastValue[iX];
654 0 : panThisY[iX] = panLastY[iX];
655 : }
656 : else
657 : {
658 0 : panThisY[iX] = nNoDataVal;
659 : }
660 : }
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* Load the last y and corresponding value from the top down pass. */
664 : /* -------------------------------------------------------------------- */
665 : eErr =
666 : GDALRasterIO( hYBand, GF_Read, 0, iY, nXSize, 1,
667 20 : panTopDownY, nXSize, 1, GDT_UInt32, 0, 0 );
668 :
669 20 : if( eErr != CE_None )
670 0 : break;
671 :
672 : eErr =
673 : GDALRasterIO( hValBand, GF_Read, 0, iY, nXSize, 1,
674 20 : pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0 );
675 :
676 20 : if( eErr != CE_None )
677 0 : break;
678 :
679 : /* -------------------------------------------------------------------- */
680 : /* Attempt to interpolate any pixels that are nodata. */
681 : /* -------------------------------------------------------------------- */
682 20 : memset( pabyFiltMask, 0, nXSize );
683 420 : for( iX = 0; iX < nXSize; iX++ )
684 : {
685 : int iStep, iQuad;
686 400 : int nThisMaxSearchDist = nMaxSearchDist;
687 :
688 : // If this was a valid target - no change.
689 400 : if( pabyMask[iX] )
690 400 : continue;
691 :
692 : // Quadrants 0:topleft, 1:bottomleft, 2:topright, 3:bottomright
693 : double adfQuadDist[4];
694 : double adfQuadValue[4];
695 :
696 0 : for( iQuad = 0; iQuad < 4; iQuad++ )
697 0 : adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;
698 :
699 : // Step left and right by one pixel searching for the closest
700 : // target value for each quadrant.
701 0 : for( iStep = 0; iStep < nThisMaxSearchDist; iStep++ )
702 : {
703 0 : int iLeftX = MAX(0,iX - iStep);
704 0 : int iRightX = MIN(nXSize-1,iX + iStep);
705 :
706 : // top left includes current line
707 0 : QUAD_CHECK(adfQuadDist[0],adfQuadValue[0],
708 : iLeftX, panTopDownY[iLeftX], iX, iY,
709 : pafTopDownValue[iLeftX] );
710 :
711 : // bottom left
712 0 : QUAD_CHECK(adfQuadDist[1],adfQuadValue[1],
713 : iLeftX, panLastY[iLeftX], iX, iY,
714 : pafLastValue[iLeftX] );
715 :
716 : // top right and bottom right do no include center pixel.
717 0 : if( iStep == 0 )
718 0 : continue;
719 :
720 : // top right includes current line
721 0 : QUAD_CHECK(adfQuadDist[2],adfQuadValue[2],
722 : iRightX, panTopDownY[iRightX], iX, iY,
723 : pafTopDownValue[iRightX] );
724 :
725 : // bottom right
726 0 : QUAD_CHECK(adfQuadDist[3],adfQuadValue[3],
727 : iRightX, panLastY[iRightX], iX, iY,
728 : pafLastValue[iRightX] );
729 :
730 : // every four steps, recompute maximum distance.
731 0 : if( (iStep & 0x3) == 0 )
732 : nThisMaxSearchDist = (int) floor(
733 0 : MAX(MAX(adfQuadDist[0],adfQuadDist[1]),
734 : MAX(adfQuadDist[2],adfQuadDist[3])) );
735 : }
736 :
737 0 : double dfWeightSum = 0.0;
738 0 : double dfValueSum = 0.0;
739 :
740 0 : for( iQuad = 0; iQuad < 4; iQuad++ )
741 : {
742 0 : if( adfQuadDist[iQuad] <= dfMaxSearchDist )
743 : {
744 0 : double dfWeight = 1.0 / adfQuadDist[iQuad];
745 :
746 0 : dfWeightSum += dfWeight;
747 0 : dfValueSum += adfQuadValue[iQuad] * dfWeight;
748 : }
749 : }
750 :
751 0 : if( dfWeightSum > 0.0 )
752 : {
753 0 : pabyMask[iX] = 255;
754 0 : pabyFiltMask[iX] = 255;
755 0 : pafScanline[iX] = dfValueSum / dfWeightSum;
756 : }
757 :
758 : }
759 :
760 : /* -------------------------------------------------------------------- */
761 : /* Write out the updated data and mask information. */
762 : /* -------------------------------------------------------------------- */
763 : eErr =
764 : GDALRasterIO( hTargetBand, GF_Write, 0, iY, nXSize, 1,
765 20 : pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
766 :
767 20 : if( eErr != CE_None )
768 0 : break;
769 :
770 : eErr =
771 : GDALRasterIO( hFiltMaskBand, GF_Write, 0, iY, nXSize, 1,
772 20 : pabyFiltMask, nXSize, 1, GDT_Byte, 0, 0 );
773 :
774 20 : if( eErr != CE_None )
775 0 : break;
776 :
777 : /* -------------------------------------------------------------------- */
778 : /* Flip this/last buffers. */
779 : /* -------------------------------------------------------------------- */
780 : {
781 20 : float *pafTmp = pafThisValue;
782 20 : pafThisValue = pafLastValue;
783 20 : pafLastValue = pafTmp;
784 :
785 20 : GUInt32 *panTmp = panThisY;
786 20 : panThisY = panLastY;
787 20 : panLastY = panTmp;
788 : }
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* report progress. */
792 : /* -------------------------------------------------------------------- */
793 20 : if( eErr == CE_None
794 : && !pfnProgress( dfProgressRatio*(0.5+0.5*(nYSize-iY) / (double)nYSize),
795 : "Filling...", pProgressArg ) )
796 : {
797 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
798 0 : eErr = CE_Failure;
799 : }
800 : }
801 :
802 : /* ==================================================================== */
803 : /* Now we will do iterative average filters over the */
804 : /* interpolated values to smooth things out and make linear */
805 : /* artifacts less obvious. */
806 : /* ==================================================================== */
807 1 : if( eErr == CE_None && nSmoothingIterations > 0 )
808 : {
809 : // force masks to be to flushed and recomputed.
810 0 : GDALFlushRasterCache( hMaskBand );
811 :
812 : void *pScaledProgress;
813 : pScaledProgress =
814 0 : GDALCreateScaledProgress( dfProgressRatio, 1.0, pfnProgress, NULL );
815 :
816 : eErr = GDALMultiFilter( hTargetBand, hMaskBand, hFiltMaskBand,
817 : nSmoothingIterations,
818 0 : GDALScaledProgress, pScaledProgress );
819 :
820 0 : GDALDestroyScaledProgress( pScaledProgress );
821 : }
822 :
823 : /* -------------------------------------------------------------------- */
824 : /* Close and clean up temporary files. Free working buffers */
825 : /* -------------------------------------------------------------------- */
826 1 : end:
827 1 : CPLFree(panLastY);
828 1 : CPLFree(panThisY);
829 1 : CPLFree(panTopDownY);
830 1 : CPLFree(pafLastValue);
831 1 : CPLFree(pafThisValue);
832 1 : CPLFree(pafTopDownValue);
833 1 : CPLFree(pafScanline);
834 1 : CPLFree(pabyMask);
835 1 : CPLFree(pabyFiltMask);
836 :
837 1 : GDALClose( hYDS );
838 1 : GDALClose( hValDS );
839 1 : GDALClose( hFiltMaskDS );
840 :
841 1 : GDALDeleteDataset( hDriver, osYTmpFile );
842 1 : GDALDeleteDataset( hDriver, osValTmpFile );
843 1 : GDALDeleteDataset( hDriver, osFiltMaskTmpFile );
844 :
845 1 : return eErr;
846 : }
|
