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