|
1 ###################################################################### |
|
2 # Predictive online PDPS for optical flow with unknown velocity field |
|
3 ###################################################################### |
|
4 |
|
5 __precompile__() |
|
6 |
|
7 module AlgorithmBothMulti |
|
8 |
|
9 identifier = "pdps_unknownmulti" |
|
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 horn_schunck_reg_prox_op!, |
|
25 mldivide_step_plus_sym2x2!, |
|
26 ConstantDisplacementHornSchunckData, |
|
27 filter_hs |
|
28 |
|
29 using ..Algorithm: step_lengths |
|
30 |
|
31 ######################### |
|
32 # Iterate initialisation |
|
33 ######################### |
|
34 |
|
35 function init_displ(xinit::Image, ::Type{DisplacementConstant}, n::Integer) |
|
36 return xinit, zeros(n, 2) |
|
37 end |
|
38 |
|
39 # function init_displ(xinit::Image, ::Type{DisplacementFull}, n::Integer) |
|
40 # return xinit, zeros(n, 2, size(xinit)...) |
|
41 # end |
|
42 |
|
43 function init_rest(x::Image, u::DisplacementT) where DisplacementT |
|
44 imdim=size(x) |
|
45 |
|
46 y = zeros(2, imdim...) |
|
47 Δx = copy(x) |
|
48 Δy = copy(y) |
|
49 x̄ = copy(x) |
|
50 |
|
51 return x, y, Δx, Δy, x̄, u |
|
52 end |
|
53 |
|
54 function init_iterates( :: Type{DisplacementT}, |
|
55 xinit::Image, |
|
56 n::Integer) where DisplacementT |
|
57 return init_rest(init_displ(copy(xinit), DisplacementT, n)...) |
|
58 end |
|
59 |
|
60 function init_iterates( :: Type{DisplacementT}, |
|
61 dim::ImageSize, |
|
62 n::Integer) where DisplacementT |
|
63 return init_rest(init_displ(zeros(dim...), DisplacementT, n)...) |
|
64 end |
|
65 |
|
66 ############ |
|
67 # Algorithm |
|
68 ############ |
|
69 |
|
70 function solve( :: Type{DisplacementT}; |
|
71 dim :: ImageSize, |
|
72 iterate = AlgTools.simple_iterate, |
|
73 params::NamedTuple) where DisplacementT |
|
74 |
|
75 ###################### |
|
76 # Initialise iterates |
|
77 ###################### |
|
78 |
|
79 n = params.displacement_count |
|
80 k = 0 # number of displacements we have already |
|
81 |
|
82 x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim, n) |
|
83 init_data = (params.init == :data) |
|
84 hs = [ConstantDisplacementHornSchunckData() for i=1:n] |
|
85 #hs = Array{ConstantDisplacementHornSchunckData}(undef, n) |
|
86 A = Array{Float64,3}(undef, n, 2, 2) |
|
87 d = Array{Float64,2}(undef, n, 2) |
|
88 |
|
89 # … for tracking cumulative movement |
|
90 ucumulbase = [0.0, 0.0] |
|
91 |
|
92 ############################################# |
|
93 # Extract parameters and set up step lengths |
|
94 ############################################# |
|
95 |
|
96 α, ρ, λ, θ = params.α, params.ρ, params.λ, params.θ |
|
97 R_K² = ∇₂_norm₂₂_est² |
|
98 γ = 1 |
|
99 τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) |
|
100 |
|
101 kernel = params.kernel |
|
102 T = params.timestep |
|
103 |
|
104 #################### |
|
105 # Run the algorithm |
|
106 #################### |
|
107 |
|
108 b_next_filt = nothing |
|
109 diffu = similar(u[1, :]) |
|
110 |
|
111 v = iterate(params) do verbose :: Function, |
|
112 b :: Image, |
|
113 🚫unused_v_known :: DisplacementT, |
|
114 b_next :: Image |
|
115 |
|
116 #################################### |
|
117 # Smooth data for Horn–Schunck term |
|
118 #################################### |
|
119 |
|
120 b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) |
|
121 |
|
122 ################################################ |
|
123 # Prediction step |
|
124 ################################################ |
|
125 |
|
126 # Predict x and y |
|
127 if k==0 |
|
128 if init_data |
|
129 x .= b |
|
130 init_data = false |
|
131 end |
|
132 else |
|
133 # Displacement from previous to this image is estimated as |
|
134 # the difference of their displacements from beginning of window. |
|
135 if k>1 |
|
136 @. @views diffu = u[k, :] - u[k-1, :] |
|
137 else |
|
138 @. @views diffu = u[k, :] |
|
139 end |
|
140 |
|
141 pdflow!(x, Δx, y, Δy, diffu, params.dual_flow) |
|
142 end |
|
143 |
|
144 # Shift stored prox matrices |
|
145 if k==n |
|
146 tmp = copy(u[1, :]) |
|
147 ucumulbase .+= tmp |
|
148 for j=1:(n-1) |
|
149 @. @views u[j, :] = u[j+1, :] - tmp |
|
150 hs[j] = hs[j+1] |
|
151 end |
|
152 # Create new struct as original contains references to objects that |
|
153 # have been moved to index n-1. |
|
154 hs[n]=ConstantDisplacementHornSchunckData() |
|
155 else |
|
156 k += 1 |
|
157 end |
|
158 |
|
159 # Predict u: zero displacement from current to next image, i.e., |
|
160 # same displacement to beginning of window. |
|
161 if k==1 |
|
162 @. @views u[k, :] = 0.0 |
|
163 else |
|
164 @. @views u[k, :] = u[k-1, :] |
|
165 end |
|
166 |
|
167 # Predictor proximals tep |
|
168 if params.prox_predict |
|
169 ∇₂!(Δy, x) |
|
170 @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) |
|
171 proj_norm₂₁ball!(y, α) |
|
172 end |
|
173 |
|
174 ################################################################################# |
|
175 # PDPS step |
|
176 # |
|
177 # For the displacements, with τ̃=τ/k, we need to solve for 2≤j<k, |
|
178 # |
|
179 # (1) (I/τ̃+M₀^j+M₀^{j+1})u^j = M₀^ju^{j-1} + M₀^{j+1}u^{j+1} |
|
180 # + ũ^j/τ̃ - z^j + z^{j+1}, |
|
181 # |
|
182 # as well as |
|
183 # |
|
184 # (2) (I/τ̃+M₀^k)u^k = M₀^k u^{k-1} + ũ^k/τ̃ - z^k |
|
185 # |
|
186 # and |
|
187 # |
|
188 # (3) (I/τ̃+M₀^1+M₀^2)u^1 = 0 + M₀^{2}u^{2} + ũ^1/τ̃ - z^1 + z^{2} |
|
189 # |
|
190 # We first construct from (2) that |
|
191 # |
|
192 # u^k = A^k u^{k-1} + d^k |
|
193 # |
|
194 # for |
|
195 # |
|
196 # A^k := (I/τ̃+M₀^k)^{-1} M₀^k |
|
197 # d_k := (I/τ̃+M₀^k)^{-1} (ũ^k/τ̃ - z^k). |
|
198 # |
|
199 # Inserting this into (1) we need |
|
200 # |
|
201 # (4) (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))u^j = M₀^ju^{j-1} + M₀^{j+1}d^{j+1} |
|
202 # + ũ^j/τ̃ - z^j + z^{j+1}. |
|
203 # |
|
204 # This is well-defined because A^{j+1} < I. It also has the same form as (1), so |
|
205 # we continue with |
|
206 # |
|
207 # (5) u^j = A^j u^{j-1} + d^j |
|
208 # |
|
209 # for |
|
210 # |
|
211 # A^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} M₀^j |
|
212 # d^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} |
|
213 # (M₀^{j+1}d^{j+1} + ũ^j/τ̃ - z^j + z^{j+1}) |
|
214 # |
|
215 # Finally from (3) with these we need |
|
216 # |
|
217 # (I/τ̃+M₀^1+M₀^2(I-A^2))u^1 = M₀^2d^2 + ũ^1/τ̃ - z^1 + z^2, |
|
218 # |
|
219 # which is of the same form as (4) with u^0=0, so by (5) u^1=d^1. |
|
220 # |
|
221 ################################################################################# |
|
222 |
|
223 ∇₂ᵀ!(Δx, y) # primal step: |
|
224 @. x̄ = x # | save old x for over-relax |
|
225 @. x = (x-τ*(Δx-b))/(1+τ) # | prox |
|
226 # | | for displacement |
|
227 # Calculate matrices for latest data; rest is stored. |
|
228 @views begin |
|
229 horn_schunck_reg_prox_op!(hs[k], b_next_filt, b_filt, θ, λ, T) |
|
230 |
|
231 τ̃=τ/k |
|
232 |
|
233 B = hs[k].M₀ |
|
234 c = u[k, :]./τ̃-hs[k].z |
|
235 mldivide_step_plus_sym2x2!(A[k, :, :], B, B, τ̃) |
|
236 mldivide_step_plus_sym2x2!(d[k, :], B, c, τ̃) |
|
237 |
|
238 for j=(k-1):-1:1 |
|
239 B = hs[j].M₀+hs[j+1].M₀*([1 0; 0 1]-A[j+1, :, :]) |
|
240 c = hs[j+1].M₀*d[j+1, :]+u[j, :]./τ̃-hs[j].z+hs[j+1].z |
|
241 mldivide_step_plus_sym2x2!(A[j, :, :], B, hs[j].M₀, τ̃) |
|
242 mldivide_step_plus_sym2x2!(d[j, :], B, c, τ̃) |
|
243 end |
|
244 |
|
245 u[1, :] .= d[1, :] |
|
246 for j=2:k |
|
247 u[j, :] .= A[j, :, :]*u[j-1, :] + d[j, :] |
|
248 end |
|
249 end |
|
250 |
|
251 @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old |
|
252 ∇₂!(Δy, x̄) # dual step: y |
|
253 @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | |
|
254 proj_norm₂₁ball!(y, α) # | prox |
|
255 |
|
256 ######################################################## |
|
257 # Give function value and cumulative movement if needed |
|
258 ######################################################## |
|
259 v = verbose() do |
|
260 ∇₂!(Δy, x) |
|
261 hs_plus_reg=0 |
|
262 for j=1:k |
|
263 v=(j==1 ? u[j, :] : u[j, :]-u[j-1, :]) |
|
264 hs_plus_reg += hs[j].cv/2 + dot(hs[j].Mv*v, v)/2+dot(hs[j].av, v) |
|
265 end |
|
266 value = (norm₂²(b-x)/2 + hs_plus_reg/k + α*γnorm₂₁(Δy, ρ)) |
|
267 |
|
268 value, x, u[k, :]+ucumulbase, u[1:k,:].+ucumulbase' |
|
269 end |
|
270 |
|
271 return v |
|
272 end |
|
273 |
|
274 return x, y, v |
|
275 end |
|
276 |
|
277 end # Module |
|
278 |
|
279 |