1 : /******************************************************************************
2 : * $Id: gdalproximity.cpp 21167 2010-11-24 15:19:51Z warmerdam $
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 21167 2010-11-24 15:19:51Z warmerdam $");
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 4 : GDALComputeProximity( GDALRasterBandH hSrcBand,
96 : GDALRasterBandH hProximityBand,
97 : char **papszOptions,
98 : GDALProgressFunc pfnProgress,
99 : void * pProgressArg )
100 :
101 : {
102 4 : int nXSize, nYSize, i, bFixedBufVal = FALSE;
103 : const char *pszOpt;
104 : double dfMaxDist;
105 4 : double dfFixedBufVal = 0.0;
106 :
107 4 : VALIDATE_POINTER1( hSrcBand, "GDALComputeProximity", CE_Failure );
108 4 : VALIDATE_POINTER1( hProximityBand, "GDALComputeProximity", CE_Failure );
109 :
110 4 : if( pfnProgress == NULL )
111 3 : pfnProgress = GDALDummyProgress;
112 :
113 : /* -------------------------------------------------------------------- */
114 : /* Are we using pixels or georeferenced coordinates for distances? */
115 : /* -------------------------------------------------------------------- */
116 4 : double dfDistMult = 1.0;
117 4 : pszOpt = CSLFetchNameValue( papszOptions, "DISTUNITS" );
118 4 : 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 4 : pszOpt = CSLFetchNameValue( papszOptions, "MAXDIST" );
147 4 : if( pszOpt )
148 2 : dfMaxDist = atof(pszOpt) / dfDistMult;
149 : else
150 2 : dfMaxDist = GDALGetRasterBandXSize(hSrcBand) + GDALGetRasterBandYSize(hSrcBand);
151 :
152 4 : CPLDebug( "GDAL", "MAXDIST=%g, DISTMULT=%g", dfMaxDist, dfDistMult );
153 :
154 : /* -------------------------------------------------------------------- */
155 : /* Verify the source and destination are compatible. */
156 : /* -------------------------------------------------------------------- */
157 4 : nXSize = GDALGetRasterBandXSize( hSrcBand );
158 4 : nYSize = GDALGetRasterBandYSize( hSrcBand );
159 4 : if( nXSize != GDALGetRasterBandXSize( hProximityBand )
160 : || nYSize != GDALGetRasterBandYSize( 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 4 : pszOpt = CSLFetchNameValue( papszOptions, "NODATA" );
172 4 : if( pszOpt != NULL )
173 2 : fNoDataValue = (float) atof(pszOpt);
174 : else
175 : {
176 : int bSuccess;
177 :
178 2 : fNoDataValue = (float) GDALGetRasterNoDataValue( hProximityBand, &bSuccess );
179 2 : if( !bSuccess )
180 2 : fNoDataValue = 65535.0;
181 : }
182 :
183 : /* -------------------------------------------------------------------- */
184 : /* Is there a fixed value we wish to force the buffer area to? */
185 : /* -------------------------------------------------------------------- */
186 4 : pszOpt = CSLFetchNameValue( papszOptions, "FIXED_BUF_VAL" );
187 4 : if( pszOpt )
188 : {
189 2 : dfFixedBufVal = atof(pszOpt);
190 2 : bFixedBufVal = TRUE;
191 : }
192 :
193 : /* -------------------------------------------------------------------- */
194 : /* Get the target value(s). */
195 : /* -------------------------------------------------------------------- */
196 4 : int *panTargetValues = NULL;
197 4 : int nTargetValues = 0;
198 :
199 4 : pszOpt = CSLFetchNameValue( papszOptions, "VALUES" );
200 4 : if( pszOpt != NULL )
201 : {
202 : char **papszValuesTokens;
203 :
204 2 : papszValuesTokens = CSLTokenizeStringComplex( pszOpt, ",", FALSE,FALSE);
205 :
206 2 : nTargetValues = CSLCount(papszValuesTokens);
207 2 : panTargetValues = (int *) CPLCalloc(sizeof(int),nTargetValues);
208 :
209 6 : for( i = 0; i < nTargetValues; i++ )
210 4 : panTargetValues[i] = atoi(papszValuesTokens[i]);
211 2 : CSLDestroy( papszValuesTokens );
212 : }
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Initialize progress counter. */
216 : /* -------------------------------------------------------------------- */
217 4 : 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 4 : GDALRasterBandH hWorkProximityBand = hProximityBand;
230 4 : GDALDatasetH hWorkProximityDS = NULL;
231 4 : GDALDataType eProxType = GDALGetRasterDataType( hProximityBand );
232 4 : int *panNearX = NULL, *panNearY = NULL;
233 4 : float *pafProximity = NULL;
234 4 : GInt32 *panSrcScanline = NULL;
235 : int iLine;
236 4 : CPLErr eErr = CE_None;
237 :
238 4 : if( eProxType == GDT_Byte
239 : || eProxType == GDT_UInt16
240 : || eProxType == GDT_UInt32 )
241 : {
242 2 : GDALDriverH hDriver = GDALGetDriverByName("GTiff");
243 2 : 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 2 : CPLString osTmpFile = CPLGenerateTempFilename( "proximity" );
251 : hWorkProximityDS =
252 : GDALCreate( hDriver, osTmpFile,
253 2 : nXSize, nYSize, 1, GDT_Float32, NULL );
254 2 : if (hWorkProximityDS == NULL)
255 : {
256 0 : eErr = CE_Failure;
257 : goto end;
258 : }
259 2 : 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 4 : pafProximity = (float *) VSIMalloc2(sizeof(float), nXSize);
267 4 : panNearX = (int *) VSIMalloc2(sizeof(int), nXSize);
268 4 : panNearY = (int *) VSIMalloc2(sizeof(int), nXSize);
269 4 : panSrcScanline = (GInt32 *) VSIMalloc2(sizeof(GInt32), nXSize);
270 :
271 4 : 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 104 : for( i = 0; i < nXSize; i++ )
287 100 : panNearX[i] = panNearY[i] = -1;
288 :
289 104 : 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 100 : panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 );
294 100 : if( eErr != CE_None )
295 0 : break;
296 :
297 2600 : for( i = 0; i < nXSize; i++ )
298 2500 : pafProximity[i] = -1.0;
299 :
300 : // Left to right
301 : ProcessProximityLine( panSrcScanline, panNearX, panNearY,
302 : TRUE, iLine, nXSize, dfMaxDist,
303 100 : pafProximity, nTargetValues, panTargetValues );
304 :
305 : // Right to Left
306 : ProcessProximityLine( panSrcScanline, panNearX, panNearY,
307 : FALSE, iLine, nXSize, dfMaxDist,
308 100 : pafProximity, nTargetValues, panTargetValues );
309 :
310 : // Write out results.
311 : eErr =
312 : GDALRasterIO( hWorkProximityBand, GF_Write, 0, iLine, nXSize, 1,
313 100 : pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
314 :
315 100 : if( eErr != CE_None )
316 0 : break;
317 :
318 100 : 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 104 : for( i = 0; i < nXSize; i++ )
330 100 : panNearX[i] = panNearY[i] = -1;
331 :
332 104 : 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 100 : pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
338 :
339 100 : if( eErr != CE_None )
340 0 : break;
341 :
342 : // Read pixel values.
343 :
344 : eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,
345 100 : panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 );
346 100 : if( eErr != CE_None )
347 0 : break;
348 :
349 : // Right to left
350 : ProcessProximityLine( panSrcScanline, panNearX, panNearY,
351 : FALSE, iLine, nXSize, dfMaxDist,
352 100 : pafProximity, nTargetValues, panTargetValues );
353 :
354 : // Left to right
355 : ProcessProximityLine( panSrcScanline, panNearX, panNearY,
356 : TRUE, iLine, nXSize, dfMaxDist,
357 100 : pafProximity, nTargetValues, panTargetValues );
358 :
359 : // Final post processing of distances.
360 2600 : for( i = 0; i < nXSize; i++ )
361 : {
362 2500 : if( pafProximity[i] < 0.0 )
363 662 : pafProximity[i] = fNoDataValue;
364 1838 : else if( pafProximity[i] > 0.0 )
365 : {
366 1594 : if( bFixedBufVal )
367 570 : pafProximity[i] = (float) dfFixedBufVal;
368 : else
369 1024 : pafProximity[i] = (float)(pafProximity[i] * dfDistMult);
370 : }
371 : }
372 :
373 : // Write out results.
374 : eErr =
375 : GDALRasterIO( hProximityBand, GF_Write, 0, iLine, nXSize, 1,
376 100 : pafProximity, nXSize, 1, GDT_Float32, 0, 0 );
377 :
378 100 : if( eErr != CE_None )
379 0 : break;
380 :
381 100 : 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 4 : CPLFree( panNearX );
394 4 : CPLFree( panNearY );
395 4 : CPLFree( panSrcScanline );
396 4 : CPLFree( pafProximity );
397 4 : CPLFree(panTargetValues);
398 :
399 4 : if( hWorkProximityDS != NULL )
400 : {
401 2 : CPLString osProxFile = GDALGetDescription( hWorkProximityDS );
402 2 : GDALClose( hWorkProximityDS );
403 2 : GDALDeleteDataset( GDALGetDriverByName( "GTiff" ), osProxFile );
404 : }
405 :
406 4 : return eErr;
407 : }
408 :
409 : /************************************************************************/
410 : /* ProcessProximityLine() */
411 : /************************************************************************/
412 :
413 : static CPLErr
414 400 : 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 400 : if( bForward )
423 : {
424 200 : iStart = 0;
425 200 : iEnd = nXSize;
426 200 : iStep = 1;
427 : }
428 : else
429 : {
430 200 : iStart = nXSize-1;
431 200 : iEnd = -1;
432 200 : iStep = -1;
433 : }
434 :
435 10400 : for( iPixel = iStart; iPixel != iEnd; iPixel += iStep )
436 : {
437 10000 : int bIsTarget = FALSE;
438 :
439 : /* -------------------------------------------------------------------- */
440 : /* Is the current pixel a target pixel? */
441 : /* -------------------------------------------------------------------- */
442 10000 : if( nTargetValues == 0 )
443 5000 : bIsTarget = (panSrcScanline[iPixel] != 0);
444 : else
445 : {
446 : int i;
447 :
448 15000 : for( i = 0; i < nTargetValues; i++ )
449 : {
450 10000 : if( panSrcScanline[iPixel] == panTargetValues[i] )
451 72 : bIsTarget = TRUE;
452 : }
453 : }
454 :
455 10000 : if( bIsTarget )
456 : {
457 976 : pafProximity[iPixel] = 0.0;
458 976 : panNearX[iPixel] = iPixel;
459 976 : panNearY[iPixel] = iLine;
460 976 : continue;
461 : }
462 :
463 : /* -------------------------------------------------------------------- */
464 : /* Are we near(er) to the closest target to the above (below) */
465 : /* pixel? */
466 : /* -------------------------------------------------------------------- */
467 9024 : float fNearDistSq = (float) (MAX(dfMaxDist,nXSize) * MAX(dfMaxDist,nXSize) * 2);
468 : float fDistSq;
469 :
470 9024 : if( panNearX[iPixel] != -1 )
471 : {
472 : fDistSq = (float)
473 12536 : ((panNearX[iPixel] - iPixel) * (panNearX[iPixel] - iPixel)
474 12536 : + (panNearY[iPixel] - iLine) * (panNearY[iPixel] - iLine));
475 :
476 6268 : if( fDistSq < fNearDistSq )
477 : {
478 6268 : 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 9024 : int iLast = iPixel-iStep;
492 :
493 9024 : if( iPixel != iStart && panNearX[iLast] != -1 )
494 : {
495 : fDistSq = (float)
496 12312 : ((panNearX[iLast] - iPixel) * (panNearX[iLast] - iPixel)
497 12312 : + (panNearY[iLast] - iLine) * (panNearY[iLast] - iLine));
498 :
499 6156 : if( fDistSq < fNearDistSq )
500 : {
501 1236 : fNearDistSq = fDistSq;
502 1236 : panNearX[iPixel] = panNearX[iLast];
503 1236 : 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 9024 : int iTR = iPixel+iStep;
512 :
513 9024 : if( iTR != iEnd && panNearX[iTR] != -1 )
514 : {
515 : fDistSq = (float)
516 11972 : ((panNearX[iTR] - iPixel) * (panNearX[iTR] - iPixel)
517 11972 : + (panNearY[iTR] - iLine) * (panNearY[iTR] - iLine));
518 :
519 5986 : if( fDistSq < fNearDistSq )
520 : {
521 52 : fNearDistSq = fDistSq;
522 52 : panNearX[iPixel] = panNearX[iTR];
523 52 : panNearY[iPixel] = panNearY[iTR];
524 : }
525 : }
526 :
527 : /* -------------------------------------------------------------------- */
528 : /* Update our proximity value. */
529 : /* -------------------------------------------------------------------- */
530 17606 : if( panNearX[iPixel] != -1
531 : && fNearDistSq <= dfMaxDist * dfMaxDist
532 5088 : && (pafProximity[iPixel] < 0
533 3494 : || fNearDistSq < pafProximity[iPixel] * pafProximity[iPixel]) )
534 2206 : pafProximity[iPixel] = sqrt(fNearDistSq);
535 : }
536 :
537 400 : return CE_None;
538 : }
|
