src/PET/OpticalFlow.jl

Sun, 21 Apr 2024 13:43:18 +0300

author
Neil Dizon <neil.dizon@helsinki.fi>
date
Sun, 21 Apr 2024 13:43:18 +0300
changeset 16
98b79c837a30
parent 8
e4ad8f7ce671
child 26
ccd22bbbb02f
permissions
-rw-r--r--

added zero dual in PET

8
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
1 ################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
2 # Code relevant to optical flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
3 ################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
4
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
5 __precompile__()
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
6
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
7 module OpticalFlow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
8
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
9 using AlgTools.Util
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
10 using ImageTools.Gradient
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
11 import ImageTools.Translate
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
12 using ImageTools.ImFilter
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
13
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
14 # using ImageTransformations
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
15 # using Images, CoordinateTransformations, Rotations, OffsetArrays
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
16 # using Interpolations
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
17
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
18 import Images: center, warp
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
19 import CoordinateTransformations: recenter
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
20 import Rotations: RotMatrix
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
21 import Interpolations: Flat
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
22
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
23 ##########
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
24 # Exports
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
25 ##########
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
26
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
27 export flow!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
28 pdflow!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
29 flow_grad!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
30 flow_interp!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
31 estimate_Λ²,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
32 estimate_linear_Λ²,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
33 pointwise_gradiprod_2d!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
34 pointwise_gradiprod_2dᵀ!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
35 horn_schunck_reg_prox!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
36 horn_schunck_reg_prox_op!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
37 mldivide_step_plus_sym2x2!,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
38 linearised_optical_flow_error,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
39 Image, AbstractImage, ImageSize,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
40 Gradient, Displacement,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
41 DisplacementFull, DisplacementConstant,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
42 HornSchunckData,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
43 filter_hs,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
44 petpdflow!
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
45
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
46 ###############################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
47 # Types (several imported from ImageTools.Translate)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
48 ###############################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
49
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
50 Image = Translate.Image
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
51 AbstractImage = AbstractArray{Float64,2}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
52 Displacement = Translate.Displacement
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
53 DisplacementFull = Translate.DisplacementFull
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
54 DisplacementConstant = Translate.DisplacementConstant
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
55 Gradient = Array{Float64,3}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
56 ImageSize = Tuple{Int64,Int64}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
57
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
58 #################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
59 # Displacement field based flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
60 #################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
61
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
62 function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage;
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
63 threads = false)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
64 tmp .= x
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
65 Translate.translate_image!(x, tmp, u; threads=threads)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
66 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
67
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
68 function flow_interp!(x::AbstractImage, u::Displacement;
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
69 threads = false)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
70 tmp = copy(x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
71 Translate.translate_image!(x, tmp, u; threads=threads)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
72 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
73
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
74 flow! = flow_interp!
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
75
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
76 function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
77 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
78 #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
79 # (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :])))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
80 @backgroundif (threads==:outer) begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
81 flow!(x, u, Δx; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
82 end begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
83 flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
84 flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
85 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
86 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
87 flow!(x, u, Δx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
88 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
89 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
90
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
91 function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
92 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
93 @backgroundif (threads==:outer) begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
94 flow!(x, u, Δx; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
95 flow!(z, u, Δz; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
96 end begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
97 flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
98 flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
99 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
100 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
101 flow!(x, u, Δx; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
102 flow!(z, u, Δz; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
103 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
104 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
105
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
106 # Additional method for Greedy
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
107 function pdflow!(x, Δx, y, Δy, u; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
108 @assert(size(u)==(2,))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
109 Δx .= x
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
110 Δy .= y
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
111 flow!(x, u; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
112 Dxx = similar(Δy)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
113 DΔx = similar(Δy)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
114 ∇₂!(Dxx, x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
115 ∇₂!(DΔx, Δx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
116 inds = abs.(Dxx) .≤ 1e-1
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
117 Dxx[inds] .= 1
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
118 DΔx[inds] .= 1
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
119 y .= y.* DΔx ./ Dxx
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
120 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
121
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
122 # Additional method for Rotation
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
123 function pdflow!(x, Δx, y, u; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
124 @assert(size(u)==(2,))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
125 Δx .= x
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
126 flow!(x, u; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
127
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
128 (m,n) = size(x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
129 dx = similar(y)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
130 dx_banana = similar(y)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
131 ∇₂!(dx, Δx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
132 ∇₂!(dx_banana, x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
133
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
134 for i=1:m
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
135 for j=1:n
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
136 ndx = @views sum(dx[:, i, j].^2)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
137 ndx_banana = @views sum(dx_banana[:, i, j].^2)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
138 if ndx > 1e-4 && ndx_banana > 1e-4
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
139 A = dx[:, i, j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
140 B = dx_banana[:, i, j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
141 theta = atan(B[1] * A[2] - B[2] * A[1], B[1] * A[1] + B[2] * A[2]) # Oriented angle from A to B
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
142 cos_theta = cos(theta)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
143 sin_theta = sin(theta)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
144 a = cos_theta * y[1, i, j] - sin_theta * y[2, i, j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
145 b = sin_theta * y[1, i, j] + cos_theta * y[2, i, j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
146 y[1, i, j] = a
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
147 y[2, i, j] = b
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
148 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
149 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
150 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
151 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
152
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
153 # Additional method for Dual Scaling
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
154 function pdflow!(x, y, u; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
155 @assert(size(u)==(2,))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
156 oldx = copy(x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
157 flow!(x, u; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
158 C = similar(y)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
159 cc = abs.(x-oldx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
160 cm = max(1e-12,maximum(cc))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
161 c = 1 .* (1 .- cc./ cm) .^(10)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
162 C[1,:,:] .= c
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
163 C[2,:,:] .= c
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
164 y .= C.*y
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
165 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
166
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
167
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
168 ##########################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
169 # PET
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
170 ##########################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
171 function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant;
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
172 threads = false)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
173 tmp .= x
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
174 center_point = center(x) .+ u
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
175 tform = recenter(RotMatrix(theta_known[1]), center_point)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
176 tmp = warp(x, tform, axes(x), fillvalue=Flat())
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
177 x .= tmp
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
178 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
179
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
180 petflow! = petflow_interp!
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
181
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
182 function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
183 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
184 @backgroundif (threads==:outer) begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
185 petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
186 end begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
187 petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
188 petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
189 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
190 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
191 petflow!(x, Δx, u, theta_known)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
192 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
193 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
194
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
195 # Method for greedy predictor
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
196 function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
197 oldx = copy(x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
198 center_point = center(x) .+ u
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
199 tform = recenter(RotMatrix(theta_known[1]), center_point)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
200 Δx = warp(x, tform, axes(x), fillvalue=Flat())
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
201 @. x = Δx
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
202 @. Δy = y
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
203 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
204 Dxx = copy(Δy)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
205 DΔx = copy(Δy)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
206 ∇₂!(Dxx, x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
207 ∇₂!(DΔx, oldx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
208 inds = abs.(Dxx) .≤ β
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
209 Dxx[inds] .= 1
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
210 DΔx[inds] .= 1
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
211 y .= y.* DΔx ./ Dxx
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
212 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
213 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
214
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
215 # Method for affine predictor
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
216 function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
217 oldx = copy(x)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
218 center_point = center(x) .+ u
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
219 tform = recenter(RotMatrix(theta_known[1]), center_point)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
220 Δx = warp(x, tform, axes(x), fillvalue=Flat())
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
221 @. x = Δx
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
222 C = similar(y)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
223 cc = abs.(x-oldx)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
224 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
225 cm = max(1e-12,maximum(cc))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
226 c = 1 .* (1 .- cc./ cm) .^(10)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
227 C[1,:,:] .= c
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
228 C[2,:,:] .= c
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
229 y .= C.*y
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
230 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
231 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
232
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
233 # Method for rotation prediction (exploiting property of inverse rotation)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
234 function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β₁, β₂; threads=:none)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
235 if dual_flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
236 @backgroundif (threads==:outer) begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
237 petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
238 end begin
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
239 petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
240 petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
241 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
242 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
243 petflow!(x, Δx, u, theta_known)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
244 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
245 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
246
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
247 ##########################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
248 # Linearised optical flow
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
249 ##########################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
250
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
251 # ⟨⟨u, ∇b⟩⟩
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
252 function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
253 u::DisplacementFull, b::Image;
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
254 add = false)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
255 ∇₂c!(vtmp, b)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
256
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
257 u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
258 vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
259 y′=reshape(y, prod(size(y)))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
260
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
261 if add
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
262 @simd for i = 1:length(y′)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
263 @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i]))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
264 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
265 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
266 @simd for i = 1:length(y′)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
267 @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i]))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
268 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
269 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
270 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
271
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
272 function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
273 u::DisplacementConstant, b::Image;
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
274 add = false)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
275 ∇₂c!(vtmp, b)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
276
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
277 vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
278 y′=reshape(y, prod(size(y)))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
279
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
280 if add
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
281 @simd for i = 1:length(y′)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
282 @inbounds y′[i] += dot(u, @view(vtmp′[:, i]))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
283 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
284 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
285 @simd for i = 1:length(y′)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
286 @inbounds y′[i] = dot(u, @view(vtmp′[:, i]))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
287 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
288 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
289 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
290
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
291 # ∇b ⋅ y
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
292 function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
293 ∇₂c!(u, b)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
294
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
295 u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
296 y′=reshape(y, prod(size(y)))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
297
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
298 @simd for i=1:length(y′)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
299 @inbounds @. u′[:, i] *= y′[i]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
300 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
301 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
302
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
303 function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
304 @assert(size(y)==size(b) && size(u)==(2,))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
305 u .= 0
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
306 ∇₂cfold!(b, nothing) do g, st, (i, j)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
307 @inbounds u .+= g.*y[i, j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
308 return st
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
309 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
310 # Reweight to be with respect to 𝟙^*𝟙 inner product.
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
311 u ./= prod(size(b))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
312 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
313
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
314 mutable struct ConstantDisplacementHornSchunckData
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
315 M₀::Array{Float64,2}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
316 z::Array{Float64,1}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
317 Mv::Array{Float64,2}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
318 av::Array{Float64,1}
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
319 cv::Float64
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
320
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
321 function ConstantDisplacementHornSchunckData()
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
322 return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
323 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
324 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
325
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
326 # For DisplacementConstant, for the simple prox step
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
327 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
328 # (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
329 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
330 # construct matrix M₀ and vector z such that we can solve u from
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
331 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
332 # (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
333 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
334 # Note that the problem
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
335 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
336 # argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
337 # + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
338 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
339 # has with respect to u the system
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
340 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
341 # (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
342 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
343 # where the primed variables correspond to (2) for (1) for uʹ in place of u:
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
344 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
345 # argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
346 #
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
347 function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
348 bnext::Image, b::Image, θ, λ, T)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
349 @assert(size(b)==size(bnext))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
350 w = prod(size(b))
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
351 z = hs.z
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
352 cv = 0
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
353 # Factors of symmetric matrix [a c; c d]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
354 a, c, d = 0.0, 0.0, 0.0
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
355 # This used to use ∇₂cfold but it is faster to allocate temporary
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
356 # storage for the full gradient due to probably better memory and SIMD
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
357 # instruction usage.
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
358 g = zeros(2, size(b)...)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
359 ∇₂c!(g, b)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
360 @inbounds for i=1:size(b, 1)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
361 for j=1:size(b, 2)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
362 δ = bnext[i,j]-b[i,j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
363 @. z += g[:,i,j]*δ
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
364 cv += δ*δ
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
365 a += g[1,i,j]*g[1,i,j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
366 c += g[1,i,j]*g[2,i,j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
367 d += g[2,i,j]*g[2,i,j]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
368 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
369 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
370 w₀ = λ
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
371 w₂ = θ/w
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
372 aʹ = w₀ + w₂*a
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
373 cʹ = w₂*c
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
374 dʹ = w₀ + w₂*d
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
375 hs.M₀ .= [aʹ cʹ; cʹ dʹ]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
376 hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
377 hs.cv = cv*θ
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
378 hs.av .= hs.z.*θ
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
379 hs.z .*= w₂/T
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
380 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
381
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
382 # Solve the 2D system (I/τ+M₀)u = z
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
383 @inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
384 a = 1/τ+M₀[1, 1]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
385 c = M₀[1, 2]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
386 d = 1/τ+M₀[2, 2]
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
387 u .= ([d -c; -c a]*z)./(a*d-c*c)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
388 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
389
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
390 function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image,
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
391 θ, λ, T, τ)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
392 hs=ConstantDisplacementHornSchunckData()
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
393 horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
394 mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
395 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
396
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
397 function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
398 if !isnothing(δ)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
399 u = δ.*u
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
400 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
401 pointwise_gradiprod_2d!(x, vtmp, u, x; add=true)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
402 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
403
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
404 # Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
405 function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
406 imdim = size(b)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
407 vtmp = zeros(2, imdim...)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
408 tmp = b-b_prev
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
409 pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
410 return tmp
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
411 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
412
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
413 ##############################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
414 # Helper to smooth data for Horn–Schunck term
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
415 ##############################################
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
416
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
417 function filter_hs(b, b_next, b_next_filt, kernel)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
418 if kernel==nothing
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
419 f = x -> x
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
420 else
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
421 f = x -> simple_imfilter(x, kernel; threads=true)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
422 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
423
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
424 # We already filtered b in the previous step (b_next in that step)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
425 b_filt = b_next_filt==nothing ? f(b) : b_next_filt
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
426 b_next_filt = f(b_next)
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
427
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
428 return b_filt, b_next_filt
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
429 end
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
430
e4ad8f7ce671 Added PET and updated README
Neil Dizon <neil.dizon@helsinki.fi>
parents:
diff changeset
431 end # Module

mercurial