src/AlgorithmBothCumul.jl

changeset 0
a55e35d20336
equal deleted inserted replaced
-1:000000000000 0:a55e35d20336
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

mercurial