Thu, 25 Apr 2024 13:05:40 -0500
Reduce code duplication.
--- a/src/Algorithm.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,148 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module Algorithm - -identifier = "pdps_known" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function step_lengths(params, γ, R_K²) - ρ̃₀, τ₀, σ₀, σ̃₀ = params.ρ̃₀, params.τ₀, params.σ₀, params.σ̃₀ - δ = params.δ - ρ = isdefined(params, :phantom_ρ) ? params.phantom_ρ : params.ρ - Λ = params.Λ - Θ = params.dual_flow ? Λ : 1 - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*min(1/(τ*R_K²), 1/max(0, τ*R_K²/((1+γ*τ-Λ)*(1-δ))-ρ)) - q = δ*(1+σ*ρ)/Θ - if 1 ≥ q - σ̃ = σ̃₀*σ/q - #ρ̃ = ρ̃₀*max(0, ((Θ*σ)/(2*δ*σ̃^2*(1+σ*ρ))+1/(2σ)-1/σ̃)) - ρ̃ = max(0, (1-q)/(2*σ)) - else - σ̃ = σ̃₀*σ/(q*(1-√(1-1/q))) - ρ̃ = 0 - end - - println("Step length parameters: τ=$(τ), σ=$(σ), σ̃=$(σ̃), ρ̃=$(ρ̃)") - - return τ, σ, σ̃, ρ̃ -end - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, params.dual_flow) - - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothMulti.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/AlgorithmBothMulti.jl Thu Apr 25 13:05:40 2024 -0500 @@ -26,7 +26,7 @@ ConstantDisplacementHornSchunckData, filter_hs -using ..Algorithm: step_lengths +using ..AlgorithmProximal: step_lengths ######################### # Iterate initialisation
--- a/src/AlgorithmBothMulti.jl.orig Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,277 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothMulti - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.ImFilter - -using ..OpticalFlow: Image, - ImageSize, - DisplacementConstant, - pdflow!, - horn_schunck_reg_prox!, - pointwise_gradiprod_2d!, - horn_schunck_reg_prox_op!, - mldivide_step_plus_sym2x2!, - ConstantDisplacementHornSchunckData, - filter_hs - -using ..Algorithm: step_lengths - -######################### -# Iterate initialisation -######################### - -function init_displ(xinit::Image, ::Type{DisplacementConstant}, n::Integer) - return xinit, zeros(n, 2) -end - -# function init_displ(xinit::Image, ::Type{DisplacementFull}, n::Integer) -# return xinit, zeros(n, 2, size(xinit)...) -# end - -function init_rest(x::Image, u::DisplacementT) where DisplacementT - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄, u -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Image, - n::Integer) where DisplacementT - return init_rest(init_displ(copy(xinit), DisplacementT, n)...) -end - -function init_iterates( :: Type{DisplacementT}, - dim::ImageSize, - n::Integer) where DisplacementT - return init_rest(init_displ(zeros(dim...), DisplacementT, n)...) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - n = params.displacement_count - k = 0 # number of displacements we have already - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim, n) - init_data = (params.init == :data) - hs = [ConstantDisplacementHornSchunckData() for i=1:n] - #hs = Array{ConstantDisplacementHornSchunckData}(undef, n) - A = Array{Float64,3}(undef, n, 2, 2) - d = Array{Float64,2}(undef, n, 2) - - # … for tracking cumulative movement - ucumulbase = [0.0, 0.0] - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ = params.α, params.ρ, params.λ, params.θ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - T = params.timestep - - #################### - # Run the algorithm - #################### - - b_next_filt = nothing - diffu = similar(u[1, :]) - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - b_next :: Image - - #################################### - # Smooth data for Horn–Schunck term - #################################### - - b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) - - ################################################ - # Prediction step - ################################################ - - # Predict x and y - if k==0 - if init_data - x .= b - init_data = false - end - else - # Displacement from previous to this image is estimated as - # the difference of their displacements from beginning of window. - if k>1 - @. @views diffu = u[k, :] - u[k-1, :] - else - @. @views diffu = u[k, :] - end - - pdflow!(x, Δx, y, Δy, diffu, params.dual_flow) - end - - # Shift stored prox matrices - if k==n - tmp = copy(u[1, :]) - ucumulbase .+= tmp - for j=1:(n-1) - @. @views u[j, :] = u[j+1, :] - tmp - hs[j] = hs[j+1] - end - # Create new struct as original contains references to objects that - # have been moved to index n-1. - hs[n]=ConstantDisplacementHornSchunckData() - else - k += 1 - end - - # Predict u: zero displacement from current to next image, i.e., - # same displacement to beginning of window. - if k==1 - @. @views u[k, :] = 0.0 - else - @. @views u[k, :] = u[k-1, :] - end - - # Predictor proximals tep - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - ################################################################################# - # PDPS step - # - # For the displacements, with τ̃=τ/k, we need to solve for 2≤j<k, - # - # (1) (I/τ̃+M₀^j+M₀^{j+1})u^j = M₀^ju^{j-1} + M₀^{j+1}u^{j+1} - # + ũ^j/τ̃ - z^j + z^{j+1}, - # - # as well as - # - # (2) (I/τ̃+M₀^k)u^k = M₀^k u^{k-1} + ũ^k/τ̃ - z^k - # - # and - # - # (3) (I/τ̃+M₀^1+M₀^2)u^1 = 0 + M₀^{2}u^{2} + ũ^1/τ̃ - z^1 + z^{2} - # - # We first construct from (2) that - # - # u^k = A^k u^{k-1} + d^k - # - # for - # - # A^k := (I/τ̃+M₀^k)^{-1} M₀^k - # d_k := (I/τ̃+M₀^k)^{-1} (ũ^k/τ̃ - z^k). - # - # Inserting this into (1) we need - # - # (4) (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))u^j = M₀^ju^{j-1} + M₀^{j+1}d^{j+1} - # + ũ^j/τ̃ - z^j + z^{j+1}. - # - # This is well-defined because A^{j+1} < I. It also has the same form as (1), so - # we continue with - # - # (5) u^j = A^j u^{j-1} + d^j - # - # for - # - # A^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} M₀^j - # d^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} - # (M₀^{j+1}d^{j+1} + ũ^j/τ̃ - z^j + z^{j+1}) - # - # Finally from (3) with these we need - # - # (I/τ̃+M₀^1+M₀^2(I-A^2))u^1 = M₀^2d^2 + ũ^1/τ̃ - z^1 + z^2, - # - # which is of the same form as (4) with u^0=0, so by (5) u^1=d^1. - # - ################################################################################# - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - # | | for displacement - # Calculate matrices for latest data; rest is stored. - @views begin - horn_schunck_reg_prox_op!(hs[k], b_next_filt, b_filt, θ, λ, T) - - τ̃=τ/k - - B = hs[k].M₀ - c = u[k, :]./τ̃-hs[k].z - mldivide_step_plus_sym2x2!(A[k, :, :], B, B, τ̃) - mldivide_step_plus_sym2x2!(d[k, :], B, c, τ̃) - - for j=(k-1):-1:1 - B = hs[j].M₀+hs[j+1].M₀*([1 0; 0 1]-A[j+1, :, :]) - c = hs[j+1].M₀*d[j+1, :]+u[j, :]./τ̃-hs[j].z+hs[j+1].z - mldivide_step_plus_sym2x2!(A[j, :, :], B, hs[j].M₀, τ̃) - mldivide_step_plus_sym2x2!(d[j, :], B, c, τ̃) - end - - u[1, :] .= d[1, :] - for j=2:k - u[j, :] .= A[j, :, :]*u[j-1, :] + d[j, :] - end - end - - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - ∇₂!(Δy, x) - hs_plus_reg=0 - for j=1:k - v=(j==1 ? u[j, :] : u[j, :]-u[j-1, :]) - hs_plus_reg += hs[j].cv/2 + dot(hs[j].Mv*v, v)/2+dot(hs[j].av, v) - end - value = (norm₂²(b-x)/2 + hs_plus_reg/k + α*γnorm₂₁(Δy, ρ)) - - value, x, u[k, :]+ucumulbase, u[1:k,:].+ucumulbase' - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmDualScaling.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,126 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmDualScaling - -identifier = "pdps_known_dualscaling" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow!, - DualScaling - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - ds = DualScaling(params.ds_exponent, params.ds_threshold) - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, ds) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmGreedy.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmGreedy - -identifier = "pdps_known_greedy" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow!, - Greedy - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, Greedy()) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmNew.jl Thu Apr 25 13:05:40 2024 -0500 @@ -0,0 +1,127 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmNew + +# Default identifier for when none is given by predictor +identifier = "pdps_known_noprediction" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient + +using ..OpticalFlow: ImageSize, + Image, + pdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + 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) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1.0 + Λ = params.Λ + τ₀, σ₀ = params.τ₀, params.σ₀ + + τ = τ₀/γ + @assert(1+γ*τ ≥ Λ) + σ = σ₀*1/(τ*R_K²) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄ = init_iterates(dim) + init_data = (params.init == :data) + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, + v_known :: DisplacementT, + 🚫unused_b_next :: Image + + ################## + # Prediction step + ################## + if init_data + x .= b + init_data = false + end + + if haskey(params, :predictor) && ~isnothing(params.predictor) + pdflow!(x, Δx, y, Δy, v_known, params.predictor) + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + @. x = (x-τ*(Δx-b))/(1+τ) # | prox + @. x̄ = 2x - x̄ # over-relax + ∇₂!(Δy, x̄) # dual step: y + @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | + proj_norm₂₁ball!(y, α) # | prox + + ################################ + # Give function value if needed + ################################ + v = verbose() do + ∇₂!(Δy, x) + value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- a/src/AlgorithmNoPrediction.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,124 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmNoPrediction - -identifier = "pdps_known_noprediction" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - # No prediction - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmPrimalOnly.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,124 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmPrimalOnly - -identifier = "pdps_known_primalonly" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, false) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmRotation.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmRotation - -identifier = "pdps_known_rotation" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow!, - Rotation - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, Rotation()) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmZeroDual.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmZeroDual - -identifier = "pdps_known_zerodual" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - 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) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, false) - y .= zeros(size(y)...) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/ImGenerate.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/ImGenerate.jl Thu Apr 25 13:05:40 2024 -0500 @@ -100,10 +100,10 @@ θ = 2π*rand(rng,Float64) r = params.shake return [r*cos(θ), r*sin(θ)] - end + end else error("Unknown shaketype $(params.shaketype)") - end + end end pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) @@ -181,15 +181,14 @@ :: Type{DisplacementConstant}, datachannel :: Channel{OnlineData{DisplacementConstant}}, params :: NamedTuple) - # Set up counter and zero factor for stabilisation @assert(params.maxiter ≥ maximum(params.stable_interval)) indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 + zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - # Restart the seed to enable comparison across predictors - Random.seed!(rng,951508) # algreadygood + # Restart the seed to enable comparison across predictors + Random.seed!(rng, if haskey(params, :seed) params.seed else 951508 end) nextv = shake(params) v_true = zero_factor.*nextv() @@ -228,17 +227,17 @@ # PETscan ######################################################################## function generate_sinogram(im, sz, - :: Type{DisplacementConstant}, - datachannel :: Channel{PetOnlineData{DisplacementConstant}}, - params :: NamedTuple) + :: Type{DisplacementConstant}, + datachannel :: Channel{PetOnlineData{DisplacementConstant}}, + params :: NamedTuple) # Set up counter and zero factor for stabilisation @assert(params.maxiter ≥ maximum(params.stable_interval)) indx = 1 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - # Restart the seed to enable comparison across predictors - Random.seed!(rng,314159) + # Restart the seed to enable comparison across predictors + Random.seed!(rng, if haskey(params, :seed) params.seed else 314159 end) nextv = shake(params) v_true = zero_factor.*nextv() @@ -282,7 +281,7 @@ # Next theta theta_true = zero_factor*rotatebytheta(params) theta_cumul += theta_true - # Next sinogram mask + # Next sinogram mask S_true = generate_radonmask(params) # Update indx and zero factor indx += 1
--- a/src/OpticalFlow.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/OpticalFlow.jl Thu Apr 25 13:05:40 2024 -0500 @@ -42,7 +42,8 @@ HornSchunckData, filter_hs, petpdflow!, - DualScaling, Greedy, Rotation + DualScaling, Greedy, Rotation, ZeroDual, PrimalOnly, + identifier ############################################### # Types (several imported from ImageTools.Translate) @@ -68,6 +69,28 @@ struct Greedy end struct Rotation end +struct ZeroDual end +struct PrimalOnly end + +function identifier(::DualScaling) + "pdps_known_dualscaling" +end + +function identifier(::Rotation) + "pdps_known_rotation" +end + +function identifier(::Greedy) + "pdps_known_greedy" +end + +function identifier(::ZeroDual) + "pdps_known_zerodual" +end + +function identifier(::PrimalOnly) + "pdps_known_primalonly" +end ################################# # Displacement field based flow @@ -102,7 +125,7 @@ end end -function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none) +function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow :: Bool; threads=:none) if dual_flow @backgroundif (threads==:outer) begin flow!(x, u, Δx; threads=(threads==:inner)) @@ -119,7 +142,7 @@ # Additional method for Greedy function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none) - @assert(size(u)==(2,)) + @assert(size(u)==(2,)) Δx .= x Δy .= y flow!(x, u; threads=(threads==:inner)) @@ -130,11 +153,11 @@ inds = abs.(Dxx) .≤ 1e-1 Dxx[inds] .= 1 DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx + y .= y.* DΔx ./ Dxx end # Additional method for Rotation -function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) @assert(size(u)==(2,)) Δx .= x flow!(x, u; threads=(threads==:inner)) @@ -173,13 +196,24 @@ c = 1 .* (1 .- cc./ cm) .^flow.exponent C[1,:,:] .= c C[2,:,:] .= c - y .= C.*y + y .= C.*y end +function pdflow!(x, Δx, y, Δy, u, :: ZeroDual; threads=:none) + @assert(size(u)==(2,)) + flow!(x, u; threads=(threads==:inner)) + y .= 0.0 +end + +function pdflow!(x, Δx, y, Δy, u, :: PrimalOnly; threads=:none) + @assert(size(u)==(2,)) + flow!(x, u; threads=(threads==:inner)) +end ########################## # PET ########################## + function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant; threads = false) tmp .= x @@ -191,7 +225,7 @@ petflow! = petflow_interp! -function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none) +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow :: Bool; threads=:none) if dual_flow @backgroundif (threads==:outer) begin petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) @@ -219,7 +253,7 @@ inds = abs.(Dxx) .≤ 1e-2 Dxx[inds] .= 1 DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx + y .= y.* DΔx ./ Dxx end # Method for dual scaling predictor @@ -248,6 +282,15 @@ end end +function petpdflow!(x, Δx, y, Δy, u, theta_known, :: ZeroDual; threads=:none) + petflow!(x, Δx, u, theta_known, threads=(threads==:inner)) + y .= 0.0 +end + +function petpdflow!(x, Δx, y, Δy, u, theta_known, :: PrimalOnly; threads=:none) + petflow!(x, Δx, u, theta_known, threads=(threads==:inner)) +end + ########################## # Linearised optical flow ########################## @@ -358,7 +401,7 @@ a, c, d = 0.0, 0.0, 0.0 # This used to use ∇₂cfold but it is faster to allocate temporary # storage for the full gradient due to probably better memory and SIMD - # instruction usage. + # instruction usage. g = zeros(2, size(b)...) ∇₂c!(g, b) @inbounds for i=1:size(b, 1) @@ -425,7 +468,7 @@ f = x -> simple_imfilter(x, kernel; threads=true) end - # We already filtered b in the previous step (b_next in that step) + # We already filtered b in the previous step (b_next in that step) b_filt = b_next_filt==nothing ? f(b) : b_next_filt b_next_filt = f(b_next)
--- a/src/OpticalFlow.jl.orig Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,280 +0,0 @@ -################################ -# Code relevant to optical flow -################################ - -__precompile__() - -module OpticalFlow - -using AlgTools.Util -using ImageTools.Gradient -import ImageTools.Translate -using ImageTools.ImFilter - -########## -# Exports -########## - -export flow!, - pdflow!, - flow_grad!, - flow_interp!, - estimate_Λ², - estimate_linear_Λ², - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - horn_schunck_reg_prox!, - horn_schunck_reg_prox_op!, - mldivide_step_plus_sym2x2!, - linearised_optical_flow_error, - Image, AbstractImage, ImageSize, - Gradient, Displacement, - DisplacementFull, DisplacementConstant, - HornSchunckData, - filter_hs - -############################################### -# Types (several imported from ImageTools.Translate) -############################################### - -Image = Translate.Image -AbstractImage = AbstractArray{Float64,2} -Displacement = Translate.Displacement -DisplacementFull = Translate.DisplacementFull -DisplacementConstant = Translate.DisplacementConstant -Gradient = Array{Float64,3} -ImageSize = Tuple{Int64,Int64} - -################################# -# Displacement field based flow -################################# - -function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage; - threads = false) - tmp .= x - Translate.translate_image!(x, tmp, u; threads=threads) -end - -function flow_interp!(x::AbstractImage, u::Displacement; - threads = false) - tmp = copy(x) - Translate.translate_image!(x, tmp, u; threads=threads) -end - -flow! = flow_interp! - -function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none) - if dual_flow - #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu, - # (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :]))) - @backgroundif (threads==:outer) begin - flow!(x, u, Δx; threads=(threads==:inner)) - end begin - flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) - flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) - end - else - flow!(x, u, Δx) - end -end - -function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none) - if dual_flow - @backgroundif (threads==:outer) begin - flow!(x, u, Δx; threads=(threads==:inner)) - flow!(z, u, Δz; threads=(threads==:inner)) - end begin - flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) - flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) - end - else - flow!(x, u, Δx; threads=(threads==:inner)) - flow!(z, u, Δz; threads=(threads==:inner)) - end -end - -########################## -# Linearised optical flow -########################## - -# ⟨⟨u, ∇b⟩⟩ -function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, - u::DisplacementFull, b::Image; - add = false) - ∇₂c!(vtmp, b) - - u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) - vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) - y′=reshape(y, prod(size(y))) - - if add - @simd for i = 1:length(y′) - @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i])) - end - else - @simd for i = 1:length(y′) - @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i])) - end - end -end - -function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, - u::DisplacementConstant, b::Image; - add = false) - ∇₂c!(vtmp, b) - - vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) - y′=reshape(y, prod(size(y))) - - if add - @simd for i = 1:length(y′) - @inbounds y′[i] += dot(u, @view(vtmp′[:, i])) - end - else - @simd for i = 1:length(y′) - @inbounds y′[i] = dot(u, @view(vtmp′[:, i])) - end - end -end - -# ∇b ⋅ y -function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image) - ∇₂c!(u, b) - - u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) - y′=reshape(y, prod(size(y))) - - @simd for i=1:length(y′) - @inbounds @. u′[:, i] *= y′[i] - end -end - -function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image) - @assert(size(y)==size(b) && size(u)==(2,)) - u .= 0 - ∇₂cfold!(b, nothing) do g, st, (i, j) - @inbounds u .+= g.*y[i, j] - return st - end - # Reweight to be with respect to 𝟙^*𝟙 inner product. - u ./= prod(size(b)) -end - -mutable struct ConstantDisplacementHornSchunckData - M₀::Array{Float64,2} - z::Array{Float64,1} - Mv::Array{Float64,2} - av::Array{Float64,1} - cv::Float64 - - function ConstantDisplacementHornSchunckData() - return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0) - end -end - -# For DisplacementConstant, for the simple prox step -# -# (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2, -# -# construct matrix M₀ and vector z such that we can solve u from -# -# (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z -# -# Note that the problem -# -# argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2 -# + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2 -# -# has with respect to u the system -# -# (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ, -# -# where the primed variables correspond to (2) for (1) for uʹ in place of u: -# -# argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2 -# -function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData, - bnext::Image, b::Image, θ, λ, T) - @assert(size(b)==size(bnext)) - w = prod(size(b)) - z = hs.z - cv = 0 - # Factors of symmetric matrix [a c; c d] - a, c, d = 0.0, 0.0, 0.0 - # This used to use ∇₂cfold but it is faster to allocate temporary - # storage for the full gradient due to probably better memory and SIMD - # instruction usage. - g = zeros(2, size(b)...) - ∇₂c!(g, b) - @inbounds for i=1:size(b, 1) - for j=1:size(b, 2) - δ = bnext[i,j]-b[i,j] - @. z += g[:,i,j]*δ - cv += δ*δ - a += g[1,i,j]*g[1,i,j] - c += g[1,i,j]*g[2,i,j] - d += g[2,i,j]*g[2,i,j] - end - end - w₀ = λ - w₂ = θ/w - aʹ = w₀ + w₂*a - cʹ = w₂*c - dʹ = w₀ + w₂*d - hs.M₀ .= [aʹ cʹ; cʹ dʹ] - hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d] - hs.cv = cv*θ - hs.av .= hs.z.*θ - hs.z .*= w₂/T -end - -# Solve the 2D system (I/τ+M₀)u = z -@inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ) - a = 1/τ+M₀[1, 1] - c = M₀[1, 2] - d = 1/τ+M₀[2, 2] - u .= ([d -c; -c a]*z)./(a*d-c*c) -end - -function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image, - θ, λ, T, τ) - hs=ConstantDisplacementHornSchunckData() - horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T) - mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ) -end - -function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing) - if !isnothing(δ) - u = δ.*u - end - pointwise_gradiprod_2d!(x, vtmp, u, x; add=true) -end - -# Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation -function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image) - imdim = size(b) - vtmp = zeros(2, imdim...) - tmp = b-b_prev - pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true) - return tmp -end - -############################################## -# Helper to smooth data for Horn–Schunck term -############################################## - -function filter_hs(b, b_next, b_next_filt, kernel) - if kernel==nothing - f = x -> x - else - f = x -> simple_imfilter(x, kernel; threads=false) - end - - # We already filtered b in the previous step (b_next in that step) - b_filt = b_next_filt==nothing ? f(b) : b_next_filt - b_next_filt = f(b_next) - - return b_filt, b_next_filt -end - -end # Module
--- a/src/PET/AlgorithmDualScaling.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmDualScaling - -identifier = "pdps_known_dualscaling" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow!, - DualScaling - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - petpdflow!(x, Δx, y, Δy, v_known, theta_known, DualScaling()) - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/PET/AlgorithmGreedy.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmGreedy - -identifier = "pdps_known_greedy" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow!, - Greedy - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - petpdflow!(x, Δx, y, Δy, v_known, theta_known, Greedy()) - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmNew.jl Thu Apr 25 13:05:40 2024 -0500 @@ -0,0 +1,183 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmNew + +# Default identifier for when none is given by predictor +identifier = "pdps_known_noprediction" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ...Radon +using ...OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + if haskey(params, :predictor) && ~isnothing(params.predictor) + petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.predictor) + end + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- a/src/PET/AlgorithmNoPrediction.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmNoPrediction - -identifier = "pdps_known_noprediction" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - # No prediction - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/PET/AlgorithmPrimalOnly.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmPrimalOnly - -identifier = "pdps_known_primalonly" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - petpdflow!(x, Δx, y, Δy, v_known, theta_known, false) - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/PET/AlgorithmProximal.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/PET/AlgorithmProximal.jl Thu Apr 25 13:05:40 2024 -0500 @@ -15,14 +15,14 @@ using ImageTools.Gradient using ImageTools.Translate -using ..Radon using ImageTransformations using Images, CoordinateTransformations, Rotations, OffsetArrays using ImageCore, Interpolations -using ..OpticalFlow: ImageSize, - Image, - petpdflow! +using ...Radon +using ...OpticalFlow: ImageSize, + Image, + petpdflow! ######################### # Iterate initialisation @@ -103,9 +103,9 @@ iterate = AlgTools.simple_iterate, params::NamedTuple) where DisplacementT - ################################ + ################################ # Extract and set up parameters - ################################ + ################################ α, ρ = params.α, params.ρ R_K² = ∇₂_norm₂₂_est² γ = 1 @@ -139,9 +139,9 @@ v_known :: DisplacementT, theta_known :: DisplacementT, b_true :: Image, - S :: Image + S :: Image - ################### + ################### # Prediction steps ################### @@ -151,14 +151,14 @@ @. oldx = x end - ############################## + ############################## # Proximal step of prediction ############################## ∇₂!(Δy, x) @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - #@. cc = y + 1000000*σ̃*Δy + #@. cc = y + 1000000*σ̃*Δy #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc - proj_norm₂₁ball!(y, α) + proj_norm₂₁ball!(y, α) ############ # PDPS step @@ -177,7 +177,7 @@ ##################### # L update if needed - ##################### + ##################### if params.L_experiment petgrad!(petgradx, x, b, c, S) petgrad!(oldpetgradx, oldx, b, c, S) @@ -186,18 +186,18 @@ println("Step length parameters: L=$(L)") τ = τ₀/L σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end + end + end ################################ # Give function value if needed ################################ - v = verbose() do + v = verbose() do ∇₂!(Δy, x) value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) value, x, [NaN, NaN], nothing - end + end v end
--- a/src/PET/AlgorithmRotation.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmRotation - -identifier = "pdps_known_rotation" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow!, - Rotation - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - petpdflow!(x, Δx, y, Δy, v_known, theta_known, Rotation()) - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/PET/AlgorithmZeroDual.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmZeroDual - -identifier = "pdps_known_zerodual" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.Translate - -using ..Radon -using ImageTransformations -using Images, CoordinateTransformations, Rotations, OffsetArrays -using ImageCore, Interpolations - -using ..OpticalFlow: ImageSize, - Image, - petpdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - radonx = copy(x) - return x, y, Δx, Δy, x̄, radonx -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -end - -######################### -# PETscan related -######################### -function petvalue(x, b, c) - tmp = similar(b) - radon!(tmp, x) - return sum(@. tmp - b*log(tmp+c)) -end - -function petgrad!(res, x, b, c, S) - tmp = similar(b) - radon!(tmp, x) - @. tmp = S .- b/(tmp+c) - backproject!(res, S.*tmp) -end - -function proj_nonneg!(y) - @inbounds @simd for i=1:length(y) - if y[i] < 0 - y[i] = 0 - end - end - return y -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ################################ - # Extract and set up parameters - ################################ - α, ρ = params.α, params.ρ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - L = params.L - τ₀, σ₀ = params.τ₀, params.σ₀ - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - λ = params.λ - c = params.c*ones(params.radondims...) - - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) - - if params.L_experiment - oldpetgradx = zeros(size(x)...) - petgradx = zeros(size(x)) - oldx = ones(size(x)) - end - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, # noisy_sinogram - v_known :: DisplacementT, - theta_known :: DisplacementT, - b_true :: Image, - S :: Image - - ################### - # Prediction steps - ################### - - petpdflow!(x, Δx, y, Δy, v_known, theta_known, false) - y .= zeros(size(y)...) - - if params.L_experiment - @. oldx = x - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term - - @. x = x-(τ*λ)*r∇-τ*Δx # | - proj_nonneg!(x) # | non-negativity constaint prox - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: - @. y = y + σ*Δy # | - proj_norm₂₁ball!(y, α) # | prox - - ##################### - # L update if needed - ##################### - if params.L_experiment - petgrad!(petgradx, x, b, c, S) - petgrad!(oldpetgradx, oldx, b, c, S) - if norm₂(x-oldx)>1e-12 - L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end - - ################################ - # Give function value if needed - ################################ - - v = verbose() do - ∇₂!(Δy, x) - value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) - value, x, [NaN, NaN], nothing, τ, σ - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/PET/ImGenerate.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,308 +0,0 @@ -################### -# Image generation -################### - -module ImGenerate - -using ColorTypes: Gray -import TestImages -# We don't really *directly* depend on QuartzImageIO. The import here is -# merely a workaround to suppress warnings when loading TestImages. -# Something is broken in those packages. -import QuartzImageIO - -using AlgTools.Util -using AlgTools.Comms -using ImageTools.Translate - -using ..OpticalFlow: Image, DisplacementConstant, DisplacementFull -using ..Radon - -# Added for reproducibility -import StableRNGs: StableRNG, Random -const rng = StableRNG(314159) - -# Added for PET -import PoissonRandom: pois_rand -import Random: shuffle -import Images: center, warp, imresize -import CoordinateTransformations: recenter -import Rotations: RotMatrix -import Interpolations: Flat -import MAT: matread - - -############## -# Our exports -############## - -export ImGen, - OnlineData, - imgen_square, - imgen_shake, - PetOnlineData, - imgen_shepplogan_radon, - imgen_brainphantom_radon - -################## -# Data structures -################## - -struct ImGen - f :: Function - dim :: Tuple{Int64,Int64} - Λ :: Float64 - dynrange :: Float64 - name :: String -end - -struct OnlineData{DisplacementT} - b_true :: Image - b_noisy :: Image - v :: DisplacementT - v_true :: DisplacementT - v_cumul_true :: DisplacementT -end - -struct PetOnlineData{DisplacementT} - b_true :: Image - sinogram_true :: Image - sinogram_noisy :: Image - v :: DisplacementT - v_true :: DisplacementT - v_cumul_true :: DisplacementT - theta :: DisplacementT # theta = thetaknown, theta_cumul - S :: Image -end - -################### -# Shake generation -################### - -function make_const_v(displ, sz) - v = zeros(2, sz...) - v[1, :, :] .= displ[1] - v[2, :, :] .= displ[2] - return v -end - -function shake(params) - if !haskey(params, :shaketype) || params.shaketype == :gaussian - return () -> params.shake.*randn(rng,2) - elseif params.shaketype == :disk - return () -> begin - θ = 2π*rand(rng,Float64) - r = params.shake*√(rand(rng,Float64)) - return [r*cos(θ), r*sin(θ)] - end - elseif params.shaketype == :circle - return () -> begin - θ = 2π*rand(rng,Float64) - r = params.shake - return [r*cos(θ), r*sin(θ)] - end - else - error("Unknown shaketype $(params.shaketype)") - end -end - -pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) - - -function rotatebytheta(params) - r = params.rotation_factor*randn(rng) - return r -end - -function generate_radonmask(params) - imdim = params.radondims - sino_sparse = params.sino_sparsity - numzero = Int64(round(sino_sparse*imdim[1]*imdim[2])) - numone = imdim[1]*imdim[2]-numzero - A = shuffle(rng,reshape([ones(numone); zeros(numzero)],(imdim[1],imdim[2]))) - return A -end - -################ -# Moving square -################ - -function generate_square(sz, - :: Type{DisplacementT}, - datachannel :: Channel{OnlineData{DisplacementT}}, - params) where DisplacementT - - if false - v₀ = make_const_v(0.1.*(-1, 1), sz) - nextv = () -> v₀ - elseif DisplacementT == DisplacementFull - nextv = pixelwise(shake(params), sz) - elseif DisplacementT == DisplacementConstant - nextv = shake(params) - else - @error "Invalid DisplacementT" - end - - # Constant linear displacement everywhere has Jacobian determinant one - # (modulo the boundaries which we ignore here) - m = round(Int, sz[1]/5) - b_orig = zeros(sz...) - b_orig[sz[1].-(2*m:3*m), 2*m:3*m] .= 1 - - v_true = nextv() - v_cumul = copy(v_true) - - while true - # Flow original data and add noise - b_true = zeros(sz...) - translate_image!(b_true, b_orig, v_cumul; threads=true) - b = b_true .+ params.noise_level.*randn(rng,sz...) - v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) - # Pass true data to iteration routine - data = OnlineData{DisplacementT}(b_true, b, v, v_true, v_cumul) - if !put_unless_closed!(datachannel, data) - return - end - # Next step shake - v_true = nextv() - v_cumul .+= v_true - end -end - -function imgen_square(sz) - return ImGen(curry(generate_square, sz), sz, 1, 1, "square$(sz[1])x$(sz[2])") -end - -################ -# Shake a photo -################ - -function generate_shake_image(im, sz, - :: Type{DisplacementConstant}, - datachannel :: Channel{OnlineData{DisplacementConstant}}, - params :: NamedTuple) - - - # Set up counter and zero factor for stabilisation - @assert(params.maxiter ≥ maximum(params.stable_interval)) - indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - - # Restart the seed to enable comparison across predictors - Random.seed!(rng,314159) - - - nextv = shake(params) - v_true = zero_factor.*nextv() - v_cumul = copy(v_true) - - while true - # Extract subwindow of original image and add noise - b_true = zeros(sz...) - extract_subimage!(b_true, im, v_cumul; threads=true) - b = b_true .+ params.noise_level.*randn(rng,sz...) - v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) - # Pass data to iteration routine - data = OnlineData{DisplacementConstant}(b_true, b, v, v_true, v_cumul) - if !put_unless_closed!(datachannel, data) - return - end - # Next step shake - v_true = zero_factor.*nextv() - v_cumul .+= v_true - # Update indx and zero factor - indx += 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - end -end - -function imgen_shake(imname, sz) - im = Float64.(Gray.(TestImages.testimage(imname))) - dynrange = maximum(im) - return ImGen(curry(generate_shake_image, im, sz), sz, 1, dynrange, - "$(imname)$(sz[1])x$(sz[2])") -end - - - -######################################################################## -# PETscan -######################################################################## -function generate_sinogram(im, sz, - :: Type{DisplacementConstant}, - datachannel :: Channel{PetOnlineData{DisplacementConstant}}, - params :: NamedTuple) - - # Set up counter and zero factor for stabilisation - @assert(params.maxiter ≥ maximum(params.stable_interval)) - indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - - # Restart the seed to enable comparison across predictors - Random.seed!(rng,314159) - - nextv = shake(params) - v_true = zero_factor.*nextv() - v_cumul = copy(v_true) - - S_true = generate_radonmask(params) - - theta_true = zero_factor*rotatebytheta(params) - theta_cumul = copy(theta_true) - - while true - # Define the transformation matrix - center_point = (sz[1]/2 + v_true[1], sz[2]/2 + v_true[2]) - tform = recenter(RotMatrix(theta_cumul), center_point) - - # Apply the transformation to the image using warp - b_true = copy(warp(im, tform, axes(im), fillvalue=Flat())) - - v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) - theta = theta_true*(1.0 + params.rotation_noise_level.*randn(rng)) - - # Generate the true and noisy sinograms - sinogram_true = zeros(params.radondims...) - sinogram_true .*= params.scale - radon!(sinogram_true, b_true) - sinogram_noisy = copy(sinogram_true) - - for i=1:params.radondims[1], j=1:params.radondims[2] - sinogram_noisy[i, j] += pois_rand(rng,params.noise_level) - end - - # Pass data to iteration routine - data = PetOnlineData{DisplacementConstant}(b_true, sinogram_true, sinogram_noisy, v, v_true, v_cumul, [theta, theta_cumul], S_true) - if !put_unless_closed!(datachannel, data) - return - end - - # Next step shake - v_true = zero_factor.*nextv() - v_cumul .+= v_true - # Next theta - theta_true = zero_factor*rotatebytheta(params) - theta_cumul += theta_true - # Next sinogram mask - S_true = generate_radonmask(params) - # Update indx and zero factor - indx += 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - end -end - -function imgen_shepplogan_radon(sz) - im = convert(Array{Float64},TestImages.shepp_logan(sz[1], highContrast=true)) - dynrange = maximum(im) - return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])") -end - -function imgen_brainphantom_radon(sz) - data = matread("src/PET/phantom_slice.mat") - im = normalise(imresize(convert(Array{Float64},data["square_data"]),sz)) - dynrange = maximum(im) - return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "brainphantom$(sz[1])x$(sz[2])") -end - -normalise = (data) -> data./maximum(data) -end # Module
--- a/src/PET/OpticalFlow.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,430 +0,0 @@ -################################ -# Code relevant to optical flow -################################ - -__precompile__() - -module OpticalFlow - -using AlgTools.Util -using ImageTools.Gradient -import ImageTools.Translate -using ImageTools.ImFilter - -# using ImageTransformations -# using Images, CoordinateTransformations, Rotations, OffsetArrays -# using Interpolations - -import Images: center, warp -import CoordinateTransformations: recenter -import Rotations: RotMatrix -import Interpolations: Flat - -########## -# Exports -########## - -export flow!, - pdflow!, - flow_grad!, - flow_interp!, - estimate_Λ², - estimate_linear_Λ², - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - horn_schunck_reg_prox!, - horn_schunck_reg_prox_op!, - mldivide_step_plus_sym2x2!, - linearised_optical_flow_error, - Image, AbstractImage, ImageSize, - Gradient, Displacement, - DisplacementFull, DisplacementConstant, - HornSchunckData, - filter_hs, - petpdflow!, - DualScaling, Greedy, Rotation - -############################################### -# Types (several imported from ImageTools.Translate) -############################################### - -Image = Translate.Image -AbstractImage = AbstractArray{Float64,2} -Displacement = Translate.Displacement -DisplacementFull = Translate.DisplacementFull -DisplacementConstant = Translate.DisplacementConstant -Gradient = Array{Float64,3} -ImageSize = Tuple{Int64,Int64} - - -################################# -# Struct for flow -################################# -struct DualScaling end -struct Greedy end -struct Rotation end - -################################# -# Displacement field based flow -################################# - -function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage; - threads = false) - tmp .= x - Translate.translate_image!(x, tmp, u; threads=threads) -end - -function flow_interp!(x::AbstractImage, u::Displacement; - threads = false) - tmp = copy(x) - Translate.translate_image!(x, tmp, u; threads=threads) -end - -flow! = flow_interp! - -function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none) - if dual_flow - #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu, - # (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :]))) - @backgroundif (threads==:outer) begin - flow!(x, u, Δx; threads=(threads==:inner)) - end begin - flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) - flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) - end - else - flow!(x, u, Δx) - end -end - -function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none) - if dual_flow - @backgroundif (threads==:outer) begin - flow!(x, u, Δx; threads=(threads==:inner)) - flow!(z, u, Δz; threads=(threads==:inner)) - end begin - flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) - flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) - end - else - flow!(x, u, Δx; threads=(threads==:inner)) - flow!(z, u, Δz; threads=(threads==:inner)) - end -end - -# Additional method for Greedy -function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none) - @assert(size(u)==(2,)) - Δx .= x - Δy .= y - flow!(x, u; threads=(threads==:inner)) - Dxx = similar(Δy) - DΔx = similar(Δy) - ∇₂!(Dxx, x) - ∇₂!(DΔx, Δx) - inds = abs.(Dxx) .≤ 1e-1 - Dxx[inds] .= 1 - DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx -end - -# Additional method for Rotation -function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) - @assert(size(u)==(2,)) - Δx .= x - flow!(x, u; threads=(threads==:inner)) - (m,n) = size(x) - dx = similar(y) - dx_banana = similar(y) - ∇₂!(dx, Δx) - ∇₂!(dx_banana, x) - for i=1:m - for j=1:n - ndx = @views sum(dx[:, i, j].^2) - ndx_banana = @views sum(dx_banana[:, i, j].^2) - if ndx > 1e-4 && ndx_banana > 1e-4 - A = dx[:, i, j] - B = dx_banana[:, i, j] - theta = atan(B[1] * A[2] - B[2] * A[1], B[1] * A[1] + B[2] * A[2]) # Oriented angle from A to B - cos_theta = cos(theta) - sin_theta = sin(theta) - a = cos_theta * y[1, i, j] - sin_theta * y[2, i, j] - b = sin_theta * y[1, i, j] + cos_theta * y[2, i, j] - y[1, i, j] = a - y[2, i, j] = b - end - end - end -end - -# Additional method for Dual Scaling -function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none) - @assert(size(u)==(2,)) - oldx = copy(x) - flow!(x, u; threads=(threads==:inner)) - C = similar(y) - cc = abs.(x-oldx) - cm = max(1e-12,maximum(cc)) - c = 1 .* (1 .- cc./ cm) .^(10) - C[1,:,:] .= c - C[2,:,:] .= c - y .= C.*y -end - - -########################## -# PET -########################## -function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant; - threads = false) - tmp .= x - center_point = center(x) .+ u - tform = recenter(RotMatrix(theta_known[1]), center_point) - tmp = warp(x, tform, axes(x), fillvalue=Flat()) - x .= tmp -end - -petflow! = petflow_interp! - -function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none) - if dual_flow - @backgroundif (threads==:outer) begin - petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) - end begin - petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, theta_known; threads=(threads==:inner)) - petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, theta_known; threads=(threads==:inner)) - end - else - petflow!(x, Δx, u, theta_known) - end -end - -# Method for greedy predictor -function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Greedy; threads=:none) - oldx = copy(x) - center_point = center(x) .+ u - tform = recenter(RotMatrix(theta_known[1]), center_point) - Δx = warp(x, tform, axes(x), fillvalue=Flat()) - @. x = Δx - @. Δy = y - Dxx = copy(Δy) - DΔx = copy(Δy) - ∇₂!(Dxx, x) - ∇₂!(DΔx, oldx) - inds = abs.(Dxx) .≤ 1e-2 - Dxx[inds] .= 1 - DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx -end - -# Method for dual scaling predictor -function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: DualScaling; threads=:none) - oldx = copy(x) - center_point = center(x) .+ u - tform = recenter(RotMatrix(theta_known[1]), center_point) - Δx = warp(x, tform, axes(x), fillvalue=Flat()) - @. x = Δx - C = similar(y) - cc = abs.(x-oldx) - cm = max(1e-12,maximum(cc)) - c = 1 .* (1 .- cc./ cm) .^(10) - C[1,:,:] .= c - C[2,:,:] .= c - y .= C.*y -end - -# Method for rotation prediction (exploiting property of inverse rotation) -function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Rotation; threads=:none) - @backgroundif (threads==:outer) begin - petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) - end begin - petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner)) - petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner)) - end -end - -########################## -# Linearised optical flow -########################## - -# ⟨⟨u, ∇b⟩⟩ -function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, - u::DisplacementFull, b::Image; - add = false) - ∇₂c!(vtmp, b) - - u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) - vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) - y′=reshape(y, prod(size(y))) - - if add - @simd for i = 1:length(y′) - @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i])) - end - else - @simd for i = 1:length(y′) - @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i])) - end - end -end - -function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, - u::DisplacementConstant, b::Image; - add = false) - ∇₂c!(vtmp, b) - - vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) - y′=reshape(y, prod(size(y))) - - if add - @simd for i = 1:length(y′) - @inbounds y′[i] += dot(u, @view(vtmp′[:, i])) - end - else - @simd for i = 1:length(y′) - @inbounds y′[i] = dot(u, @view(vtmp′[:, i])) - end - end -end - -# ∇b ⋅ y -function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image) - ∇₂c!(u, b) - - u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) - y′=reshape(y, prod(size(y))) - - @simd for i=1:length(y′) - @inbounds @. u′[:, i] *= y′[i] - end -end - -function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image) - @assert(size(y)==size(b) && size(u)==(2,)) - u .= 0 - ∇₂cfold!(b, nothing) do g, st, (i, j) - @inbounds u .+= g.*y[i, j] - return st - end - # Reweight to be with respect to 𝟙^*𝟙 inner product. - u ./= prod(size(b)) -end - -mutable struct ConstantDisplacementHornSchunckData - M₀::Array{Float64,2} - z::Array{Float64,1} - Mv::Array{Float64,2} - av::Array{Float64,1} - cv::Float64 - - function ConstantDisplacementHornSchunckData() - return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0) - end -end - -# For DisplacementConstant, for the simple prox step -# -# (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2, -# -# construct matrix M₀ and vector z such that we can solve u from -# -# (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z -# -# Note that the problem -# -# argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2 -# + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2 -# -# has with respect to u the system -# -# (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ, -# -# where the primed variables correspond to (2) for (1) for uʹ in place of u: -# -# argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2 -# -function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData, - bnext::Image, b::Image, θ, λ, T) - @assert(size(b)==size(bnext)) - w = prod(size(b)) - z = hs.z - cv = 0 - # Factors of symmetric matrix [a c; c d] - a, c, d = 0.0, 0.0, 0.0 - # This used to use ∇₂cfold but it is faster to allocate temporary - # storage for the full gradient due to probably better memory and SIMD - # instruction usage. - g = zeros(2, size(b)...) - ∇₂c!(g, b) - @inbounds for i=1:size(b, 1) - for j=1:size(b, 2) - δ = bnext[i,j]-b[i,j] - @. z += g[:,i,j]*δ - cv += δ*δ - a += g[1,i,j]*g[1,i,j] - c += g[1,i,j]*g[2,i,j] - d += g[2,i,j]*g[2,i,j] - end - end - w₀ = λ - w₂ = θ/w - aʹ = w₀ + w₂*a - cʹ = w₂*c - dʹ = w₀ + w₂*d - hs.M₀ .= [aʹ cʹ; cʹ dʹ] - hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d] - hs.cv = cv*θ - hs.av .= hs.z.*θ - hs.z .*= w₂/T -end - -# Solve the 2D system (I/τ+M₀)u = z -@inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ) - a = 1/τ+M₀[1, 1] - c = M₀[1, 2] - d = 1/τ+M₀[2, 2] - u .= ([d -c; -c a]*z)./(a*d-c*c) -end - -function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image, - θ, λ, T, τ) - hs=ConstantDisplacementHornSchunckData() - horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T) - mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ) -end - -function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing) - if !isnothing(δ) - u = δ.*u - end - pointwise_gradiprod_2d!(x, vtmp, u, x; add=true) -end - -# Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation -function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image) - imdim = size(b) - vtmp = zeros(2, imdim...) - tmp = b-b_prev - pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true) - return tmp -end - -############################################## -# Helper to smooth data for Horn–Schunck term -############################################## - -function filter_hs(b, b_next, b_next_filt, kernel) - if kernel==nothing - f = x -> x - else - f = x -> simple_imfilter(x, kernel; threads=true) - end - - # We already filtered b in the previous step (b_next in that step) - b_filt = b_next_filt==nothing ? f(b) : b_next_filt - b_next_filt = f(b_next) - - return b_filt, b_next_filt -end - -end # Module
--- a/src/PET/PET.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/PET/PET.jl Thu Apr 25 13:05:40 2024 -0500 @@ -18,7 +18,7 @@ using AlgTools.StructTools using AlgTools.LinkedLists using AlgTools.Comms -using ImageTools.Visualise: secs_ns, grayimg, do_visualise +using ImageTools.Visualise: secs_ns, grayimg, do_visualise using ImageTools.ImFilter: gaussian # For PET @@ -26,30 +26,18 @@ ##################### # Load local modules -##################### -include("OpticalFlow.jl") -include("Radon.jl") -include("ImGenerate.jl") -include("AlgorithmDualScaling.jl") -include("AlgorithmGreedy.jl") -include("AlgorithmNoPrediction.jl") -include("AlgorithmPrimalOnly.jl") +#####################a +include("AlgorithmNew.jl") include("AlgorithmProximal.jl") -include("AlgorithmRotation.jl") -include("AlgorithmZeroDual.jl") #include("PlotResults.jl") -import .AlgorithmDualScaling -import .AlgorithmGreedy -import .AlgorithmNoPrediction -import .AlgorithmPrimalOnly +import .AlgorithmNew import .AlgorithmProximal -import .AlgorithmRotation -import .AlgorithmZeroDual -using .Radon: backproject! -using .ImGenerate -using .OpticalFlow: DisplacementFull, DisplacementConstant +using ..Radon: backproject! +using ..ImGenerate +using ..OpticalFlow +using ..Run #using .PlotResults @@ -57,9 +45,9 @@ # Our exports ############## -export demo_petS1, demo_petS2, demo_petS3, +export demo_petS1, demo_petS2, demo_petS3, demo_petS4, demo_petS5, demo_petS6, demo_petS7, - demo_petB1, demo_petB2, demo_petB3, + demo_petB1, demo_petB2, demo_petB3, demo_petB4, demo_petB5, demo_petB6, demo_petB7, batchrun_shepplogan, batchrun_brainphantom, batchrun_pet #plot_pet @@ -68,24 +56,6 @@ # Parameterisation and experiments ################################### -struct Experiment - mod :: Module - DisplacementT :: Type - imgen :: ImGen - params :: NamedTuple -end - -function Base.show(io::IO, e::Experiment) - displacementname(::Type{DisplacementFull}) = "DisplacementFull" - displacementname(::Type{DisplacementConstant}) = "DisplacementConstant" - print(io, " - mod: $(e.mod) - DisplacementT: $(displacementname(e.DisplacementT)) - imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2]) - params: $(e.params ⬿ (kernel = "(not shown)",)) - ") -end - const default_save_prefix="img/" const default_params = ( @@ -107,66 +77,68 @@ stable_interval = Set(0), ) -const p_known₀_pet = ( +const p_known₀_pet = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.1, - shake = 1.0, - rotation_factor = 0.075, - rotation_noise_level = 0.0075, + shake = 1.0, + rotation_factor = 0.075, + rotation_noise_level = 0.0075, α = 0.15, ρ̃₀ = 1.0, σ̃₀ = 1.0, δ = 0.9, σ₀ = 1.0, τ₀ = 0.9, - λ = 1, + λ = 1, radondims = [128,64], sz = (256,256), scale = 1, c = 1.0, - sino_sparsity = 0.5, + sino_sparsity = 0.5, L = 300.0, L_experiment = false, #stable_interval = Set(0), stable_interval = union(Set(1000:2000),Set(3500:4000)), ) +const p_known₀_petb = p_known₀_pet# ⬿ ( seed = 313159, ) + const shepplogan = imgen_shepplogan_radon(p_known₀_pet.sz) const brainphantom = imgen_brainphantom_radon(p_known₀_pet.sz) const shepplogan_experiments_pdps_known = ( - Experiment(AlgorithmDualScaling, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmGreedy, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmNoPrediction, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, shepplogan, - p_known₀_pet), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, shepplogan, p_known₀_pet ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmZeroDual, DisplacementConstant, shepplogan, - p_known₀_pet), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pet ⬿ (predictor=ZeroDual(),)), ) const brainphantom_experiments_pdps_known = ( - Experiment(AlgorithmDualScaling, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmGreedy, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmNoPrediction, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, brainphantom, - p_known₀_pet), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, brainphantom, - p_known₀_pet ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmZeroDual, DisplacementConstant, brainphantom, - p_known₀_pet), + p_known₀_petb ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=ZeroDual(),)), ) @@ -178,167 +150,15 @@ brainphantom_experiments_pdps_known, )) -################ -# Log -################ - -struct LogEntry <: IterableStruct - iter :: Int - time :: Float64 - function_value :: Float64 - #v_cumul_true_y :: Float64 - #v_cumul_true_x :: Float64 - #v_cumul_est_y :: Float64 - #v_cumul_est_x :: Float64 - psnr :: Float64 - ssim :: Float64 - #psnr_data :: Float64 - #ssim_data :: Float64 -end - -struct LogEntryHiFi <: IterableStruct - iter :: Int - v_cumul_true_y :: Float64 - v_cumul_true_x :: Float64 -end - -############### -# Main routine -############### - -struct State - vis :: Union{Channel,Bool,Nothing} - start_time :: Union{Real,Nothing} - wasted_time :: Real - log :: LinkedList{LogEntry} - log_hifi :: LinkedList{LogEntryHiFi} - aborted :: Bool -end - -function name(e::Experiment, p) - ig = e.imgen - # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))" - return "$(ig.name)_$(e.mod.identifier)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))" -end - -function write_tex(texfile, e_params) - open(texfile, "w") do io - wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}") - wf = (n, s) -> if isdefined(e_params, s) - wp(n, getfield(e_params, s)) - end - wf("alpha", :α) - wf("sigmazero", :σ₀) - wf("tauzero", :τ₀) - wf("tildetauzero", :τ̃₀) - wf("delta", :δ) - wf("lambda", :λ) - wf("theta", :θ) - wf("maxiter", :maxiter) - wf("noiselevel", :noise_level) - wf("shakenoiselevel", :shake_noise_level) - wf("shake", :shake) - wf("timestep", :timestep) - wf("displacementcount", :displacementcount) - wf("phantomrho", :phantom_ρ) - if isdefined(e_params, :σ₀) - wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}") - end - end -end - -function run_experiments(;visualise=true, - recalculate=true, - experiments, - save_prefix=default_save_prefix, - fullscreen=false, - kwargs...) - - # Create visualisation - if visualise - rc = Channel(1) - visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen) - bind(rc, visproc) - vis = rc - else - vis = false - end - - # Run all experiments - for e ∈ experiments - - # Parameters for this experiment - e_params = default_params ⬿ e.params ⬿ kwargs - ename = name(e, e_params) - e_params = e_params ⬿ (save_prefix = save_prefix * ename, - dynrange = e.imgen.dynrange, - Λ = e.imgen.Λ) - - if recalculate || !isfile(e_params.save_prefix * ".txt") - println("Running experiment \"$(ename)\"") - - # Start data generation task - datachannel = Channel{PetOnlineData{e.DisplacementT}}(2) - gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params) - bind(datachannel, gentask) - - # Run algorithm - iterate = curry(iterate_visualise, datachannel, - State(vis, nothing, 0.0, nothing, nothing, false)) - - x, y, st = e.mod.solve(e.DisplacementT; - dim=e.imgen.dim, - iterate=iterate, - params=e_params) - - # Clear non-saveable things - st = @set st.vis = nothing - - println("Wasted_time: $(st.wasted_time)s") - - if e_params.save_results - println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)") - - perffile = e_params.save_prefix * ".txt" - hififile = e_params.save_prefix * "_hifi.txt" - texfile = e_params.save_prefix * "_params.tex" - # datafile = e_params.save_prefix * ".jld2" - - write_log(perffile, st.log, "# params = $(e_params)\n") - #write_log(hififile, st.log_hifi, "# params = $(e_params)\n") - #write_tex(texfile, e_params) - # @save datafile x y st params - end - - close(datachannel) - wait(gentask) - - if st.aborted - break - end - else - println("Skipping already computed experiment \"$(ename)\"") - # texfile = e_params.save_prefix * "_params.tex" - # write_tex(texfile, e_params) - end - end - - if visualise - # Tell subprocess to finish, and wait - put!(rc, nothing) - close(rc) - wait(visproc) - end - - return -end - ####################### # Demos and batch runs ####################### function demo(experiment; kwargs...) run_experiments(;experiments=(experiment,), + save_prefix=default_save_prefix, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, save_results=false, save_images=false, visualise=true, @@ -367,6 +187,9 @@ function batchrun_shepplogan(;kwargs...) run_experiments(;experiments=shepplogan_experiments_all, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -376,6 +199,9 @@ function batchrun_brainphantom(;kwargs...) run_experiments(;experiments=brainphantom_experiments_all, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -395,21 +221,21 @@ function rescale(arr, new_range) old_min = minimum(arr) old_max = maximum(arr) - scale_factor = (new_range[2] - new_range[1]) / (old_max - old_min) + scale_factor = (new_range[2] - new_range[1]) / (old_max - old_min) scaled_arr = new_range[1] .+ (arr .- old_min) * scale_factor return scaled_arr end -function iterate_visualise(datachannel::Channel{PetOnlineData{DisplacementT}}, - st :: State, - step :: Function, - params :: NamedTuple) where DisplacementT +function iterate_visualise_pet(datachannel::Channel{PetOnlineData{DisplacementT}}, + st :: State, + step :: Function, + params :: NamedTuple) where DisplacementT try sc = nothing d = take!(datachannel) - for iter=1:params.maxiter + for iter=1:params.maxiter dnext = take!(datachannel) st = step(d.sinogram_noisy, d.v, d.theta, d.b_true, d.S) do calc_objective stn = st @@ -435,7 +261,7 @@ entry = LogEntry(iter, tm, value, #sc*d.v_cumul_true[1], #sc*d.v_cumul_true[2], - #sc*v[1], sc*v[2], + #sc*v[1], sc*v[2], assess_psnr(x, d.b_true), assess_ssim(x, d.b_true), #assess_psnr(d.b_noisy, d.b_true), @@ -509,39 +335,6 @@ return st end -function bg_visualise_enhanced(rc; fullscreen=false) - process_channel(rc) do d - imgs, plot_movement, log, vhist = d - do_visualise(imgs, refresh=false, fullscreen=fullscreen) - # Overlay movement - GR.settextcolorind(5) - GR.setcharheight(0.015) - GR.settextpath(GR.TEXT_PATH_RIGHT) - tx, ty = GR.wctondc(0, 1) - GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr) - if plot_movement - sc=1.0 - p=unfold_linked_list(log) - x=map(e -> 1.5+sc*e.v_cumul_true_x, p) - y=map(e -> 0.5+sc*e.v_cumul_true_y, p) - GR.setlinewidth(2) - GR.setlinecolorind(2) - GR.polyline(x, y) - x=map(e -> 1.5+sc*e.v_cumul_est_x, p) - y=map(e -> 0.5+sc*e.v_cumul_est_y, p) - GR.setlinecolorind(3) - GR.polyline(x, y) - if vhist != nothing - GR.setlinecolorind(4) - x=map(v -> 1.5+sc*v, vhist[:,2]) - y=map(v -> 0.5+sc*v, vhist[:,1]) - GR.polyline(x, y) - end - end - GR.updatews() - end -end - # Clip image values to allowed range clip = x -> min(max(x, 0.0), 1.0) @@ -566,13 +359,4 @@ # fv_plot("brainphantom") #end - -############### -# Precompiling -############### - -# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}}) -# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool}) -# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}}) - end # Module
--- a/src/PET/Radon.jl Thu Apr 25 11:14:41 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,85 +0,0 @@ -# From https://github.com/LMescheder/Radon.jl/blob/master/src/Radon.jl -# Updated to modern Julia. Fixed radon! to zero out `radonim`. - -__precompile__() - - -module Radon -export radon!, backproject! -# package code goes here - -function radon(im::AbstractArray{T, 2}) where T<:AbstractFloat - Nx, Ny = size(im) - radonim = similar(im, min(Nx, Ny), max(Nx, Ny)) - radon!(radonim, im) - radonim -end - -function radon!(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2}) where T<:AbstractFloat - Nφ, Ns = size(radonim) - Nx, Ny = size(im) - Ldiag = hypot(Nx, Ny) - K = round(Int, Ldiag) - - dLz = Ldiag/(K-1) - dLs = Ldiag/(Ns-1) - dφ = π/Nφ - - @inbounds for is = 1:Ns, iφ = 1:Nφ - φ = (iφ-1)*dφ - s = -Ldiag/2 + (is-1)*dLs - tx = cos(φ) - ty = sin(φ) - - radonim[iφ, is] = zero(T) - for k=1:K - z = -Ldiag/2 + (k-1)*dLz - x = round(Int, Nx/2 + z*ty + s*tx) - y = round(Int, Ny/2 - z*tx + s*ty) - - if 1 <= x <= Nx && 1 <= y <= Ny - radonim[iφ, is] += im[x, y]*dLz - end - end - end - radonim -end - -function backproject(radonim::AbstractArray{T, 2}) where T<:AbstractFloat - Nφ, Ns = size(radonim) - im = similar(radonim, Ns, Ns) - backproject!(im, radonim) - im -end - -function backproject!(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2}) where T<:AbstractFloat - Nφ, Ns = size(radonim) - Nx, Ny = size(im) - Ldiag = hypot(Nx, Ny) - K = round(Int, Ldiag) - - dLz = Ldiag/(K-1) - dLs = Ldiag/(Ns-1) - dφ = π/Nφ - - im .= 0. - @inbounds for is = 1:Ns, iφ = 1:Nφ - φ = (iφ-1)*dφ - s = -Ldiag/2 + (is-1)*dLs - tx = cos(φ) - ty = sin(φ) - - for k=1:K - z = -Ldiag/2 + (k-1)*dLz - x = round(Int, Nx/2 + z*ty + s*tx) - y = round(Int, Ny/2 - z*tx + s*ty) - - if 1 <= x <= Nx && 1 <= y <= Ny - im[x, y] += radonim[iφ, is]*dLz - end - end - end - im -end - -end # module
--- a/src/PredictPDPS.jl Thu Apr 25 11:14:41 2024 -0500 +++ b/src/PredictPDPS.jl Thu Apr 25 13:05:40 2024 -0500 @@ -10,20 +10,8 @@ # Load external modules ######################## -using Printf -using FileIO -#using JLD2 -using Setfield -using ImageQualityIndexes: assess_psnr, assess_ssim -using DelimitedFiles -import GR - +using ImageTools.ImFilter: gaussian using AlgTools.Util -using AlgTools.StructTools -using AlgTools.LinkedLists -using AlgTools.Comms -using ImageTools.Visualise: secs_ns, grayimg, do_visualise -using ImageTools.ImFilter: gaussian ##################### # Load local modules @@ -32,50 +20,29 @@ include("OpticalFlow.jl") include("Radon.jl") include("ImGenerate.jl") -include("Algorithm.jl") -include("AlgorithmBoth.jl") -include("AlgorithmBothGreedyV.jl") -include("AlgorithmBothCumul.jl") +include("Run.jl") +include("AlgorithmProximal.jl") include("AlgorithmBothMulti.jl") -include("AlgorithmBothNL.jl") include("AlgorithmFB.jl") include("AlgorithmFBDual.jl") +include("AlgorithmNew.jl") include("Stats.jl") #include("PlotResults.jl") - - -# Additional -include("AlgorithmProximal.jl") -include("AlgorithmGreedy.jl") -include("AlgorithmRotation.jl") -include("AlgorithmNoPrediction.jl") -include("AlgorithmPrimalOnly.jl") -include("AlgorithmDualScaling.jl") -include("AlgorithmZeroDual.jl") include("PET/PET.jl") -import .Algorithm, - .AlgorithmBoth, - .AlgorithmBothGreedyV, - .AlgorithmBothCumul, - .AlgorithmBothMulti, - .AlgorithmBothNL, +import .AlgorithmBothMulti, .AlgorithmFB, .AlgorithmFBDual, .AlgorithmProximal, - .AlgorithmGreedy, - .AlgorithmRotation, - .AlgorithmNoPrediction, - .AlgorithmPrimalOnly, - .AlgorithmDualScaling, - .AlgorithmZeroDual + .AlgorithmNew using .ImGenerate -using .OpticalFlow: DisplacementFull, DisplacementConstant +using .OpticalFlow using .Stats #using .PlotResults using .PET +using .Run ############## # Our exports @@ -101,24 +68,6 @@ # Parameterisation and experiments ################################### -struct Experiment - mod :: Module - DisplacementT :: Type - imgen :: ImGen - params :: NamedTuple -end - -function Base.show(io::IO, e::Experiment) - displacementname(::Type{DisplacementFull}) = "DisplacementFull" - displacementname(::Type{DisplacementConstant}) = "DisplacementConstant" - print(io, " - mod: $(e.mod) - DisplacementT: $(displacementname(e.DisplacementT)) - imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2]) - params: $(e.params ⬿ (kernel = "(not shown)",)) - ") -end - const default_save_prefix="img/" const default_params = ( @@ -133,20 +82,17 @@ 1000, 2000, 2500, 3000, 4000, 5000, 6000, 7000, 7500, 8000, 9000, 10000, 8700]), pixelwise_displacement=false, - dual_flow = true, - prox_predict = true, + dual_flow = true, # For AlgorithmProximalfrom 2019 paper handle_interrupt = true, init = :zero, plot_movement = false, stable_interval = Set(0), - ds_exponent = 10, - ds_threshold = 10e-12 ) const square = imgen_square((200, 300)) const lighthouse = imgen_shake("lighthouse", (200, 300)) -const p_known₀_denoising = ( +const p_known₀_denoising = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.025, shake = 2.0, @@ -160,7 +106,7 @@ stable_interval = union(Set(2500:5000),Set(8700:10000)), ) -const p_known₀ = ( +const p_known₀ = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.05, shake = 2, @@ -172,7 +118,7 @@ τ₀ = 0.01, ) -const p_unknown₀ = ( +const p_unknown₀ = default_params ⬿ ( noise_level = 0.3, shake_noise_level = 0.05, shake = 2, @@ -190,12 +136,14 @@ ) +# Experiments for 2019 paper + const experiments_pdps_known = ( - Experiment(Algorithm, DisplacementConstant, lighthouse, + Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀_denoising ⬿ (phantom_ρ = 0,)), - Experiment(Algorithm, DisplacementConstant, lighthouse, + Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀ ⬿ (phantom_ρ = 100,)), - Experiment(Algorithm, DisplacementConstant, square, + Experiment(AlgorithmProximal, DisplacementConstant, square, p_known₀ ⬿ (phantom_ρ = 0,)) ) @@ -219,182 +167,29 @@ experiments_fb_known )) +# Image stabilisation experiments for 2024 paper. PET experiments are in PET/PET.jl + const denoising_experiments_pdps_known = ( - Experiment(AlgorithmDualScaling, DisplacementConstant, lighthouse, - p_known₀_denoising), - Experiment(AlgorithmGreedy, DisplacementConstant, lighthouse, - p_known₀_denoising), - Experiment(AlgorithmNoPrediction, DisplacementConstant, lighthouse, - p_known₀_denoising), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, lighthouse, - p_known₀_denoising), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀_denoising ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, lighthouse, - p_known₀_denoising), - Experiment(AlgorithmZeroDual, DisplacementConstant, lighthouse, - p_known₀_denoising), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=ZeroDual(),)), ) const denoising_experiments_all = Iterators.flatten(( denoising_experiments_pdps_known, )) -################ -# Log -################ - -struct LogEntry <: IterableStruct - iter :: Int - time :: Float64 - function_value :: Float64 - #v_cumul_true_y :: Float64 - #v_cumul_true_x :: Float64 - #v_cumul_est_y :: Float64 - #v_cumul_est_x :: Float64 - psnr :: Float64 - ssim :: Float64 - #psnr_data :: Float64 - #ssim_data :: Float64 -end - -struct LogEntryHiFi <: IterableStruct - iter :: Int - v_cumul_true_y :: Float64 - v_cumul_true_x :: Float64 -end - -############### -# Main routine -############### - -struct State - vis :: Union{Channel,Bool,Nothing} - start_time :: Union{Real,Nothing} - wasted_time :: Real - log :: LinkedList{LogEntry} - log_hifi :: LinkedList{LogEntryHiFi} - aborted :: Bool -end - -function name(e::Experiment, p) - ig = e.imgen - # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))" - return "$(ig.name)_$(e.mod.identifier)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))" -end - -function write_tex(texfile, e_params) - open(texfile, "w") do io - wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}") - wf = (n, s) -> if isdefined(e_params, s) - wp(n, getfield(e_params, s)) - end - wf("alpha", :α) - wf("sigmazero", :σ₀) - wf("tauzero", :τ₀) - wf("tildetauzero", :τ̃₀) - wf("delta", :δ) - wf("lambda", :λ) - wf("theta", :θ) - wf("maxiter", :maxiter) - wf("noiselevel", :noise_level) - wf("shakenoiselevel", :shake_noise_level) - wf("shake", :shake) - wf("timestep", :timestep) - wf("displacementcount", :displacementcount) - wf("phantomrho", :phantom_ρ) - if isdefined(e_params, :σ₀) - wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}") - end - end -end - -function run_experiments(;visualise=true, - recalculate=true, - experiments, - save_prefix=default_save_prefix, - fullscreen=false, - kwargs...) - - # Create visualisation - if visualise - rc = Channel(1) - visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen) - bind(rc, visproc) - vis = rc - else - vis = false - end - - # Run all experiments - for e ∈ experiments - - # Parameters for this experiment - e_params = default_params ⬿ e.params ⬿ kwargs - ename = name(e, e_params) - e_params = e_params ⬿ (save_prefix = save_prefix * ename, - dynrange = e.imgen.dynrange, - Λ = e.imgen.Λ) - - if recalculate || !isfile(e_params.save_prefix * ".txt") - println("Running experiment \"$(ename)\"") - - # Start data generation task - datachannel = Channel{OnlineData{e.DisplacementT}}(2) - gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params) - bind(datachannel, gentask) - - # Run algorithm - iterate = curry(iterate_visualise, datachannel, - State(vis, nothing, 0.0, nothing, nothing, false)) - - x, y, st = e.mod.solve(e.DisplacementT; - dim=e.imgen.dim, - iterate=iterate, - params=e_params) - - # Clear non-saveable things - st = @set st.vis = nothing - - println("Wasted_time: $(st.wasted_time)s") - - if e_params.save_results - println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)") - - perffile = e_params.save_prefix * ".txt" - hififile = e_params.save_prefix * "_hifi.txt" - texfile = e_params.save_prefix * "_params.tex" - # datafile = e_params.save_prefix * ".jld2" - - write_log(perffile, st.log, "# params = $(e_params)\n") - #write_log(hififile, st.log_hifi, "# params = $(e_params)\n") - #write_tex(texfile, e_params) - # @save datafile x y st params - end - - close(datachannel) - wait(gentask) - - if st.aborted - break - end - else - println("Skipping already computed experiment \"$(ename)\"") - # texfile = e_params.save_prefix * "_params.tex" - # write_tex(texfile, e_params) - end - end - - if visualise - # Tell subprocess to finish, and wait - put!(rc, nothing) - close(rc) - wait(visproc) - end - - return -end - ####################### # Demos and batch runs ####################### @@ -403,6 +198,7 @@ run_experiments(;experiments=(experiment,), save_results=false, save_images=false, + save_prefix=default_save_prefix, visualise=true, recalculate=true, verbose_iter=50, @@ -428,6 +224,7 @@ function batchrun_article(kwargs...) run_experiments(;experiments=experiments_all, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -437,6 +234,7 @@ function batchrun_denoising(;kwargs...) run_experiments(;experiments=denoising_experiments_all, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -450,150 +248,6 @@ batchrun_pet(;kwargs...) end -###################################################### -# Iterator that does visualisation and log collection -###################################################### - -function iterate_visualise(datachannel::Channel{OnlineData{DisplacementT}}, - st :: State, - step :: Function, - params :: NamedTuple) where DisplacementT - try - sc = nothing - - d = take!(datachannel) - - for iter=1:params.maxiter - dnext = take!(datachannel) - st = step(d.b_noisy, d.v, dnext.b_noisy) do calc_objective - stn = st - - if isnothing(stn.start_time) - # The Julia precompiler is a miserable joke, apparently not crossing module - # boundaries, so only start timing after the first iteration. - stn = @set stn.start_time=secs_ns() - end - - verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0 - - # Normalise movement to image dimensions so - # our TikZ plotting code doesn't need to know - # the image pixel size. - sc = 1.0./maximum(size(d.b_true)) - - if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0) - verb_start = secs_ns() - tm = verb_start - stn.start_time - stn.wasted_time - value, x, v, vhist = calc_objective() - - entry = LogEntry(iter, tm, value, - #sc*d.v_cumul_true[1], - #sc*d.v_cumul_true[2], - #sc*v[1], sc*v[2], - assess_psnr(x, d.b_true), - assess_ssim(x, d.b_true), - #assess_psnr(d.b_noisy, d.b_true), - #assess_ssim(d.b_noisy, d.b_true) - ) - - # (**) Collect a singly-linked list of log to avoid array resizing - # while iterating - stn = @set stn.log=LinkedListEntry(entry, stn.log) - - if !isnothing(vhist) - vhist=vhist.*sc - end - - if verb - @printf("%d/%d J=%f, PSNR=%f, SSIM=%f, avg. FPS=%f\n", - iter, params.maxiter, value, entry.psnr, - entry.ssim, entry.iter/entry.time) - if isa(stn.vis, Channel) - put_onlylatest!(stn.vis, ((d.b_noisy, x), - params.plot_movement, - stn.log, vhist)) - - end - end - - if params.save_images && (!haskey(params, :save_images_iters) - || iter ∈ params.save_images_iters) - fn = (t, ext) -> "$(params.save_prefix)_$(t)_frame$(iter).$(ext)" - # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true)) - # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy)) - save(File(format"PNG", fn("reco", "png")), grayimg(x)) - if !isnothing(vhist) - open(fn("movement", "txt"), "w") do io - writedlm(io, ["est_y" "est_x"]) - writedlm(io, vhist) - end - end - end - - stn = @set stn.wasted_time += (secs_ns() - verb_start) - - return stn - end - - hifientry = LogEntryHiFi(iter, sc*d.v_cumul_true[1], sc*d.v_cumul_true[2]) - st = @set st.log_hifi=LinkedListEntry(hifientry, st.log_hifi) - - return st - end - d=dnext - end - catch ex - if params.handle_interrupt && 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) - st = @set st.aborted = true - else - rethrow(ex) - end - end - - return st -end - -function bg_visualise_enhanced(rc; fullscreen=false) - process_channel(rc) do d - imgs, plot_movement, log, vhist = d - do_visualise(imgs, refresh=false, fullscreen=fullscreen) - # Overlay movement - GR.settextcolorind(5) - GR.setcharheight(0.015) - GR.settextpath(GR.TEXT_PATH_RIGHT) - tx, ty = GR.wctondc(0, 1) - GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr) - if plot_movement - sc=1.0 - p=unfold_linked_list(log) - x=map(e -> 1.5+sc*e.v_cumul_true_x, p) - y=map(e -> 0.5+sc*e.v_cumul_true_y, p) - GR.setlinewidth(2) - GR.setlinecolorind(2) - GR.polyline(x, y) - x=map(e -> 1.5+sc*e.v_cumul_est_x, p) - y=map(e -> 0.5+sc*e.v_cumul_est_y, p) - GR.setlinecolorind(3) - GR.polyline(x, y) - if vhist != nothing - GR.setlinecolorind(4) - x=map(v -> 1.5+sc*v, vhist[:,2]) - y=map(v -> 0.5+sc*v, vhist[:,1]) - GR.polyline(x, y) - end - end - GR.updatews() - end -end - ######################### # Plotting SSIM and PSNR ######################### @@ -604,13 +258,4 @@ # fv_plot("lighthouse") #end - -############### -# Precompiling -############### - -# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}}) -# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool}) -# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}}) - end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Run.jl Thu Apr 25 13:05:40 2024 -0500 @@ -0,0 +1,360 @@ +################## +# Experiment running and interactive visualisation +################## + +module Run + +import GR +using Setfield +using Printf +using ImageQualityIndexes: assess_psnr, assess_ssim +using FileIO +#using JLD2 +using DelimitedFiles + +using AlgTools.Util +using AlgTools.StructTools +using AlgTools.LinkedLists +using AlgTools.Comms + +using ImageTools.Visualise: secs_ns, grayimg, do_visualise + +using ..ImGenerate +using ..OpticalFlow: identifier, DisplacementFull, DisplacementConstant + +export run_experiments, + Experiment, + State, + LogEntry, + LogEntryHiFi + +################ +# Experiment +################ + +struct Experiment + mod :: Module + DisplacementT :: Type + imgen :: ImGen + params :: NamedTuple +end + +function Base.show(io::IO, e::Experiment) + displacementname(::Type{DisplacementFull}) = "DisplacementFull" + displacementname(::Type{DisplacementConstant}) = "DisplacementConstant" + print(io, " + mod: $(e.mod) + DisplacementT: $(displacementname(e.DisplacementT)) + imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2]) + params: $(e.params ⬿ (kernel = "(not shown)",)) + ") +end + +################ +# Log +################ + +struct LogEntry <: IterableStruct + iter :: Int + time :: Float64 + function_value :: Float64 + #v_cumul_true_y :: Float64 + #v_cumul_true_x :: Float64 + #v_cumul_est_y :: Float64 + #v_cumul_est_x :: Float64 + psnr :: Float64 + ssim :: Float64 + #psnr_data :: Float64 + #ssim_data :: Float64 +end + +struct LogEntryHiFi <: IterableStruct + iter :: Int + v_cumul_true_y :: Float64 + v_cumul_true_x :: Float64 +end + +############### +# Main routine +############### + +struct State + vis :: Union{Channel,Bool,Nothing} + start_time :: Union{Real,Nothing} + wasted_time :: Real + log :: LinkedList{LogEntry} + log_hifi :: LinkedList{LogEntryHiFi} + aborted :: Bool +end + +function name(e::Experiment, p) + ig = e.imgen + id = if haskey(p, :predictor) && ~isnothing(p.predictor) + identifier(p.predictor) + else + e.mod.identifier + end + # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))" + return "$(ig.name)_$(id)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))" +end + +function write_tex(texfile, e_params) + open(texfile, "w") do io + wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}") + wf = (n, s) -> if isdefined(e_params, s) + wp(n, getfield(e_params, s)) + end + wf("alpha", :α) + wf("sigmazero", :σ₀) + wf("tauzero", :τ₀) + wf("tildetauzero", :τ̃₀) + wf("delta", :δ) + wf("lambda", :λ) + wf("theta", :θ) + wf("maxiter", :maxiter) + wf("noiselevel", :noise_level) + wf("shakenoiselevel", :shake_noise_level) + wf("shake", :shake) + wf("timestep", :timestep) + wf("displacementcount", :displacementcount) + wf("phantomrho", :phantom_ρ) + if isdefined(e_params, :σ₀) + wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}") + end + end +end + +function run_experiments(;visualise=true, + visfn = iterate_visualise, + datatype = OnlineData, + recalculate=true, + experiments, + save_prefix="", + fullscreen=false, + kwargs...) + + # Create visualisation + if visualise + rc = Channel(1) + visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen) + bind(rc, visproc) + vis = rc + else + vis = false + end + + # Run all experiments + for e ∈ experiments + + # Parameters for this experiment + e_params = e.params ⬿ kwargs + ename = name(e, e_params) + e_params = e_params ⬿ (save_prefix = save_prefix * ename, + dynrange = e.imgen.dynrange, + Λ = e.imgen.Λ) + + if recalculate || !isfile(e_params.save_prefix * ".txt") + println("Running experiment \"$(ename)\"") + + # Start data generation task + datachannel = Channel{datatype{e.DisplacementT}}(2) + gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params) + bind(datachannel, gentask) + + # Run algorithm + iterate = curry(visfn, datachannel, + State(vis, nothing, 0.0, nothing, nothing, false)) + + x, y, st = e.mod.solve(e.DisplacementT; + dim=e.imgen.dim, + iterate=iterate, + params=e_params) + + # Clear non-saveable things + st = @set st.vis = nothing + + println("Wasted_time: $(st.wasted_time)s") + + if e_params.save_results + println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)") + + perffile = e_params.save_prefix * ".txt" + hififile = e_params.save_prefix * "_hifi.txt" + texfile = e_params.save_prefix * "_params.tex" + # datafile = e_params.save_prefix * ".jld2" + + write_log(perffile, st.log, "# params = $(e_params)\n") + #write_log(hififile, st.log_hifi, "# params = $(e_params)\n") + #write_tex(texfile, e_params) + # @save datafile x y st params + end + + close(datachannel) + wait(gentask) + + if st.aborted + break + end + else + println("Skipping already computed experiment \"$(ename)\"") + # texfile = e_params.save_prefix * "_params.tex" + # write_tex(texfile, e_params) + end + end + + if visualise + # Tell subprocess to finish, and wait + put!(rc, nothing) + close(rc) + wait(visproc) + end + + return +end + + +###################################################### +# Iterator that does visualisation and log collection +###################################################### + +function iterate_visualise(datachannel::Channel{OnlineData{DisplacementT}}, + st :: State, + step :: Function, + params :: NamedTuple) where DisplacementT + try + sc = nothing + + d = take!(datachannel) + + for iter=1:params.maxiter + dnext = take!(datachannel) + st = step(d.b_noisy, d.v, dnext.b_noisy) do calc_objective + stn = st + + if isnothing(stn.start_time) + # The Julia precompiler is a miserable joke, apparently not crossing module + # boundaries, so only start timing after the first iteration. + stn = @set stn.start_time=secs_ns() + end + + verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0 + + # Normalise movement to image dimensions so + # our TikZ plotting code doesn't need to know + # the image pixel size. + sc = 1.0./maximum(size(d.b_true)) + + if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0) + verb_start = secs_ns() + tm = verb_start - stn.start_time - stn.wasted_time + value, x, v, vhist = calc_objective() + + entry = LogEntry(iter, tm, value, + #sc*d.v_cumul_true[1], + #sc*d.v_cumul_true[2], + #sc*v[1], sc*v[2], + assess_psnr(x, d.b_true), + assess_ssim(x, d.b_true), + #assess_psnr(d.b_noisy, d.b_true), + #assess_ssim(d.b_noisy, d.b_true) + ) + + # (**) Collect a singly-linked list of log to avoid array resizing + # while iterating + stn = @set stn.log=LinkedListEntry(entry, stn.log) + + if !isnothing(vhist) + vhist=vhist.*sc + end + + if verb + @printf("%d/%d J=%f, PSNR=%f, SSIM=%f, avg. FPS=%f\n", + iter, params.maxiter, value, entry.psnr, + entry.ssim, entry.iter/entry.time) + if isa(stn.vis, Channel) + put_onlylatest!(stn.vis, ((d.b_noisy, x), + params.plot_movement, + stn.log, vhist)) + + end + end + + if params.save_images && (!haskey(params, :save_images_iters) + || iter ∈ params.save_images_iters) + fn = (t, ext) -> "$(params.save_prefix)_$(t)_frame$(iter).$(ext)" + # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true)) + # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy)) + save(File(format"PNG", fn("reco", "png")), grayimg(x)) + if !isnothing(vhist) + open(fn("movement", "txt"), "w") do io + writedlm(io, ["est_y" "est_x"]) + writedlm(io, vhist) + end + end + end + + stn = @set stn.wasted_time += (secs_ns() - verb_start) + + return stn + end + + hifientry = LogEntryHiFi(iter, sc*d.v_cumul_true[1], sc*d.v_cumul_true[2]) + st = @set st.log_hifi=LinkedListEntry(hifientry, st.log_hifi) + + return st + end + d=dnext + end + catch ex + if params.handle_interrupt && 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) + st = @set st.aborted = true + else + rethrow(ex) + end + end + + return st +end + +function bg_visualise_enhanced(rc; fullscreen=false) + process_channel(rc) do d + imgs, plot_movement, log, vhist = d + do_visualise(imgs, refresh=false, fullscreen=fullscreen) + # Overlay movement + GR.settextcolorind(5) + GR.setcharheight(0.015) + GR.settextpath(GR.TEXT_PATH_RIGHT) + tx, ty = GR.wctondc(0, 1) + GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr) + if plot_movement + sc=1.0 + p=unfold_linked_list(log) + x=map(e -> 1.5+sc*e.v_cumul_true_x, p) + y=map(e -> 0.5+sc*e.v_cumul_true_y, p) + GR.setlinewidth(2) + GR.setlinecolorind(2) + GR.polyline(x, y) + x=map(e -> 1.5+sc*e.v_cumul_est_x, p) + y=map(e -> 0.5+sc*e.v_cumul_est_y, p) + GR.setlinecolorind(3) + GR.polyline(x, y) + if vhist != nothing + GR.setlinecolorind(4) + x=map(v -> 1.5+sc*v, vhist[:,2]) + y=map(v -> 0.5+sc*v, vhist[:,1]) + GR.polyline(x, y) + end + end + GR.updatews() + end +end + +end # module