From 59a567a9dd433c78da23c86b7595fee33caeb9f0 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 2 Nov 2022 17:38:31 +0100 Subject: [PATCH 1/2] Use LocalFilters in terrain functions. Fix skb. Use approximate circular in PMF. --- Project.toml | 2 +- src/skew.jl | 5 +- src/terrain.jl | 299 +++++++++++++++++++++++++++++++++++-------------- src/utils.jl | 38 ++++++- 4 files changed, 251 insertions(+), 93 deletions(-) diff --git a/Project.toml b/Project.toml index b9764fe..ffa6551 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoArrayOps" uuid = "e6005643-8f00-58df-85c5-7d93a3520d70" authors = ["Maarten Pronk ", "Deltares"] -version = "0.4.0" +version = "0.5.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" diff --git a/src/skew.jl b/src/skew.jl index 65af3dd..8fc2cbd 100644 --- a/src/skew.jl +++ b/src/skew.jl @@ -9,7 +9,7 @@ Improved the performance by applying a binary search to find the threshold value [^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(iA::AbstractArray{T}; mean::T) where {T<:Real} +function skb(iA::AbstractArray{T}; mean::T=mean(iA)) where {T<:Real} m = .!isfinite.(iA) if sum(m) > 0 A = copy(iA) @@ -45,9 +45,6 @@ function skb(iA::AbstractArray{T}; mean::T) where {T<:Real} return m end -function skb(A::AbstractArray{T}) where {T<:Real} - return skb(A; mean=mean(A)) -end """ mask = skbr(A; iterations=10) diff --git a/src/terrain.jl b/src/terrain.jl index 4ec8104..f195435 100644 --- a/src/terrain.jl +++ b/src/terrain.jl @@ -2,6 +2,8 @@ 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] const neib_8_mask = @SMatrix Bool[1 1 1; 1 0 1; 1 1 1] +const nbkernel = LocalFilters.Kernel{Int8,2}(reshape(1:9, 3, 3)) + abstract type DerivativeMethod end """Second order finite difference estimator using all 4 neighbors (Zevenbergen and Thorne, 1987).""" @@ -13,40 +15,6 @@ struct Horn <: DerivativeMethod end """Maximum Downward Gradient""" struct MDG <: DerivativeMethod end -function noise(factor=0.01) - (rand() - 0.5) * factor -end - -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 - # 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] - # Set corners to mirror corners of center - oA[begin, begin] = A[begin, begin] + noise() - oA[begin, end] = A[begin, end] + noise() - oA[end, begin] = A[end, begin] + noise() - oA[end, end] = A[end, end] + noise() - return oA -end - -function terrain_kernel(dem::AbstractMatrix{T}, f::Function, t=T) where {T<:Real} - dembuffer = buffer_array(dem) - output = similar(dem, t) - Δ = CartesianIndex(1, 1) - - @inbounds for I ∈ CartesianIndices(size(output)) - patch = I-Δ:I+Δ - dempatch = SMatrix{3,3}(view(dembuffer, patch)) - output[I] = f(dempatch) - end - output -end - """ roughness(dem::Matrix{<:Real}) @@ -77,7 +45,6 @@ function TPI(dem::AbstractMatrix{<:Real}) return localfilter!(dst, dem, 3, initial, update, store!) end - """ TRI(dem::Matrix{<:Real}) @@ -124,49 +91,147 @@ end Slope is the rate of change between a cell and its neighbors as defined in Burrough, P. A., and McDonell, R. A., (1998, Principles of Geographical Information Systems). """ function slope(dem::AbstractMatrix{<:Real}; cellsize=1.0, method::DerivativeMethod=Horn()) - terrain_kernel(dem, f -> slope_kernel(f, cellsize, method)) + dst = copy(dem) + slope(method, dst, dem, cellsize) end -function slope_kernel(A, cellsize=1.0, method::DerivativeMethod=Horn()) - δzδx, δzδy = derivative(method, A, cellsize) - return atand(√(δzδx^2 + δzδy^2)) + +function slope(::Horn, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + store!(d, i, v) = @inbounds d[i] = atand( + √( + ((v[1] - v[2]) / (8 * v[5]))^2 + ((v[3] - v[4]) / (8 * v[5]))^2 + )) + return localfilter!(dst, dem, nbkernel, initial, horn, store!) +end + +function slope(::ZevenbergenThorne, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), cellsize) + store!(d, i, v) = @inbounds d[i] = atand( + √( + (v[1] / (2 * v[3]))^2 + (v[2] / (2 * v[3]))^2 + )) + return localfilter!(dst, dem, nbkernel, initial, zevenbergenthorne, store!) end +function slope(::MDG, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + m = max( + abs(v[1]) / (v[5] * sqrt2), + abs(v[2]) / v[5], + abs(v[3]) / (v[5] * sqrt2), + abs(v[4]) / v[5], + ) / 2 + d[i] = atand(√(m^2 + m^2)) + end + return localfilter!(dst, dem, nbkernel, initial, mdg, store!) +end + + """ aspect(dem::Matrix{<:Real}, method=Horn()) Aspect is direction of [`slope`](@ref), as defined in Burrough, P. A., and McDonell, R. A., (1998, Principles of Geographical Information Systems). """ function aspect(dem::AbstractMatrix{<:Real}; method::DerivativeMethod=Horn()) - terrain_kernel(dem, f -> aspect_kernel(f, method)) + dst = copy(dem) + aspect(method, dst, dem, 1) +end + +# Useless, as there's no x/y component. +function aspect(::MDG, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + δzδx = max( + abs(v[1]) / (v[5] * sqrt2), + abs(v[2]) / v[5], + abs(v[3]) / (v[5] * sqrt2), + abs(v[4]) / v[5], + ) / 2 + d[i] = compass(atand(δzδx, -δzδx)) + end + return localfilter!(dst, dem, nbkernel, initial, mdg, store!) end -function aspect_kernel(A, method::DerivativeMethod=Horn()) - δzδx, δzδy = derivative(method, A, 1.0) - return compass(atand(δzδy, -δzδx)) +function aspect(::Horn, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + δzδx = (v[1] - v[2]) / (8 * v[5]) + δzδy = (v[3] - v[4]) / (8 * v[5]) + # @info δzδx, δzδy + d[i] = compass(atand(-δzδx, δzδy)) + end + return localfilter!(dst, dem, nbkernel, initial, horn, store!) end -function derivative(::ZevenbergenThorne, A, cellsize=1.0) - δzδx = (A[1, 2] - A[3, 2]) / 2 * cellsize - δzδy = (A[2, 1] - A[2, 3]) / 2 * cellsize - return δzδx, δzδy +function aspect(::ZevenbergenThorne, dst, dem::AbstractMatrix{<:Real}, cellsize) + initial(A) = (zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + δzδx = v[1] / 2 + δzδy = v[2] / 2 + # @info δzδx, δzδy + d[i] = compass(atand(δzδy, -δzδx)) + end + return localfilter!(dst, dem, nbkernel, initial, zevenbergenthorne, store!) end -function derivative(::Horn, A, cellsize=1.0) - δzδx = ((A[1, 3] + 2A[2, 3] + A[3, 3]) - (A[1, 1] + 2A[2, 1] + A[3, 1])) / (8 * cellsize) - δzδy = ((A[3, 1] + 2A[3, 2] + A[3, 3]) - (A[1, 1] + 2A[1, 2] + A[1, 3])) / (8 * cellsize) - return δzδx, δzδy +@inline @inbounds function zevenbergenthorne(v, a, b) + if b == 2 + return (v[1], v[2] + a, v[3]) + elseif b == 4 + return (v[1] + a, v[2], v[3]) + elseif b == 6 + return (v[1] - a, v[2], v[3]) + elseif b == 8 + return (v[1], v[2] - a, v[3]) + else + return v + end end -function derivative(::MDG, A, cellsize=1.0) - δzδx = max( - abs(A[1, 1] - A[3, 3]) / (cellsize * sqrt2), - abs(A[1, 2] - A[3, 2]) / cellsize, - abs(A[1, 3] - A[3, 1]) / (cellsize * sqrt2), - abs(A[2, 1] - A[2, 3]) / cellsize, - ) / 2 +@inline @inbounds function horn(v, a, b) + if b == 1 + return (v[1], v[2] + a, v[3], v[4] + a, v[5]) + elseif b == 2 + return (v[1], v[2] + 2a, v[3], v[4], v[5]) + elseif b == 3 + return (v[1], v[2] + a, v[3] + a, v[4], v[5]) + elseif b == 4 + return (v[1], v[2], v[3], v[4] + 2a, v[5]) + elseif b == 5 + return v + elseif b == 6 + return (v[1], v[2], v[3] + 2a, v[4], v[5]) + elseif b == 7 + return (v[1] + a, v[2], v[3], v[4] + a, v[5]) + elseif b == 8 + return (v[1] + 2a, v[2], v[3], v[4], v[5]) + elseif b == 9 + return (v[1] + a, v[2], v[3] + a, v[4], v[5]) + end +end - return δzδx, δzδx +@inbounds function mdg(v, a, b) + if b == 1 + return (v[1] + a, v[2], v[3], v[4], v[5]) + elseif b == 2 + return (v[1], v[2], v[3], v[4] + a, v[5]) + elseif b == 3 + return (v[1], v[2], v[3] - a, v[4], v[5]) + elseif b == 4 + return (v[1], v[2] + a, v[3], v[4], v[5]) + elseif b == 5 + return v + elseif b == 6 + return (v[1], v[2] - a, v[3], v[4], v[5]) + elseif b == 7 + return (v[1], v[2], v[3] + a, v[4], v[5]) + elseif b == 8 + return (v[1], v[2], v[3], v[4] - a, v[5]) + elseif b == 9 + return (v[1] - a, v[2], v[3], v[4], v[5]) + end end function compass(aspect) @@ -176,7 +241,7 @@ function compass(aspect) end function aspect(compass::Real) - return (360 - compass + 90) % 360 + return (450 - compass) % 360 end """ @@ -185,44 +250,106 @@ end Curvature is derivative of [`slope`](@ref), so the second derivative of the DEM. """ function curvature(dem::AbstractMatrix{<:Real}; cellsize=1.0) - return terrain_kernel(dem, f -> curvature_kernel(f, cellsize)) -end - -function curvature_kernel(A, cellsize=1.0) - δzδx = ((A[1, 2] + A[3, 2]) / 2 - A[2, 2]) / cellsize^2 - δzδy = ((A[2, 1] + A[2, 3]) / 2 - A[2, 2]) / cellsize^2 + dst = copy(dem) - return -2(δzδx + δzδy) * 100 + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function curvature(v, a, b) + if b == 2 + return (v[1], v[2] + a, v[3], v[4]) + elseif b == 4 + return (v[1] + a, v[2], v[3], v[4]) + elseif b == 5 + return (v[1], v[2], v[3] + a, v[4]) + elseif b == 6 + return (v[1] + a, v[2], v[3], v[4]) + elseif b == 8 + return (v[1], v[2] + a, v[3], v[4]) + else + return v + end + end + function store!(d, i, v) + δzδx = (v[1] / 2 - v[3]) / v[4]^2 + δzδy = (v[2] / 2 - v[3]) / v[4]^2 + d[i] = -2(δzδx + δzδy) * 100 + end + return localfilter!(dst, dem, nbkernel, initial, curvature, store!) end + """ - hillshade(dem::Matrix{<:Real}; azimuth=315.0, zenith=45.0, cellsize=1.0, method=Horn()) + hillshade(dem::Matrix{<:Real}; azimuth=315.0, zenith=45.0, cellsize=1.0) hillshade is the simulated illumination of a surface based on its [`slope`](@ref) and [`aspect`](@ref) given a light source with azimuth and zenith angles in °, , as defined in Burrough, P. A., and McDonell, R. A., (1998, Principles of Geographical Information Systems). """ -function hillshade(dem::AbstractMatrix{<:Real}; azimuth=315.0, zenith=45.0, cellsize=1.0, method::DerivativeMethod=Horn()) - return terrain_kernel(dem, f -> hillshade_kernel(f, cellsize, azimuth, zenith, method), UInt8) +function hillshade(dem::AbstractMatrix{<:Real}; azimuth=315.0, zenith=45.0, cellsize=1.0) + dst = similar(dem, UInt8) + zenithr = deg2rad(zenith) + azimuthr = deg2rad(aspect(azimuth)) + + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + δzδx, δzδy = (v[1] - v[2]) / (8 * v[5]), (v[3] - v[4]) / (8 * v[5]) + if δzδx != 0 + a = atan(-δzδx, δzδy) + if a < 0 + a += 2π + end + else + a = π / 2 + if δzδy < 0 + a += 2π + end + end + slope = atan(√(δzδx^2 + δzδy^2)) + d[i] = round(UInt8, max(0, 255 * ((cos(zenithr) * cos(slope)) + (sin(zenithr) * sin(slope) * cos(azimuthr - a))))) + end + return localfilter!(dst, dem, nbkernel, initial, horn, store!) end -function hillshade_kernel(A, cellsize=1.0, azimuth=315.0, zenith=45.0, method::DerivativeMethod=Horn()) - z = deg2rad(zenith) - c = deg2rad(aspect(azimuth)) +""" + multihillshade(dem::Matrix{<:Real}; cellsize=1.0) - δzδx, δzδy = derivative(method, A, cellsize) - if δzδx != 0 - a = atan(δzδy, -δzδx) - if a < 0 - a += 2π - end - else - a = π / 2 - if δzδy > 0 - a += 2π +multihillshade is the simulated illumination of a surface based on its [`slope`](@ref) and +[`aspect`](@ref). Like [`hillshade`](@ref), but now using multiple sources as defined in +https://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf, similar to GDALs -multidirectional. +""" +function multihillshade(dem::AbstractMatrix{<:Real}; cellsize=1.0) + dst = similar(dem, UInt8) + zenithr = deg2rad(60) + + initial(A) = (zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), zero(eltype(A)), cellsize) + function store!(d, i, v) + δzδx, δzδy = (v[1] - v[2]) / (8 * v[5]), (v[3] - v[4]) / (8 * v[5]) + if δzδx != 0 + a = atan(-δzδx, δzδy) + if a < 0 + a += 2π + end + else + a = π / 2 + if δzδy < 0 + a += 2π + end end + slope = atan(√(δzδx^2 + δzδy^2)) + + w225 = 0.5 * (1 - cos(2(a - deg2rad(aspect(225))))) + w270 = 0.5 * (1 - cos(2(a - deg2rad(aspect(270))))) + w315 = 0.5 * (1 - cos(2(a - deg2rad(aspect(315))))) + w360 = 0.5 * (1 - cos(2(a - deg2rad(aspect(360))))) + + α = cos(zenithr) * cos(slope) + β = sin(zenithr) * sin(slope) + something = ( + w225 * (α + β * cos(deg2rad(aspect(225)) - a)) + + w270 * (α + β * cos(deg2rad(aspect(270)) - a)) + + w315 * (α + β * cos(deg2rad(aspect(315)) - a)) + + w360 * (α + β * cos(deg2rad(aspect(360)) - a))) / 2 + + d[i] = round(UInt8, max(0, 255 * something)) end - - s = atan(√(δzδx^2 + δzδy^2)) - return round(UInt8, 255 * ((cos(z) * cos(s)) + (sin(z) * sin(s) * cos(c - a)))) + return localfilter!(dst, dem, nbkernel, initial, horn, store!) end diff --git a/src/utils.jl b/src/utils.jl index 438e544..387b847 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,8 +14,8 @@ 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 + mapwindowcirc_approx!(minimum_mask, A, ω, out, Inf) # erosion + mapwindowcirc_approx!(maximum_mask, out, ω, A, -Inf) # dilation A end @@ -120,6 +120,40 @@ function mapwindowcirc_approx!(f, img, window, out, fill=Inf) out end +# function opening_circ_approx2!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real} +# iterations = ω:-2:3 + +# B = circmask(1) +# for _ in iterations +# localfilter!(out, A, B, +# (a) -> typemax(a), +# (v, a, b) -> b ? min(v, a) : v, +# (d, i, v) -> @inbounds(d[i] = v)) +# A, out = out, A +# end +# for _ in iterations +# localfilter!(out, A, B, +# (a) -> typemin(a), +# (v, a, b) -> b ? max(v, a) : v, +# (d, i, v) -> @inbounds(d[i] = v)) +# A, out = out, A +# end +# A +# end + +# function opening_circ2!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real} +# B = circmask(ω >> 1) +# localfilter!(out, A, B, +# (a) -> typemax(a), +# (v, a, b) -> b ? min(v, a) : v, +# (d, i, v) -> @inbounds(d[i] = v)) +# localfilter!(A, out, B, +# (a) -> typemin(a), +# (v, a, b) -> b ? max(v, a) : v, +# (d, i, v) -> @inbounds(d[i] = v)) +# A +# end + @inline function maximum_mask(x, m) o = -Inf @inbounds for I in eachindex(x) From 3bb64c7f65335426efbacbd029db9b4b7fe39617 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 2 Nov 2022 18:10:55 +0100 Subject: [PATCH 2/2] Small cleanup. --- src/GeoArrayOps.jl | 2 +- src/terrain.jl | 2 -- src/utils.jl | 1 + 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl index cb7b645..f88038c 100644 --- a/src/GeoArrayOps.jl +++ b/src/GeoArrayOps.jl @@ -25,7 +25,7 @@ export pmf, smf, psf export pssm export pitremoval export spread, spread2 -export roughness, TRI, TPI, slope, aspect, curvature, hillshade +export roughness, TRI, TPI, slope, aspect, curvature, hillshade, multihillshade export skb, skbr end # module diff --git a/src/terrain.jl b/src/terrain.jl index f195435..e076cc8 100644 --- a/src/terrain.jl +++ b/src/terrain.jl @@ -159,7 +159,6 @@ function aspect(::Horn, dst, dem::AbstractMatrix{<:Real}, cellsize) function store!(d, i, v) δzδx = (v[1] - v[2]) / (8 * v[5]) δzδy = (v[3] - v[4]) / (8 * v[5]) - # @info δzδx, δzδy d[i] = compass(atand(-δzδx, δzδy)) end return localfilter!(dst, dem, nbkernel, initial, horn, store!) @@ -170,7 +169,6 @@ function aspect(::ZevenbergenThorne, dst, dem::AbstractMatrix{<:Real}, cellsize) function store!(d, i, v) δzδx = v[1] / 2 δzδy = v[2] / 2 - # @info δzδx, δzδy d[i] = compass(atand(δzδy, -δzδx)) end return localfilter!(dst, dem, nbkernel, initial, zevenbergenthorne, store!) diff --git a/src/utils.jl b/src/utils.jl index 387b847..cf34e67 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -120,6 +120,7 @@ function mapwindowcirc_approx!(f, img, window, out, fill=Inf) out end +# Functions for future changes, based on LocalFiltering # function opening_circ_approx2!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real} # iterations = ω:-2:3