graydist

Gray-weighted distance transform of grayscale image

Description

T = graydist(I, mask)

Computes the gray-weighted distance transform of the grayscale image I. Locations where mask is true are seed locations.

T = graydist(I, C, R)

Uses vectors C and R to specify the column and row coordinates of seed locations.

T = graydist(I, ind)

Specifies the linear indices of seed locations using the vector ind.

T = graydist(I, ___, method)

Specifies an alternate distance metric. The method determines the chamfer weights that are assigned to the local neighborhood during outward propagation. Each pixel's contribution to the geodesic time is based on the chamfer weight in a particular direction multiplied by the pixel intensity.

Input Arguments

I — GMTimage/GMTgrid or Matrix grayscale image

Grayscale image, specified as a numeric or logical array. Can be of any numeric type and must be nonsparse.

Data Types: Float64, Float32, Int, UInt8, Bool

mask — Binary mask that specifies seed locations

Binary mask that specifies seed locations, specified as a logical array the same size as I. Locations where mask is true are treated as seed locations.

Data Types: Bool

C, R — Column and row coordinates of seed locations

Column and row coordinates of seed locations, specified as vectors of positive integers. Coordinate values must be valid subscripts in I.

Data Types: Vector{Int}

ind — Linear indices of seed locations

Linear indices of seed locations, specified as a vector of CartesianIndex or integer values that represent valid linear indices in I.

Data Types: Vector{Int} or Vector{CartesianIndex}

method — Distance metric

Distance metric, specified as one of the following values:

The method determines the chamfer weights assigned to the local neighborhood during outward propagation.

Output Arguments

T — Gray-weighted distance transform

Gray-weighted distance transform, returned as an array of the same size as I. Each element represents the gray-weighted geodesic distance from the corresponding pixel to the nearest seed location.

If the input numeric type of I is Float64, the output T is Float64. For any other numeric input type, the output T is Float32.

Data Types: Float64, Float32

Algorithm

graydist uses the geodesic time algorithm [1]. The basic equation for geodesic time along a path is:

τ_f(P) = f(p_0)/2 + f(p_l)/2 + Σ(i=1 to l-1) f(p_i)

where:

• τ_f(P) is the geodesic time along path P

• f(p_i) is the intensity value at point p_i

• The path P consists of points p_0, p_1, ..., p_l

The method parameter determines the chamfer weights assigned to the local neighborhood during outward propagation. Each pixel's contribution to the geodesic time is based on:

• The chamfer weight in a particular direction

• Multiplied by the pixel intensity

This approach finds minimum-cost paths through the image where the cost is weighted by pixel intensities.

Examples

Compute Minimum Path Through Magic Square

Create a magic square. Matrices generated by the magic function have equal row, column, and diagonal sums. The minimum path between the upper-left and lower-right corner is along the diagonal:

A = magic(3)

Output:

3×3 Matrix{Int64}:
 8  1  6
 3  5  7
 4  9  2

Calculate the gray-weighted distance transform, specifying the upper left corner and the lower right corner of the square as seed locations:

T1 = graydist(A, 1, 1)
T2 = graydist(A, 3, 3)

Output for T1:

3×3 Matrix{Float64}:
 0.0   4.5   8.0
 5.5   6.5   8.5
 9.0  11.5  10.0

Sum the two transforms to find the minimum path between the seed locations. As expected, there is a constant-value minimum path along the diagonal:

T = T1 + T2

Output:

3×3 Matrix{Float64}:
 13.0  14.5  18.0
 15.5  13.0  15.5
 18.0  14.5  13.0

Find Least-Cost Path in Terrain

This example demonstrates finding a path through steep terrain that minimizes the cumulative gradient:

# Load a DEM (Digital Elevation Model)
DEM = gmtread("terrain.tif")
imshow(DEM)

# Calculate gradient as weights for distance calculation
G = gradient(DEM)

# Define start and end points
start_mask = falses(size(DEM))
start_mask[100, 100] = true
end_mask = falses(size(DEM))
end_mask[500, 500] = true

# Compute gray-weighted distance from both points
T1 = graydist(G, start_mask)
T2 = graydist(G, end_mask)

# Find minimum path
T_total = T1 + T2

Using Column and Row Coordinates

Specify seed locations using coordinate vectors:

# Load grayscale image
I = gmtread("sample.png")

# Define multiple seed points
C = [10, 50, 100]  # Column coordinates
R = [20, 60, 120]  # Row coordinates

# Compute distance transform
T = graydist(I, C, R)
imshow(T)

Using Linear Indices

Specify seed locations using linear indices:

I = gmtread("sample.png")

# Find dark pixels as seed locations
dark_pixels = findall(I .< 50)
ind = [CartesianIndex(p) for p in dark_pixels]

# Compute distance transform from dark regions
T = graydist(I, ind)

References

[1] Soille, P. "Generalized geodesy via geodesic time." Pattern Recognition Letters. Vol. 15, December 1994, pp. 1235–1240.

Tips

• Use graydist to find least-cost paths through images where pixel values represent costs (e.g., elevation gradients, edge strengths)

• For routing through flat areas or topographic depressions in digital elevation models

• To create distance transforms that account for the underlying image structure

• Combine distance transforms from multiple seed points to find optimal paths between locations

See Also

bwdist, watershed