diff -r c5aabb2c41d9 -r 1cffd3d07fe2 src/Denoise.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Denoise.jl	Thu Nov 21 18:44:27 2019 -0500
@@ -0,0 +1,116 @@
+########################################################
+# Basic TV denoising via primal–dual proximal splitting
+########################################################
+
+__precompile__()
+
+module Denoise
+
+using AlgTools.Util
+import AlgTools.Iterate
+using ImageTools.Gradient
+
+##############
+# Our exports
+##############
+
+export denoise_pdps
+
+#############
+# Data types
+#############
+
+ImageSize = Tuple{Integer,Integer}
+Image = Array{Float64,2}
+Primal = Image
+Dual = Array{Float64,3}
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_rest(x::Primal)
+    imdim=size(x)
+
+    y = zeros(2, imdim...)
+    Δx = copy(x)
+    Δy = copy(y)
+    x̄ = copy(x)
+
+    return x, y, Δx, Δy, x̄
+end
+
+function init_iterates(xinit::Image, b)
+    return init_rest(copy(xinit))
+end
+
+function init_iterates(xinit::Nothing, b :: Image)
+    return init_rest(zeros(size(b)...))
+end
+
+############
+# Algorithm
+############
+
+function denoise_pdps(b :: Image;
+                      xinit :: Union{Image,Nothing} = nothing,
+                      iterate = AlgTools.simple_iterate,
+                      params::NamedTuple) where DisplacementT
+
+    ################################                                        
+    # Extract and set up parameters
+    ################################                    
+
+    α, ρ = params.α, params.ρ
+    τ₀, σ₀ =  params.τ₀, params.σ₀
+
+    R_K = ∇₂_norm₂₂_est
+    γ = 1
+
+    @assert(τ₀*σ₀ < 1)
+    σ = σ₀/R_K
+    τ = τ₀/R_K
+    
+    ######################
+    # Initialise iterates
+    ######################
+
+    x, y, Δx, Δy, x̄ = init_iterates(xinit, b)
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    v = iterate(params) do verbose :: Function
+        ω = params.accel ? 1/√(1+2*γ*τ) : 1
+        
+        ∇₂ᵀ!(Δx, y)                    # primal step:
+        @. x̄ = x                       # |  save old x for over-relax
+        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
+        @. x̄ = (1+ω)*x - ω*x̄           # over-relax: x̄ = 2x-x_old
+        ∇₂!(Δy, x̄)                     # dual step: y
+        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
+        proj_norm₂₁ball!(y, α)         # |  prox
+
+        if params.accel
+            τ, σ = τ*ω, σ/ω
+        end
+                
+        ################################
+        # Give function value if needed
+        ################################
+        v = verbose() do            
+            ∇₂!(Δy, x)
+            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
+            value, x
+        end
+
+        v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+