Mercurial > repos > PredictPDPS / changeset
changeset
Remove unused old alg attempts
Thu, 25 Apr 2024 13:25:48 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Thu, 25 Apr 2024 13:25:48 -0500
- changeset 37
- bba159cf1636
- parent 36
- e4a8f662a1ac
- child 39
- b180a4f3b9bd
Remove unused old alg attempts
| src/AlgorithmBoth.jl | file | annotate | diff | comparison | revisions | |
| src/AlgorithmBothCumul.jl | file | annotate | diff | comparison | revisions | |
| src/AlgorithmBothGreedyV.jl | file | annotate | diff | comparison | revisions | |
| src/AlgorithmBothNL.jl | file | annotate | diff | comparison | revisions |
--- a/src/AlgorithmBoth.jl Thu Apr 25 13:05:40 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,239 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBoth - -identifier = "pdps_unknown_basic" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - Gradient, - DisplacementConstant, - DisplacementFull, - pdflow!, - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - filter_hs - -using ..Algorithm: step_lengths - -############# -# Data types -############# - -struct Primal{DisplacementT} - x :: Image - u :: DisplacementT -end - -function Base.similar(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.similar(x.x), Base.similar(x.u)) -end - -function Base.copy(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.copy(x.x), Base.copy(x.u)) - end - -struct Dual - tv :: Gradient - flow :: Image -end - -function Base.similar(y::Dual) - return Dual(Base.similar(y.tv), Base.similar(y.flow)) -end - -function Base.copy(y::Dual) - return Dual(Base.copy(y.tv), Base.copy(y.flow)) - end - -######################### -# Iterate initialisation -######################### - -function init_primal(xinit::Image, ::Type{DisplacementConstant}) - return Primal{DisplacementConstant}(xinit, zeros(2)) -end - -function init_primal(xinit::Image, ::Type{DisplacementFull}) - return Primal{DisplacementFull}(xinit, zeros(2, size(xinit)...)) -end - -function init_rest(x::Primal{DisplacementT}) where DisplacementT - imdim=size(x.x) - - y = Dual(zeros(2, imdim...), zeros(imdim)) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Primal{DisplacementT}) where DisplacementT - return init_rest(copy(xinit)) -end - -function init_iterates( :: Type{DisplacementT}, xinit::Image) where DisplacementT - return init_rest(init_primal(copy(xinit), DisplacementT)) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_primal(zeros(dim...), DisplacementT)) -end - -############################################## -# Weighting for different displacements types -############################################## - -norm²weight( :: Type{DisplacementConstant}, sz ) = prod(sz) -norm²weight( :: Type{DisplacementFull}, sz ) = 1 - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = max(∇₂_norm₂₂_est², ∇₂_norm₂∞_est²*params.dynrange^2) - γ = min(1, λ*norm²weight(DisplacementT, size(x.x))) - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - 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) - - ############################ - # Construct K for this step - ############################ - - K! = (yʹ, xʹ) -> begin - # Optical flow part: ⟨⟨u, ∇b_k⟩⟩. - # Use y.tv as temporary gradient storage. - pointwise_gradiprod_2d!(yʹ.flow, yʹ.tv, xʹ.u, b_filt) - #@. yʹ.flow = -yʹ.flow - # TV part - ∇₂!(yʹ.tv, xʹ.x) - end - Kᵀ! = (xʹ, yʹ) -> begin - # Optical flow part: ∇b_k ⋅ y - pointwise_gradiprod_2dᵀ!(xʹ.u, yʹ.flow, b_filt) - #@. xʹ.u = -xʹ.u - # TV part - ∇₂ᵀ!(xʹ.x, yʹ.tv) - end - - ################## - # Prediction step - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x.x, Δx.x, y.tv, Δy.tv, y.flow, Δy.flow, x.u, params.dual_flow) - - # Predict zero displacement - x.u .= 0 - if params.prox_predict - K!(Δy, x) - @. y.tv = (y.tv + σ̃*Δy.tv)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y.tv, α) - @. y.flow = (y.flow+σ̃*((b_next_filt-b_filt)/T+Δy.flow))/(1+σ̃*(ρ̃+1/θ)) - end - - ############ - # PDPS step - # - # NOTE: For DisplacementConstant, the x.u update is supposed to be with - # respect to the 𝟙^*𝟙 norm/inner product that makes the norm equivalent - # to full-space norm when restricted to constant displacements. Since - # `OpticalFlow.pointwise_gradiprod_2dᵀ!` already uses this inner product, - # and the λ-weighted term in the problem is with respect to this norm, - # all the norm weights disappear in this update. - ############ - - Kᵀ!(Δx, y) # primal step: - @. x̄.x = x.x # | save old x for over-relax - @. x̄.u = x.u # | - @. x.x = (x.x-τ*(Δx.x-b))/(1+τ) # | prox - @. x.u = (x.u-τ*Δx.u)/(1+τ*λ) # | - @. x̄.x = 2x.x - x̄.x # over-relax - @. x̄.u = 2x.u - x̄.u # | - K!(Δy, x̄) # dual step: y - @. y.tv = (y.tv + σ*Δy.tv)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y.tv, α) # | prox - @. y.flow = (y.flow+σ*((b_next_filt-b_filt)/T+Δy.flow))/(1+σ/θ) - - if DisplacementT == DisplacementConstant - ucumul .+= x.u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - K!(Δy, x) - value = (norm₂²(b-x.x)/2 + θ*norm₂²((b_next_filt-b_filt)./T+Δy.flow) - + λ*norm₂²(x.u)/2 + α*γnorm₂₁(Δy.tv, ρ)) - - value, x.x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothCumul.jl Thu Apr 25 13:05:40 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,139 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothCumul - -identifier = "pdps_unknown_cumul" - -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! - -using ..AlgorithmBothGreedyV: init_iterates -using ..Algorithm: step_lengths - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = zeros(size(u)...) - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b₀=nothing - b₀_filt=nothing - u_prev=zeros(size(u)) - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ######################################################### - # Smoothen data for Horn–Schunck term; zero initial data - ######################################################### - - b_filt = (kernel==nothing ? b : simple_imfilter(b, kernel)) - - if b₀ == nothing - b₀ = b - b₀_filt = b_filt - end - - ################################################ - # Prediction step - # We leave u as-is in this cumulative version - ################################################ - - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, u-u_prev, params.dual_flow) - - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - # Store current cumulative displacement before updating in next step. - u_prev .= u - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - horn_schunck_reg_prox!(u, b_filt, b₀_filt, τ, θ, λ, T) - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δ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) - tmp = zeros(size(b_filt)) - pointwise_gradiprod_2d!(tmp, Δy, u, b₀_filt) - value = (norm₂²(b-x)/2 + θ*norm₂²((b_filt-b₀_filt)./T+tmp) - + λ*norm₂²(u)/2 + α*γnorm₂₁(Δy, ρ)) - - value, x, u, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothGreedyV.jl Thu Apr 25 13:05:40 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,167 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothGreedyV - -identifier = "pdps_unknown_greedyv" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: Image, - ImageSize, - DisplacementConstant, - DisplacementFull, - pdflow!, - horn_schunck_reg_prox!, - pointwise_gradiprod_2d!, - filter_hs - -using ..Algorithm: step_lengths - -######################### -# Iterate initialisation -######################### - -function init_displ(xinit::Image, ::Type{DisplacementConstant}) - return xinit, zeros(2) -end - -function init_displ(xinit::Image, ::Type{DisplacementFull}) - return xinit, zeros(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) where DisplacementT - return init_rest(init_displ(copy(xinit), DisplacementT)...) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_displ(zeros(dim...), DisplacementT)...) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - 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 - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, u, params.dual_flow) - - # Predict zero displacement - u .= 0 - 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 - horn_schunck_reg_prox!(u, b_next_filt, b_filt, θ, λ, T, τ) - @. x̄ = 2x - x̄ # over-relax - ∇₂!(y, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - if DisplacementT == DisplacementConstant - ucumul .+= u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - ∇₂!(Δy, x) - tmp = zeros(size(b_filt)) - pointwise_gradiprod_2d!(tmp, Δy, u, b_filt) - value = (norm₂²(b-x)/2 + θ*norm₂²((b_next_filt-b_filt)./T+tmp) - + λ*norm₂²(u)/2 + α*γnorm₂₁(Δy, ρ)) - - value, x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothNL.jl Thu Apr 25 13:05:40 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,242 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothNL - -identifier = "pdps_unknown_nl" - - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - Gradient, - DisplacementConstant, - DisplacementFull, - pdflow!, - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - filter_hs - -using ..Algorithm: step_lengths - -############# -# Data types -############# - -struct Primal{DisplacementT} - x :: Image - u :: DisplacementT -end - -function Base.similar(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.similar(x.x), Base.similar(x.u)) -end - -function Base.copy(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.copy(x.x), Base.copy(x.u)) - end - -struct Dual - tv :: Gradient - flow :: Image -end - -function Base.similar(y::Dual) - return Dual(Base.similar(y.tv), Base.similar(y.flow)) -end - -function Base.copy(y::Dual) - return Dual(Base.copy(y.tv), Base.copy(y.flow)) - end - -######################### -# Iterate initialisation -######################### - -function init_primal(xinit::Image, ::Type{DisplacementConstant}) - return Primal{DisplacementConstant}(xinit, zeros(2)) -end - -function init_primal(xinit::Image, ::Type{DisplacementFull}) - return Primal{DisplacementFull}(xinit, zeros(2, size(xinit)...)) -end - -function init_rest(x::Primal{DisplacementT}) where DisplacementT - imdim=size(x.x) - - y = Dual(zeros(2, imdim...), zeros(imdim)) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Primal{DisplacementT}) where DisplacementT - return init_rest(copy(xinit)) -end - -function init_iterates( :: Type{DisplacementT}, xinit::Image) where DisplacementT - return init_rest(init_primal(copy(xinit), DisplacementT)) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_primal(zeros(dim...), DisplacementT)) -end - -############################################## -# Weighting for different displacements types -############################################## - -norm²weight( :: Type{DisplacementConstant}, sz ) = prod(sz) -norm²weight( :: Type{DisplacementFull}, sz ) = 1 - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = max(∇₂_norm₂₂_est², ∇₂_norm₂∞_est²*params.dynrange^2) - γ = min(1, λ*norm²weight(DisplacementT, size(x.x))) - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - 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) - - ############################ - # Construct K for this step - ############################ - - K! = (yʹ, xʹ) -> begin - # Optical flow part - @. yʹ.flow = b_filt - flow!(yʹ.flow, Δx.x, xʹ.u) - @. yʹ.flow = yʹ.flow - b_next_filt - # TV part - ∇₂!(yʹ.tv, xʹ.x) - end - Kᵀ! = (xʹ, yʹ) -> begin - # Optical flow part: ∇b_k ⋅ y - # - # TODO: This really should depend x.u, but x.u is zero. - # - pointwise_gradiprod_2dᵀ!(xʹ.u, yʹ.flow, b_filt) - # TV part - ∇₂ᵀ!(xʹ.x, yʹ.tv) - end - - ################## - # Prediction step - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x.x, Δx.x, y.tv, Δy.tv, y.flow, Δy.flow, x.u, params.dual_flow) - - # Predict zero displacement - x.u .= 0 - if params.prox_predict - K!(Δy, x) - @. y.tv = (y.tv + σ̃*Δy.tv)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y.tv, α) - @. y.flow = (y.flow+σ̃*Δy.flow)/(1+σ̃*(ρ̃+1/θ)) - end - - ############ - # PDPS step - # - # NOTE: For DisplacementConstant, the x.u update is supposed to be with - # respect to the 𝟙^*𝟙 norm/inner product that makes the norm equivalent - # to full-space norm when restricted to constant displacements. Since - # `OpticalFlow.pointwise_gradiprod_2dᵀ!` already uses this inner product, - # and the λ-weighted term in the problem is with respect to this norm, - # all the norm weights disappear in this update. - ############ - - Kᵀ!(Δx, y) # primal step: - @. x̄.x = x.x # | save old x for over-relax - @. x̄.u = x.u # | - @. x.x = (x.x-τ*(Δx.x-b))/(1+τ) # | prox - @. x.u = (x.u-τ*Δx.u)/(1+τ*λ) # | - @. x̄.x = 2x.x - x̄.x # over-relax - @. x̄.u = 2x.u - x̄.u # | - K!(Δy, x̄) # dual step: y - @. y.tv = (y.tv + σ*Δy.tv)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y.tv, α) # | prox - @. y.flow = (y.flow+σ*Δy.flow)/(1+σ/θ) - - if DisplacementT == DisplacementConstant - ucumul .+= x.u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - K!(Δy, x) - value = (norm₂²(b-x.x)/2 + θ*norm₂²(Δy.flow) - + λ*norm₂²(x.u)/2 + α*γnorm₂₁(Δy.tv, ρ)) - - value, x.x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
