Mercurial > repos > ImageTools / file comparison
comparison: src/TVRecon.jl
src/TVRecon.jl
- changeset 47
- b32ced049cbd
- child 53
- f8a3bc920f6a
equal
deleted
inserted
replaced
| 46:46053b3af251 | 47:b32ced049cbd |
|---|---|
| 1 ######################################################## | |
| 2 # TV reconstruction via primal–dual proximal splitting | |
| 3 ######################################################## | |
| 4 | |
| 5 __precompile__() | |
| 6 | |
| 7 module TVRecon | |
| 8 | |
| 9 using AlgTools.Util | |
| 10 using AlgTools.LinOps | |
| 11 import AlgTools.Iterate | |
| 12 using ImageTools.Gradient | |
| 13 | |
| 14 ############## | |
| 15 # Our exports | |
| 16 ############## | |
| 17 | |
| 18 export recon_pdps | |
| 19 | |
| 20 ############# | |
| 21 # Data types | |
| 22 ############# | |
| 23 | |
| 24 ImageSize = Tuple{Integer,Integer} | |
| 25 Image = Array{Float64,2} | |
| 26 Primal = Image | |
| 27 Dual = Array{Float64,3} | |
| 28 | |
| 29 ######################### | |
| 30 # Iterate initialisation | |
| 31 ######################### | |
| 32 | |
| 33 function init_rest(x::Primal, b::Array{Float64, N}) where N | |
| 34 imdim=size(x) | |
| 35 datadim=size(b) | |
| 36 | |
| 37 y = zeros(2, imdim...) | |
| 38 λ = zeros(datadim...) | |
| 39 Δx₁ = copy(x) | |
| 40 Δx₂ = copy(x) | |
| 41 Δy = copy(y) | |
| 42 Δλ = copy(λ) | |
| 43 x̄ = copy(x) | |
| 44 | |
| 45 return x, y, λ, Δx₁, Δx₂, Δy, Δλ, x̄ | |
| 46 end | |
| 47 | |
| 48 function init_primal(xinit::Image) | |
| 49 return copy(xinit) | |
| 50 end | |
| 51 | |
| 52 ############ | |
| 53 # Algorithm | |
| 54 ############ | |
| 55 | |
| 56 function recon_pdps(b :: Data, op :: LinOp{Image, Data}; | |
| 57 xinit :: Union{Image, Nothing} = nothing, | |
| 58 iterate = Iterate.simple_iterate, | |
| 59 params::NamedTuple) where Data | |
| 60 | |
| 61 ################################ | |
| 62 # Extract and set up parameters | |
| 63 ################################ | |
| 64 | |
| 65 α, ρ = params.α, params.ρ | |
| 66 τ₀, σ₀ = params.τ₀, params.σ₀ | |
| 67 | |
| 68 R_K = √(∇₂_norm₂₂_est^2+opnorm_estimate(op)^2) | |
| 69 γ = 1 | |
| 70 | |
| 71 @assert(τ₀*σ₀ < 1) | |
| 72 σ = σ₀/R_K | |
| 73 τ = τ₀/R_K | |
| 74 | |
| 75 ###################### | |
| 76 # Initialise iterates | |
| 77 ###################### | |
| 78 | |
| 79 x, y, λ, Δx₁, Δx₂, Δy, Δλ, x̄ = init_rest(copy(xinit), b) | |
| 80 | |
| 81 #################### | |
| 82 # Run the algorithm | |
| 83 #################### | |
| 84 | |
| 85 v = iterate(params) do verbose :: Function | |
| 86 ω = params.accel ? 1/√(1+2*γ*τ) : 1 | |
| 87 | |
| 88 ∇₂ᵀ!(Δx₁, y) # primal step: x | |
| 89 inplace!(Δx₂, op', λ) # | | |
| 90 @. x̄ = x # | | |
| 91 @. x = x-τ*(Δx₁+Δx₂) # | | |
| 92 @. x̄ = (1+ω)*x - ω*x̄ # over-relax: x̄ = 2x-x_old | |
| 93 ∇₂!(Δy, x̄) # dual step: y, λ | |
| 94 inplace!(Δλ, op, x̄) # | | |
| 95 @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | | |
| 96 proj_norm₂₁ball!(y, α) # | | |
| 97 @. λ = (λ + σ*(Δλ-b))/(1 + σ) # | | |
| 98 | |
| 99 if params.accel | |
| 100 τ, σ = τ*ω, σ/ω | |
| 101 end | |
| 102 | |
| 103 ################################ | |
| 104 # Give function value if needed | |
| 105 ################################ | |
| 106 v = verbose() do | |
| 107 ∇₂!(Δy, x) | |
| 108 inplace!(Δλ, op, x) | |
| 109 value = norm₂²(Δλ-b)/2 + params.α*γnorm₂₁(Δy, params.ρ) | |
| 110 value, x | |
| 111 end | |
| 112 | |
| 113 v | |
| 114 end | |
| 115 | |
| 116 return x, y, v | |
| 117 end | |
| 118 | |
| 119 end # Module | |
| 120 | |
| 121 |
