User Guide

This guide covers both basic and advanced usage of DistanceTransforms.jl.

Getting Started

Installation

using Pkg
Pkg.add("DistanceTransforms")

Basic Usage

The primary function in DistanceTransforms.jl is transform. This function processes an array of 0s and 1s, converting each background element (0) into a value representing its squared Euclidean distance to the nearest foreground element (1).

using DistanceTransforms: transform, boolean_indicator
using CairoMakie: Figure, Axis, heatmap!, hidedecorations!

# Create a random binary array
arr = rand([0, 1], 10, 10)

# Apply distance transform
transformed = transform(boolean_indicator(arr))

# Create visualization
fig = Figure(size = (800, 400))
ax1 = Axis(fig[1, 1], title = "Original")
ax2 = Axis(fig[1, 2], title = "Distance Transform")
heatmap!(ax1, arr, colormap = :grays)
heatmap!(ax2, transformed, colormap = :grays)
fig

Real-World Example

Let’s apply a distance transform to a real image:

using Images

# Download and load example image
img = load(download("http://docs.opencv.org/3.1.0/water_coins.jpg"))

# Convert to binary image
img_bw = Gray.(img) .> 0.5

# Apply distance transform
img_tfm = transform(boolean_indicator(img_bw))

# Visualize
fig = Figure(size = (900, 300))
ax1 = Axis(fig[1, 1], title = "Original Image")
ax2 = Axis(fig[1, 2], title = "Segmented Image")
ax3 = Axis(fig[1, 3], title = "Distance Transform")
heatmap!(ax1, rotr90(img), colormap = :grays)
heatmap!(ax2, rotr90(img_bw), colormap = :grays)
heatmap!(ax3, rotr90(img_tfm), colormap = :grays)
hidedecorations!.([ax1, ax2, ax3])
fig

Understanding Euclidean Distance

The library, by default, returns the squared Euclidean distance. If you need the true Euclidean distance, you can take the square root of each element:

# Create sample binary array
array2 = [
    0 1 1 0 1
    0 0 0 1 0
    1 1 0 0 0
]

# Convert to boolean indicator
array2_bool = boolean_indicator(array2)

# Apply squared Euclidean distance transform
sq_euc_transform = transform(array2_bool)

# Convert to true Euclidean distance
euc_transform = sqrt.(sq_euc_transform)

# Display results
println("Squared Euclidean Distance:")
display(sq_euc_transform)
println("\nEuclidean Distance:")
display(euc_transform)

Squared Euclidean Distance:

Euclidean Distance:

3×5 Matrix{Float32}:
 1.0  0.0  0.0  1.0  0.0
 1.0  1.0  1.0  0.0  1.0
 0.0  0.0  1.0  1.0  2.0

3×5 Matrix{Float32}:
 1.0  0.0  0.0  1.0  0.0
 1.0  1.0  1.0  0.0  1.0
 0.0  0.0  1.0  1.0  1.41421

Comparison with ImageMorphology.jl

using ImageMorphology: distance_transform, feature_transform

# Apply ImageMorphology distance transform
euc_transform2 = distance_transform(feature_transform(Bool.(array2)))

# Compare results
println("ImageMorphology result:")
display(euc_transform2)
println("\nAre the results approximately equal?")
isapprox(euc_transform2, euc_transform; rtol = 1e-2)

ImageMorphology result:

Are the results approximately equal?

3×5 Matrix{Float64}:
 1.0  0.0  0.0  1.0  0.0
 1.0  1.0  1.0  0.0  1.0
 0.0  0.0  1.0  1.0  1.41421

true

Advanced Features

Multi-threading

DistanceTransforms.jl efficiently utilizes multi-threading, particularly in its Felzenszwalb distance transform algorithm.

⚠️ Julia might only load with 1 thread which makes the actual benchmark below meaningless and potentially confusing. Use this for understanding how to use the threaded kwarg

using BenchmarkTools
using Base.Threads: nthreads

# Create a random binary array
x = boolean_indicator(rand([0, 1], 1000, 1000))

# Single-threaded benchmark
single_threaded = @benchmark transform($x; threaded = false)

# Multi-threaded benchmark
multi_threaded = @benchmark transform($x; threaded = true)

# Display results
println("Number of threads: $(nthreads())")
println("Single-threaded minimum time: $(minimum(single_threaded).time / 1e6) ms")
println("Multi-threaded minimum time: $(minimum(multi_threaded).time / 1e6) ms")
println("Speedup factor: $(minimum(single_threaded).time / minimum(multi_threaded).time)")

Number of threads: 1
Single-threaded minimum time: 17.071291 ms
Multi-threaded minimum time: 17.319625 ms
Speedup factor: 0.9856616987954416

GPU Acceleration

DistanceTransforms.jl extends its performance capabilities with GPU acceleration. The library uses Julia’s multiple dispatch to automatically leverage GPU resources when available.

CUDA Example

using CUDA
using DistanceTransforms

# Create a random array on GPU
x_gpu = CUDA.rand([0, 1], 1000, 1000)
x_gpu = boolean_indicator(x_gpu)

# The transform function automatically uses GPU
result_gpu = transform(x_gpu)

# Transfer result back to CPU if needed
result_cpu = Array(result_gpu)

Metal Example

using Metal
using DistanceTransforms

# Create a random array on GPU
x_gpu = Metal.rand([0, 1], 1000, 1000)
x_gpu = boolean_indicator(x_gpu)

# The transform function automatically uses GPU
result_gpu = transform(x_gpu)

# Transfer result back to CPU if needed
result_cpu = Array(result_gpu)

Performance Benchmarks

Performance comparison across different implementations for 2D and 3D distance transforms:

As shown in the graph, the GPU implementation demonstrates superior performance, especially for larger arrays.

Best Practices

For small arrays (<100x100): CPU with multi-threading is often sufficient
For medium arrays: Multi-threaded CPU may be faster than GPU due to lower overhead
For large arrays (>1000x1000): GPU acceleration provides the best performance
For 3D data: GPU acceleration is strongly recommended due to the computational complexity

Algorithm Details

On the CPU, DistanceTransforms.jl uses the squared Euclidean distance transform algorithm by Felzenszwalb and Huttenlocher, known for its accuracy and efficiency. On the GPU, thanks to the amazing AcceleratedKernels.jl and KernelAbstraction.jl packages, the algorithm is essentially identical.

--- title: "User Guide" sidebar: julia format: html: toc: true execute: engine: julia julia: threads: 4 --- This guide covers both basic and advanced usage of DistanceTransforms.jl. ## Getting Started ### Installation ```julia using Pkg Pkg.add("DistanceTransforms") ``` ### Basic Usage The primary function in DistanceTransforms.jl is `transform`. This function processes an array of 0s and 1s, converting each background element (0) into a value representing its squared Euclidean distance to the nearest foreground element (1). ```{julia} using DistanceTransforms: transform, boolean_indicator using CairoMakie: Figure, Axis, heatmap!, hidedecorations! # Create a random binary array arr = rand([0, 1], 10, 10) # Apply distance transform transformed = transform(boolean_indicator(arr)) # Create visualization fig = Figure(size = (800, 400)) ax1 = Axis(fig[1, 1], title = "Original") ax2 = Axis(fig[1, 2], title = "Distance Transform") heatmap!(ax1, arr, colormap = :grays) heatmap!(ax2, transformed, colormap = :grays) fig ``` ### Real-World Example Let's apply a distance transform to a real image: ```{julia} using Images # Download and load example image img = load(download("http://docs.opencv.org/3.1.0/water_coins.jpg")) # Convert to binary image img_bw = Gray.(img) .> 0.5 # Apply distance transform img_tfm = transform(boolean_indicator(img_bw)) # Visualize fig = Figure(size = (900, 300)) ax1 = Axis(fig[1, 1], title = "Original Image") ax2 = Axis(fig[1, 2], title = "Segmented Image") ax3 = Axis(fig[1, 3], title = "Distance Transform") heatmap!(ax1, rotr90(img), colormap = :grays) heatmap!(ax2, rotr90(img_bw), colormap = :grays) heatmap!(ax3, rotr90(img_tfm), colormap = :grays) hidedecorations!.([ax1, ax2, ax3]) fig ``` ### Understanding Euclidean Distance The library, by default, returns the squared Euclidean distance. If you need the true Euclidean distance, you can take the square root of each element: ```{julia} # Create sample binary array array2 = [ 0 1 1 0 1 0 0 0 1 0 1 1 0 0 0 ] # Convert to boolean indicator array2_bool = boolean_indicator(array2) # Apply squared Euclidean distance transform sq_euc_transform = transform(array2_bool) # Convert to true Euclidean distance euc_transform = sqrt.(sq_euc_transform) # Display results println("Squared Euclidean Distance:") display(sq_euc_transform) println("\nEuclidean Distance:") display(euc_transform) ``` ### Comparison with ImageMorphology.jl ```{julia} using ImageMorphology: distance_transform, feature_transform # Apply ImageMorphology distance transform euc_transform2 = distance_transform(feature_transform(Bool.(array2))) # Compare results println("ImageMorphology result:") display(euc_transform2) println("\nAre the results approximately equal?") isapprox(euc_transform2, euc_transform; rtol = 1e-2) ``` ## Advanced Features ### Multi-threading DistanceTransforms.jl efficiently utilizes multi-threading, particularly in its Felzenszwalb distance transform algorithm. ⚠️ Julia might only load with `1 thread` which makes the actual benchmark below meaningless and potentially confusing. Use this for understanding how to use the `threaded` kwarg ```{julia} using BenchmarkTools using Base.Threads: nthreads # Create a random binary array x = boolean_indicator(rand([0, 1], 1000, 1000)) # Single-threaded benchmark single_threaded = @benchmark transform($x; threaded = false) # Multi-threaded benchmark multi_threaded = @benchmark transform($x; threaded = true) # Display results println("Number of threads: $(nthreads())") println("Single-threaded minimum time: $(minimum(single_threaded).time / 1e6) ms") println("Multi-threaded minimum time: $(minimum(multi_threaded).time / 1e6) ms") println("Speedup factor: $(minimum(single_threaded).time / minimum(multi_threaded).time)") ``` ### GPU Acceleration DistanceTransforms.jl extends its performance capabilities with GPU acceleration. The library uses Julia's multiple dispatch to automatically leverage GPU resources when available. #### CUDA Example ```julia using CUDA using DistanceTransforms # Create a random array on GPU x_gpu = CUDA.rand([0, 1], 1000, 1000) x_gpu = boolean_indicator(x_gpu) # The transform function automatically uses GPU result_gpu = transform(x_gpu) # Transfer result back to CPU if needed result_cpu = Array(result_gpu) ``` #### Metal Example ```julia using Metal using DistanceTransforms # Create a random array on GPU x_gpu = Metal.rand([0, 1], 1000, 1000) x_gpu = boolean_indicator(x_gpu) # The transform function automatically uses GPU result_gpu = transform(x_gpu) # Transfer result back to CPU if needed result_cpu = Array(result_gpu) ``` ### Performance Benchmarks Performance comparison across different implementations for 2D and 3D distance transforms: ![Performance Comparison](../images/julia_distance_transforms.png) As shown in the graph, the GPU implementation demonstrates superior performance, especially for larger arrays. ### Best Practices 1. **For small arrays (<100x100)**: CPU with multi-threading is often sufficient 2. **For medium arrays**: Multi-threaded CPU may be faster than GPU due to lower overhead 3. **For large arrays (>1000x1000)**: GPU acceleration provides the best performance 4. **For 3D data**: GPU acceleration is strongly recommended due to the computational complexity ## Algorithm Details On the CPU, DistanceTransforms.jl uses the squared Euclidean distance transform algorithm by [Felzenszwalb and Huttenlocher](https://theoryofcomputing.org/articles/v008a019/), known for its accuracy and efficiency. On the GPU, thanks to the amazing [AcceleratedKernels.jl](https://github.com/JuliaGPU/AcceleratedKernels.jl) and [KernelAbstraction.jl](https://github.com/JuliaGPU/KernelAbstractions.jl) packages, the algorithm is essentially identical.