# HG changeset patch # User Tuomo Valkonen # Date 1574379867 18000 # Node ID 1cffd3d07fe2c97fb630ed6414d6b3656742e295 # Parent c5aabb2c41d99d0a6995a5126db61df8f09281c3 Denoising routine for testing + visualisation tools diff -r c5aabb2c41d9 -r 1cffd3d07fe2 Project.toml --- a/Project.toml Tue Nov 19 10:15:38 2019 -0500 +++ b/Project.toml Thu Nov 21 18:44:27 2019 -0500 @@ -2,3 +2,13 @@ uuid = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9" authors = ["Tuomo Valkonen "] version = "0.1.0" + +[deps] +AlgTools = "c46e2e78-5339-41fd-a966-983ff60ab8e7" +Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" 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 + + diff -r c5aabb2c41d9 -r 1cffd3d07fe2 src/ImageTools.jl --- a/src/ImageTools.jl Tue Nov 19 10:15:38 2019 -0500 +++ b/src/ImageTools.jl Thu Nov 21 18:44:27 2019 -0500 @@ -6,6 +6,8 @@ include("Gradient.jl") include("Translate.jl") +include("Denoise.jl") +include("Visualise.jl") end # module diff -r c5aabb2c41d9 -r 1cffd3d07fe2 src/Visualise.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Visualise.jl Thu Nov 21 18:44:27 2019 -0500 @@ -0,0 +1,192 @@ +################################################## +# Visualising and data-collecting iteration tools +################################################## + +module Visualise + +using Printf +using Distributed +using FileIO +using Setfield +using Images, Plots, Measures + +using AlgTools.Util +using AlgTools.StructTools +using AlgTools.LinkedLists + +############## +# Our exports +############## + +export LogEntry, + bg_visualise, + visualise, + clip, + grayimg, + secs_ns, + iterate_visualise, + initialise_visualisation, + finalise_visualisation + +################## +# Data structures +################## + +struct LogEntry <: IterableStruct + iter :: Int + time :: Float64 + function_value :: Float64 +end + +struct State + vis :: Union{Distributed.RemoteChannel,Bool,Nothing} + visproc :: Union{Nothing,Future} + start_time :: Union{Real,Nothing} + wasted_time :: Real + log :: LinkedList{LogEntry} +end + +################## +# Helper routines +################## + +@inline function secs_ns() + return convert(Float64, time_ns())*1e-9 +end + +clip = x -> min(max(x, 0.0), 1.0) +grayimg = im -> Gray.(clip.(im)) + +################ +# Visualisation +################ + +function bg_visualise(rc) + while true + imgs=take!(rc) + # Take only the latest image to visualise + while isready(rc) + imgs=take!(rc) + end + # We're done if we were fed an empty image + if isnothing(imgs) + break + end + do_visualise(imgs) + end + return +end + +function do_visualise(imgs) + plt = im -> plot(grayimg(im), showaxis=false, grid=false, aspect_ratio=:equal, margin=2mm) + display(plot([plt(imgs[i]) for i =1:length(imgs)]..., reuse=true, margin=0mm)) +end + +function visualise(rc_or_vis, imgs) + if isa(rc_or_vis, RemoteChannel) + rc = rc_or_vis + while isready(rc) + take!(rc) + end + put!(rc, imgs) + elseif isa(rc_or_vis, Bool) && rc_or_vis + do_visualise(imgs) + end +end + +###################################################### +# Iterator that does visualisation and log collection +###################################################### + +function iterate_visualise(st :: State, + step :: Function, + params :: NamedTuple) where DisplacementT + try + for iter=1:params.maxiter + st = step() do calc_objective + if isnothing(st.start_time) + # The Julia precompiler is a miserable joke, apparently not crossing module + # boundaries, so only start timing after the first iteration. + st = @set st.start_time=secs_ns() + end + + verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0 + + if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0) + verb_start = secs_ns() + tm = verb_start - st.start_time - st.wasted_time + value, x = calc_objective() + + entry = LogEntry(iter, tm, value) + + # (**) Collect a singly-linked list of log to avoid array resizing + # while iterating + st = @set st.log=LinkedListEntry(entry, st.log) + + if verb + @printf("%d/%d J=%f\n", iter, params.maxiter, value) + visualise(st.vis, (x,)) + end + + if params.save_iterations + fn = t -> "$(params.save_prefix)_$(t)_iter$(iter).png" + save(File(format"PNG", fn("reco")), grayimg(x)) + end + + st = @set st.wasted_time += (secs_ns() - verb_start) + end + + return st + end + end + catch ex + if isa(ex, InterruptException) + # If SIGINT is received (user pressed ^C), terminate computations, + # returning current status. Effectively, we do not call `step()` again, + # ending the iterations, but letting the algorithm finish up. + # Assuming (**) above occurs atomically, `st.log` should be valid, but + # any results returned by the algorithm itself may be partial, as for + # reasons of efficiency we do *not* store results of an iteration until + # the next iteration is finished. + printstyled("\rUser interrupt—finishing up.\n", bold=true, color=202) + else + throw(ex) + end + end + + return st +end + +#################### +# Launcher routines +#################### + +function initialise_visualisation(visualise; iterator=iterate_visualise) + # Create visualisation + if visualise + rc = RemoteChannel() + visproc = @spawn bg_visualise(rc) + vis =rc + #vis = true + + sleep(0.01) + else + vis = false + visproc = nothing + end + + st = State(vis, visproc, nothing, 0.0, nothing) + iterate = curry(iterate_visualise, st) + + return st, iterate +end + +function finalise_visualisation(st) + if isa(st.rc, RemoteChannel) + # Tell subprocess to finish, and wait + put!(st.rc, nothing) + wait(st.visproc) + end +end + +end # Module \ No newline at end of file diff -r c5aabb2c41d9 -r 1cffd3d07fe2 test/denoise.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/denoise.jl Thu Nov 21 18:44:27 2019 -0500 @@ -0,0 +1,64 @@ +################## +# Denoise testing +################## + +__precompile__() + +using Printf +using Images +import TestImages + +using AlgTools.Util +using AlgTools.LinkedLists +using ImageTools.Denoise +using ImageTools.Visualise + +const default_save_prefix="denoise_result_" + +const default_params = ( + α = 1, + τ₀ = 5, + σ₀ = 0.99/5, + ρ = 0, + accel = true, + noise_level = 0.5, + verbose_iter = 10, + maxiter = 1000, + save_results = false, + save_iterations = false, + image_name = "lighthouse", +) + +####################### +# Main testing routine +####################### + +function test_denoise(; + visualise=true, + save_prefix=default_save_prefix, + kwargs...) + + + # Parameters for this experiment + params = default_params ⬿ kwargs + params = params ⬿ (save_prefix = save_prefix * params.image_name,) + + # Load image and add noise + b = Float64.(Gray.(TestImages.testimage(params.image_name))) + b_noisy = b .+ params.noise_level.*randn(size(b)...) + + # Launch (background) visualiser + st, iterate = initialise_visualisation(visualise) + + # Run algorithm + x, y, st = denoise_pdps(b_noisy; iterate=iterate, params=params) + + if params.save_results + perffile = params.save_prefix * ".txt" + println("Saving " * perffile) + write_log(perffile, st.log, "# params = $(params)\n") + end + + # Exit background visualiser + finish_visualisation(st) +end