From 89b855dcb13053be298dcbc6bf847dabd80aa4cc Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Fri, 17 Dec 2021 17:36:41 +0100 Subject: [PATCH 1/2] Add slope filter. --- Project.toml | 1 + src/GeoArrayOps.jl | 4 ++- src/psf.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++ src/utils.jl | 4 +-- 4 files changed, 81 insertions(+), 3 deletions(-) create mode 100644 src/psf.jl diff --git a/Project.toml b/Project.toml index fd56951..ac890b5 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ 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" diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl index 64ebb1e..fae1925 100644 --- a/src/GeoArrayOps.jl +++ b/src/GeoArrayOps.jl @@ -4,11 +4,12 @@ module GeoArrayOps include("utils.jl") include("pmf.jl") include("smf.jl") +include("psf.jl") include("plot.jl") include("spread.jl") include("terrain.jl") -export pmf, smf +export pmf, smf, psf export pssm export spread, spread2 export roughness, TRI, TPI @@ -26,6 +27,7 @@ export roughness, TRI, TPI # end precompile(pmf, (Matrix{Float64},)) +precompile(psf, (Matrix{Float64},)) precompile(smf, (Matrix{Float64},)) precompile(pssm, (Matrix{Float64},)) precompile(spread, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64})) diff --git a/src/psf.jl b/src/psf.jl new file mode 100644 index 0000000..749392e --- /dev/null +++ b/src/psf.jl @@ -0,0 +1,75 @@ +using ProgressMeter + +""" +``` +B, flags = psf(A; ωₘ, slope, dhₘ, dh₀, cellsize) +``` +Applies a progressive slope filter to `A`. + +# Output +- `B::Array{T,2}` Maximum allowable values based A + slope, dhₘ, dh₀ +- `flags::Array{Float64,2}` A sized array with window sizes if filtered, zero if not filtered. + +Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .== 0.`. + +# Arguments +- `A::Array{T,2}` Input Array +- `ωₘ::Float64=20.` Maximum window size [m] +- `slope::Float64=0.01` Terrain slope [m/m] +- `dhₘ::Float64=2.5` Maximum elevation threshold [m] +- `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} + + # Compute windowsizes and thresholds + ωₘ = round(Int, ωₘ / cellsize) + κ_max = floor(Int, log2(ωₘ - 1)) # determine # iterations based on exp growth + windowsizes = Int.(exp2.(1:κ_max)) .+ 1 + + # 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" + + # Set up arrays + Af = copy(A) # array to be morphed + nan_mask = isnan.(Af) + Af[nan_mask] .= Inf # Replace NaN with Inf, as to always filter these + + B = copy(A) # max_elevation raster + + flags = zeros(size(A)) # 0 = ground, other values indicate window size + flags[nan_mask] .= NaN + + mask = falses(size(A)) + + # Iterate over window sizes and height tresholds + p = Progress(sum(windowsizes .^ 2)) + progress = 0 + for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds) + if circular + mapwindowcirc!(minimum_mask, A, ωₖ, Af, Inf) + else + mapwindow!(minimum, A, ωₖ, Af) + end + mask .= (A .- Af) .> dhₜ + for I in eachindex(flags) + if mask[I] && flags[I] == 0 + flags[I] = ωₖ + end + end + B .= min.(B, Af .+ dhₜ) + progress += ωₖ^2 + ProgressMeter.update!(p, progress) + end + + B, flags +end diff --git a/src/utils.jl b/src/utils.jl index 929bd75..bb08854 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -71,7 +71,7 @@ function mapwindowcirc!(f, img, window, out, fill = Inf) end function maximum_mask(x, m) - o = x[1] + o = -Inf for I in eachindex(x) if m[I] o = max(o, x[I]) @@ -81,7 +81,7 @@ function maximum_mask(x, m) end function minimum_mask(x, m) - o = x[1] + o = Inf for I in eachindex(x) if m[I] o = min(o, x[I]) From 8ae2a3be410b50bbed9a372002b907840ce0289a Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Fri, 24 Dec 2021 13:26:56 +0100 Subject: [PATCH 2/2] Add pit filter. --- src/GeoArrayOps.jl | 3 +++ src/psf.jl | 1 - src/terrain.jl | 34 ++++++++++++++++++++++++++++++++++ src/utils.jl | 20 ++++++++------------ test/runtests.jl | 6 ++++++ 5 files changed, 51 insertions(+), 13 deletions(-) diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl index fae1925..9050794 100644 --- a/src/GeoArrayOps.jl +++ b/src/GeoArrayOps.jl @@ -1,5 +1,6 @@ __precompile__() module GeoArrayOps +using ProgressMeter include("utils.jl") include("pmf.jl") @@ -11,6 +12,7 @@ include("terrain.jl") export pmf, smf, psf export pssm +export pitremoval export spread, spread2 export roughness, TRI, TPI @@ -28,6 +30,7 @@ export roughness, TRI, TPI 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})) diff --git a/src/psf.jl b/src/psf.jl index 749392e..9f5fab1 100644 --- a/src/psf.jl +++ b/src/psf.jl @@ -1,4 +1,3 @@ -using ProgressMeter """ ``` diff --git a/src/terrain.jl b/src/terrain.jl index 568f254..06b47ae 100644 --- a/src/terrain.jl +++ b/src/terrain.jl @@ -1,4 +1,5 @@ const neib_8 = @SMatrix [1.0 1 1; 1 0 1; 1 1 1] +const neib_8_inf = @SMatrix [1.0 1 1; 1 Inf 1; 1 1 1] function buffer_array(A::AbstractMatrix{<:Real}) oA = OffsetMatrix(fill(zero(eltype(A)), size(A) .+ 2), UnitRange.(0, size(A) .+ 1)) @@ -48,6 +49,7 @@ function TPI(dem::AbstractMatrix{<:Real}) ex_dem = buffer_array(dem) tpi = similar(dem) x = @MMatrix zeros(3, 3) + Δ = CartesianIndex(1, 1) @inbounds for I ∈ CartesianIndices(size(tpi)) patch = I-Δ:I+Δ @@ -71,6 +73,7 @@ function TRI(dem::AbstractMatrix{<:Real}) ex_dem = buffer_array(dem) tri = similar(dem) x = @MMatrix zeros(3, 3) + Δ = CartesianIndex(1, 1) @inbounds for I ∈ CartesianIndices(size(tri)) patch = I-Δ:I+Δ @@ -81,3 +84,34 @@ function TRI(dem::AbstractMatrix{<:Real}) end return tri end + + +""" +pitremoval(dem::Matrix{<:Real}) + +Remove pits from a DEM Array if the center cell of a 3x3 patch is `limit` lower or than the surrounding cells. +""" +function pitremoval(dem::AbstractMatrix{<:Real}, limit = 2.0) + + ex_dem = buffer_array(dem) + tri = similar(dem) + x = @MMatrix zeros(3, 3) + Δ = CartesianIndex(1, 1) + + @inbounds for I ∈ CartesianIndices(size(tri)) + patch = I-Δ:I+Δ + x .= view(ex_dem, patch) + tri[I] = _pitremoval(x, limit) + end + return tri +end + +@inline function _pitremoval(x, limit) + A = vec(x) + order = sortperm(A) + ifelse( + (order[1] == 5) && + ((A[order[2]] - A[order[1]]) >= limit), + Inf, + x[2, 2]) +end diff --git a/src/utils.jl b/src/utils.jl index bb08854..b4b63b7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -47,7 +47,7 @@ function mapwindow!(f, img, window, out) R = CartesianIndices(img) I_first, I_last = first(R), last(R) Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img))) - for I in R + Threads.@threads for I in R patch = max(I_first, I - Δ):min(I_last, I + Δ) out[I] = f(view(img, patch)) end @@ -63,29 +63,25 @@ function mapwindowcirc!(f, img, window, out, fill = Inf) patch = (-Δ):(Δ) m = euclidean.(Tuple.(patch), Ref((0, 0))) .<= Δ[1] - for I in R + Threads.@threads for I in R patch = (I-Δ):(I+Δ) out[I] = f(view(A, patch), m) end out end -function maximum_mask(x, m) +@inline function maximum_mask(x, m) o = -Inf - for I in eachindex(x) - if m[I] - o = max(o, x[I]) - end + @inbounds for I in eachindex(x) + o = ifelse(m[I], max(o, x[I]), o) end o end -function minimum_mask(x, m) +@inline function minimum_mask(x, m) o = Inf - for I in eachindex(x) - if m[I] - o = min(o, x[I]) - end + @inbounds for I in eachindex(x) + o = ifelse(m[I], min(o, x[I]), o) end o end diff --git a/test/runtests.jl b/test/runtests.jl index 93c7ca7..744cec9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,12 @@ end @testset "pssm" begin B = pssm(rand(25, 25)) end +@testset "psf" begin + B = psf(rand(25, 25)) +end +@testset "pitremoval" begin + B = pitremoval(rand(25, 25)) +end @testset "spread" begin points = [0.0 0 0 0 2; 0 0 0 0 0; 0 0 0 0 0; 0 1 0 0 0; 0 0 0 0 0] initial = [8.0 8 8 8 4; 8 8 8 8 8; 8 8 8 8 8; 0 0 8 8 8; 0 0 8 8 8]