Add skewness filter (#18)

.github/workflows/CI.yml

@@ -9,8 +9,9 @@ jobs:
  fail-fast: false
  matrix:
  version:
- - '1.6'
- - 'nightly'
+ - "1.6"
+ - "1"
+ - "nightly"
  os:
  - ubuntu-latest
  - macOS-latest

.gitignore

@@ -4,3 +4,7 @@
 .DS_Store
 Manifest.toml
 docs/build/
+*.tif
+*.gif
+test/benchmark.jl
+.vscode

Project.toml

@@ -1,7 +1,7 @@
 name = "GeoArrayOps"
 uuid = "e6005643-8f00-58df-85c5-7d93a3520d70"
 authors = ["Maarten Pronk <[email protected]>", "Deltares"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -11,21 +11,21 @@ ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
-ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 DataStructures = "0.18"
 Distances = "0.10"
 FillArrays = "0.12, 0.13"
 ImageCore = "0.9"
-ImageFiltering = "0.7"
+ImageFiltering = "0.6, 0.7"
 OffsetArrays = "1.10"
-ProgressMeter = "1.7"
 PaddedViews = "0.5"
 StaticArrays = "1"
-julia = "1.5, 1.6, 1.7"
+StatsBase = "0.33"
+julia = "1.5, 1.6, 1.7, 1.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -4,10 +4,8 @@
 # GeoArrayOps
 Geospatial operations, cost and filtering algorithms as used in for elevation rasters.
 
-*This is a work in progress*
-
 ## Functionality
-- Terrain filters, such as Progressive Morphological Filters (PMF, SMF)
+- Terrain filters, such as Progressive Morphological Filters (PMF, SMF) and Skewness balancing
 - Geospatial cost (friction) operations that mimic PCRaster. These functions should however be more Julian, extensible and scale better.
 - Visualization, such as Perceptually Shaded Slope Map (PSSM)
 - Terrain analysis functions, such as roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI).
@@ -17,5 +15,5 @@ The package can be installed with the Julia package manager.
 From the Julia REPL, type `]` to enter the Pkg REPL mode and run:
 
 ```
-pkg> add https://github.com/Deltares/GeoArrayOps.jl
+pkg> add GeoArrayOps
 ```

docs/src/index.md

@@ -4,10 +4,8 @@
 # GeoArrayOps
 Geospatial operations, cost and filtering algorithms as used in for elevation rasters.
 
-*This is a work in progress*
-
 ## Functionality
-- Terrain filters, such as Progressive Morphological Filters (PMF, SMF)
+- Terrain filters, such as Progressive Morphological Filters (PMF, SMF) and Skewness balancing
 - Geospatial cost (friction) operations that mimic PCRaster. These functions should however be more Julian, extensible and scale better.
 - Visualization, such as Perceptually Shaded Slope Map (PSSM)
 - Terrain analysis functions, such as roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI).
@@ -17,7 +15,7 @@ The package can be installed with the Julia package manager.
 From the Julia REPL, type `]` to enter the Pkg REPL mode and run:
 
 ```
-pkg> add https://github.com/Deltares/GeoArrayOps.jl
+pkg> add GeoArrayOps
 ```
 
 ## Index

src/GeoArrayOps.jl

@@ -1,6 +1,14 @@
-__precompile__()
 module GeoArrayOps
-using ProgressMeter
+using StatsBase: skewness
+using OffsetArrays: OffsetMatrix
+using ImageFiltering: mapwindow, sobel, imgradients
+using Distances: Euclidean, euclidean, evaluate
+using PaddedViews: PaddedView
+using FillArrays: fill
+using StaticArrays: @SMatrix, @MMatrix, SMatrix, MMatrix
+using DataStructures: Deque
+using Statistics: median, mean
+using ImageCore: scaleminmax, Gray
 
 include("utils.jl")
 include("pmf.jl")
@@ -9,36 +17,13 @@ include("psf.jl")
 include("plot.jl")
 include("spread.jl")
 include("terrain.jl")
+include("skew.jl")
 
 export pmf, smf, psf
 export pssm
 export pitremoval
 export spread, spread2
 export roughness, TRI, TPI
-
-# function __init__()
-# A = rand(Float64, 3, 3)
-# pssm(A)
-# spread(A, A, A)
-# spread2(A, A, A)
-# spread(A, A, 1.)
-# spread(A, 1., 1.)
-# roughness(A)
-# TRI(A)
-# TPI(A)
-# end
-
-precompile(pmf, (Matrix{Float64},))
-precompile(psf, (Matrix{Float64},))
-precompile(pitremoval, (Matrix{Float64},))
-precompile(smf, (Matrix{Float64},))
-precompile(pssm, (Matrix{Float64},))
-precompile(spread, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}))
-precompile(spread2, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}))
-precompile(spread, (Matrix{Float64}, Matrix{Float64}, Float64))
-precompile(spread, (Matrix{Float64}, Float64, Float64))
-precompile(roughness, (Matrix{Float64},))
-precompile(TRI, (Matrix{Float64},))
-precompile(TPI, (Matrix{Float64},))
+export skb
 
 end # module

src/plot.jl

@@ -1,7 +1,4 @@
 # using RecipesBase
-using ImageFiltering
-using ImageCore
-
 """
 ```
 image = pssm(A; exaggeration, resolution)
@@ -18,13 +15,16 @@ Perceptually Shaded Slope Map by *Pingel, Clarke. 2014* [^pingel2014].
 
 [^pingel2014]: Pingel, Thomas, and Clarke, Keith. 2014. ‘Perceptually Shaded Slope Maps for the Visualization of Digital Surface Models’. Cartographica: The International Journal for Geographic Information and Geovisualization 49 (4): 225–40. <https://doi.org/10/ggnthv>.
 """
-function pssm(A::AbstractMatrix{<:Real}; exaggeration=2.3, resolution=1.)
- x, y = imgradients(A * exaggeration, ImageFiltering.sobel)
- G = sqrt.(x.^2 .+ y.^2)
-
- G /= resolution # account for horizontal resolution
-
+function pssm(A::AbstractMatrix{<:Real}; exaggeration=2.3, resolution=1.0)
+ G = slope(A, exaggeration, resolution)
  f = scaleminmax(0, 1)
  clamped = f.(G)
  Gray.(1 .- clamped)
 end
+
+function slope(A, exaggeration=1.0, resolution=1.0)
+ x, y = imgradients(A .* exaggeration, sobel)
+ G = sqrt.(x .^ 2 .+ y .^ 2)
+ G /= resolution # account for horizontal resolution
+ G
+end

src/pmf.jl

@@ -21,12 +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.0,
- slope::Float64 = 0.01,
- dhₘ::Float64 = 2.5,
- dh₀::Float64 = 0.2,
- cellsize::Float64 = 1.0,
- circular = false) 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)
@@ -37,7 +37,7 @@ function pmf(A::AbstractMatrix{T};
  dwindows = vcat(windowsizes[1], windowsizes) # prepend first element so we get 0 as diff
  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"
+ # @info "Using the following thresholds: $height_tresholds for the following windows: $windowsizes"
 
  # Set up arrays
  Af = copy(A) # array to be morphed

src/psf.jl

@@ -1,4 +1,3 @@
-
 """
 ```
 B, flags = psf(A; ωₘ, slope, dhₘ, dh₀, cellsize)
@@ -19,13 +18,13 @@ Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .
 - `dh₀::Float64=0.2` Initial elevation threshold [m]
 - `cellsize::Float64=1.` Cellsize in [m]
 """
-function psf(A::AbstractMatrix{T};
- ωₘ::Float64 = 20.0,
- slope::Float64 = 0.01,
- dhₘ::Float64 = 2.5,
- dh₀::Float64 = 0.2,
- cellsize::Float64 = 1.0,
- circular = false) where {T<:Real}
+function apsf(A::AbstractMatrix{T};
+ ωₘ::Float64=20.0,
+ slope::Matrix{Float64},
+ dhₘ::Float64=2.5,
+ dh₀::Float64=0.2,
+ cellsize::Float64=1.0,
+ circular=false) where {T<:Real}
 
  # Compute windowsizes and thresholds
  ωₘ = round(Int, ωₘ / cellsize)
@@ -35,8 +34,8 @@ function psf(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 = 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"
+ # height_tresholds = [min(dhₘ, slope[1] * window_diff * cellsize + dh₀) for window_diff in window_diffs]
+ # @debug "Using the following thresholds: $height_tresholds for the following windows: $windowsizes"
 
  # Set up arrays
  Af = copy(A) # array to be morphed
@@ -49,15 +48,15 @@ function psf(A::AbstractMatrix{T};
  flags[nan_mask] .= NaN
 
  mask = falses(size(A))
+ dhₜ = min.(dhₘ, slope * window_diffs[1] * cellsize .+ dh₀)
 
  # Iterate over window sizes and height tresholds
- p = Progress(sum(windowsizes .^ 2))
- progress = 0
- for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds)
+ for (i, ωₖ) in enumerate(windowsizes)
+ dhₜ = min.(dhₘ, slope * window_diffs[i] * cellsize .+ dh₀)
  if circular
- mapwindowcirc!(minimum_mask, A, ωₖ, Af, Inf)
+ mapwindowcirc_approx!(minimum_mask, A, ωₖ, Af, Inf)
  else
- mapwindow!(minimum, A, ωₖ, Af)
+ mapwindow_stack!(minimum, A, ωₖ, Af)
  end
  mask .= (A .- Af) .> dhₜ
  for I in eachindex(flags)
@@ -66,9 +65,25 @@ function psf(A::AbstractMatrix{T};
  end
  end
  B .= min.(B, Af .+ dhₜ)
- progress += ωₖ^2
- ProgressMeter.update!(p, progress)
  end
 
- B, flags
+ B, flags, dhₜ
+end
+
+
+
+function psf(A::AbstractMatrix{T};
+ ωₘ::Float64=20.0,
+ slope::Float64=0.01,
+ dhₘ::Float64=2.5,
+ dh₀::Float64=0.2,
+ cellsize::Float64=1.0,
+ circular=false) where {T<:Real}
+ return apsf(A;
+ ωₘ=ωₘ,
+ slope=fill(slope, size(A)),
+ dhₘ=dhₘ,
+ dh₀=dh₀,
+ cellsize=cellsize,
+ circular=circular)
 end

src/skew.jl

@@ -0,0 +1,34 @@
+# Skewness balancing as by Bartels (2006)
+"""
+```
+mask = skb(A, mean)
+mask = skb(A)
+```
+Applies skewness balancing by *Bartels e.a (2006)* [^bartels2006] to `A`.
+
+# Output
+- `mask::BitMatrix` Maximum allowable values
+
+Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .== 0.`.
+
+[^bartels2006]: Bartels, M., Hong Wei, and D.C. Mason. 2006. “DTM Generation from LIDAR Data Using Skewness Balancing.” In 18th International Conference on Pattern Recognition (ICPR’06), 1:566–69. https://doi.org/10/cwk4v2.
+"""
+function skb(A::AbstractArray{T}, mean::T) where {T<:Real}
+ I = sortperm(vec(A))
+ AA = A[I]
+ s = 1
+ for i ∈ length(A):-1:1
+ X = skewness(AA[begin:i], mean)
+ if X <= 0
+ s = i + 1
+ break
+ end
+ end
+ m = trues(size(A))
+ m[I[s+1:end]] .= false
+ return m
+end
+
+function skb(A::AbstractArray{T}) where {T<:Real}
+ return skb(A, mean(A))
+end

src/smf.jl

@@ -16,9 +16,9 @@ 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.0,
- 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))