|
1 #################################################################### |
|
2 # Predictive online PDPS for optical flow with known velocity field |
|
3 #################################################################### |
|
4 |
|
5 __precompile__() |
|
6 |
|
7 module AlgorithmFBDual |
|
8 |
|
9 using Printf |
|
10 |
|
11 using AlgTools.Util |
|
12 import AlgTools.Iterate |
|
13 using ImageTools.Gradient |
|
14 |
|
15 using ..OpticalFlow: Image, |
|
16 ImageSize, |
|
17 flow! |
|
18 |
|
19 ######################### |
|
20 # Iterate initialisation |
|
21 ######################### |
|
22 |
|
23 function init_rest(x::Image) |
|
24 imdim=size(x) |
|
25 |
|
26 y = zeros(2, imdim...) |
|
27 Δx = copy(x) |
|
28 Δy = copy(y) |
|
29 |
|
30 return x, y, Δx, Δy |
|
31 end |
|
32 |
|
33 function init_iterates(xinit::Image) |
|
34 return init_rest(copy(xinit)) |
|
35 end |
|
36 |
|
37 function init_iterates(dim::ImageSize) |
|
38 return init_rest(zeros(dim...)) |
|
39 end |
|
40 |
|
41 ############ |
|
42 # Algorithm |
|
43 ############ |
|
44 |
|
45 function solve( :: Type{DisplacementT}; |
|
46 dim :: ImageSize, |
|
47 iterate = AlgTools.simple_iterate, |
|
48 params::NamedTuple) where DisplacementT |
|
49 |
|
50 ################################ |
|
51 # Extract and set up parameters |
|
52 ################################ |
|
53 |
|
54 α, ρ = params.α, params.ρ |
|
55 τ₀, τ̃₀ = params.τ₀, params.τ̃₀ |
|
56 |
|
57 R_K² = ∇₂_norm₂₂_est² |
|
58 τ = τ₀/R_K² |
|
59 |
|
60 ###################### |
|
61 # Initialise iterates |
|
62 ###################### |
|
63 |
|
64 x, y, Δx, Δy = init_iterates(dim) |
|
65 |
|
66 #################### |
|
67 # Run the algorithm |
|
68 #################### |
|
69 |
|
70 v = iterate(params) do verbose :: Function, |
|
71 b :: Image, |
|
72 v_known :: DisplacementT, |
|
73 🚫unused_b_next :: Image |
|
74 |
|
75 ################## |
|
76 # Prediction step |
|
77 ################## |
|
78 |
|
79 # Δx is a temporary storage variable of correct dimensions |
|
80 flow!(@view(y[1,:, :]), Δx, v_known) |
|
81 flow!(@view(y[2,:, :]), Δx, v_known) |
|
82 |
|
83 ∇₂ᵀ!(Δx, y) |
|
84 @. x = b - Δx |
|
85 ∇₂!(Δy, x) |
|
86 @. y = (y - τ*Δy)/(1 + τ*ρ/α) |
|
87 proj_norm₂₁ball!(y, α) |
|
88 |
|
89 ################################ |
|
90 # Give function value if needed |
|
91 ################################ |
|
92 v = verbose() do |
|
93 ∇₂ᵀ!(Δx, y) |
|
94 @. x = b - Δx |
|
95 ∇₂!(Δy, x) |
|
96 value = norm₂²(b-x)/2 + α*γnorm₂₁(Δy, ρ) |
|
97 value, x, [NaN, NaN], nothing |
|
98 end |
|
99 |
|
100 v |
|
101 end |
|
102 |
|
103 @warn "Using old x value. Better approach unimplemented as this algorithm is not good." |
|
104 # ∇₂ᵀ!(Δx, y) |
|
105 # @. x = b - Δx |
|
106 |
|
107 return x, y, v |
|
108 end |
|
109 |
|
110 end # Module |
|
111 |
|
112 |