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