Add circular PMF. (#14)

Project.toml

@@ -1,7 +1,7 @@
 name = "GeoArrayOps"
 uuid = "e6005643-8f00-58df-85c5-7d93a3520d70"
 authors = ["Maarten Pronk <[email protected]>", "Deltares"]
-version = "0.1.1"
+version = "0.2.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -10,13 +10,19 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+DataStructures = "0.18"
+Distances = "0.10"
+FillArrays = "0.12"
 ImageCore = "^0.9"
 ImageFiltering = "^0.7"
 OffsetArrays = "^1.10"
+PaddedViews = "0.5"
+StaticArrays = "1"
 julia = "1.5, 1.6"
 
 [extras]

src/pmf.jl

@@ -21,11 +21,12 @@ Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .
 [^zhang2003]: Zhang, Keqi, Shu-Ching Chen, Dean Whitman, Mei-Ling Shyu, Jianhua Yan, and Chengcui Zhang. “A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LIDAR Data.” IEEE Transactions on Geoscience and Remote Sensing 41, no. 4 (2003): 872–82. <https://doi.org/10.1109/TGRS.2003.810682>.
 """
 function pmf(A::AbstractMatrix{T};
- ωₘ::Float64=20.,
- slope::Float64=0.01,
- dhₘ::Float64=2.5,
- dh₀::Float64=0.2,
- cellsize::Float64=1.0) where T <: Real
+ ωₘ::Float64 = 20.0,
+ slope::Float64 = 0.01,
+ dhₘ::Float64 = 2.5,
+ dh₀::Float64 = 0.2,
+ cellsize::Float64 = 1.0,
+ circular = false) where {T<:Real}
 
  # Compute windowsizes and thresholds
  ωₘ = round(Int, ωₘ / cellsize)
@@ -34,8 +35,9 @@ function pmf(A::AbstractMatrix{T};
 
  # Compute tresholds
  dwindows = vcat(windowsizes[1], windowsizes) # prepend first element so we get 0 as diff
- window_diffs = [dwindows[i] - dwindows[i - 1] for i in 2:length(dwindows)]
+ window_diffs = [dwindows[i] - dwindows[i-1] for i = 2:length(dwindows)]
  height_tresholds = [min(dhₘ, slope * window_diff * cellsize + dh₀) for window_diff in window_diffs]
+ @info "Using the following thresholds: $height_tresholds for the following windows: $windowsizes"
 
  # Set up arrays
  Af = copy(A) # array to be morphed
@@ -52,7 +54,11 @@ function pmf(A::AbstractMatrix{T};
 
  # Iterate over window sizes and height tresholds
  for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds)
- opening!(Af, ωₖ, out)
+ if circular
+ opening_circ!(Af, ωₖ, out)
+ else
+ opening!(Af, ωₖ, out)
+ end
  mask .= (A .- Af) .> dhₜ
  for I in eachindex(flags)
  if mask[I] && flags[I] == 0

src/smf.jl

@@ -16,13 +16,14 @@ Applies the simple morphological filter by *Pingel et al. (2013)* [^pingel2013]
 [^pingel2013]: Pingel, Thomas J., Keith C. Clarke, and William A. McBride. 2013. ‘An Improved Simple Morphological Filter for the Terrain Classification of Airborne LIDAR Data’. ISPRS Journal of Photogrammetry and Remote Sensing 77 (March): 21–30. <https://doi.org/10.1016/j.isprsjprs.2012.12.002>.
 """
 function smf(A::AbstractMatrix{T};
- ω::Real=17.,
- slope::Real=0.01,
- cellsize::Real=1.0) where T <: Real
+ ω::Real = 17.0,
+ slope::Real = 0.01,
+ cellsize::Real = 1.0) where {T<:Real}
 
  lastsurface = -copy(A)
  is_low = falses(size(A))
  radii = 3:2:round(Int, ω / cellsize)
+
  for radius ∈ radii
  threshold = slope * radius * cellsize
  thissurface = opening_circ(lastsurface, radius)
@@ -39,6 +40,6 @@ function smf(A::AbstractMatrix{T};
  lastsurface = thissurface
  end
  out = Array{Union{Missing,T}}(A)
- out[is_low .| is_obj] .= missing
+ out[is_low.|is_obj] .= missing
  return out
 end

src/terrain.jl

@@ -1,14 +1,14 @@
-const neib_8 = @SMatrix[1. 1 1; 1 0 1; 1 1 1]
+const neib_8 = @SMatrix [1.0 1 1; 1 0 1; 1 1 1]
 
-function buffer_array(A::Array)
- oA = OffsetMatrix(fill(Inf, size(A) .+ 2), UnitRange.(0, size(A) .+ 1))
+function buffer_array(A::AbstractMatrix{<:Real})
+ oA = OffsetMatrix(fill(zero(eltype(A)), size(A) .+ 2), UnitRange.(0, size(A) .+ 1))
  # Update center
- oA[begin + 1:end - 1,begin + 1:end - 1] .= A
+ oA[begin+1:end-1, begin+1:end-1] .= A
  # Set edges to mirror center
- oA[begin, begin + 1:end - 1] .= A[begin, :]
- oA[end, begin + 1:end - 1] .= A[end, :]
- oA[begin + 1:end - 1, begin] .= A[:, begin]
- oA[begin + 1:end - 1, end] .= A[:, end]
+ oA[begin, begin+1:end-1] .= A[begin, :]
+ oA[end, begin+1:end-1] .= A[end, :]
+ oA[begin+1:end-1, begin] .= A[:, begin]
+ oA[begin+1:end-1, end] .= A[:, end]
  # Set corners to mirror corners of center
  oA[begin, begin] = A[begin, begin]
  oA[begin, end] = A[begin, end]
@@ -29,10 +29,10 @@ function roughness(dem::AbstractMatrix{<:Real})
  x = @MMatrix zeros(3, 3)
 
  @inbounds for I ∈ CartesianIndices(size(roughness))
- patch = I - Δ:I + Δ
+ patch = I-Δ:I+Δ
  rdata = view(ex_dem, patch)
 
- x .= rdata .- rdata[2,2]
+ x .= rdata .- rdata[2, 2]
  roughness[I] = maximum(abs.(x))
  end
  roughness
@@ -50,11 +50,11 @@ function TPI(dem::AbstractMatrix{<:Real})
  x = @MMatrix zeros(3, 3)
 
  @inbounds for I ∈ CartesianIndices(size(tpi))
- patch = I - Δ:I + Δ
+ patch = I-Δ:I+Δ
  rdata = view(ex_dem, patch)
 
  x .= rdata .* neib_8
- tpi[I] = rdata[2,2] - mean(x)
+ tpi[I] = rdata[2, 2] - mean(x)
  end
  return tpi
 end
@@ -73,10 +73,10 @@ function TRI(dem::AbstractMatrix{<:Real})
  x = @MMatrix zeros(3, 3)
 
  @inbounds for I ∈ CartesianIndices(size(tri))
- patch = I - Δ:I + Δ
+ patch = I-Δ:I+Δ
  rdata = view(ex_dem, patch)
 
- x .= (rdata .- rdata[2,2]).^2
+ x .= (rdata .- rdata[2, 2]) .^ 2
  tri[I] = sqrt(sum(x))
  end
  return tri

src/utils.jl

@@ -1,20 +1,29 @@
 using OffsetArrays
 using ImageFiltering
+using Distances
+using PaddedViews
 
 """Apply the opening operation to `A` with window size `ω`."""
-function opening(A::Array{T,2}, ω::Integer) where T <: Real
+function opening(A::Array{T,2}, ω::Integer) where {T<:Real}
  A = mapwindow(minimum, A, (ω, ω)) # erosion
  A = mapwindow(maximum, A, (ω, ω)) # dilation
  A
 end
 
 """Apply the opening operation to `A` with window size `ω`."""
-function opening!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where T <: Real
+function opening!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real}
  mapwindow!(minimum, A, ω, out) # erosion
  mapwindow!(maximum, out, ω, A) # dilation
  A
 end
 
+"""Apply the opening operation to `A` with window size `ω`."""
+function opening_circ!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real}
+ mapwindowcirc!(minimum_mask, A, ω, out, Inf) # erosion
+ mapwindowcirc!(maximum_mask, out, ω, A, -Inf) # dilation
+ A
+end
+
 function circmask(n::Integer)
  # TODO This could be precomputed for first N integers
  kern = falses(-n:n, -n:n)
@@ -25,7 +34,7 @@ function circmask(n::Integer)
  return kern.parent
 end
 
-function opening_circ(A::Array{T,2}, ω::Integer) where T <: Real
+function opening_circ(A::Array{T,2}, ω::Integer) where {T<:Real}
  m = circmask(ω >> 1)
  A = mapwindow(x -> minimum(x[m]), A, (ω, ω)) # erosion
  A = mapwindow(x -> maximum(x[m]), A, (ω, ω)) # dilation
@@ -38,9 +47,45 @@ function mapwindow!(f, img, window, out)
  R = CartesianIndices(img)
  I_first, I_last = first(R), last(R)
  Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))
- @inbounds @simd for I in R
+ for I in R
  patch = max(I_first, I - Δ):min(I_last, I + Δ)
  out[I] = f(view(img, patch))
  end
  out
 end
+
+function mapwindowcirc!(f, img, window, out, fill = Inf)
+ R = CartesianIndices(img)
+ Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))
+
+ w, h = size(img)
+ A = PaddedView(fill, img, (-Δ[1]+1:w+Δ[1], -Δ[2]+1:h+Δ[2]))
+
+ patch = (-Δ):(Δ)
+ m = euclidean.(Tuple.(patch), Ref((0, 0))) .<= Δ[1]
+ for I in R
+ patch = (I-Δ):(I+Δ)
+ out[I] = f(view(A, patch), m)
+ end
+ out
+end
+
+function maximum_mask(x, m)
+ o = x[1]
+ for I in eachindex(x)
+ if m[I]
+ o = max(o, x[I])
+ end
+ end
+ o
+end
+
+function minimum_mask(x, m)
+ o = x[1]
+ for I in eachindex(x)
+ if m[I]
+ o = min(o, x[I])
+ end
+ end
+ o
+end

test/runtests.jl

@@ -4,14 +4,16 @@ using Test
 @testset "pmf" begin
  # Write your own tests here.
  A = rand(25, 25)
- A[2,2] = NaN
+ A[2, 2] = NaN
  B, flags = pmf(A)
- @test (A .<= B) == (flags .== 0.)
+ @test (A .<= B) == (flags .== 0.0)
+ B, flags = pmf(A; circular = true)
+ @test (A .<= B) == (flags .== 0.0)
 end
 @testset "smf" begin
  # Write your own tests here.
  A = rand(25, 25)
- A[2,2] = NaN
+ A[2, 2] = NaN
  B = smf(A)
 end
 @testset "pssm" begin