Add skewness filter (#18)

Author: evetion

.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