diff --git a/Project.toml b/Project.toml
index 0e50687..b9764fe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GeoArrayOps"
 uuid = "e6005643-8f00-58df-85c5-7d93a3520d70"
 authors = ["Maarten Pronk <git@evetion.nl>", "Deltares"]
-version = "0.3.2"
+version = "0.4.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -9,6 +9,7 @@ Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
+LocalFilters = "085fde7c-5f94-55e4-8448-8bbb5db6dde9"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -21,11 +22,12 @@ Distances = "0.10"
 FillArrays = "0.12, 0.13"
 ImageCore = "0.9"
 ImageFiltering = "0.6, 0.7"
+LocalFilters = "1.2"
 OffsetArrays = "1.10"
 PaddedViews = "0.5"
 StaticArrays = "1"
 StatsBase = "0.33"
-julia = "1.5, 1.6, 1.7, 1.8"
+julia = "1.6, 1.7, 1.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index ebb3fb3..39cd96b 100644
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@ Geospatial operations, cost and filtering algorithms as used in for elevation ra
 - 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 slope, aspect, roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI).
+- Terrain analysis functions, such as slope, aspect, roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI), curvature and hillslope.
 
 ## Installation
 The package can be installed with the Julia package manager.
diff --git a/docs/src/index.md b/docs/src/index.md
index 6696ed8..d5e9f95 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -8,7 +8,8 @@ Geospatial operations, cost and filtering algorithms as used in for elevation ra
 - 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 slope, aspect, roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI).
+- Terrain analysis functions, such as slope, aspect, roughness, Topographic Position Index (TPI), Terrain Ruggedness Index (TRI), curvature and hillslope.
+
 ## Installation
 The package can be installed with the Julia package manager.
 From the Julia REPL, type `]` to enter the Pkg REPL mode and run:
diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl
index a5fcc3c..cb7b645 100644
--- a/src/GeoArrayOps.jl
+++ b/src/GeoArrayOps.jl
@@ -9,6 +9,7 @@ using StaticArrays: @SMatrix, @MMatrix, SMatrix, MMatrix, MVector
 using DataStructures: Deque
 using Statistics: median, mean
 using ImageCore: scaleminmax, Gray
+using LocalFilters
 
 include("utils.jl")
 include("pmf.jl")
@@ -19,11 +20,12 @@ include("spread.jl")
 include("terrain.jl")
 include("skew.jl")
 
+export ZevenbergenThorne, Horn, MDG
 export pmf, smf, psf
 export pssm
 export pitremoval
 export spread, spread2
-export roughness, TRI, TPI, slope, aspect
-export skb
+export roughness, TRI, TPI, slope, aspect, curvature, hillshade
+export skb, skbr
 
 end # module
diff --git a/src/pmf.jl b/src/pmf.jl
index 6f6a7e4..77f0a3e 100644
--- a/src/pmf.jl
+++ b/src/pmf.jl
@@ -36,7 +36,6 @@ 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"
 
     # Set up arrays
     Af = copy(A)  # array to be morphed
@@ -56,7 +55,7 @@ function pmf(A::AbstractMatrix{T};
         if circular
             opening_circ!(Af, ωₖ, out)
         else
-            opening!(Af, ωₖ, out)
+            LocalFilters.opening!(Af, out, Af, ωₖ)
         end
         mask .= (A .- Af) .> dhₜ
         for I in eachindex(flags)
diff --git a/src/skew.jl b/src/skew.jl
index b4dd46d..65af3dd 100644
--- a/src/skew.jl
+++ b/src/skew.jl
@@ -1,6 +1,5 @@
 """
-    mask = skb(A, mean)
-    mask = skb(A)
+    mask = skb(A; mean=mean(A))
 
 Applies skewness balancing by *Bartels e.a (2006)* [^bartels2006] to `A`.
 Improved the performance by applying a binary search to find the threshold value.
@@ -10,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) where {T<:Real}
     m = .!isfinite.(iA)
     if sum(m) > 0
         A = copy(iA)
@@ -47,5 +46,29 @@ function skb(iA::AbstractArray{T}, mean::T) where {T<:Real}
 end
 
 function skb(A::AbstractArray{T}) where {T<:Real}
-    return skb(A, mean(A))
+    return skb(A; mean=mean(A))
+end
+
+"""
+    mask = skbr(A; iterations=10)
+
+Applies recursive skewness balancing by *Bartels e.a (2006)* [^bartels2006] to `A`.
+Applies `skb` `iterations` times to the object (non-terrain) mask, as to include
+more (sloped) terrain.
+
+# Output
+- `mask::BitMatrix` Mask of allowed values
+
+[^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 skbr(A; iterations=10)
+    terrain_mask = skb(A)
+    object_mask = .!terrain_mask
+    while iterations > 1 && sum(object_mask) > 0
+        # @info "$(round(Int, sum(object_mask) / length(object_mask) * 100))% objects..."
+        terrain_mask[object_mask] = terrain_mask[object_mask] .| skb(A[object_mask])
+        object_mask .= .!terrain_mask
+        iterations -= 1
+    end
+    terrain_mask
 end
diff --git a/src/terrain.jl b/src/terrain.jl
index 6c6993a..4ec8104 100644
--- a/src/terrain.jl
+++ b/src/terrain.jl
@@ -2,6 +2,21 @@ 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]
 
+abstract type DerivativeMethod end
+
+"""Second order finite difference estimator using all 4 neighbors (Zevenbergen and Thorne, 1987)."""
+struct ZevenbergenThorne <: DerivativeMethod end
+
+"""Third order finite difference estimator using all 8 neighbors (Horn, 1981)."""
+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
@@ -12,16 +27,16 @@ function buffer_array(A::AbstractMatrix{<:Real})
     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]
-    oA[end, begin] = A[end, begin]
-    oA[end, end] = A[end, end]
+    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{<:Real}, f::Function)
+function terrain_kernel(dem::AbstractMatrix{T}, f::Function, t=T) where {T<:Real}
     dembuffer = buffer_array(dem)
-    output = similar(dem)
+    output = similar(dem, t)
     Δ = CartesianIndex(1, 1)
 
     @inbounds for I ∈ CartesianIndices(size(output))
@@ -38,12 +53,13 @@ end
 Roughness is the largest inter-cell difference of a central pixel and its surrounding cell, as defined in Wilson et al (2007, Marine Geodesy 30:3-35).
 """
 function roughness(dem::AbstractMatrix{<:Real})
-    terrain_kernel(dem, roughness_kernel)
-end
+    dst = copy(dem)
+
+    initial(a) = (zero(a), a)
+    update(v, a, _) = (max(v[1], abs(a - v[2])), v[2])
+    store!(d, i, v) = @inbounds d[i] = v[1]
 
-function roughness_kernel(A)
-    x = A .- A[2, 2]
-    return maximum(abs.(x))
+    return localfilter!(dst, dem, 3, initial, update, store!)
 end
 
 """
@@ -52,14 +68,16 @@ end
 TPI stands for Topographic Position Index, which is defined as the difference between a central pixel and the mean of its surrounding cells (see Wilson et al 2007, Marine Geodesy 30:3-35).
 """
 function TPI(dem::AbstractMatrix{<:Real})
-    terrain_kernel(dem, TPI_kernel)
-end
+    dst = copy(dem)
 
-function TPI_kernel(A)
-    x = A .* neib_8
-    return A[2, 2] - sum(x) / 8
+    initial(a) = zero(a)
+    update(v, a, b) = v + a
+    store!(d, i, v) = @inbounds d[i] = d[i] - (v - d[i]) / 8
+
+    return localfilter!(dst, dem, 3, initial, update, store!)
 end
 
+
 """
     TRI(dem::Matrix{<:Real})
 
@@ -68,65 +86,143 @@ This algorithm uses the square root of the sum of the square of the difference b
 This is recommended for terrestrial use cases.
 """
 function TRI(dem::AbstractMatrix{<:Real})
-    terrain_kernel(dem, TRI_kernel)
-end
+    dst = copy(dem)
 
-function TRI_kernel(A)
-    x = (A .- A[2, 2]) .^ 2
-    return sqrt(sum(x))
-end
+    initial(a) = (zero(a), a)
+    update(v, a, _) = (v[1] + (a - v[2])^2, v[2])
+    store!(d, i, v) = @inbounds d[i] = sqrt(v[1])
 
+    return localfilter!(dst, dem, 3, initial, update, store!)
+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)
-    terrain_kernel(dem, f -> _pitremoval(f, limit))
-end
+function pitremoval(dem::AbstractMatrix{<:Real}; limit=0.0)
+    dst = copy(dem)
+
+    initial(a) = (true, typemax(a), a, limit)  # (is_pit, min, center, limit)
+    @inbounds function update(v, a, _)
+        if v[3] == a
+            return v
+        elseif (a - v[3]) > v[4]
+            return (v[1] & true, min(a, v[2]), v[3], v[4])
+        else
+            (false, v[2], v[3], v[4])
+        end
+    end
+    store!(d, i, v) = @inbounds d[i] = v[1] ? v[2] : v[3]
 
-function _pitremoval(A, limit)
-    order = sortperm(vec(A))
-    return (order[1] == 5) && ((A[order[2]] - A[order[1]]) >= limit) ? Inf : A[2, 2]
+    return localfilter!(dst, dem, 3, initial, update, store!)
 end
 
 """
-    slope(dem::Matrix{<:Real})
+    slope(dem::Matrix{<:Real}, cellsize=1.0, method=Horn())
 
 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)
-    terrain_kernel(dem, f -> slope_kernel(f, cellsize))
+function slope(dem::AbstractMatrix{<:Real}; cellsize=1.0, method::DerivativeMethod=Horn())
+    terrain_kernel(dem, f -> slope_kernel(f, cellsize, method))
+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))
 end
 
 """
-    aspect(dem::Matrix{<:Real})
+    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})
-    terrain_kernel(dem, aspect_kernel)
+function aspect(dem::AbstractMatrix{<:Real}; method::DerivativeMethod=Horn())
+    terrain_kernel(dem, f -> aspect_kernel(f, method))
 end
 
-function slope_kernel(A, cellsize=1.0)
-    δzδx, δzδy = derivative(A, cellsize)
-    return atand(√(δzδx^2 + δzδy^2))
+function aspect_kernel(A, method::DerivativeMethod=Horn())
+    δzδx, δzδy = derivative(method, A, 1.0)
+    return compass(atand(δzδy, -δzδx))
 end
 
-function aspect_kernel(A)
-    δzδx, δzδy = derivative(A)
-    return compass(atand(δzδy, -δzδx))
+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
 end
 
-function derivative(A, cellsize=1.0)
+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
 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
+
+    return δzδx, δzδx
+end
+
 function compass(aspect)
     a = 90 - aspect
     a < 0 && (a += 360)
     return a
 end
+
+function aspect(compass::Real)
+    return (360 - compass + 90) % 360
+end
+
+"""
+    curvature(dem::Matrix{<:Real})
+
+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
+
+    return -2(δzδx + δzδy) * 100
+end
+
+"""
+    hillshade(dem::Matrix{<:Real}; azimuth=315.0, zenith=45.0, cellsize=1.0, method=Horn())
+
+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)
+end
+
+function hillshade_kernel(A, cellsize=1.0, azimuth=315.0, zenith=45.0, method::DerivativeMethod=Horn())
+    z = deg2rad(zenith)
+    c = deg2rad(aspect(azimuth))
+
+    δ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π
+        end
+    end
+
+    s = atan(√(δzδx^2 + δzδy^2))
+    return round(UInt8, 255 * ((cos(z) * cos(s)) + (sin(z) * sin(s) * cos(c - a))))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d651d7c..136671c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,9 +23,14 @@ using Test
     @testset "psf" begin
         B = psf(rand(25, 25))
     end
+    @testset "skb" begin
+        B = skb(rand(25, 25))
+        B = skb(rand(25, 25); mean=0.25)
+        B = skbr(rand(25, 25); iterations=5)
+    end
     @testset "pitremoval" begin
         B = pitremoval(rand(25, 25))
-        B = pitremoval(rand(25, 25), 0.1)
+        B = pitremoval(rand(25, 25); limit=0.1)
     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]
@@ -38,7 +43,12 @@ using Test
         TRI(A)
         TPI(A)
         roughness(A)
-        slope(A)
+        slope(A; method=Horn())
+        slope(A; cellsize=5, method=ZevenbergenThorne())
+        slope(A; cellsize=10, method=MDG())
         aspect(A)
+        aspect(A, method=MDG())
+        curvature(A)
+        hillshade(A)
     end
 end