|
1 ###################################################################### |
|
2 # Predictive online PDPS for optical flow with unknown velocity field |
|
3 ###################################################################### |
|
4 |
|
5 __precompile__() |
|
6 |
|
7 module AlgorithmBothCumul |
|
8 |
|
9 identifier = "pdps_unknown_cumul" |
|
10 |
|
11 using Printf |
|
12 |
|
13 using AlgTools.Util |
|
14 import AlgTools.Iterate |
|
15 using ImageTools.Gradient |
|
16 using ImageTools.ImFilter |
|
17 |
|
18 using ..OpticalFlow: Image, |
|
19 ImageSize, |
|
20 DisplacementConstant, |
|
21 pdflow!, |
|
22 horn_schunck_reg_prox!, |
|
23 pointwise_gradiprod_2d! |
|
24 |
|
25 using ..AlgorithmBothGreedyV: init_iterates |
|
26 using ..Algorithm: step_lengths |
|
27 |
|
28 ############ |
|
29 # Algorithm |
|
30 ############ |
|
31 |
|
32 function solve( :: Type{DisplacementT}; |
|
33 dim :: ImageSize, |
|
34 iterate = AlgTools.simple_iterate, |
|
35 params::NamedTuple) where DisplacementT |
|
36 |
|
37 ###################### |
|
38 # Initialise iterates |
|
39 ###################### |
|
40 |
|
41 x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim) |
|
42 init_data = (params.init == :data) |
|
43 |
|
44 # … for tracking cumulative movement |
|
45 if DisplacementT == DisplacementConstant |
|
46 ucumul = zeros(size(u)...) |
|
47 end |
|
48 |
|
49 ############################################# |
|
50 # Extract parameters and set up step lengths |
|
51 ############################################# |
|
52 |
|
53 α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep |
|
54 R_K² = ∇₂_norm₂₂_est² |
|
55 γ = 1 |
|
56 τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) |
|
57 |
|
58 kernel = params.kernel |
|
59 |
|
60 #################### |
|
61 # Run the algorithm |
|
62 #################### |
|
63 |
|
64 b₀=nothing |
|
65 b₀_filt=nothing |
|
66 u_prev=zeros(size(u)) |
|
67 |
|
68 v = iterate(params) do verbose :: Function, |
|
69 b :: Image, |
|
70 🚫unused_v_known :: DisplacementT, |
|
71 🚫unused_b_next :: Image |
|
72 |
|
73 ######################################################### |
|
74 # Smoothen data for Horn–Schunck term; zero initial data |
|
75 ######################################################### |
|
76 |
|
77 b_filt = (kernel==nothing ? b : simple_imfilter(b, kernel)) |
|
78 |
|
79 if b₀ == nothing |
|
80 b₀ = b |
|
81 b₀_filt = b_filt |
|
82 end |
|
83 |
|
84 ################################################ |
|
85 # Prediction step |
|
86 # We leave u as-is in this cumulative version |
|
87 ################################################ |
|
88 |
|
89 if init_data |
|
90 x .= b |
|
91 init_data = false |
|
92 end |
|
93 |
|
94 pdflow!(x, Δx, y, Δy, u-u_prev, params.dual_flow) |
|
95 |
|
96 if params.prox_predict |
|
97 ∇₂!(Δy, x) |
|
98 @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) |
|
99 proj_norm₂₁ball!(y, α) |
|
100 end |
|
101 |
|
102 # Store current cumulative displacement before updating in next step. |
|
103 u_prev .= u |
|
104 |
|
105 ############ |
|
106 # PDPS step |
|
107 ############ |
|
108 |
|
109 ∇₂ᵀ!(Δx, y) # primal step: |
|
110 @. x̄ = x # | save old x for over-relax |
|
111 @. x = (x-τ*(Δx-b))/(1+τ) # | prox |
|
112 horn_schunck_reg_prox!(u, b_filt, b₀_filt, τ, θ, λ, T) |
|
113 @. x̄ = 2x - x̄ # over-relax |
|
114 ∇₂!(Δy, x̄) # dual step: y |
|
115 @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | |
|
116 proj_norm₂₁ball!(y, α) # | prox |
|
117 |
|
118 ######################################################## |
|
119 # Give function value and cumulative movement if needed |
|
120 ######################################################## |
|
121 v = verbose() do |
|
122 ∇₂!(Δy, x) |
|
123 tmp = zeros(size(b_filt)) |
|
124 pointwise_gradiprod_2d!(tmp, Δy, u, b₀_filt) |
|
125 value = (norm₂²(b-x)/2 + θ*norm₂²((b_filt-b₀_filt)./T+tmp) |
|
126 + λ*norm₂²(u)/2 + α*γnorm₂₁(Δy, ρ)) |
|
127 |
|
128 value, x, u, nothing |
|
129 end |
|
130 |
|
131 return v |
|
132 end |
|
133 |
|
134 return x, y, v |
|
135 end |
|
136 |
|
137 end # Module |
|
138 |
|
139 |