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