# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1610084087 18000
# Node ID b32ced049cbd1365e27f969d4ce272e5a3080b24
# Parent  46053b3af25169f313c3fc21189831047fb3e504
Add TV reconstruction code for abstract linear operators

diff -r 46053b3af251 -r b32ced049cbd Manifest.toml
--- a/Manifest.toml	Thu Jan 07 17:52:42 2021 -0500
+++ b/Manifest.toml	Fri Jan 08 00:34:47 2021 -0500
@@ -181,10 +181,10 @@
 version = "1.0.2"
 
 [[MbedTLS_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "c0b1286883cac4e2b617539de41111e0776d02e8"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "0eef589dd1c26a3ac9d753fe1a8bcad63f956fa6"
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.16.8+0"
+version = "2.16.8+1"
 
 [[Missings]]
 deps = ["DataAPI"]
diff -r 46053b3af251 -r b32ced049cbd src/ImageTools.jl
--- a/src/ImageTools.jl	Thu Jan 07 17:52:42 2021 -0500
+++ b/src/ImageTools.jl	Fri Jan 08 00:34:47 2021 -0500
@@ -9,6 +9,7 @@
 include("Denoise.jl")
 include("Visualise.jl")
 include("ImFilter.jl")
+include("TVRecon.jl")
 
 end # module
 
diff -r 46053b3af251 -r b32ced049cbd src/TVRecon.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/TVRecon.jl	Fri Jan 08 00:34:47 2021 -0500
@@ -0,0 +1,121 @@
+########################################################
+# TV reconstruction via primal–dual proximal splitting
+########################################################
+
+__precompile__()
+
+module TVRecon
+
+using AlgTools.Util
+using AlgTools.LinOps
+import AlgTools.Iterate
+using ImageTools.Gradient
+
+##############
+# Our exports
+##############
+
+export recon_pdps
+
+#############
+# Data types
+#############
+
+ImageSize = Tuple{Integer,Integer}
+Image = Array{Float64,2}
+Primal = Image
+Dual = Array{Float64,3}
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_rest(x::Primal, b::Array{Float64, N}) where N
+    imdim=size(x)
+    datadim=size(b)
+
+    y = zeros(2, imdim...)
+    λ = zeros(datadim...)
+    Δx₁ = copy(x)
+    Δx₂ = copy(x)
+    Δy = copy(y)
+    Δλ = copy(λ)
+    x̄ = copy(x)
+
+    return x, y, λ, Δx₁, Δx₂, Δy, Δλ, x̄
+end
+
+function init_primal(xinit::Image)
+    return copy(xinit)
+end
+
+############
+# Algorithm
+############
+
+function recon_pdps(b :: Data, op :: LinOp{Image, Data};
+                    xinit :: Union{Image, Nothing} = nothing,
+                    iterate = Iterate.simple_iterate,
+                    params::NamedTuple) where Data
+
+    ################################
+    # Extract and set up parameters
+    ################################
+
+    α, ρ = params.α, params.ρ
+    τ₀, σ₀ =  params.τ₀, params.σ₀
+
+    R_K = √(∇₂_norm₂₂_est^2+opnorm_estimate(op)^2)
+    γ = 1
+
+    @assert(τ₀*σ₀ < 1)
+    σ = σ₀/R_K
+    τ = τ₀/R_K
+    
+    ######################
+    # Initialise iterates
+    ######################
+
+    x, y, λ, Δx₁, Δx₂, Δy, Δλ, x̄ = init_rest(copy(xinit), b)
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    v = iterate(params) do verbose :: Function
+        ω = params.accel ? 1/√(1+2*γ*τ) : 1
+        
+        ∇₂ᵀ!(Δx₁, y)                      # primal step: x
+        inplace!(Δx₂, op', λ)             # |
+        @. x̄ = x                          # |
+        @. x = x-τ*(Δx₁+Δx₂)              # |
+        @. x̄ = (1+ω)*x - ω*x̄              # over-relax: x̄ = 2x-x_old
+        ∇₂!(Δy, x̄)                        # dual step: y, λ
+        inplace!(Δλ, op, x̄)               # |
+        @. y = (y + σ*Δy)/(1 + σ*ρ/α)     # |
+        proj_norm₂₁ball!(y, α)            # |
+        @. λ = (λ + σ*(Δλ-b))/(1 + σ)     # |
+
+        if params.accel
+            τ, σ = τ*ω, σ/ω
+        end
+                
+        ################################
+        # Give function value if needed
+        ################################
+        v = verbose() do
+            ∇₂!(Δy, x)
+            inplace!(Δλ, op, x)
+            value = norm₂²(Δλ-b)/2 + params.α*γnorm₂₁(Δy, params.ρ)
+            value, x
+        end
+
+        v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+