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