Mercurial > repos > PredictPDPS / changeset

--- a/src/AlgorithmDualScaling.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/AlgorithmDualScaling.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -16,7 +16,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     pdflow!
+                     pdflow!,
+                     DualScaling
 
 #########################
 # Iterate initialisation
@@ -90,7 +91,7 @@
             init_data = false
         end
 
-        pdflow!(x, y, v_known)
+        pdflow!(x, Δx, y, Δy, v_known, DualScaling())
 
         ############
         # PDPS step

--- a/src/AlgorithmGreedy.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/AlgorithmGreedy.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -16,7 +16,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     pdflow!
+                     pdflow!,
+                     Greedy
 
 #########################
 # Iterate initialisation
@@ -90,7 +91,7 @@
             init_data = false
         end
 
-        pdflow!(x, Δx, y, Δy, v_known)
+        pdflow!(x, Δx, y, Δy, v_known, Greedy())
 
         ############
         # PDPS step

--- a/src/AlgorithmNoPrediction.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/AlgorithmNoPrediction.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -90,7 +90,7 @@
             init_data = false
         end
 
-        
+        # No prediction
 
         ############
         # PDPS step

--- a/src/AlgorithmRotation.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/AlgorithmRotation.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -16,7 +16,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     pdflow!
+                     pdflow!,
+                     Rotation
 
 #########################
 # Iterate initialisation
@@ -90,7 +91,7 @@
             init_data = false
         end
 
-        pdflow!(x, Δx, y, v_known)
+        pdflow!(x, Δx, y, Δy, v_known, Rotation())
 
         ############
         # PDPS step

--- a/src/OpticalFlow.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/OpticalFlow.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -41,7 +41,8 @@
        DisplacementFull, DisplacementConstant,
        HornSchunckData,
        filter_hs,
-       petpdflow!
+       petpdflow!,
+       DualScaling, Greedy, Rotation
 
 ###############################################
 # Types (several imported from ImageTools.Translate)
@@ -55,6 +56,14 @@
 Gradient = Array{Float64,3}
 ImageSize = Tuple{Int64,Int64}
 
+
+#################################
+# Struct for flow
+#################################
+struct DualScaling end
+struct Greedy end
+struct Rotation end
+
 #################################
 # Displacement field based flow
 #################################
@@ -104,7 +113,7 @@
 end
 
 # Additional method for Greedy
-function pdflow!(x, Δx, y, Δy, u; threads=:none)
+function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none)
     @assert(size(u)==(2,))  
     Δx .= x
     Δy .= y
@@ -120,17 +129,15 @@
 end
 
 # Additional method for Rotation
-function pdflow!(x, Δx, y, u; threads=:none) 
+function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) 
     @assert(size(u)==(2,))
     Δx .= x
     flow!(x, u; threads=(threads==:inner))
-
     (m,n) = size(x)
     dx = similar(y)
     dx_banana = similar(y)
     ∇₂!(dx, Δx)
     ∇₂!(dx_banana, x)
-
     for i=1:m
         for j=1:n
             ndx = @views sum(dx[:, i, j].^2)
@@ -151,7 +158,7 @@
 end
 
 # Additional method for Dual Scaling
-function pdflow!(x, y, u; threads=:none)
+function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none)
     @assert(size(u)==(2,))
     oldx = copy(x)
     flow!(x, u; threads=(threads==:inner))
@@ -193,27 +200,25 @@
 end
 
 # Method for greedy predictor
-function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none)
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Greedy; threads=:none)
     oldx = copy(x)
     center_point = center(x) .+ u
     tform = recenter(RotMatrix(theta_known[1]), center_point)
     Δx = warp(x, tform, axes(x), fillvalue=Flat())
     @. x = Δx
     @. Δy = y
-    if dual_flow
-        Dxx = copy(Δy)
-        DΔx = copy(Δy)
-        ∇₂!(Dxx, x)
-        ∇₂!(DΔx, oldx)
-        inds = abs.(Dxx) .≤ β
-        Dxx[inds] .= 1
-        DΔx[inds] .= 1
-        y .= y.* DΔx ./ Dxx            
-    end
+    Dxx = copy(Δy)
+    DΔx = copy(Δy)
+    ∇₂!(Dxx, x)
+    ∇₂!(DΔx, oldx)
+    inds = abs.(Dxx) .≤ 1e-2
+    Dxx[inds] .= 1
+    DΔx[inds] .= 1
+    y .= y.* DΔx ./ Dxx            
 end
 
-# Method for affine predictor
-function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none)
+# Method for dual scaling predictor
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: DualScaling; threads=:none)
     oldx = copy(x)
     center_point = center(x) .+ u
     tform = recenter(RotMatrix(theta_known[1]), center_point)
@@ -221,26 +226,20 @@
     @. x = Δx
     C = similar(y)
     cc = abs.(x-oldx)
-    if dual_flow
-        cm = max(1e-12,maximum(cc))
-        c = 1 .* (1 .- cc./ cm) .^(10)
-        C[1,:,:] .= c
-        C[2,:,:] .= c
-        y .= C.*y 
-    end
+    cm = max(1e-12,maximum(cc))
+    c = 1 .* (1 .- cc./ cm) .^(10)
+    C[1,:,:] .= c
+    C[2,:,:] .= c
+    y .= C.*y
 end
 
 # Method for rotation prediction (exploiting property of inverse rotation)
-function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β₁, β₂; threads=:none)
-    if dual_flow
-        @backgroundif (threads==:outer) begin
-            petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
-        end begin
-            petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner))
-            petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner))
-        end
-    else
-        petflow!(x, Δx, u, theta_known)
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Rotation; threads=:none)
+    @backgroundif (threads==:outer) begin
+        petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
+    end begin
+        petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner))
+        petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner))
     end
 end
 

--- a/src/PET/AlgorithmDualScaling.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/PET/AlgorithmDualScaling.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -22,7 +22,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     petpdflow!
+                     petpdflow!,
+                     DualScaling
 
 #########################
 # Iterate initialisation
@@ -124,7 +125,7 @@
         # Prediction steps
         ###################
     
-        petpdflow!(x, Δx, y, v_known, theta_known, true)  
+        petpdflow!(x, Δx, y, Δy, v_known, theta_known, DualScaling())  
         
         if params.L_experiment
             @. oldx = x

--- a/src/PET/AlgorithmGreedy.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/PET/AlgorithmGreedy.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -22,7 +22,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     petpdflow!
+                     petpdflow!,
+                     Greedy
 
 #########################
 # Iterate initialisation
@@ -124,7 +125,7 @@
         # Prediction steps
         ###################
     
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, true, 1e-2)   
+        petpdflow!(x, Δx, y, Δy, v_known, theta_known, Greedy())   
         
         if params.L_experiment
             @. oldx = x

--- a/src/PET/AlgorithmProximal.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/PET/AlgorithmProximal.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -151,13 +151,14 @@
             @. oldx = x
         end
 
-        if params.prox_predict
-            ∇₂!(Δy, x)
-            @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
-            #@. cc = y + 1000000*σ̃*Δy 
-            #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc
-            proj_norm₂₁ball!(y, α) 
-        end
+        ##############################    
+        # Proximal step of prediction
+        ##############################
+        ∇₂!(Δy, x)
+        @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
+        #@. cc = y + 1000000*σ̃*Δy 
+        #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc
+        proj_norm₂₁ball!(y, α) 
 
         ############
         # PDPS step

--- a/src/PET/AlgorithmRotation.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/PET/AlgorithmRotation.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -22,7 +22,8 @@
 
 using ..OpticalFlow: ImageSize,
                      Image,
-                     petpdflow!
+                     petpdflow!,
+                     Rotation
 
 #########################
 # Iterate initialisation
@@ -124,7 +125,7 @@
         # Prediction steps
         ###################
     
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow, 0,0)   
+        petpdflow!(x, Δx, y, Δy, v_known, theta_known, Rotation())   
         
         if params.L_experiment
             @. oldx = x

--- a/src/PET/OpticalFlow.jl	Sun Apr 21 22:35:03 2024 +0300
+++ b/src/PET/OpticalFlow.jl	Mon Apr 22 09:26:19 2024 +0300
@@ -41,7 +41,8 @@
        DisplacementFull, DisplacementConstant,
        HornSchunckData,
        filter_hs,
-       petpdflow!
+       petpdflow!,
+       DualScaling, Greedy, Rotation
 
 ###############################################
 # Types (several imported from ImageTools.Translate)
@@ -55,6 +56,14 @@
 Gradient = Array{Float64,3}
 ImageSize = Tuple{Int64,Int64}
 
+
+#################################
+# Struct for flow
+#################################
+struct DualScaling end
+struct Greedy end
+struct Rotation end
+
 #################################
 # Displacement field based flow
 #################################
@@ -104,7 +113,7 @@
 end
 
 # Additional method for Greedy
-function pdflow!(x, Δx, y, Δy, u; threads=:none)
+function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none)
     @assert(size(u)==(2,))  
     Δx .= x
     Δy .= y
@@ -120,17 +129,15 @@
 end
 
 # Additional method for Rotation
-function pdflow!(x, Δx, y, u; threads=:none) 
+function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) 
     @assert(size(u)==(2,))
     Δx .= x
     flow!(x, u; threads=(threads==:inner))
-
     (m,n) = size(x)
     dx = similar(y)
     dx_banana = similar(y)
     ∇₂!(dx, Δx)
     ∇₂!(dx_banana, x)
-
     for i=1:m
         for j=1:n
             ndx = @views sum(dx[:, i, j].^2)
@@ -151,7 +158,7 @@
 end
 
 # Additional method for Dual Scaling
-function pdflow!(x, y, u; threads=:none)
+function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none)
     @assert(size(u)==(2,))
     oldx = copy(x)
     flow!(x, u; threads=(threads==:inner))
@@ -193,27 +200,25 @@
 end
 
 # Method for greedy predictor
-function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none)
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Greedy; threads=:none)
     oldx = copy(x)
     center_point = center(x) .+ u
     tform = recenter(RotMatrix(theta_known[1]), center_point)
     Δx = warp(x, tform, axes(x), fillvalue=Flat())
     @. x = Δx
     @. Δy = y
-    if dual_flow
-        Dxx = copy(Δy)
-        DΔx = copy(Δy)
-        ∇₂!(Dxx, x)
-        ∇₂!(DΔx, oldx)
-        inds = abs.(Dxx) .≤ β
-        Dxx[inds] .= 1
-        DΔx[inds] .= 1
-        y .= y.* DΔx ./ Dxx            
-    end
+    Dxx = copy(Δy)
+    DΔx = copy(Δy)
+    ∇₂!(Dxx, x)
+    ∇₂!(DΔx, oldx)
+    inds = abs.(Dxx) .≤ 1e-2
+    Dxx[inds] .= 1
+    DΔx[inds] .= 1
+    y .= y.* DΔx ./ Dxx            
 end
 
-# Method for affine predictor
-function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none)
+# Method for dual scaling predictor
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: DualScaling; threads=:none)
     oldx = copy(x)
     center_point = center(x) .+ u
     tform = recenter(RotMatrix(theta_known[1]), center_point)
@@ -221,26 +226,20 @@
     @. x = Δx
     C = similar(y)
     cc = abs.(x-oldx)
-    if dual_flow
-        cm = max(1e-12,maximum(cc))
-        c = 1 .* (1 .- cc./ cm) .^(10)
-        C[1,:,:] .= c
-        C[2,:,:] .= c
-        y .= C.*y 
-    end
+    cm = max(1e-12,maximum(cc))
+    c = 1 .* (1 .- cc./ cm) .^(10)
+    C[1,:,:] .= c
+    C[2,:,:] .= c
+    y .= C.*y
 end
 
 # Method for rotation prediction (exploiting property of inverse rotation)
-function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β₁, β₂; threads=:none)
-    if dual_flow
-        @backgroundif (threads==:outer) begin
-            petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
-        end begin
-            petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner))
-            petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner))
-        end
-    else
-        petflow!(x, Δx, u, theta_known)
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Rotation; threads=:none)
+    @backgroundif (threads==:outer) begin
+        petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
+    end begin
+        petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner))
+        petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner))
     end
 end