diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 9d798ea..df412de 100644
--- a/.github/workflows/CI.yml
+++ b/.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
diff --git a/.gitignore b/.gitignore
index c531dd1..3176d51 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,7 @@
 .DS_Store
 Manifest.toml
 docs/build/
+*.tif
+*.gif
+test/benchmark.jl
+.vscode
diff --git a/Project.toml b/Project.toml
index ca0628e..49604dc 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.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"
diff --git a/README.md b/README.md
index c7fecdc..98e32d2 100644
--- a/README.md
+++ b/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
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 69db735..3c94b9d 100644
--- a/docs/src/index.md
+++ b/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
diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl
index 9050794..55396b7 100644
--- a/src/GeoArrayOps.jl
+++ b/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
diff --git a/src/plot.jl b/src/plot.jl
index 6668bdc..002777b 100644
--- a/src/plot.jl
+++ b/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
diff --git a/src/pmf.jl b/src/pmf.jl
index 6cf27d1..290ab2b 100644
--- a/src/pmf.jl
+++ b/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
diff --git a/src/psf.jl b/src/psf.jl
index 9f5fab1..91e29e6 100644
--- a/src/psf.jl
+++ b/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
diff --git a/src/skew.jl b/src/skew.jl
new file mode 100644
index 0000000..8d9041c
--- /dev/null
+++ b/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
diff --git a/src/smf.jl b/src/smf.jl
index dfd794a..3bbeadb 100644
--- a/src/smf.jl
+++ b/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))
diff --git a/src/spread.jl b/src/spread.jl
index 18b7da9..55ee74f 100644
--- a/src/spread.jl
+++ b/src/spread.jl
@@ -1,11 +1,3 @@
-
-using FillArrays
-using StaticArrays
-using Distances
-using OffsetArrays
-using DataStructures
-using Statistics
-
 const sqrt2 = sqrt(2.0)
 const distance_8 = @SMatrix[sqrt2 1 sqrt2; 1 Inf 1; sqrt2 1 sqrt2]
 const distance_4 = @SMatrix[Inf 1 Inf; 1 Inf 1; Inf 1 Inf]
@@ -33,16 +25,16 @@ This is also the method implemented by [PCRaster](https://pcraster.geo.uu.nl/pcr
 function spread(points::AbstractMatrix{<:Real}, initial::AbstractMatrix{<:Real}, friction::AbstractMatrix{<:Real}; res=1, limit=Inf)
 
     ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
-    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction
+    ofriction[begin+1:end-1, begin+1:end-1] .= friction
 
     result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
-    r = @view result[1:end - 1,1:end - 1]
+    r = @view result[1:end-1, 1:end-1]
     locations = points .> 0
     r[locations] .= initial[locations]
 
     # Construct stack for locations
     mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
-    mask[begin + 1:end - 1,begin + 1:end - 1] .= false
+    mask[begin+1:end-1, begin+1:end-1] .= false
 
     II = CartesianIndices(size(points))
     stack = Deque{CartesianIndex}()
@@ -67,14 +59,14 @@ end
 function spread!(stack, mask, result, ofriction, sdata, mcell, res)
     I = popfirst!(stack)
     mask[I] = true
-    patch = I - Δ:I + Δ
+    patch = I-Δ:I+Δ
 
     rdata = view(result, patch)
     fdata = view(ofriction, patch)
 
     # New distance is cell_distance + average friction values
     for i ∈ eachindex(sdata)
-        sdata[i] = muladd(fdata[i] + fdata[2,2], res / 2 * distance_8[i], rdata[2,2])
+        sdata[i] = muladd(fdata[i] + fdata[2, 2], res / 2 * distance_8[i], rdata[2, 2])
         mcell[i] = sdata[i] < rdata[i]  # cells where new distance is lower
     end
     rdata[mcell] .= sdata[mcell]
@@ -111,31 +103,31 @@ This method should scale much better (linearly) than the [^tomlin1983] method, b
 function spread2(points::AbstractMatrix{<:Real}, initial::AbstractMatrix{<:Real}, friction::AbstractMatrix{<:Real}; res=1, limit=Inf, iterations=3)
 
     ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
-    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction
+    ofriction[begin+1:end-1, begin+1:end-1] .= friction
 
     result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
-    r = @view result[1:end - 1,1:end - 1]
+    r = @view result[1:end-1, 1:end-1]
     locations = points .> 0
     r[locations] .= initial[locations]
 
     mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
-    mask[begin + 1:end - 1,begin + 1:end - 1] .= false
+    mask[begin+1:end-1, begin+1:end-1] .= false
 
-    minval, minidx = [0.], [CartesianIndex(1, 1)]
+    minval, minidx = [0.0], [CartesianIndex(1, 1)]
     x = @MMatrix zeros(3, 3)
 
     indices = CartesianIndices(size(points))
     for i in 1:iterations
         II = (i % 2 == 1) ? indices : reverse(indices)
         for I ∈ II
-            patch = I - Δ:I + Δ
+            patch = I-Δ:I+Δ
 
             rdata = view(result, patch)
             fdata = view(ofriction, patch)
 
-            x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
+            x .= (fdata .+ fdata[2, 2]) .* res ./ 2 .* distance_8 .+ rdata
             findmin!(minval, minidx, x)
-            rdata[2,2] = min(rdata[2,2], minval[1])
+            rdata[2, 2] = min(rdata[2, 2], minval[1])
         end
     end
     r
diff --git a/src/terrain.jl b/src/terrain.jl
index 06b47ae..11f96f9 100644
--- a/src/terrain.jl
+++ b/src/terrain.jl
@@ -91,7 +91,7 @@ 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)
+function pitremoval(dem::AbstractMatrix{<:Real}, limit=2.0)
 
     ex_dem = buffer_array(dem)
     tri = similar(dem)
diff --git a/src/utils.jl b/src/utils.jl
index 2c80d83..438e544 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,8 +1,3 @@
-using OffsetArrays
-using ImageFiltering
-using Distances
-using PaddedViews
-
 """Apply the opening operation to `A` with window size `ω`."""
 function opening(A::Array{T,2}, ω::Integer) where {T<:Real}
     A = mapwindow(minimum, A, (ω, ω))  # erosion
@@ -70,7 +65,7 @@ function mapwindow_stack!(f, img, window, out)
     out
 end
 
-function mapwindow_sep!(f, img, window, out, fill = Inf)
+function mapwindow_sep!(f, img, window, out, fill=Inf)
     Δ = window ÷ 2
 
     w, h = size(img)
@@ -89,7 +84,7 @@ function mapwindow_sep!(f, img, window, out, fill = Inf)
 end
 
 
-function mapwindowcirc!(f, img, window, out, fill = Inf)
+function mapwindowcirc!(f, img, window, out, fill=Inf)
     R = CartesianIndices(img)
     Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))
 
@@ -105,7 +100,7 @@ function mapwindowcirc!(f, img, window, out, fill = Inf)
 end
 
 
-function mapwindowcirc_approx!(f, img, window, out, fill = Inf)
+function mapwindowcirc_approx!(f, img, window, out, fill=Inf)
     R = CartesianIndices(img)
     Δ = CartesianIndex(1, 1)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 744cec9..b164f6c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,39 +1,41 @@
 using GeoArrayOps
 using Test
 
-@testset "pmf" begin
-    # Write your own tests here.
-    A = rand(25, 25)
-    A[2, 2] = NaN
-    B, flags = pmf(A)
-    @test (A .<= B) == (flags .== 0.0)
-    B, flags = pmf(A; circular = true)
-    @test (A .<= B) == (flags .== 0.0)
-end
-@testset "smf" begin
-    # Write your own tests here.
-    A = rand(25, 25)
-    A[2, 2] = NaN
-    B = smf(A)
-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]
-    friction = [1.0 200 1 1 1; 200 1 1 4 4; 1 1 4 4 4; 1 1 3 200 200; 1 Inf 3 200 4]
-    @test spread(points, initial, friction) == spread2(points, initial, friction)
-end
-@testset "terrain" begin
-    A = rand(25, 25)
-    TRI(A)
-    TPI(A)
-    roughness(A)
+@testset "GeoArrayOps" begin
+    @testset "pmf" begin
+        # Write your own tests here.
+        A = rand(25, 25)
+        A[2, 2] = NaN
+        B, flags = pmf(A)
+        @test (A .<= B) == (flags .== 0.0)
+        B, flags = pmf(A; circular=true)
+        @test (A .<= B) == (flags .== 0.0)
+    end
+    @testset "smf" begin
+        # Write your own tests here.
+        A = rand(25, 25)
+        A[2, 2] = NaN
+        B = smf(A)
+    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]
+        friction = [1.0 200 1 1 1; 200 1 1 4 4; 1 1 4 4 4; 1 1 3 200 200; 1 Inf 3 200 4]
+        @test spread(points, initial, friction) == spread2(points, initial, friction)
+    end
+    @testset "terrain" begin
+        A = rand(25, 25)
+        TRI(A)
+        TPI(A)
+        roughness(A)
+    end
 end