# HG changeset patch # User Neil Dizon # Date 1713767179 -10800 # Node ID ccd22bbbb02f0d312c8a708355f9d283208df845 # Parent c9b06736a47722284c785e524a8198f1a5a1dcfa added dummy structs for different flows diff -r c9b06736a477 -r ccd22bbbb02f src/AlgorithmDualScaling.jl --- a/src/AlgorithmDualScaling.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/AlgorithmDualScaling.jl Mon Apr 22 09:26:19 2024 +0300 @@ -16,7 +16,8 @@ using ..OpticalFlow: ImageSize, Image, - pdflow! + pdflow!, + DualScaling ######################### # Iterate initialisation @@ -90,7 +91,7 @@ init_data = false end - pdflow!(x, y, v_known) + pdflow!(x, Δx, y, Δy, v_known, DualScaling()) ############ # PDPS step diff -r c9b06736a477 -r ccd22bbbb02f src/AlgorithmGreedy.jl --- a/src/AlgorithmGreedy.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/AlgorithmGreedy.jl Mon Apr 22 09:26:19 2024 +0300 @@ -16,7 +16,8 @@ using ..OpticalFlow: ImageSize, Image, - pdflow! + pdflow!, + Greedy ######################### # Iterate initialisation @@ -90,7 +91,7 @@ init_data = false end - pdflow!(x, Δx, y, Δy, v_known) + pdflow!(x, Δx, y, Δy, v_known, Greedy()) ############ # PDPS step diff -r c9b06736a477 -r ccd22bbbb02f src/AlgorithmNoPrediction.jl --- a/src/AlgorithmNoPrediction.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/AlgorithmNoPrediction.jl Mon Apr 22 09:26:19 2024 +0300 @@ -90,7 +90,7 @@ init_data = false end - + # No prediction ############ # PDPS step diff -r c9b06736a477 -r ccd22bbbb02f src/AlgorithmRotation.jl --- a/src/AlgorithmRotation.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/AlgorithmRotation.jl Mon Apr 22 09:26:19 2024 +0300 @@ -16,7 +16,8 @@ using ..OpticalFlow: ImageSize, Image, - pdflow! + pdflow!, + Rotation ######################### # Iterate initialisation @@ -90,7 +91,7 @@ init_data = false end - pdflow!(x, Δx, y, v_known) + pdflow!(x, Δx, y, Δy, v_known, Rotation()) ############ # PDPS step diff -r c9b06736a477 -r ccd22bbbb02f src/OpticalFlow.jl --- a/src/OpticalFlow.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/OpticalFlow.jl Mon Apr 22 09:26:19 2024 +0300 @@ -41,7 +41,8 @@ DisplacementFull, DisplacementConstant, HornSchunckData, filter_hs, - petpdflow! + petpdflow!, + DualScaling, Greedy, Rotation ############################################### # Types (several imported from ImageTools.Translate) @@ -55,6 +56,14 @@ Gradient = Array{Float64,3} ImageSize = Tuple{Int64,Int64} + +################################# +# Struct for flow +################################# +struct DualScaling end +struct Greedy end +struct Rotation end + ################################# # Displacement field based flow ################################# @@ -104,7 +113,7 @@ end # Additional method for Greedy -function pdflow!(x, Δx, y, Δy, u; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none) @assert(size(u)==(2,)) Δx .= x Δy .= y @@ -120,17 +129,15 @@ end # Additional method for Rotation -function pdflow!(x, Δx, y, u; 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)) - (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) @@ -151,7 +158,7 @@ end # Additional method for Dual Scaling -function pdflow!(x, y, u; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none) @assert(size(u)==(2,)) oldx = copy(x) flow!(x, u; threads=(threads==:inner)) @@ -193,27 +200,25 @@ end # Method for greedy predictor -function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none) +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 - if dual_flow - Dxx = copy(Δy) - DΔx = copy(Δy) - ∇₂!(Dxx, x) - ∇₂!(DΔx, oldx) - inds = abs.(Dxx) .≤ β - Dxx[inds] .= 1 - DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx - end + 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 affine predictor -function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none) +# 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) @@ -221,26 +226,20 @@ @. x = Δx C = similar(y) cc = abs.(x-oldx) - if dual_flow - cm = max(1e-12,maximum(cc)) - c = 1 .* (1 .- cc./ cm) .^(10) - C[1,:,:] .= c - C[2,:,:] .= c - y .= C.*y - end + 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, 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) +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 diff -r c9b06736a477 -r ccd22bbbb02f src/PET/AlgorithmDualScaling.jl --- a/src/PET/AlgorithmDualScaling.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/PET/AlgorithmDualScaling.jl Mon Apr 22 09:26:19 2024 +0300 @@ -22,7 +22,8 @@ using ..OpticalFlow: ImageSize, Image, - petpdflow! + petpdflow!, + DualScaling ######################### # Iterate initialisation @@ -124,7 +125,7 @@ # Prediction steps ################### - petpdflow!(x, Δx, y, v_known, theta_known, true) + petpdflow!(x, Δx, y, Δy, v_known, theta_known, DualScaling()) if params.L_experiment @. oldx = x diff -r c9b06736a477 -r ccd22bbbb02f src/PET/AlgorithmGreedy.jl --- a/src/PET/AlgorithmGreedy.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/PET/AlgorithmGreedy.jl Mon Apr 22 09:26:19 2024 +0300 @@ -22,7 +22,8 @@ using ..OpticalFlow: ImageSize, Image, - petpdflow! + petpdflow!, + Greedy ######################### # Iterate initialisation @@ -124,7 +125,7 @@ # Prediction steps ################### - petpdflow!(x, Δx, y, Δy, v_known, theta_known, true, 1e-2) + petpdflow!(x, Δx, y, Δy, v_known, theta_known, Greedy()) if params.L_experiment @. oldx = x diff -r c9b06736a477 -r ccd22bbbb02f src/PET/AlgorithmProximal.jl --- a/src/PET/AlgorithmProximal.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/PET/AlgorithmProximal.jl Mon Apr 22 09:26:19 2024 +0300 @@ -151,13 +151,14 @@ @. oldx = x end - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - #@. cc = y + 1000000*σ̃*Δy - #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc - proj_norm₂₁ball!(y, α) - end + ############################## + # Proximal step of prediction + ############################## + ∇₂!(Δy, x) + @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + #@. cc = y + 1000000*σ̃*Δy + #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc + proj_norm₂₁ball!(y, α) ############ # PDPS step diff -r c9b06736a477 -r ccd22bbbb02f src/PET/AlgorithmRotation.jl --- a/src/PET/AlgorithmRotation.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/PET/AlgorithmRotation.jl Mon Apr 22 09:26:19 2024 +0300 @@ -22,7 +22,8 @@ using ..OpticalFlow: ImageSize, Image, - petpdflow! + petpdflow!, + Rotation ######################### # Iterate initialisation @@ -124,7 +125,7 @@ # Prediction steps ################### - petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow, 0,0) + petpdflow!(x, Δx, y, Δy, v_known, theta_known, Rotation()) if params.L_experiment @. oldx = x diff -r c9b06736a477 -r ccd22bbbb02f src/PET/OpticalFlow.jl --- a/src/PET/OpticalFlow.jl Sun Apr 21 22:35:03 2024 +0300 +++ b/src/PET/OpticalFlow.jl Mon Apr 22 09:26:19 2024 +0300 @@ -41,7 +41,8 @@ DisplacementFull, DisplacementConstant, HornSchunckData, filter_hs, - petpdflow! + petpdflow!, + DualScaling, Greedy, Rotation ############################################### # Types (several imported from ImageTools.Translate) @@ -55,6 +56,14 @@ Gradient = Array{Float64,3} ImageSize = Tuple{Int64,Int64} + +################################# +# Struct for flow +################################# +struct DualScaling end +struct Greedy end +struct Rotation end + ################################# # Displacement field based flow ################################# @@ -104,7 +113,7 @@ end # Additional method for Greedy -function pdflow!(x, Δx, y, Δy, u; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none) @assert(size(u)==(2,)) Δx .= x Δy .= y @@ -120,17 +129,15 @@ end # Additional method for Rotation -function pdflow!(x, Δx, y, u; 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)) - (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) @@ -151,7 +158,7 @@ end # Additional method for Dual Scaling -function pdflow!(x, y, u; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none) @assert(size(u)==(2,)) oldx = copy(x) flow!(x, u; threads=(threads==:inner)) @@ -193,27 +200,25 @@ end # Method for greedy predictor -function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none) +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 - if dual_flow - Dxx = copy(Δy) - DΔx = copy(Δy) - ∇₂!(Dxx, x) - ∇₂!(DΔx, oldx) - inds = abs.(Dxx) .≤ β - Dxx[inds] .= 1 - DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx - end + 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 affine predictor -function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none) +# 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) @@ -221,26 +226,20 @@ @. x = Δx C = similar(y) cc = abs.(x-oldx) - if dual_flow - cm = max(1e-12,maximum(cc)) - c = 1 .* (1 .- cc./ cm) .^(10) - C[1,:,:] .= c - C[2,:,:] .= c - y .= C.*y - end + 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, 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) +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