1 /******************************************************************************
2 *
3 * Project: GDAL
4 * Purpose: Compute each pixel's proximity to a set of target pixels.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2008, Frank Warmerdam
9 * Copyright (c) 2009-2010, Even Rouault <even dot rouault at spatialys.com>
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 "cpl_port.h"
31 #include "gdal_alg.h"
32
33 #include <cmath>
34 #include <cstdlib>
35
36 #include <algorithm>
37
38 #include "cpl_conv.h"
39 #include "cpl_error.h"
40 #include "cpl_progress.h"
41 #include "cpl_string.h"
42 #include "cpl_vsi.h"
43 #include "gdal.h"
44
45 CPL_CVSID("$Id: gdalproximity.cpp b1c9c12ad373e40b955162b45d704070d4ebf7b0 2019-06-19 16:50:15 +0200 Even Rouault $")
46
47 static CPLErr
48 ProcessProximityLine( GInt32 *panSrcScanline, int *panNearX, int *panNearY,
49 int bForward, int iLine, int nXSize, double nMaxDist,
50 float *pafProximity, double *pdfSrcNoDataValue,
51 int nTargetValues, int *panTargetValues );
52
53 /************************************************************************/
54 /* GDALComputeProximity() */
55 /************************************************************************/
56
57 /**
58 Compute the proximity of all pixels in the image to a set of pixels in
59 the source image.
60
61 This function attempts to compute the proximity of all pixels in
62 the image to a set of pixels in the source image. The following
63 options are used to define the behavior of the function. By
64 default all non-zero pixels in hSrcBand will be considered the
65 "target", and all proximities will be computed in pixels. Note
66 that target pixels are set to the value corresponding to a distance
67 of zero.
68
69 The progress function args may be NULL or a valid progress reporting function
70 such as GDALTermProgress/NULL.
71
72 Options:
73
74 VALUES=n[,n]*
75
76 A list of target pixel values to measure the distance from. If this
77 option is not provided proximity will be computed from non-zero
78 pixel values. Currently pixel values are internally processed as
79 integers.
80
81 DISTUNITS=[PIXEL]/GEO
82
83 Indicates whether distances will be computed in pixel units or
84 in georeferenced units. The default is pixel units. This also
85 determines the interpretation of MAXDIST.
86
87 MAXDIST=n
88
89 The maximum distance to search. Proximity distances greater than
90 this value will not be computed. Instead output pixels will be
91 set to a nodata value.
92
93 NODATA=n
94
95 The NODATA value to use on the output band for pixels that are
96 beyond MAXDIST. If not provided, the hProximityBand will be
97 queried for a nodata value. If one is not found, 65535 will be used.
98
99 USE_INPUT_NODATA=YES/NO
100
101 If this option is set, the input data set no-data value will be
102 respected. Leaving no data pixels in the input as no data pixels in
103 the proximity output.
104
105 FIXED_BUF_VAL=n
106
107 If this option is set, all pixels within the MAXDIST threadhold are
108 set to this fixed value instead of to a proximity distance.
109 */
110
111 CPLErr CPL_STDCALL
GDALComputeProximity(GDALRasterBandH hSrcBand,GDALRasterBandH hProximityBand,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)112 GDALComputeProximity( GDALRasterBandH hSrcBand,
113 GDALRasterBandH hProximityBand,
114 char **papszOptions,
115 GDALProgressFunc pfnProgress,
116 void * pProgressArg )
117
118 {
119 VALIDATE_POINTER1( hSrcBand, "GDALComputeProximity", CE_Failure );
120 VALIDATE_POINTER1( hProximityBand, "GDALComputeProximity", CE_Failure );
121
122 if( pfnProgress == nullptr )
123 pfnProgress = GDALDummyProgress;
124
125 /* -------------------------------------------------------------------- */
126 /* Are we using pixels or georeferenced coordinates for distances? */
127 /* -------------------------------------------------------------------- */
128 double dfDistMult = 1.0;
129 const char *pszOpt = CSLFetchNameValue( papszOptions, "DISTUNITS" );
130 if( pszOpt )
131 {
132 if( EQUAL(pszOpt, "GEO") )
133 {
134 GDALDatasetH hSrcDS = GDALGetBandDataset( hSrcBand );
135 if( hSrcDS )
136 {
137 double adfGeoTransform[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
138
139 GDALGetGeoTransform( hSrcDS, adfGeoTransform );
140 if( std::abs(adfGeoTransform[1]) !=
141 std::abs(adfGeoTransform[5]) )
142 CPLError(
143 CE_Warning, CPLE_AppDefined,
144 "Pixels not square, distances will be inaccurate." );
145 dfDistMult = std::abs(adfGeoTransform[1]);
146 }
147 }
148 else if( !EQUAL(pszOpt, "PIXEL") )
149 {
150 CPLError(
151 CE_Failure, CPLE_AppDefined,
152 "Unrecognized DISTUNITS value '%s', should be GEO or PIXEL.",
153 pszOpt );
154 return CE_Failure;
155 }
156 }
157
158 /* -------------------------------------------------------------------- */
159 /* What is our maxdist value? */
160 /* -------------------------------------------------------------------- */
161 pszOpt = CSLFetchNameValue( papszOptions, "MAXDIST" );
162 const double dfMaxDist = pszOpt ?
163 CPLAtof(pszOpt) / dfDistMult :
164 GDALGetRasterBandXSize(hSrcBand) + GDALGetRasterBandYSize(hSrcBand);
165
166 CPLDebug( "GDAL", "MAXDIST=%g, DISTMULT=%g", dfMaxDist, dfDistMult );
167
168 /* -------------------------------------------------------------------- */
169 /* Verify the source and destination are compatible. */
170 /* -------------------------------------------------------------------- */
171 const int nXSize = GDALGetRasterBandXSize( hSrcBand );
172 const int nYSize = GDALGetRasterBandYSize( hSrcBand );
173 if( nXSize != GDALGetRasterBandXSize( hProximityBand )
174 || nYSize != GDALGetRasterBandYSize( hProximityBand ))
175 {
176 CPLError( CE_Failure, CPLE_AppDefined,
177 "Source and proximity bands are not the same size." );
178 return CE_Failure;
179 }
180
181 /* -------------------------------------------------------------------- */
182 /* Get input NODATA value. */
183 /* -------------------------------------------------------------------- */
184 double dfSrcNoDataValue = 0.0;
185 double *pdfSrcNoData = nullptr;
186 if( CPLFetchBool( papszOptions, "USE_INPUT_NODATA", false ) )
187 {
188 int bSrcHasNoData = 0;
189 dfSrcNoDataValue = GDALGetRasterNoDataValue( hSrcBand, &bSrcHasNoData );
190 if( bSrcHasNoData )
191 pdfSrcNoData = &dfSrcNoDataValue;
192 }
193
194 /* -------------------------------------------------------------------- */
195 /* Get output NODATA value. */
196 /* -------------------------------------------------------------------- */
197 float fNoDataValue = 0.0f;
198 pszOpt = CSLFetchNameValue( papszOptions, "NODATA" );
199 if( pszOpt != nullptr )
200 {
201 fNoDataValue = static_cast<float>(CPLAtof(pszOpt));
202 }
203 else
204 {
205 int bSuccess = FALSE;
206
207 fNoDataValue = static_cast<float>(
208 GDALGetRasterNoDataValue( hProximityBand, &bSuccess ) );
209 if( !bSuccess )
210 fNoDataValue = 65535.0;
211 }
212
213 /* -------------------------------------------------------------------- */
214 /* Is there a fixed value we wish to force the buffer area to? */
215 /* -------------------------------------------------------------------- */
216 double dfFixedBufVal = 0.0;
217 bool bFixedBufVal = false;
218 pszOpt = CSLFetchNameValue( papszOptions, "FIXED_BUF_VAL" );
219 if( pszOpt )
220 {
221 dfFixedBufVal = CPLAtof(pszOpt);
222 bFixedBufVal = true;
223 }
224
225 /* -------------------------------------------------------------------- */
226 /* Get the target value(s). */
227 /* -------------------------------------------------------------------- */
228 int *panTargetValues = nullptr;
229 int nTargetValues = 0;
230
231 pszOpt = CSLFetchNameValue( papszOptions, "VALUES" );
232 if( pszOpt != nullptr )
233 {
234 char **papszValuesTokens =
235 CSLTokenizeStringComplex( pszOpt, ",", FALSE, FALSE);
236
237 nTargetValues = CSLCount(papszValuesTokens);
238 panTargetValues = static_cast<int *>(
239 CPLCalloc(sizeof(int), nTargetValues) );
240
241 for( int i = 0; i < nTargetValues; i++ )
242 panTargetValues[i] = atoi(papszValuesTokens[i]);
243 CSLDestroy( papszValuesTokens );
244 }
245
246 /* -------------------------------------------------------------------- */
247 /* Initialize progress counter. */
248 /* -------------------------------------------------------------------- */
249 if( !pfnProgress( 0.0, "", pProgressArg ) )
250 {
251 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
252 CPLFree(panTargetValues);
253 return CE_Failure;
254 }
255
256 /* -------------------------------------------------------------------- */
257 /* We need a signed type for the working proximity values kept */
258 /* on disk. If our proximity band is not signed, then create a */
259 /* temporary file for this purpose. */
260 /* -------------------------------------------------------------------- */
261 GDALRasterBandH hWorkProximityBand = hProximityBand;
262 GDALDatasetH hWorkProximityDS = nullptr;
263 const GDALDataType eProxType = GDALGetRasterDataType(hProximityBand);
264 CPLErr eErr = CE_None;
265
266 // TODO(schwehr): Localize after removing gotos.
267 float *pafProximity = nullptr;
268 int *panNearX = nullptr;
269 int *panNearY = nullptr;
270 GInt32 *panSrcScanline = nullptr;
271 bool bTempFileAlreadyDeleted = false;
272
273 if( eProxType == GDT_Byte
274 || eProxType == GDT_UInt16
275 || eProxType == GDT_UInt32 )
276 {
277 GDALDriverH hDriver = GDALGetDriverByName("GTiff");
278 if( hDriver == nullptr )
279 {
280 CPLError( CE_Failure, CPLE_AppDefined,
281 "GDALComputeProximity needs GTiff driver" );
282 eErr = CE_Failure;
283 goto end;
284 }
285 CPLString osTmpFile = CPLGenerateTempFilename( "proximity" );
286 hWorkProximityDS =
287 GDALCreate( hDriver, osTmpFile,
288 nXSize, nYSize, 1, GDT_Float32, nullptr );
289 if( hWorkProximityDS == nullptr )
290 {
291 eErr = CE_Failure;
292 goto end;
293 }
294 // On Unix, attempt at deleting the temporary file now, so that
295 // if the process gets interrupted, it is automatically destroyed
296 // by the operating system.
297 bTempFileAlreadyDeleted = VSIUnlink( osTmpFile ) == 0;
298 hWorkProximityBand = GDALGetRasterBand( hWorkProximityDS, 1 );
299 }
300
301 /* -------------------------------------------------------------------- */
302 /* Allocate buffer for two scanlines of distances as floats */
303 /* (the current and last line). */
304 /* -------------------------------------------------------------------- */
305 pafProximity =
306 static_cast<float *>(VSI_MALLOC2_VERBOSE(sizeof(float), nXSize));
307 panNearX =
308 static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nXSize));
309 panNearY =
310 static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nXSize));
311 panSrcScanline =
312 static_cast<GInt32 *>(VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
313
314 if( pafProximity == nullptr
315 || panNearX == nullptr
316 || panNearY == nullptr
317 || panSrcScanline == nullptr)
318 {
319 eErr = CE_Failure;
320 goto end;
321 }
322
323 /* -------------------------------------------------------------------- */
324 /* Loop from top to bottom of the image. */
325 /* -------------------------------------------------------------------- */
326
327 for( int i = 0; i < nXSize; i++ )
328 {
329 panNearX[i] = -1;
330 panNearY[i] = -1;
331 }
332
333 for( int iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ )
334 {
335 // Read for target values.
336 eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,
337 panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 );
338 if( eErr != CE_None )
339 break;
340
341 for( int i = 0; i < nXSize; i++ )
342 pafProximity[i] = -1.0;
343
344 // Left to right.
345 ProcessProximityLine( panSrcScanline, panNearX, panNearY,
346 TRUE, iLine, nXSize, dfMaxDist, pafProximity,
347 pdfSrcNoData, nTargetValues, panTargetValues );
348
349 // Right to Left.
350 ProcessProximityLine( panSrcScanline, panNearX, panNearY,
351 FALSE, iLine, nXSize, dfMaxDist, pafProximity,
352 pdfSrcNoData, nTargetValues, panTargetValues );
353
354 // Write out results.
355 eErr =
356 GDALRasterIO( hWorkProximityBand, GF_Write, 0, iLine, nXSize, 1,
357 pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
358
359 if( eErr != CE_None )
360 break;
361
362 if( !pfnProgress( 0.5 * (iLine+1) / static_cast<double>(nYSize),
363 "", pProgressArg ) )
364 {
365 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
366 eErr = CE_Failure;
367 }
368 }
369
370 /* -------------------------------------------------------------------- */
371 /* Loop from bottom to top of the image. */
372 /* -------------------------------------------------------------------- */
373 for( int i = 0; i < nXSize; i++ )
374 {
375 panNearX[i] = -1;
376 panNearY[i] = -1;
377 }
378
379 for( int iLine = nYSize-1; eErr == CE_None && iLine >= 0; iLine-- )
380 {
381 // Read first pass proximity.
382 eErr =
383 GDALRasterIO( hWorkProximityBand, GF_Read, 0, iLine, nXSize, 1,
384 pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
385
386 if( eErr != CE_None )
387 break;
388
389 // Read pixel values.
390
391 eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,
392 panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 );
393 if( eErr != CE_None )
394 break;
395
396 // Right to left.
397 ProcessProximityLine( panSrcScanline, panNearX, panNearY,
398 FALSE, iLine, nXSize, dfMaxDist, pafProximity,
399 pdfSrcNoData, nTargetValues, panTargetValues );
400
401 // Left to right.
402 ProcessProximityLine( panSrcScanline, panNearX, panNearY,
403 TRUE, iLine, nXSize, dfMaxDist, pafProximity,
404 pdfSrcNoData, nTargetValues, panTargetValues );
405
406 // Final post processing of distances.
407 for( int i = 0; i < nXSize; i++ )
408 {
409 if( pafProximity[i] < 0.0 )
410 pafProximity[i] = fNoDataValue;
411 else if( pafProximity[i] > 0.0 )
412 {
413 if( bFixedBufVal )
414 pafProximity[i] = static_cast<float>( dfFixedBufVal );
415 else
416 pafProximity[i] =
417 static_cast<float>(pafProximity[i] * dfDistMult);
418 }
419 }
420
421 // Write out results.
422 eErr =
423 GDALRasterIO( hProximityBand, GF_Write, 0, iLine, nXSize, 1,
424 pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
425
426 if( eErr != CE_None )
427 break;
428
429 if( !pfnProgress( 0.5 +
430 0.5 * (nYSize-iLine) / static_cast<double>( nYSize ),
431 "", pProgressArg ) )
432 {
433 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
434 eErr = CE_Failure;
435 }
436 }
437
438 /* -------------------------------------------------------------------- */
439 /* Cleanup */
440 /* -------------------------------------------------------------------- */
441 end:
442 CPLFree( panNearX );
443 CPLFree( panNearY );
444 CPLFree( panSrcScanline );
445 CPLFree( pafProximity );
446 CPLFree( panTargetValues );
447
448 if( hWorkProximityDS != nullptr )
449 {
450 CPLString osProxFile = GDALGetDescription( hWorkProximityDS );
451 GDALClose( hWorkProximityDS );
452 if( !bTempFileAlreadyDeleted )
453 {
454 GDALDeleteDataset( GDALGetDriverByName( "GTiff" ), osProxFile );
455 }
456 }
457
458 return eErr;
459 }
460
461 /************************************************************************/
462 /* SquareDistance() */
463 /************************************************************************/
464
SquareDistance(double dfX1,double dfX2,double dfY1,double dfY2)465 static double SquareDistance(double dfX1, double dfX2,
466 double dfY1, double dfY2)
467 {
468 const double dfDX = dfX1 - dfX2;
469 const double dfDY = dfY1 - dfY2;
470 return dfDX * dfDX + dfDY * dfDY;
471 }
472
473 /************************************************************************/
474 /* ProcessProximityLine() */
475 /************************************************************************/
476
477 static CPLErr
ProcessProximityLine(GInt32 * panSrcScanline,int * panNearX,int * panNearY,int bForward,int iLine,int nXSize,double dfMaxDist,float * pafProximity,double * pdfSrcNoDataValue,int nTargetValues,int * panTargetValues)478 ProcessProximityLine( GInt32 *panSrcScanline, int *panNearX, int *panNearY,
479 int bForward, int iLine, int nXSize, double dfMaxDist,
480 float *pafProximity, double *pdfSrcNoDataValue,
481 int nTargetValues, int *panTargetValues )
482
483 {
484 const int iStart = bForward ? 0 : nXSize - 1;
485 const int iEnd = bForward ? nXSize : -1;
486 const int iStep = bForward ? 1 : -1;
487
488 for( int iPixel = iStart; iPixel != iEnd; iPixel += iStep )
489 {
490 bool bIsTarget = false;
491
492 /* -------------------------------------------------------------------- */
493 /* Is the current pixel a target pixel? */
494 /* -------------------------------------------------------------------- */
495 if( nTargetValues == 0 )
496 {
497 bIsTarget = panSrcScanline[iPixel] != 0;
498 }
499 else
500 {
501 for( int i = 0; i < nTargetValues; i++ )
502 {
503 if( panSrcScanline[iPixel] == panTargetValues[i] )
504 bIsTarget = TRUE;
505 }
506 }
507
508 if( bIsTarget )
509 {
510 pafProximity[iPixel] = 0.0;
511 panNearX[iPixel] = iPixel;
512 panNearY[iPixel] = iLine;
513 continue;
514 }
515
516 /* -------------------------------------------------------------------- */
517 /* Are we near(er) to the closest target to the above (below) */
518 /* pixel? */
519 /* -------------------------------------------------------------------- */
520 double dfNearDistSq =
521 std::max(dfMaxDist, static_cast<double>(nXSize)) *
522 std::max(dfMaxDist, static_cast<double>(nXSize)) * 2.0;
523
524 if( panNearX[iPixel] != -1 )
525 {
526 const double dfDistSq =
527 SquareDistance(panNearX[iPixel], iPixel,
528 panNearY[iPixel], iLine);
529
530 if( dfDistSq < dfNearDistSq )
531 {
532 dfNearDistSq = dfDistSq;
533 }
534 else
535 {
536 panNearX[iPixel] = -1;
537 panNearY[iPixel] = -1;
538 }
539 }
540
541 /* -------------------------------------------------------------------- */
542 /* Are we near(er) to the closest target to the left (right) */
543 /* pixel? */
544 /* -------------------------------------------------------------------- */
545 const int iLast = iPixel - iStep;
546
547 if( iPixel != iStart && panNearX[iLast] != -1 )
548 {
549 const double dfDistSq =
550 SquareDistance(panNearX[iLast], iPixel,
551 panNearY[iLast], iLine);
552
553 if( dfDistSq < dfNearDistSq )
554 {
555 dfNearDistSq = dfDistSq;
556 panNearX[iPixel] = panNearX[iLast];
557 panNearY[iPixel] = panNearY[iLast];
558 }
559 }
560
561 /* -------------------------------------------------------------------- */
562 /* Are we near(er) to the closest target to the topright */
563 /* (bottom left) pixel? */
564 /* -------------------------------------------------------------------- */
565 const int iTR = iPixel + iStep;
566
567 if( iTR != iEnd && panNearX[iTR] != -1 )
568 {
569 const double dfDistSq =
570 SquareDistance(panNearX[iTR], iPixel,
571 panNearY[iTR], iLine);
572
573 if( dfDistSq < dfNearDistSq )
574 {
575 dfNearDistSq = dfDistSq;
576 panNearX[iPixel] = panNearX[iTR];
577 panNearY[iPixel] = panNearY[iTR];
578 }
579 }
580
581 /* -------------------------------------------------------------------- */
582 /* Update our proximity value. */
583 /* -------------------------------------------------------------------- */
584 if( panNearX[iPixel] != -1
585 && (pdfSrcNoDataValue == nullptr
586 || panSrcScanline[iPixel] != *pdfSrcNoDataValue)
587 && dfNearDistSq <= dfMaxDist * dfMaxDist
588 && (pafProximity[iPixel] < 0
589 || dfNearDistSq < pafProximity[iPixel] * pafProximity[iPixel]) )
590 pafProximity[iPixel] = static_cast<float>(sqrt(dfNearDistSq));
591 }
592
593 return CE_None;
594 }
595