From d3be5b678ed8a8cd8e398c346b742c51ad27df89 Mon Sep 17 00:00:00 2001
From: Maarten Pronk <git@evetion.nl>
Date: Fri, 18 Feb 2022 14:40:30 +0100
Subject: [PATCH] First version of CSF

---
 Project.toml       |   5 +-
 src/GeoArrayOps.jl |   4 +-
 src/csf.jl         | 120 +++++++++++++++++++++++++++++++++++++++++++++
 test/Project.toml  |   6 +++
 4 files changed, 132 insertions(+), 3 deletions(-)
 create mode 100644 src/csf.jl
 create mode 100644 test/Project.toml

diff --git a/Project.toml b/Project.toml
index f54e005..17d4dd8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,13 +22,14 @@ FillArrays = "0.12"
 ImageCore = "0.9"
 ImageFiltering = "0.7"
 OffsetArrays = "1.10"
-ProgressMeter = "1.7"
 PaddedViews = "0.5"
+ProgressMeter = "1.7"
 StaticArrays = "1"
 julia = "1.5, 1.6, 1.7"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Aqua"]
diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl
index 9050794..b4ba4f8 100644
--- a/src/GeoArrayOps.jl
+++ b/src/GeoArrayOps.jl
@@ -5,12 +5,13 @@ using ProgressMeter
 include("utils.jl")
 include("pmf.jl")
 include("smf.jl")
+include("csf.jl")
 include("psf.jl")
 include("plot.jl")
 include("spread.jl")
 include("terrain.jl")
 
-export pmf, smf, psf
+export pmf, smf, psf, csf
 export pssm
 export pitremoval
 export spread, spread2
@@ -30,6 +31,7 @@ export roughness, TRI, TPI
 
 precompile(pmf, (Matrix{Float64},))
 precompile(psf, (Matrix{Float64},))
+precompile(csf, (Matrix{Float64},))
 precompile(pitremoval, (Matrix{Float64},))
 precompile(smf, (Matrix{Float64},))
 precompile(pssm, (Matrix{Float64},))
diff --git a/src/csf.jl b/src/csf.jl
new file mode 100644
index 0000000..2fa6608
--- /dev/null
+++ b/src/csf.jl
@@ -0,0 +1,120 @@
+# Cloth Simulation Filter (Zhang et al., 2016)"
+# using Plots
+
+const gravity = 9.85
+const smoothThreshold = 0.3
+CI = CartesianIndex
+
+const nb = (
+    (CI(0, 0), CI(1, 0)),
+    (CI(0, 0), CI(0, 1)),
+    (CI(0, 0), CI(1, 1)),
+    (CI(0, 0), CI(-1, 1)),
+    (CI(0, 0), CI(1, -1)),
+    (CI(0, 0), CI(2, 0)),
+    (CI(0, 0), CI(0, 2)),
+    (CI(0, 0), CI(2, 2)),
+    (CI(0, 0), CI(-2, 2)),
+    (CI(0, 0), CI(2, -2)),
+)
+
+function csf(A; iterations = 500, time_step = 0.25, rigidness = 3, postprocess = true, threshold = 0.05, damping = 0.01)
+
+    min_A, max_A = extrema(skipmissing(A))
+    cloth = copy(A)
+
+    nodata = ismissing.(A)
+    moveable = .~nodata
+    cloth[moveable] .= min_A
+    prev_cloth = copy(cloth)
+    acceleration = gravity * time_step * time_step
+    temp = zero(eltype(A))
+
+    rigidness_factor = (1 - 0.5^rigidness)
+
+    R = CartesianIndices(cloth)
+    w, h = size(A)
+
+    max_diff = zero(eltype(A))
+    prev_max_diff = typemax(max_diff)
+
+    @showprogress 1 "Cloth simulating..." for i ∈ 1:iterations
+
+        # Drop cloth each iteration (gravity)
+        @inbounds for I in R
+            if moveable[I]
+                temp = cloth[I]
+                move = (cloth[I] - prev_cloth[I]) * (1.0 - damping) + acceleration * time_step * time_step
+                cloth[I] += move
+                prev_cloth[I] = temp
+                # println("$move $(cloth[I] - prev_cloth[I]) $(acceleration * time_step * time_step)")
+            end
+        end
+        # plot((
+        #         heatmap(A; colorbar = false, clims = (min_A, max_A), c = :viridis),
+        #         heatmap(cloth; colorbar = false, clims = (min_A, max_A), c = :viridis),
+        #         heatmap(moveable; colorbar = false, clims = (0, 1), c = :grays)
+        #     )...,
+        #     layout = (1, 3), ratio = :equal)
+
+
+        # If we've reached the ground, set node to unmoveable
+        max_diff = zero(eltype(A))
+        @inbounds for I in R
+            if moveable[I]
+                diff = abs(prev_cloth[I] - cloth[I])
+                max_diff = max(max_diff, diff)
+
+                if cloth[I] > A[I]
+                    cloth[I] = A[I]
+                    moveable[I] = false
+                end
+            end
+        end
+        # println("$max_diff $prev_max_diff")
+
+        # Stop if the cloth stopped moving
+        if (max_diff < (threshold / 100)) || (abs(prev_max_diff - max_diff) < (threshold / 100))
+            @info "Done in $i iterations."
+            break
+        end
+        prev_max_diff = max_diff
+
+        # Move each cloth node based on spring interaction with neighbors
+        @inbounds for I in CartesianIndices((3:w-2, 3:h-2))
+
+            # Only get graph, so
+            for (a, b) in nb
+                ai, bi = I + a, I + b
+
+                # vertical difference between nodes
+                correction = cloth[bi] - cloth[ai]
+                correction *= rigidness_factor
+
+                if (moveable[ai] && moveable[bi])
+                    cloth[ai] += correction * 0.5
+                    cloth[bi] -= correction * 0.5
+                elseif (moveable[ai] && !moveable[bi] && !nodata[bi])
+                    cloth[ai] += correction
+                elseif (!moveable[ai] && !nodata[ai] && moveable[bi])
+                    cloth[bi] -= correction
+                end
+            end
+        end
+
+
+
+        # plot((
+        #         heatmap(A; colorbar = false, clims = (min_A, max_A), c = :viridis),
+        #         heatmap(cloth; colorbar = false, clims = (min_A, max_A), c = :viridis),
+        #         heatmap(moveable; colorbar = false, clims = (0, 1), c = :grays)
+        #     )...,
+        #     layout = (1, 3), ratio = :equal)
+
+        # if (params.bSloopSmooth) {
+        #     cloth.movableFilter();
+        # }
+    end
+
+    return cloth, moveable
+end
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..a48142a
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,6 @@
+[deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+Aqua = "0.5"