# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1714068340 18000
# Node ID e4a8f662a1ac1de914be6e1760524c12aad2f74a
# Parent  74b1a9f0c35ea7d4c99eea263399c3a4f6fa471e
Reduce code duplication.

diff -r 74b1a9f0c35e -r e4a8f662a1ac src/Algorithm.jl
--- a/src/Algorithm.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,148 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module Algorithm
-
-identifier = "pdps_known"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function step_lengths(params, γ, R_K²)
-    ρ̃₀, τ₀, σ₀, σ̃₀ =  params.ρ̃₀, params.τ₀, params.σ₀, params.σ̃₀
-    δ = params.δ
-    ρ = isdefined(params, :phantom_ρ) ? params.phantom_ρ : params.ρ
-    Λ = params.Λ
-    Θ = params.dual_flow ? Λ : 1
-
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*min(1/(τ*R_K²), 1/max(0, τ*R_K²/((1+γ*τ-Λ)*(1-δ))-ρ))
-    q = δ*(1+σ*ρ)/Θ
-    if 1 ≥ q
-        σ̃ = σ̃₀*σ/q
-        #ρ̃ = ρ̃₀*max(0, ((Θ*σ)/(2*δ*σ̃^2*(1+σ*ρ))+1/(2σ)-1/σ̃))
-        ρ̃ = max(0, (1-q)/(2*σ))
-    else
-        σ̃ = σ̃₀*σ/(q*(1-√(1-1/q)))
-        ρ̃ = 0
-    end
-    
-    println("Step length parameters: τ=$(τ), σ=$(σ), σ̃=$(σ̃), ρ̃=$(ρ̃)")
-
-    return τ, σ, σ̃, ρ̃
-end
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²)
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, params.dual_flow)
-
-        if params.prox_predict
-            ∇₂!(Δy, x)
-            @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
-            proj_norm₂₁ball!(y, α) 
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmBothMulti.jl
--- a/src/AlgorithmBothMulti.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/AlgorithmBothMulti.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -26,7 +26,7 @@
                      ConstantDisplacementHornSchunckData,
                      filter_hs
 
-using ..Algorithm: step_lengths
+using ..AlgorithmProximal: step_lengths
 
 #########################
 # Iterate initialisation
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmBothMulti.jl.orig
--- a/src/AlgorithmBothMulti.jl.orig	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,277 +0,0 @@
-######################################################################
-# Predictive online PDPS for optical flow with unknown velocity field
-######################################################################
-
-__precompile__()
-
-module AlgorithmBothMulti
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.ImFilter
-
-using ..OpticalFlow: Image,
-                     ImageSize,
-                     DisplacementConstant,
-                     pdflow!,
-                     horn_schunck_reg_prox!,
-                     pointwise_gradiprod_2d!,
-                     horn_schunck_reg_prox_op!,
-                     mldivide_step_plus_sym2x2!,
-                     ConstantDisplacementHornSchunckData,
-                     filter_hs
-
-using ..Algorithm: step_lengths
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_displ(xinit::Image, ::Type{DisplacementConstant}, n::Integer)
-    return xinit, zeros(n, 2)
-end
-
-# function init_displ(xinit::Image, ::Type{DisplacementFull}, n::Integer)
-#     return xinit, zeros(n, 2, size(xinit)...)
-# end
-
-function init_rest(x::Image, u::DisplacementT) where DisplacementT
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄, u
-end
-
-function init_iterates( :: Type{DisplacementT},
-                       xinit::Image,
-                       n::Integer) where DisplacementT    
-    return init_rest(init_displ(copy(xinit), DisplacementT, n)...)
-end
-
-function init_iterates( :: Type{DisplacementT},
-                       dim::ImageSize,
-                       n::Integer) where DisplacementT
-    return init_rest(init_displ(zeros(dim...), DisplacementT, n)...)
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate, 
-               params::NamedTuple) where DisplacementT
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    n = params.displacement_count
-    k = 0 # number of displacements we have already
-
-    x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim, n)
-    init_data = (params.init == :data)
-    hs = [ConstantDisplacementHornSchunckData() for i=1:n]
-    #hs = Array{ConstantDisplacementHornSchunckData}(undef, n)
-    A = Array{Float64,3}(undef, n, 2, 2)
-    d = Array{Float64,2}(undef, n, 2)
-
-    # … for tracking cumulative movement
-    ucumulbase = [0.0, 0.0]
-    
-    #############################################
-    # Extract parameters and set up step lengths
-    #############################################
-
-    α, ρ, λ, θ = params.α, params.ρ, params.λ, params.θ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²)
-
-    kernel = params.kernel
-    T = params.timestep 
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    b_next_filt = nothing
-    diffu = similar(u[1, :])
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           🚫unused_v_known :: DisplacementT,
-                           b_next :: Image
-
-        ####################################
-        # Smooth data for Horn–Schunck term
-        ####################################
-
-        b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel)
-
-        ################################################
-        # Prediction step
-        ################################################
-
-        # Predict x and y
-        if k==0
-            if init_data
-                x .= b
-                init_data = false
-            end 
-        else
-            # Displacement from previous to this image is estimated as
-            # the difference of their displacements from beginning of window.
-            if k>1
-                @. @views diffu = u[k, :] - u[k-1, :]
-            else
-                @. @views diffu = u[k, :]
-            end
-
-            pdflow!(x, Δx, y, Δy, diffu, params.dual_flow)
-        end
-
-        # Shift stored prox matrices
-        if k==n
-            tmp = copy(u[1, :])
-            ucumulbase .+= tmp
-            for j=1:(n-1)
-                @. @views u[j, :] = u[j+1, :] - tmp
-                hs[j] = hs[j+1]
-            end
-            # Create new struct as original contains references to objects that
-            # have been moved to index n-1.
-            hs[n]=ConstantDisplacementHornSchunckData()
-        else
-            k += 1
-        end
-
-        # Predict u: zero displacement from current to next image, i.e.,
-        # same displacement to beginning of window.
-        if k==1        
-            @. @views u[k, :] = 0.0
-        else
-            @. @views u[k, :] = u[k-1, :]
-        end
-
-        # Predictor proximals tep
-        if params.prox_predict
-            ∇₂!(Δy, x) 
-            @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
-            proj_norm₂₁ball!(y, α) 
-        end
-
-        #################################################################################
-        # PDPS step
-        #
-        # For the displacements, with τ̃=τ/k, we need to solve for 2≤j<k,
-        #
-        # (1) (I/τ̃+M₀^j+M₀^{j+1})u^j = M₀^ju^{j-1} + M₀^{j+1}u^{j+1}
-        #                              + ũ^j/τ̃ - z^j + z^{j+1}, 
-        #
-        # as well as
-        #
-        # (2) (I/τ̃+M₀^k)u^k = M₀^k u^{k-1} + ũ^k/τ̃ - z^k
-        #
-        # and
-        #
-        # (3) (I/τ̃+M₀^1+M₀^2)u^1 = 0 + M₀^{2}u^{2} + ũ^1/τ̃ - z^1 + z^{2}
-        #
-        # We first construct from (2) that
-        #
-        #     u^k = A^k u^{k-1} + d^k
-        #
-        # for
-        #
-        #     A^k := (I/τ̃+M₀^k)^{-1} M₀^k
-        #     d_k := (I/τ̃+M₀^k)^{-1} (ũ^k/τ̃ - z^k).
-        #
-        # Inserting this into (1) we need
-        #
-        # (4)  (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))u^j = M₀^ju^{j-1} + M₀^{j+1}d^{j+1}
-        #                                          + ũ^j/τ̃ - z^j + z^{j+1}.
-        # 
-        # This is well-defined because A^{j+1} < I. It also has the same form as (1), so 
-        # we continue with
-        #
-        # (5)   u^j = A^j u^{j-1} + d^j
-        #
-        # for
-        #
-        #      A^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} M₀^j
-        #      d^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1}
-        #               (M₀^{j+1}d^{j+1} + ũ^j/τ̃ - z^j + z^{j+1})
-        # 
-        # Finally from (3) with these we need
-        #
-        #      (I/τ̃+M₀^1+M₀^2(I-A^2))u^1 = M₀^2d^2 + ũ^1/τ̃ - z^1 + z^2,
-        #
-        # which is of the same form as (4) with u^0=0, so by (5) u^1=d^1.
-        #
-        #################################################################################
-
-        ∇₂ᵀ!(Δx, y)                         # primal step:
-        @. x̄ = x                            # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)           # |  prox
-                                            # |  | for displacement
-        # Calculate matrices for latest data; rest is stored. 
-        @views begin
-            horn_schunck_reg_prox_op!(hs[k], b_next_filt, b_filt, θ, λ, T)
-
-            τ̃=τ/k
-
-            B = hs[k].M₀
-            c = u[k, :]./τ̃-hs[k].z
-            mldivide_step_plus_sym2x2!(A[k, :, :], B, B, τ̃)
-            mldivide_step_plus_sym2x2!(d[k, :], B, c, τ̃)
-
-            for j=(k-1):-1:1
-                B = hs[j].M₀+hs[j+1].M₀*([1 0; 0 1]-A[j+1, :, :])
-                c = hs[j+1].M₀*d[j+1, :]+u[j, :]./τ̃-hs[j].z+hs[j+1].z
-                mldivide_step_plus_sym2x2!(A[j, :, :], B, hs[j].M₀, τ̃)
-                mldivide_step_plus_sym2x2!(d[j, :], B, c, τ̃)
-            end
-
-            u[1, :] .= d[1, :]
-            for j=2:k
-                u[j, :] .= A[j, :, :]*u[j-1, :] + d[j, :]
-            end
-        end
-
-        @. x̄ = 2x - x̄                       # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                          # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)       # |
-        proj_norm₂₁ball!(y, α)              # |  prox
-
-        ########################################################
-        # Give function value and cumulative movement if needed
-        ########################################################
-        v = verbose() do
-            ∇₂!(Δy, x)
-            hs_plus_reg=0            
-            for j=1:k
-                v=(j==1 ? u[j, :] : u[j, :]-u[j-1, :])                
-                hs_plus_reg += hs[j].cv/2 + dot(hs[j].Mv*v, v)/2+dot(hs[j].av, v)
-            end
-            value = (norm₂²(b-x)/2 + hs_plus_reg/k + α*γnorm₂₁(Δy, ρ))
-
-            value, x, u[k, :]+ucumulbase, u[1:k,:].+ucumulbase'
-        end
-
-        return v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmDualScaling.jl
--- a/src/AlgorithmDualScaling.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,126 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmDualScaling
-
-identifier = "pdps_known_dualscaling"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!,
-                     DualScaling
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################
-    # Extract and set up parameters
-    ################################
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    ds = DualScaling(params.ds_exponent, params.ds_threshold)
-
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, ds)
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmGreedy.jl
--- a/src/AlgorithmGreedy.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,125 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmGreedy
-
-identifier = "pdps_known_greedy"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!,
-                     Greedy
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, Greedy())
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmNew.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/AlgorithmNew.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -0,0 +1,127 @@
+####################################################################
+# Predictive online PDPS for optical flow with known velocity field
+####################################################################
+
+__precompile__()
+
+module AlgorithmNew
+
+# Default identifier for when none is given by predictor
+identifier = "pdps_known_noprediction"
+
+using Printf
+
+using AlgTools.Util
+import AlgTools.Iterate
+using ImageTools.Gradient
+
+using ..OpticalFlow: ImageSize,
+                     Image,
+                     pdflow!
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_rest(x::Image)
+    imdim=size(x)
+
+    y = zeros(2, imdim...)
+    Δx = copy(x)
+    Δy = copy(y)
+    x̄ = copy(x)
+
+    return x, y, Δx, Δy, x̄
+end
+
+function init_iterates(xinit::Image)
+    return init_rest(copy(xinit))
+end
+
+function init_iterates(dim::ImageSize)
+    return init_rest(zeros(dim...))
+end
+
+############
+# Algorithm
+############
+
+function solve( :: Type{DisplacementT};
+               dim :: ImageSize,
+               iterate = AlgTools.simple_iterate,
+               params::NamedTuple) where DisplacementT
+
+    ################################
+    # Extract and set up parameters
+    ################################
+
+    α, ρ = params.α, params.ρ
+    R_K² = ∇₂_norm₂₂_est²
+    γ = 1.0
+    Λ = params.Λ
+    τ₀, σ₀ = params.τ₀, params.σ₀
+
+    τ = τ₀/γ
+    @assert(1+γ*τ ≥ Λ)
+    σ = σ₀*1/(τ*R_K²)
+
+    println("Step length parameters: τ=$(τ), σ=$(σ)")
+
+    ######################
+    # Initialise iterates
+    ######################
+
+    x, y, Δx, Δy, x̄ = init_iterates(dim)
+    init_data = (params.init == :data)
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    v = iterate(params) do verbose :: Function,
+                           b :: Image,
+                           v_known :: DisplacementT,
+                           🚫unused_b_next :: Image
+
+        ##################
+        # Prediction step
+        ##################
+        if init_data
+            x .= b
+            init_data = false
+        end
+
+        if haskey(params, :predictor) && ~isnothing(params.predictor)
+            pdflow!(x, Δx, y, Δy, v_known, params.predictor)
+        end
+
+        ############
+        # PDPS step
+        ############
+
+        ∇₂ᵀ!(Δx, y)                    # primal step:
+        @. x̄ = x                       # |  save old x for over-relax
+        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
+        @. x̄ = 2x - x̄                  # over-relax
+        ∇₂!(Δy, x̄)                     # dual step: y
+        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
+        proj_norm₂₁ball!(y, α)         # |  prox
+
+        ################################
+        # Give function value if needed
+        ################################
+        v = verbose() do
+            ∇₂!(Δy, x)
+            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
+            value, x, [NaN, NaN], nothing
+        end
+
+        v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmNoPrediction.jl
--- a/src/AlgorithmNoPrediction.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,124 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmNoPrediction
-
-identifier = "pdps_known_noprediction"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        # No prediction
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmPrimalOnly.jl
--- a/src/AlgorithmPrimalOnly.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,124 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmPrimalOnly
-
-identifier = "pdps_known_primalonly"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, false)
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmRotation.jl
--- a/src/AlgorithmRotation.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,125 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmRotation
-
-identifier = "pdps_known_rotation"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!,
-                     Rotation
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, Rotation())
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/AlgorithmZeroDual.jl
--- a/src/AlgorithmZeroDual.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,125 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmZeroDual
-
-identifier = "pdps_known_zerodual"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     pdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-
-    return x, y, Δx, Δy, x̄
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1.0
-    Λ = params.Λ
-    τ₀, σ₀ = params.τ₀, params.σ₀
-
-    τ = τ₀/γ
-    @assert(1+γ*τ ≥ Λ)
-    σ = σ₀*1/(τ*R_K²)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄ = init_iterates(dim)
-    init_data = (params.init == :data)
-
-    ####################
-    # Run the algorithm
-    ####################
-
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,
-                           v_known :: DisplacementT,
-                           🚫unused_b_next :: Image
-
-        ##################
-        # Prediction step
-        ##################
-        if init_data
-            x .= b
-            init_data = false
-        end
-
-        pdflow!(x, Δx, y, Δy, v_known, false)
-        y .= zeros(size(y)...)
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # |  save old x for over-relax
-        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
-        @. x̄ = 2x - x̄                  # over-relax
-        ∇₂!(Δy, x̄)                     # dual step: y
-        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        ################################
-        # Give function value if needed
-        ################################
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
-            value, x, [NaN, NaN], nothing
-        end
-
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/ImGenerate.jl
--- a/src/ImGenerate.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/ImGenerate.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -100,10 +100,10 @@
             θ = 2π*rand(rng,Float64)
             r = params.shake
             return [r*cos(θ), r*sin(θ)]
-        end        
+        end
     else
         error("Unknown shaketype $(params.shaketype)")
-    end 
+    end
 end
 
 pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz)
@@ -181,15 +181,14 @@
                               :: Type{DisplacementConstant},
                               datachannel :: Channel{OnlineData{DisplacementConstant}},
                               params :: NamedTuple)
-    
 
     # Set up counter and zero factor for stabilisation
     @assert(params.maxiter ≥ maximum(params.stable_interval))
     indx = 1
-    zero_factor = indx in params.stable_interval ? 0.0 : 1.0                          
+    zero_factor = indx in params.stable_interval ? 0.0 : 1.0
 
-    # Restart the seed to enable comparison across predictors                       
-    Random.seed!(rng,951508)   # algreadygood
+    # Restart the seed to enable comparison across predictors
+    Random.seed!(rng, if haskey(params, :seed) params.seed else 951508 end)
 
     nextv = shake(params)
     v_true = zero_factor.*nextv()
@@ -228,17 +227,17 @@
 # PETscan
 ########################################################################
 function generate_sinogram(im, sz,
-                        :: Type{DisplacementConstant},
-                        datachannel :: Channel{PetOnlineData{DisplacementConstant}},
-                        params :: NamedTuple)
+                           :: Type{DisplacementConstant},
+                            datachannel :: Channel{PetOnlineData{DisplacementConstant}},
+                            params :: NamedTuple)
 
     # Set up counter and zero factor for stabilisation
     @assert(params.maxiter ≥ maximum(params.stable_interval))
     indx = 1
     zero_factor = indx in params.stable_interval ? 0.0 : 1.0
 
-    # Restart the seed to enable comparison across predictors                       
-    Random.seed!(rng,314159)
+    # Restart the seed to enable comparison across predictors
+    Random.seed!(rng, if haskey(params, :seed) params.seed else 314159 end)
 
     nextv = shake(params)
     v_true = zero_factor.*nextv()
@@ -282,7 +281,7 @@
         # Next theta
         theta_true = zero_factor*rotatebytheta(params)
         theta_cumul += theta_true
-        # Next sinogram mask 
+        # Next sinogram mask
         S_true = generate_radonmask(params)
         # Update indx and zero factor
         indx += 1
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/OpticalFlow.jl
--- a/src/OpticalFlow.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/OpticalFlow.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -42,7 +42,8 @@
        HornSchunckData,
        filter_hs,
        petpdflow!,
-       DualScaling, Greedy, Rotation
+       DualScaling, Greedy, Rotation, ZeroDual, PrimalOnly,
+       identifier
 
 ###############################################
 # Types (several imported from ImageTools.Translate)
@@ -68,6 +69,28 @@
 
 struct Greedy end
 struct Rotation end
+struct ZeroDual end
+struct PrimalOnly end
+
+function identifier(::DualScaling)
+    "pdps_known_dualscaling"
+end
+
+function identifier(::Rotation)
+    "pdps_known_rotation"
+end
+
+function identifier(::Greedy)
+    "pdps_known_greedy"
+end
+
+function identifier(::ZeroDual)
+    "pdps_known_zerodual"
+end
+
+function identifier(::PrimalOnly)
+    "pdps_known_primalonly"
+end
 
 #################################
 # Displacement field based flow
@@ -102,7 +125,7 @@
     end
 end
 
-function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none)
+function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow :: Bool; threads=:none)
     if dual_flow
         @backgroundif (threads==:outer) begin
             flow!(x, u, Δx; threads=(threads==:inner))
@@ -119,7 +142,7 @@
 
 # Additional method for Greedy
 function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none)
-    @assert(size(u)==(2,))  
+    @assert(size(u)==(2,))
     Δx .= x
     Δy .= y
     flow!(x, u; threads=(threads==:inner))
@@ -130,11 +153,11 @@
     inds = abs.(Dxx) .≤ 1e-1
     Dxx[inds] .= 1
     DΔx[inds] .= 1
-    y .= y.* DΔx ./ Dxx  
+    y .= y.* DΔx ./ Dxx
 end
 
 # Additional method for Rotation
-function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; 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))
@@ -173,13 +196,24 @@
     c = 1 .* (1 .- cc./ cm) .^flow.exponent
     C[1,:,:] .= c
     C[2,:,:] .= c
-    y .= C.*y 
+    y .= C.*y
 end
 
+function pdflow!(x, Δx, y, Δy, u, :: ZeroDual; threads=:none)
+    @assert(size(u)==(2,))
+    flow!(x, u; threads=(threads==:inner))
+    y .= 0.0
+end
+
+function pdflow!(x, Δx, y, Δy, u, :: PrimalOnly; threads=:none)
+    @assert(size(u)==(2,))
+    flow!(x, u; threads=(threads==:inner))
+end
 
 ##########################
 # PET
 ##########################
+
 function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant;
     threads = false)
     tmp .= x
@@ -191,7 +225,7 @@
 
 petflow! = petflow_interp!
 
-function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none)
+function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow :: Bool; threads=:none)
     if dual_flow
         @backgroundif (threads==:outer) begin
         petflow!(x, Δx, u, theta_known; threads=(threads==:inner))
@@ -219,7 +253,7 @@
     inds = abs.(Dxx) .≤ 1e-2
     Dxx[inds] .= 1
     DΔx[inds] .= 1
-    y .= y.* DΔx ./ Dxx            
+    y .= y.* DΔx ./ Dxx
 end
 
 # Method for dual scaling predictor
@@ -248,6 +282,15 @@
     end
 end
 
+function petpdflow!(x, Δx, y, Δy, u, theta_known, :: ZeroDual; threads=:none)
+    petflow!(x, Δx, u, theta_known, threads=(threads==:inner))
+    y .= 0.0
+end
+
+function petpdflow!(x, Δx, y, Δy, u, theta_known, :: PrimalOnly; threads=:none)
+    petflow!(x, Δx, u, theta_known, threads=(threads==:inner))
+end
+
 ##########################
 # Linearised optical flow
 ##########################
@@ -358,7 +401,7 @@
     a, c, d = 0.0, 0.0, 0.0
     # This used to use  ∇₂cfold but it is faster to allocate temporary
     # storage for the full gradient due to probably better memory and SIMD
-    # instruction usage. 
+    # instruction usage.
     g = zeros(2, size(b)...)
     ∇₂c!(g, b)
     @inbounds for i=1:size(b, 1)
@@ -425,7 +468,7 @@
         f = x -> simple_imfilter(x, kernel; threads=true)
     end
 
-    # We already filtered b in the previous step (b_next in that step) 
+    # We already filtered b in the previous step (b_next in that step)
     b_filt = b_next_filt==nothing ? f(b) : b_next_filt
     b_next_filt = f(b_next)
 
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/OpticalFlow.jl.orig
--- a/src/OpticalFlow.jl.orig	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,280 +0,0 @@
-################################
-# Code relevant to optical flow
-################################
-
-__precompile__()
-
-module OpticalFlow
-
-using AlgTools.Util
-using ImageTools.Gradient
-import ImageTools.Translate
-using ImageTools.ImFilter
-
-##########
-# Exports
-##########
-
-export flow!,
-       pdflow!,
-       flow_grad!,
-       flow_interp!,
-       estimate_Λ²,
-       estimate_linear_Λ²,
-       pointwise_gradiprod_2d!,
-       pointwise_gradiprod_2dᵀ!,
-       horn_schunck_reg_prox!,
-       horn_schunck_reg_prox_op!,
-       mldivide_step_plus_sym2x2!,
-       linearised_optical_flow_error,
-       Image, AbstractImage, ImageSize,
-       Gradient, Displacement,
-       DisplacementFull, DisplacementConstant,
-       HornSchunckData,
-       filter_hs
-
-###############################################
-# Types (several imported from ImageTools.Translate)
-###############################################
-
-Image = Translate.Image
-AbstractImage = AbstractArray{Float64,2}
-Displacement = Translate.Displacement
-DisplacementFull = Translate.DisplacementFull
-DisplacementConstant = Translate.DisplacementConstant
-Gradient = Array{Float64,3}
-ImageSize = Tuple{Int64,Int64}
-
-#################################
-# Displacement field based flow
-#################################
-
-function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage;
-                      threads = false)
-    tmp .= x
-    Translate.translate_image!(x, tmp, u; threads=threads)
-end
-
-function flow_interp!(x::AbstractImage, u::Displacement;
-                      threads = false)
-    tmp = copy(x)
-    Translate.translate_image!(x, tmp, u; threads=threads)
-end
-
-flow! = flow_interp!
-
-function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none)
-    if dual_flow
-        #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu,
-        #      (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :])))
-        @backgroundif (threads==:outer) begin
-            flow!(x, u, Δx; threads=(threads==:inner))
-        end begin
-            flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
-            flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
-        end
-    else
-        flow!(x, u, Δx)
-    end
-end
-
-function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none)
-    if dual_flow
-        @backgroundif (threads==:outer) begin
-            flow!(x, u, Δx; threads=(threads==:inner))
-            flow!(z, u, Δz; threads=(threads==:inner))
-        end begin
-            flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
-            flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
-        end
-    else
-        flow!(x, u, Δx; threads=(threads==:inner))
-        flow!(z, u, Δz; threads=(threads==:inner))
-    end
-end
-
-##########################
-# Linearised optical flow
-##########################
-
-# ⟨⟨u, ∇b⟩⟩
-function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
-                                 u::DisplacementFull, b::Image;
-                                 add = false)
-    ∇₂c!(vtmp, b)
-
-    u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
-    vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    if add
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i]))
-        end
-    else
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i]))
-        end
-    end
-end
-
-function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
-                                 u::DisplacementConstant, b::Image;
-                                 add = false)
-    ∇₂c!(vtmp, b)
-
-    vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    if add
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] += dot(u, @view(vtmp′[:, i]))
-        end
-    else
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] = dot(u, @view(vtmp′[:, i]))
-        end
-    end
-end
-
-# ∇b ⋅ y
-function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image)
-    ∇₂c!(u, b)
-
-    u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    @simd for i=1:length(y′)
-        @inbounds @. u′[:, i] *= y′[i]
-    end
-end
-
-function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image)
-    @assert(size(y)==size(b) && size(u)==(2,))
-    u .= 0
-    ∇₂cfold!(b, nothing) do g, st, (i, j)
-        @inbounds u .+= g.*y[i, j]
-        return st
-    end
-    # Reweight to be with respect to 𝟙^*𝟙 inner product.
-    u ./= prod(size(b))
-end
-
-mutable struct ConstantDisplacementHornSchunckData
-    M₀::Array{Float64,2}
-    z::Array{Float64,1}
-    Mv::Array{Float64,2}
-    av::Array{Float64,1}
-    cv::Float64
-
-    function ConstantDisplacementHornSchunckData()
-        return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0)
-    end
-end
-
-# For DisplacementConstant, for the simple prox step
-#
-# (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2,
-#
-# construct matrix M₀ and vector z such that we can solve u from
-#
-# (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z
-#
-# Note that the problem
-#
-#    argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2
-#                           + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2
-#
-# has with respect to u the system
-#
-#     (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ,
-#
-# where the primed variables correspond to (2) for (1) for uʹ in place of u:
-#
-#    argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2
-#
-function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData,
-                                   bnext::Image, b::Image, θ, λ, T)
-    @assert(size(b)==size(bnext))
-    w = prod(size(b))
-    z = hs.z
-    cv = 0
-    # Factors of symmetric matrix [a c; c d]
-    a, c, d = 0.0, 0.0, 0.0
-    # This used to use  ∇₂cfold but it is faster to allocate temporary
-    # storage for the full gradient due to probably better memory and SIMD
-    # instruction usage. 
-    g = zeros(2, size(b)...)
-    ∇₂c!(g, b)
-    @inbounds for i=1:size(b, 1)
-        for j=1:size(b, 2)
-            δ = bnext[i,j]-b[i,j]
-            @. z += g[:,i,j]*δ
-            cv += δ*δ
-            a += g[1,i,j]*g[1,i,j]
-            c += g[1,i,j]*g[2,i,j]
-            d += g[2,i,j]*g[2,i,j]
-        end
-    end
-    w₀ = λ
-    w₂ = θ/w
-    aʹ = w₀ + w₂*a
-    cʹ = w₂*c
-    dʹ = w₀ + w₂*d
-    hs.M₀ .= [aʹ cʹ; cʹ dʹ]
-    hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d]
-    hs.cv = cv*θ
-    hs.av .= hs.z.*θ
-    hs.z .*= w₂/T
-end
-
-# Solve the 2D system (I/τ+M₀)u = z
-@inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ)
-    a = 1/τ+M₀[1, 1]
-    c = M₀[1, 2]
-    d = 1/τ+M₀[2, 2]
-    u .= ([d -c; -c a]*z)./(a*d-c*c)
-end
-
-function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image,
-                                θ, λ, T, τ)
-    hs=ConstantDisplacementHornSchunckData()
-    horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T)
-    mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ)
-end
-
-function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing)
-    if !isnothing(δ)
-        u = δ.*u
-    end
-    pointwise_gradiprod_2d!(x, vtmp, u, x; add=true)
-end
-
-# Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation
-function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image)
-    imdim = size(b)
-    vtmp = zeros(2, imdim...)
-    tmp = b-b_prev
-    pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true)
-    return tmp
-end
-
-##############################################
-# Helper to smooth data for Horn–Schunck term
-##############################################
-
-function filter_hs(b, b_next, b_next_filt, kernel)
-    if kernel==nothing
-        f = x -> x
-    else
-        f = x -> simple_imfilter(x, kernel; threads=false)
-    end
-
-    # We already filtered b in the previous step (b_next in that step) 
-    b_filt = b_next_filt==nothing ? f(b) : b_next_filt
-    b_next_filt = f(b_next)
-
-    return b_filt, b_next_filt
-end
-
-end # Module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmDualScaling.jl
--- a/src/PET/AlgorithmDualScaling.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,181 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmDualScaling
-
-identifier = "pdps_known_dualscaling"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!,
-                     DualScaling
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################
-    
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, DualScaling())  
-        
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        #####################        
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end       
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmGreedy.jl
--- a/src/PET/AlgorithmGreedy.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,181 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmGreedy
-
-identifier = "pdps_known_greedy"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!,
-                     Greedy
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################
-    
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, Greedy())   
-        
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        ##################### 
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end       
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmNew.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/PET/AlgorithmNew.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -0,0 +1,183 @@
+####################################################################
+# Predictive online PDPS for optical flow with known velocity field
+####################################################################
+
+__precompile__()
+
+module AlgorithmNew
+
+# Default identifier for when none is given by predictor
+identifier = "pdps_known_noprediction"
+
+using Printf
+
+using AlgTools.Util
+import AlgTools.Iterate
+using ImageTools.Gradient
+using ImageTools.Translate
+  
+using ImageTransformations
+using Images, CoordinateTransformations, Rotations, OffsetArrays
+using ImageCore, Interpolations
+
+using ...Radon
+using ...OpticalFlow: ImageSize,
+                      Image,
+                      petpdflow!
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_rest(x::Image)
+    imdim=size(x)
+    y = zeros(2, imdim...)
+    Δx = copy(x)
+    Δy = copy(y)
+    x̄ = copy(x)
+    radonx = copy(x)
+    return x, y, Δx, Δy, x̄, radonx
+end
+
+function init_iterates(xinit::Image)
+    return init_rest(copy(xinit))
+end
+
+function init_iterates(dim::ImageSize)
+    return init_rest(zeros(dim...))
+end
+
+#########################
+# PETscan related
+#########################
+function petvalue(x, b, c)
+    tmp = similar(b)
+    radon!(tmp, x)
+    return sum(@. tmp - b*log(tmp+c))
+end
+
+function petgrad!(res, x, b, c, S)
+    tmp = similar(b)
+    radon!(tmp, x)
+    @. tmp = S .- b/(tmp+c)
+    backproject!(res, S.*tmp)
+end
+
+function proj_nonneg!(y)
+    @inbounds @simd for i=1:length(y)
+        if y[i] < 0
+            y[i] = 0
+        end
+    end
+    return y
+end
+
+############
+# Algorithm
+############
+
+function solve( :: Type{DisplacementT};
+               dim :: ImageSize,
+               iterate = AlgTools.simple_iterate,
+               params::NamedTuple) where DisplacementT
+
+    ################################
+    # Extract and set up parameters
+    ################################
+    α, ρ = params.α, params.ρ
+    R_K² = ∇₂_norm₂₂_est²
+    γ = 1
+    L = params.L
+    τ₀, σ₀ = params.τ₀, params.σ₀
+    τ = τ₀/L
+    σ = σ₀*(1-τ₀)/(R_K²*τ)
+
+    println("Step length parameters: τ=$(τ), σ=$(σ)")
+
+    λ = params.λ
+    c = params.c*ones(params.radondims...)
+
+    
+    ######################
+    # Initialise iterates
+    ######################
+
+    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
+    
+    if params.L_experiment
+        oldpetgradx = zeros(size(x)...)
+        petgradx = zeros(size(x))
+        oldx = ones(size(x))
+    end
+
+    ####################
+    # Run the algorithm
+    ####################
+                        
+    v = iterate(params) do verbose :: Function,
+                           b :: Image,                   # noisy_sinogram
+                           v_known :: DisplacementT,
+                           theta_known :: DisplacementT,
+                           b_true :: Image,
+                           S :: Image
+
+        ###################
+        # Prediction steps
+        ###################
+    
+        if haskey(params, :predictor) && ~isnothing(params.predictor)
+            petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.predictor)
+        end
+        
+        if params.L_experiment
+            @. oldx = x
+        end
+
+        ############
+        # PDPS step
+        ############
+
+        ∇₂ᵀ!(Δx, y)                    # primal step:
+        @. x̄ = x                       # | save old x for over-relax
+        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
+
+        @. x = x-(τ*λ)*r∇-τ*Δx         # |
+        proj_nonneg!(x)                # | non-negativity constaint prox
+        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
+        ∇₂!(Δy, x̄)                     # dual step:
+        @. y = y + σ*Δy                # |
+        proj_norm₂₁ball!(y, α)         # |  prox
+
+        #####################
+        # L update if needed
+        #####################
+        if params.L_experiment
+            petgrad!(petgradx, x, b, c, S)
+            petgrad!(oldpetgradx, oldx, b, c, S)
+            if norm₂(x-oldx)>1e-12
+                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
+                println("Step length parameters: L=$(L)")
+                τ = τ₀/L
+                σ = σ₀*(1-τ₀)/(R_K²*τ)
+            end
+        end
+       
+        ################################
+        # Give function value if needed
+        ################################
+        
+        v = verbose() do
+            ∇₂!(Δy, x)
+            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
+            value, x, [NaN, NaN], nothing
+        end
+        
+        v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmNoPrediction.jl
--- a/src/PET/AlgorithmNoPrediction.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,180 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmNoPrediction
-
-identifier = "pdps_known_noprediction"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################  
-        
-        # No prediction
-
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        ##################### 
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end      
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmPrimalOnly.jl
--- a/src/PET/AlgorithmPrimalOnly.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,180 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmPrimalOnly
-
-identifier = "pdps_known_primalonly"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################
-    
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, false)   
-        
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        ##################### 
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end       
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmProximal.jl
--- a/src/PET/AlgorithmProximal.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/PET/AlgorithmProximal.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -15,14 +15,14 @@
 using ImageTools.Gradient
 using ImageTools.Translate
   
-using ..Radon   
 using ImageTransformations
 using Images, CoordinateTransformations, Rotations, OffsetArrays
 using ImageCore, Interpolations
 
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!
+using ...Radon
+using ...OpticalFlow: ImageSize,
+                      Image,
+                      petpdflow!
 
 #########################
 # Iterate initialisation
@@ -103,9 +103,9 @@
                iterate = AlgTools.simple_iterate,
                params::NamedTuple) where DisplacementT
 
-    ################################                                        
+    ################################
     # Extract and set up parameters
-    ################################                    
+    ################################
     α, ρ = params.α, params.ρ
     R_K² = ∇₂_norm₂₂_est²
     γ = 1
@@ -139,9 +139,9 @@
                            v_known :: DisplacementT,
                            theta_known :: DisplacementT,
                            b_true :: Image,
-                           S :: Image    
+                           S :: Image
 
-        ###################    
+        ###################
         # Prediction steps
         ###################
 
@@ -151,14 +151,14 @@
             @. oldx = x
         end
 
-        ##############################    
+        ##############################
         # Proximal step of prediction
         ##############################
         ∇₂!(Δy, x)
         @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
-        #@. cc = y + 1000000*σ̃*Δy 
+        #@. cc = y + 1000000*σ̃*Δy
         #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc
-        proj_norm₂₁ball!(y, α) 
+        proj_norm₂₁ball!(y, α)
 
         ############
         # PDPS step
@@ -177,7 +177,7 @@
 
         #####################
         # L update if needed
-        #####################        
+        #####################
         if params.L_experiment
             petgrad!(petgradx, x, b, c, S)
             petgrad!(oldpetgradx, oldx, b, c, S)
@@ -186,18 +186,18 @@
                 println("Step length parameters: L=$(L)")
                 τ = τ₀/L
                 σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end       
+            end
+        end
        
         ################################
         # Give function value if needed
         ################################
         
-        v = verbose() do            
+        v = verbose() do
             ∇₂!(Δy, x)
             value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
             value, x, [NaN, NaN], nothing
-        end 
+        end
         
         v
     end
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmRotation.jl
--- a/src/PET/AlgorithmRotation.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,181 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmRotation
-
-identifier = "pdps_known_rotation"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!,
-                     Rotation
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################
-    
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, Rotation())   
-        
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        ##################### 
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end     
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/AlgorithmZeroDual.jl
--- a/src/PET/AlgorithmZeroDual.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,181 +0,0 @@
-####################################################################
-# Predictive online PDPS for optical flow with known velocity field
-####################################################################
-
-__precompile__()
-
-module AlgorithmZeroDual
-
-identifier = "pdps_known_zerodual"
-
-using Printf
-
-using AlgTools.Util
-import AlgTools.Iterate
-using ImageTools.Gradient
-using ImageTools.Translate
-  
-using ..Radon   
-using ImageTransformations
-using Images, CoordinateTransformations, Rotations, OffsetArrays
-using ImageCore, Interpolations
-
-using ..OpticalFlow: ImageSize,
-                     Image,
-                     petpdflow!
-
-#########################
-# Iterate initialisation
-#########################
-
-function init_rest(x::Image)
-    imdim=size(x)
-    y = zeros(2, imdim...)
-    Δx = copy(x)
-    Δy = copy(y)
-    x̄ = copy(x)
-    radonx = copy(x)
-    return x, y, Δx, Δy, x̄, radonx
-end
-
-function init_iterates(xinit::Image)
-    return init_rest(copy(xinit))
-end
-
-function init_iterates(dim::ImageSize)
-    return init_rest(zeros(dim...))
-end
-
-#########################
-# PETscan related
-#########################
-function petvalue(x, b, c)
-    tmp = similar(b)
-    radon!(tmp, x)
-    return sum(@. tmp - b*log(tmp+c))
-end
-
-function petgrad!(res, x, b, c, S)
-    tmp = similar(b)
-    radon!(tmp, x)
-    @. tmp = S .- b/(tmp+c)
-    backproject!(res, S.*tmp)
-end
-
-function proj_nonneg!(y)
-    @inbounds @simd for i=1:length(y)
-        if y[i] < 0
-            y[i] = 0
-        end
-    end
-    return y
-end
-
-############
-# Algorithm
-############
-
-function solve( :: Type{DisplacementT};
-               dim :: ImageSize,
-               iterate = AlgTools.simple_iterate,
-               params::NamedTuple) where DisplacementT
-
-    ################################                                        
-    # Extract and set up parameters
-    ################################                    
-    α, ρ = params.α, params.ρ
-    R_K² = ∇₂_norm₂₂_est²
-    γ = 1
-    L = params.L
-    τ₀, σ₀ = params.τ₀, params.σ₀
-    τ = τ₀/L
-    σ = σ₀*(1-τ₀)/(R_K²*τ)
-
-    println("Step length parameters: τ=$(τ), σ=$(σ)")
-
-    λ = params.λ
-    c = params.c*ones(params.radondims...)
-
-    
-    ######################
-    # Initialise iterates
-    ######################
-
-    x, y, Δx, Δy, x̄, r∇ = init_iterates(dim)
-    
-    if params.L_experiment
-        oldpetgradx = zeros(size(x)...)
-        petgradx = zeros(size(x))
-        oldx = ones(size(x))
-    end
-
-    ####################
-    # Run the algorithm
-    ####################
-                        
-    v = iterate(params) do verbose :: Function,
-                           b :: Image,                   # noisy_sinogram
-                           v_known :: DisplacementT,
-                           theta_known :: DisplacementT,
-                           b_true :: Image,
-                           S :: Image    
-
-        ###################    
-        # Prediction steps
-        ###################
-    
-        petpdflow!(x, Δx, y, Δy, v_known, theta_known, false)   
-        y .= zeros(size(y)...)
-
-        if params.L_experiment
-            @. oldx = x
-        end
-
-        ############
-        # PDPS step
-        ############
-
-        ∇₂ᵀ!(Δx, y)                    # primal step:
-        @. x̄ = x                       # | save old x for over-relax
-        petgrad!(r∇, x, b, c, S)       # | Calculate gradient of fidelity term
-
-        @. x = x-(τ*λ)*r∇-τ*Δx         # |
-        proj_nonneg!(x)                # | non-negativity constaint prox
-        @. x̄ = 2x - x̄                  # over-relax: x̄ = 2x-x_old
-        ∇₂!(Δy, x̄)                     # dual step:
-        @. y = y + σ*Δy                # |
-        proj_norm₂₁ball!(y, α)         # |  prox
-
-        #####################
-        # L update if needed
-        ##################### 
-        if params.L_experiment
-            petgrad!(petgradx, x, b, c, S)
-            petgrad!(oldpetgradx, oldx, b, c, S)
-            if norm₂(x-oldx)>1e-12
-                L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L)
-                println("Step length parameters: L=$(L)")
-                τ = τ₀/L
-                σ = σ₀*(1-τ₀)/(R_K²*τ)
-            end  
-        end       
-       
-        ################################
-        # Give function value if needed
-        ################################
-        
-        v = verbose() do            
-            ∇₂!(Δy, x)
-            value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy)
-            value, x, [NaN, NaN], nothing, τ, σ
-        end 
-        
-        v
-    end
-
-    return x, y, v
-end
-
-end # Module
-
-
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/ImGenerate.jl
--- a/src/PET/ImGenerate.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,308 +0,0 @@
-###################
-# Image generation
-###################
-
-module ImGenerate
-
-using ColorTypes: Gray
-import TestImages
-# We don't really *directly* depend on QuartzImageIO. The import here is
-# merely a workaround to suppress warnings when  loading TestImages.
-# Something is broken in those packages.
-import QuartzImageIO
-
-using AlgTools.Util
-using AlgTools.Comms
-using ImageTools.Translate
-
-using ..OpticalFlow: Image, DisplacementConstant, DisplacementFull
-using ..Radon
-
-# Added for reproducibility
-import StableRNGs: StableRNG, Random
-const rng = StableRNG(314159)
-
-# Added for PET
-import PoissonRandom: pois_rand
-import Random: shuffle
-import Images: center, warp, imresize
-import CoordinateTransformations: recenter
-import Rotations: RotMatrix
-import Interpolations: Flat
-import MAT: matread
-
-
-##############
-# Our exports
-##############
-
-export ImGen,
-       OnlineData,
-       imgen_square,
-       imgen_shake,
-       PetOnlineData,
-       imgen_shepplogan_radon,
-       imgen_brainphantom_radon
-
-##################
-# Data structures
-##################
-
-struct ImGen
-    f :: Function
-    dim :: Tuple{Int64,Int64}
-    Λ :: Float64
-    dynrange :: Float64
-    name :: String
-end
-
-struct OnlineData{DisplacementT}
-    b_true :: Image
-    b_noisy :: Image
-    v :: DisplacementT
-    v_true :: DisplacementT
-    v_cumul_true :: DisplacementT
-end
-
-struct PetOnlineData{DisplacementT}
-    b_true :: Image
-    sinogram_true :: Image
-    sinogram_noisy :: Image
-    v :: DisplacementT
-    v_true :: DisplacementT
-    v_cumul_true :: DisplacementT
-    theta :: DisplacementT              # theta = thetaknown, theta_cumul
-    S :: Image
-end
-
-###################
-# Shake generation
-###################
-
-function make_const_v(displ, sz)
-    v = zeros(2, sz...)
-    v[1, :, :] .= displ[1]
-    v[2, :, :] .= displ[2]
-    return v
-end
-
-function shake(params)
-    if !haskey(params, :shaketype) || params.shaketype == :gaussian
-        return () -> params.shake.*randn(rng,2)
-    elseif params.shaketype == :disk
-        return () -> begin
-            θ = 2π*rand(rng,Float64)
-            r = params.shake*√(rand(rng,Float64))
-            return [r*cos(θ), r*sin(θ)]
-        end
-    elseif params.shaketype == :circle
-        return () -> begin
-            θ = 2π*rand(rng,Float64)
-            r = params.shake
-            return [r*cos(θ), r*sin(θ)]
-        end        
-    else
-        error("Unknown shaketype $(params.shaketype)")
-    end 
-end
-
-pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz)
-
-
-function rotatebytheta(params)
-    r = params.rotation_factor*randn(rng)
-    return r
-end
-
-function generate_radonmask(params)
-    imdim = params.radondims
-    sino_sparse = params.sino_sparsity
-    numzero = Int64(round(sino_sparse*imdim[1]*imdim[2]))
-    numone = imdim[1]*imdim[2]-numzero
-    A = shuffle(rng,reshape([ones(numone); zeros(numzero)],(imdim[1],imdim[2])))
-    return A
-end
-
-################
-# Moving square
-################
-
-function generate_square(sz,
-                         :: Type{DisplacementT},
-                         datachannel :: Channel{OnlineData{DisplacementT}},
-                         params) where DisplacementT
-
-    if false
-        v₀ = make_const_v(0.1.*(-1, 1), sz)
-        nextv = () -> v₀
-    elseif DisplacementT == DisplacementFull
-        nextv = pixelwise(shake(params), sz)
-    elseif DisplacementT == DisplacementConstant
-        nextv = shake(params)
-    else
-        @error "Invalid DisplacementT"
-    end
-
-    # Constant linear displacement everywhere has Jacobian determinant one
-    # (modulo the boundaries which we ignore here)
-    m = round(Int, sz[1]/5)
-    b_orig = zeros(sz...)
-    b_orig[sz[1].-(2*m:3*m), 2*m:3*m] .= 1
-
-    v_true = nextv()
-    v_cumul = copy(v_true)
-
-    while true
-        # Flow original data and add noise
-        b_true = zeros(sz...)
-        translate_image!(b_true, b_orig, v_cumul; threads=true)
-        b = b_true .+ params.noise_level.*randn(rng,sz...)
-        v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...))
-        # Pass true data to iteration routine
-        data = OnlineData{DisplacementT}(b_true, b, v, v_true, v_cumul)
-        if !put_unless_closed!(datachannel, data)
-            return
-        end
-        # Next step shake
-        v_true = nextv()
-        v_cumul .+= v_true
-    end
-end
-
-function imgen_square(sz)
-    return ImGen(curry(generate_square, sz), sz, 1, 1, "square$(sz[1])x$(sz[2])")
-end
-
-################
-# Shake a photo
-################
-
-function generate_shake_image(im, sz,
-                              :: Type{DisplacementConstant},
-                              datachannel :: Channel{OnlineData{DisplacementConstant}},
-                              params :: NamedTuple)
-    
-
-    # Set up counter and zero factor for stabilisation
-    @assert(params.maxiter ≥ maximum(params.stable_interval))
-    indx = 1
-    zero_factor = indx in params.stable_interval ? 0.0 : 1.0                          
-
-    # Restart the seed to enable comparison across predictors                       
-    Random.seed!(rng,314159)   
-
-
-    nextv = shake(params)
-    v_true = zero_factor.*nextv()
-    v_cumul = copy(v_true)
-
-    while true
-        # Extract subwindow of original image and add noise
-        b_true = zeros(sz...)
-        extract_subimage!(b_true, im, v_cumul; threads=true)
-        b = b_true .+ params.noise_level.*randn(rng,sz...)
-        v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...))
-        # Pass data to iteration routine
-        data = OnlineData{DisplacementConstant}(b_true, b, v, v_true, v_cumul)
-        if !put_unless_closed!(datachannel, data)
-            return
-        end
-        # Next step shake
-        v_true = zero_factor.*nextv()
-        v_cumul .+= v_true
-        # Update indx and zero factor
-        indx += 1
-        zero_factor = indx in params.stable_interval ? 0.0 : 1.0
-    end
-end
-
-function imgen_shake(imname, sz)
-    im = Float64.(Gray.(TestImages.testimage(imname)))
-    dynrange = maximum(im)
-    return ImGen(curry(generate_shake_image, im, sz), sz, 1, dynrange,
-                 "$(imname)$(sz[1])x$(sz[2])")
-end
-
-
-
-########################################################################
-# PETscan
-########################################################################
-function generate_sinogram(im, sz,
-                        :: Type{DisplacementConstant},
-                        datachannel :: Channel{PetOnlineData{DisplacementConstant}},
-                        params :: NamedTuple)
-
-    # Set up counter and zero factor for stabilisation
-    @assert(params.maxiter ≥ maximum(params.stable_interval))
-    indx = 1
-    zero_factor = indx in params.stable_interval ? 0.0 : 1.0
-
-    # Restart the seed to enable comparison across predictors                       
-    Random.seed!(rng,314159)
-
-    nextv = shake(params)
-    v_true = zero_factor.*nextv()
-    v_cumul = copy(v_true)
-
-    S_true = generate_radonmask(params)
-
-    theta_true = zero_factor*rotatebytheta(params)
-    theta_cumul = copy(theta_true)
-
-    while true
-        # Define the transformation matrix
-        center_point = (sz[1]/2 + v_true[1], sz[2]/2 + v_true[2])
-        tform = recenter(RotMatrix(theta_cumul), center_point)
-
-        # Apply the transformation to the image using warp
-        b_true = copy(warp(im, tform, axes(im), fillvalue=Flat()))
-
-        v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...))
-        theta = theta_true*(1.0 + params.rotation_noise_level.*randn(rng))
-
-        # Generate the true and noisy sinograms
-        sinogram_true = zeros(params.radondims...)
-        sinogram_true .*= params.scale
-        radon!(sinogram_true, b_true)
-        sinogram_noisy = copy(sinogram_true)
-
-        for i=1:params.radondims[1], j=1:params.radondims[2]
-            sinogram_noisy[i, j] += pois_rand(rng,params.noise_level)
-        end
-
-        # Pass data to iteration routine
-        data = PetOnlineData{DisplacementConstant}(b_true, sinogram_true, sinogram_noisy, v, v_true, v_cumul, [theta, theta_cumul], S_true)
-        if !put_unless_closed!(datachannel, data)
-        return
-        end
-
-        # Next step shake
-        v_true = zero_factor.*nextv()
-        v_cumul .+= v_true
-        # Next theta
-        theta_true = zero_factor*rotatebytheta(params)
-        theta_cumul += theta_true
-        # Next sinogram mask 
-        S_true = generate_radonmask(params)
-        # Update indx and zero factor
-        indx += 1
-        zero_factor = indx in params.stable_interval ? 0.0 : 1.0
-    end
-end
-
-function imgen_shepplogan_radon(sz)
-    im = convert(Array{Float64},TestImages.shepp_logan(sz[1], highContrast=true))
-    dynrange = maximum(im)
-    return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])")
-end
-
-function imgen_brainphantom_radon(sz)
-    data = matread("src/PET/phantom_slice.mat")
-    im = normalise(imresize(convert(Array{Float64},data["square_data"]),sz))
-    dynrange = maximum(im)
-    return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "brainphantom$(sz[1])x$(sz[2])")
-end
-
-normalise = (data) -> data./maximum(data)
-end # Module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/OpticalFlow.jl
--- a/src/PET/OpticalFlow.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,430 +0,0 @@
-################################
-# Code relevant to optical flow
-################################
-
-__precompile__()
-
-module OpticalFlow
-
-using AlgTools.Util
-using ImageTools.Gradient
-import ImageTools.Translate
-using ImageTools.ImFilter
-
-# using ImageTransformations
-# using Images, CoordinateTransformations, Rotations, OffsetArrays
-# using Interpolations
-
-import Images: center, warp
-import CoordinateTransformations: recenter
-import Rotations: RotMatrix
-import Interpolations: Flat
-
-##########
-# Exports
-##########
-
-export flow!,
-       pdflow!,
-       flow_grad!,
-       flow_interp!,
-       estimate_Λ²,
-       estimate_linear_Λ²,
-       pointwise_gradiprod_2d!,
-       pointwise_gradiprod_2dᵀ!,
-       horn_schunck_reg_prox!,
-       horn_schunck_reg_prox_op!,
-       mldivide_step_plus_sym2x2!,
-       linearised_optical_flow_error,
-       Image, AbstractImage, ImageSize,
-       Gradient, Displacement,
-       DisplacementFull, DisplacementConstant,
-       HornSchunckData,
-       filter_hs,
-       petpdflow!,
-       DualScaling, Greedy, Rotation
-
-###############################################
-# Types (several imported from ImageTools.Translate)
-###############################################
-
-Image = Translate.Image
-AbstractImage = AbstractArray{Float64,2}
-Displacement = Translate.Displacement
-DisplacementFull = Translate.DisplacementFull
-DisplacementConstant = Translate.DisplacementConstant
-Gradient = Array{Float64,3}
-ImageSize = Tuple{Int64,Int64}
-
-
-#################################
-# Struct for flow
-#################################
-struct DualScaling end
-struct Greedy end
-struct Rotation end
-
-#################################
-# Displacement field based flow
-#################################
-
-function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage;
-                      threads = false)
-    tmp .= x
-    Translate.translate_image!(x, tmp, u; threads=threads)
-end
-
-function flow_interp!(x::AbstractImage, u::Displacement;
-                      threads = false)
-    tmp = copy(x)
-    Translate.translate_image!(x, tmp, u; threads=threads)
-end
-
-flow! = flow_interp!
-
-function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none)
-    if dual_flow
-        #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu,
-        #      (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :])))
-        @backgroundif (threads==:outer) begin
-            flow!(x, u, Δx; threads=(threads==:inner))
-        end begin
-            flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
-            flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
-        end
-    else
-        flow!(x, u, Δx)
-    end
-end
-
-function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none)
-    if dual_flow
-        @backgroundif (threads==:outer) begin
-            flow!(x, u, Δx; threads=(threads==:inner))
-            flow!(z, u, Δz; threads=(threads==:inner))
-        end begin
-            flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner))
-            flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner))
-        end
-    else
-        flow!(x, u, Δx; threads=(threads==:inner))
-        flow!(z, u, Δz; threads=(threads==:inner))
-    end
-end
-
-# Additional method for Greedy
-function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none)
-    @assert(size(u)==(2,))  
-    Δx .= x
-    Δy .= y
-    flow!(x, u; threads=(threads==:inner))
-    Dxx = similar(Δy)
-    DΔx = similar(Δy)
-    ∇₂!(Dxx, x)
-    ∇₂!(DΔx, Δx)
-    inds = abs.(Dxx) .≤ 1e-1
-    Dxx[inds] .= 1
-    DΔx[inds] .= 1
-    y .= y.* DΔx ./ Dxx  
-end
-
-# Additional method for Rotation
-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)
-            ndx_banana = @views sum(dx_banana[:, i, j].^2)
-            if ndx > 1e-4 && ndx_banana > 1e-4
-                A = dx[:, i, j]
-                B = dx_banana[:, i, j]
-                theta = atan(B[1] * A[2] - B[2] * A[1], B[1] * A[1] + B[2] * A[2]) # Oriented angle from A to B
-                cos_theta = cos(theta)
-                sin_theta = sin(theta)
-                a = cos_theta * y[1, i, j] - sin_theta * y[2, i, j]
-                b = sin_theta * y[1, i, j] + cos_theta * y[2, i, j]
-                y[1, i, j] = a
-                y[2, i, j] = b
-            end
-        end
-    end
-end
-
-# Additional method for Dual Scaling
-function pdflow!(x, Δx, y, Δy, u, flow :: DualScaling; threads=:none)
-    @assert(size(u)==(2,))
-    oldx = copy(x)
-    flow!(x, u; threads=(threads==:inner))
-    C = similar(y)
-    cc = abs.(x-oldx)
-    cm = max(1e-12,maximum(cc))
-    c = 1 .* (1 .- cc./ cm) .^(10)
-    C[1,:,:] .= c
-    C[2,:,:] .= c
-    y .= C.*y 
-end
-
-
-##########################
-# PET
-##########################
-function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant;
-    threads = false)
-    tmp .= x
-    center_point = center(x) .+ u
-    tform = recenter(RotMatrix(theta_known[1]), center_point)
-    tmp = warp(x, tform, axes(x), fillvalue=Flat())
-    x .= tmp
-end
-
-petflow! = petflow_interp!
-
-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)
-    end
-end
-
-# Method for greedy predictor
-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
-    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 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)
-    Δx = warp(x, tform, axes(x), fillvalue=Flat())
-    @. x = Δx
-    C = similar(y)
-    cc = abs.(x-oldx)
-    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, 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
-
-##########################
-# Linearised optical flow
-##########################
-
-# ⟨⟨u, ∇b⟩⟩
-function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
-                                 u::DisplacementFull, b::Image;
-                                 add = false)
-    ∇₂c!(vtmp, b)
-
-    u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
-    vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    if add
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i]))
-        end
-    else
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i]))
-        end
-    end
-end
-
-function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient,
-                                 u::DisplacementConstant, b::Image;
-                                 add = false)
-    ∇₂c!(vtmp, b)
-
-    vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    if add
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] += dot(u, @view(vtmp′[:, i]))
-        end
-    else
-        @simd for i = 1:length(y′)
-            @inbounds y′[i] = dot(u, @view(vtmp′[:, i]))
-        end
-    end
-end
-
-# ∇b ⋅ y
-function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image)
-    ∇₂c!(u, b)
-
-    u′=reshape(u, (size(u, 1), prod(size(u)[2:end])))
-    y′=reshape(y, prod(size(y)))
-
-    @simd for i=1:length(y′)
-        @inbounds @. u′[:, i] *= y′[i]
-    end
-end
-
-function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image)
-    @assert(size(y)==size(b) && size(u)==(2,))
-    u .= 0
-    ∇₂cfold!(b, nothing) do g, st, (i, j)
-        @inbounds u .+= g.*y[i, j]
-        return st
-    end
-    # Reweight to be with respect to 𝟙^*𝟙 inner product.
-    u ./= prod(size(b))
-end
-
-mutable struct ConstantDisplacementHornSchunckData
-    M₀::Array{Float64,2}
-    z::Array{Float64,1}
-    Mv::Array{Float64,2}
-    av::Array{Float64,1}
-    cv::Float64
-
-    function ConstantDisplacementHornSchunckData()
-        return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0)
-    end
-end
-
-# For DisplacementConstant, for the simple prox step
-#
-# (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2,
-#
-# construct matrix M₀ and vector z such that we can solve u from
-#
-# (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z
-#
-# Note that the problem
-#
-#    argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2
-#                           + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2
-#
-# has with respect to u the system
-#
-#     (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ,
-#
-# where the primed variables correspond to (2) for (1) for uʹ in place of u:
-#
-#    argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2
-#
-function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData,
-                                   bnext::Image, b::Image, θ, λ, T)
-    @assert(size(b)==size(bnext))
-    w = prod(size(b))
-    z = hs.z
-    cv = 0
-    # Factors of symmetric matrix [a c; c d]
-    a, c, d = 0.0, 0.0, 0.0
-    # This used to use  ∇₂cfold but it is faster to allocate temporary
-    # storage for the full gradient due to probably better memory and SIMD
-    # instruction usage. 
-    g = zeros(2, size(b)...)
-    ∇₂c!(g, b)
-    @inbounds for i=1:size(b, 1)
-        for j=1:size(b, 2)
-            δ = bnext[i,j]-b[i,j]
-            @. z += g[:,i,j]*δ
-            cv += δ*δ
-            a += g[1,i,j]*g[1,i,j]
-            c += g[1,i,j]*g[2,i,j]
-            d += g[2,i,j]*g[2,i,j]
-        end
-    end
-    w₀ = λ
-    w₂ = θ/w
-    aʹ = w₀ + w₂*a
-    cʹ = w₂*c
-    dʹ = w₀ + w₂*d
-    hs.M₀ .= [aʹ cʹ; cʹ dʹ]
-    hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d]
-    hs.cv = cv*θ
-    hs.av .= hs.z.*θ
-    hs.z .*= w₂/T
-end
-
-# Solve the 2D system (I/τ+M₀)u = z
-@inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ)
-    a = 1/τ+M₀[1, 1]
-    c = M₀[1, 2]
-    d = 1/τ+M₀[2, 2]
-    u .= ([d -c; -c a]*z)./(a*d-c*c)
-end
-
-function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image,
-                                θ, λ, T, τ)
-    hs=ConstantDisplacementHornSchunckData()
-    horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T)
-    mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ)
-end
-
-function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing)
-    if !isnothing(δ)
-        u = δ.*u
-    end
-    pointwise_gradiprod_2d!(x, vtmp, u, x; add=true)
-end
-
-# Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation
-function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image)
-    imdim = size(b)
-    vtmp = zeros(2, imdim...)
-    tmp = b-b_prev
-    pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true)
-    return tmp
-end
-
-##############################################
-# Helper to smooth data for Horn–Schunck term
-##############################################
-
-function filter_hs(b, b_next, b_next_filt, kernel)
-    if kernel==nothing
-        f = x -> x
-    else
-        f = x -> simple_imfilter(x, kernel; threads=true)
-    end
-
-    # We already filtered b in the previous step (b_next in that step) 
-    b_filt = b_next_filt==nothing ? f(b) : b_next_filt
-    b_next_filt = f(b_next)
-
-    return b_filt, b_next_filt
-end
-
-end # Module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/PET.jl
--- a/src/PET/PET.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/PET/PET.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -18,7 +18,7 @@
 using AlgTools.StructTools
 using AlgTools.LinkedLists
 using AlgTools.Comms
-using ImageTools.Visualise: secs_ns, grayimg, do_visualise 
+using ImageTools.Visualise: secs_ns, grayimg, do_visualise
 using ImageTools.ImFilter: gaussian
 
 # For PET
@@ -26,30 +26,18 @@
 
 #####################
 # Load local modules
-#####################
-include("OpticalFlow.jl")
-include("Radon.jl")
-include("ImGenerate.jl")
-include("AlgorithmDualScaling.jl")
-include("AlgorithmGreedy.jl")
-include("AlgorithmNoPrediction.jl")
-include("AlgorithmPrimalOnly.jl")
+#####################a
+include("AlgorithmNew.jl")
 include("AlgorithmProximal.jl")
-include("AlgorithmRotation.jl")
-include("AlgorithmZeroDual.jl")
 #include("PlotResults.jl")
 
-import .AlgorithmDualScaling
-import .AlgorithmGreedy
-import .AlgorithmNoPrediction
-import .AlgorithmPrimalOnly
+import .AlgorithmNew
 import .AlgorithmProximal
-import .AlgorithmRotation
-import .AlgorithmZeroDual
 
-using .Radon: backproject!
-using .ImGenerate
-using .OpticalFlow: DisplacementFull, DisplacementConstant
+using ..Radon: backproject!
+using ..ImGenerate
+using ..OpticalFlow
+using ..Run
 #using .PlotResults
 
 
@@ -57,9 +45,9 @@
 # Our exports
 ##############
 
-export demo_petS1, demo_petS2, demo_petS3, 
+export demo_petS1, demo_petS2, demo_petS3,
        demo_petS4, demo_petS5, demo_petS6, demo_petS7,
-       demo_petB1, demo_petB2, demo_petB3, 
+       demo_petB1, demo_petB2, demo_petB3,
        demo_petB4, demo_petB5, demo_petB6, demo_petB7,
        batchrun_shepplogan, batchrun_brainphantom, batchrun_pet
        #plot_pet
@@ -68,24 +56,6 @@
 # Parameterisation and experiments
 ###################################
 
-struct Experiment
-    mod :: Module
-    DisplacementT :: Type
-    imgen :: ImGen
-    params :: NamedTuple
-end
-
-function Base.show(io::IO, e::Experiment)
-    displacementname(::Type{DisplacementFull}) = "DisplacementFull"
-    displacementname(::Type{DisplacementConstant}) = "DisplacementConstant"
-    print(io, "
-    mod: $(e.mod)
-    DisplacementT: $(displacementname(e.DisplacementT))
-    imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2])
-    params: $(e.params ⬿ (kernel = "(not shown)",))
-    ")
-end
-
 const default_save_prefix="img/"
 
 const default_params = (
@@ -107,66 +77,68 @@
     stable_interval = Set(0),
 )
 
-const p_known₀_pet = (
+const p_known₀_pet = default_params ⬿ (
     noise_level = 0.5,
     shake_noise_level = 0.1,
-    shake = 1.0, 
-    rotation_factor = 0.075,                  
-    rotation_noise_level = 0.0075,         
+    shake = 1.0,
+    rotation_factor = 0.075,
+    rotation_noise_level = 0.0075,
     α = 0.15,
     ρ̃₀ = 1.0,
     σ̃₀ = 1.0,
     δ = 0.9,
     σ₀ = 1.0,
     τ₀ = 0.9,
-    λ = 1,            
+    λ = 1,
     radondims = [128,64],
     sz = (256,256),
     scale = 1,
     c = 1.0,
-    sino_sparsity = 0.5,            
+    sino_sparsity = 0.5,
     L = 300.0,
     L_experiment = false,
     #stable_interval = Set(0),
     stable_interval = union(Set(1000:2000),Set(3500:4000)),
 )
 
+const p_known₀_petb = p_known₀_pet# ⬿ ( seed = 313159, )
+
 const shepplogan = imgen_shepplogan_radon(p_known₀_pet.sz)
 
 const brainphantom = imgen_brainphantom_radon(p_known₀_pet.sz)
 
 const shepplogan_experiments_pdps_known = (
-    Experiment(AlgorithmDualScaling, DisplacementConstant, shepplogan,
-               p_known₀_pet),
-    Experiment(AlgorithmGreedy, DisplacementConstant, shepplogan,
-               p_known₀_pet),    
-    Experiment(AlgorithmNoPrediction, DisplacementConstant, shepplogan,
-               p_known₀_pet),          
-    Experiment(AlgorithmPrimalOnly, DisplacementConstant, shepplogan,
-               p_known₀_pet), 
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=DualScaling(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=Greedy(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=nothing,),),
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=PrimalOnly(),)),
     Experiment(AlgorithmProximal, DisplacementConstant, shepplogan,
                p_known₀_pet ⬿ (phantom_ρ = 100,)),
-    Experiment(AlgorithmRotation, DisplacementConstant, shepplogan,
-               p_known₀_pet),   
-    Experiment(AlgorithmZeroDual, DisplacementConstant, shepplogan,
-               p_known₀_pet),                       
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=Rotation(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (predictor=ZeroDual(),)),
 )
 
 const brainphantom_experiments_pdps_known = (
-    Experiment(AlgorithmDualScaling, DisplacementConstant, brainphantom,
-               p_known₀_pet),
-    Experiment(AlgorithmGreedy, DisplacementConstant, brainphantom,
-               p_known₀_pet),    
-    Experiment(AlgorithmNoPrediction, DisplacementConstant, brainphantom,
-               p_known₀_pet),          
-    Experiment(AlgorithmPrimalOnly, DisplacementConstant, brainphantom,
-               p_known₀_pet), 
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=DualScaling(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=Greedy(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=nothing,),),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=PrimalOnly(),)),
     Experiment(AlgorithmProximal, DisplacementConstant, brainphantom,
-               p_known₀_pet ⬿ (phantom_ρ = 100,)),
-    Experiment(AlgorithmRotation, DisplacementConstant, brainphantom,
-               p_known₀_pet),   
-    Experiment(AlgorithmZeroDual, DisplacementConstant, brainphantom,
-               p_known₀_pet),                       
+               p_known₀_petb ⬿ (phantom_ρ = 100,)),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=Rotation(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=ZeroDual(),)),
 )
 
 
@@ -178,167 +150,15 @@
     brainphantom_experiments_pdps_known,
 ))
 
-################
-# Log
-################
-
-struct LogEntry <: IterableStruct
-    iter :: Int
-    time :: Float64
-    function_value :: Float64
-    #v_cumul_true_y :: Float64
-    #v_cumul_true_x :: Float64
-    #v_cumul_est_y :: Float64
-    #v_cumul_est_x :: Float64
-    psnr :: Float64
-    ssim :: Float64
-    #psnr_data :: Float64
-    #ssim_data :: Float64
-end
-
-struct LogEntryHiFi <: IterableStruct
-    iter :: Int
-    v_cumul_true_y :: Float64
-    v_cumul_true_x :: Float64
-end
-
-###############
-# Main routine
-###############
-
-struct State
-    vis :: Union{Channel,Bool,Nothing}
-    start_time :: Union{Real,Nothing}
-    wasted_time :: Real
-    log :: LinkedList{LogEntry}
-    log_hifi :: LinkedList{LogEntryHiFi}
-    aborted :: Bool
-end
-
-function name(e::Experiment, p)
-    ig = e.imgen
-    # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))"
-    return "$(ig.name)_$(e.mod.identifier)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))"
-end
-
-function write_tex(texfile, e_params)
-    open(texfile, "w") do io
-        wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}")
-        wf = (n, s) -> if isdefined(e_params, s)
-                            wp(n, getfield(e_params, s))
-                        end
-        wf("alpha", :α)
-        wf("sigmazero", :σ₀)
-        wf("tauzero", :τ₀)
-        wf("tildetauzero", :τ̃₀)
-        wf("delta", :δ)
-        wf("lambda", :λ)
-        wf("theta", :θ)
-        wf("maxiter", :maxiter)
-        wf("noiselevel", :noise_level)
-        wf("shakenoiselevel", :shake_noise_level)
-        wf("shake", :shake)
-        wf("timestep", :timestep)
-        wf("displacementcount", :displacementcount)
-        wf("phantomrho", :phantom_ρ)
-        if isdefined(e_params, :σ₀)
-            wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}")
-        end
-    end
-end                
-
-function run_experiments(;visualise=true,
-                          recalculate=true,
-                          experiments,
-                          save_prefix=default_save_prefix,
-                          fullscreen=false,
-                          kwargs...)
-
-    # Create visualisation
-    if visualise
-        rc = Channel(1)
-        visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen)
-        bind(rc, visproc)
-        vis = rc
-    else
-        vis = false
-    end
-
-    # Run all experiments
-    for e ∈ experiments
-
-        # Parameters for this experiment
-        e_params = default_params ⬿ e.params ⬿ kwargs
-        ename = name(e, e_params)
-        e_params = e_params ⬿ (save_prefix = save_prefix * ename,
-                                dynrange = e.imgen.dynrange,
-                                Λ = e.imgen.Λ)
-
-        if recalculate || !isfile(e_params.save_prefix * ".txt")
-            println("Running experiment \"$(ename)\"")
-
-            # Start data generation task
-            datachannel = Channel{PetOnlineData{e.DisplacementT}}(2)
-            gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params)
-            bind(datachannel, gentask)
-
-            # Run algorithm
-            iterate = curry(iterate_visualise, datachannel,
-                            State(vis, nothing, 0.0, nothing, nothing, false))
-                            
-            x, y, st = e.mod.solve(e.DisplacementT;
-                                   dim=e.imgen.dim,
-                                   iterate=iterate,
-                                   params=e_params)
-
-            # Clear non-saveable things
-            st = @set st.vis = nothing
-
-            println("Wasted_time: $(st.wasted_time)s")
-
-            if e_params.save_results
-                println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)")
-
-                perffile = e_params.save_prefix * ".txt"
-                hififile = e_params.save_prefix * "_hifi.txt"
-                texfile = e_params.save_prefix * "_params.tex"
-                # datafile = e_params.save_prefix * ".jld2"
-
-                write_log(perffile, st.log, "# params = $(e_params)\n")
-                #write_log(hififile, st.log_hifi, "# params = $(e_params)\n")
-                #write_tex(texfile, e_params)
-                # @save datafile x y st params
-            end
-
-            close(datachannel)
-            wait(gentask)
-
-            if st.aborted
-                break
-            end
-        else
-            println("Skipping already computed experiment \"$(ename)\"")
-            # texfile = e_params.save_prefix * "_params.tex"
-            # write_tex(texfile, e_params)
-        end
-    end
-
-    if visualise
-        # Tell subprocess to finish, and wait
-        put!(rc, nothing)
-        close(rc)
-        wait(visproc)
-    end
-
-    return
-end
-
 #######################
 # Demos and batch runs
 #######################
 
 function demo(experiment; kwargs...)
     run_experiments(;experiments=(experiment,),
+                    save_prefix=default_save_prefix,
+                    visfn=iterate_visualise_pet,
+                    datatype=PetOnlineData,
                     save_results=false,
                     save_images=false,
                     visualise=true,
@@ -367,6 +187,9 @@
 
 function batchrun_shepplogan(;kwargs...)
     run_experiments(;experiments=shepplogan_experiments_all,
+                    visfn=iterate_visualise_pet,
+                    datatype=PetOnlineData,
+                    save_prefix=default_save_prefix,
                     save_results=true,
                     save_images=true,
                     visualise=false,
@@ -376,6 +199,9 @@
 
 function batchrun_brainphantom(;kwargs...)
     run_experiments(;experiments=brainphantom_experiments_all,
+                    visfn=iterate_visualise_pet,
+                    datatype=PetOnlineData,
+                    save_prefix=default_save_prefix,
                     save_results=true,
                     save_images=true,
                     visualise=false,
@@ -395,21 +221,21 @@
 function rescale(arr, new_range)
     old_min = minimum(arr)
     old_max = maximum(arr)
-    scale_factor = (new_range[2] - new_range[1]) / (old_max - old_min) 
+    scale_factor = (new_range[2] - new_range[1]) / (old_max - old_min)
     scaled_arr = new_range[1] .+ (arr .- old_min) * scale_factor
     return scaled_arr
 end
 
-function iterate_visualise(datachannel::Channel{PetOnlineData{DisplacementT}},
-                           st :: State,
-                           step :: Function,
-                           params :: NamedTuple) where DisplacementT
+function iterate_visualise_pet(datachannel::Channel{PetOnlineData{DisplacementT}},
+                               st :: State,
+                               step :: Function,
+                               params :: NamedTuple) where DisplacementT
     try
         sc = nothing
 
         d = take!(datachannel)
 
-        for iter=1:params.maxiter 
+        for iter=1:params.maxiter
             dnext = take!(datachannel)
             st = step(d.sinogram_noisy, d.v, d.theta, d.b_true, d.S) do calc_objective
                 stn = st
@@ -435,7 +261,7 @@
                     entry = LogEntry(iter, tm, value,
                                      #sc*d.v_cumul_true[1],
                                      #sc*d.v_cumul_true[2],
-                                     #sc*v[1], sc*v[2], 
+                                     #sc*v[1], sc*v[2],
                                      assess_psnr(x, d.b_true),
                                      assess_ssim(x, d.b_true),
                                      #assess_psnr(d.b_noisy, d.b_true),
@@ -509,39 +335,6 @@
     return st
 end
 
-function bg_visualise_enhanced(rc; fullscreen=false)
-    process_channel(rc) do d
-        imgs, plot_movement, log, vhist = d
-        do_visualise(imgs, refresh=false, fullscreen=fullscreen)
-        # Overlay movement
-        GR.settextcolorind(5)
-        GR.setcharheight(0.015)
-        GR.settextpath(GR.TEXT_PATH_RIGHT)
-        tx, ty = GR.wctondc(0, 1)
-        GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr)
-        if plot_movement
-            sc=1.0
-            p=unfold_linked_list(log)
-            x=map(e -> 1.5+sc*e.v_cumul_true_x, p)
-            y=map(e -> 0.5+sc*e.v_cumul_true_y, p)
-            GR.setlinewidth(2)
-            GR.setlinecolorind(2)
-            GR.polyline(x, y)
-            x=map(e -> 1.5+sc*e.v_cumul_est_x, p)
-            y=map(e -> 0.5+sc*e.v_cumul_est_y, p)
-            GR.setlinecolorind(3)
-            GR.polyline(x, y)
-            if vhist != nothing
-                GR.setlinecolorind(4)
-                x=map(v -> 1.5+sc*v, vhist[:,2])
-                y=map(v -> 0.5+sc*v, vhist[:,1])
-                GR.polyline(x, y)
-            end
-        end
-        GR.updatews() 
-    end
-end
-
 # Clip image values to allowed range
 clip = x -> min(max(x, 0.0), 1.0)
 
@@ -566,13 +359,4 @@
 #    fv_plot("brainphantom")
 #end
 
-
-###############
-# Precompiling
-###############
-
-# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}})
-# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool})
-# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}})
-
 end # Module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PET/Radon.jl
--- a/src/PET/Radon.jl	Thu Apr 25 11:14:41 2024 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,85 +0,0 @@
-# From https://github.com/LMescheder/Radon.jl/blob/master/src/Radon.jl
-# Updated to modern Julia. Fixed radon! to zero out `radonim`.
-
-__precompile__()
-
-
-module Radon
-export radon!, backproject!
-# package code goes here
-
-function radon(im::AbstractArray{T, 2}) where T<:AbstractFloat
-    Nx, Ny = size(im)
-    radonim = similar(im, min(Nx, Ny), max(Nx, Ny))
-    radon!(radonim, im)
-    radonim
-end
-
-function radon!(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2}) where T<:AbstractFloat
-    Nφ, Ns = size(radonim)
-    Nx, Ny = size(im)
-    Ldiag = hypot(Nx, Ny)
-    K = round(Int, Ldiag)
-
-    dLz = Ldiag/(K-1)
-    dLs = Ldiag/(Ns-1)
-    dφ = π/Nφ
-
-    @inbounds for is = 1:Ns, iφ = 1:Nφ
-        φ = (iφ-1)*dφ
-        s = -Ldiag/2 + (is-1)*dLs
-        tx = cos(φ)
-        ty = sin(φ)
-
-        radonim[iφ, is] = zero(T)
-        for k=1:K
-            z = -Ldiag/2 + (k-1)*dLz
-            x = round(Int, Nx/2 + z*ty + s*tx)
-            y = round(Int, Ny/2 - z*tx + s*ty)
-
-            if 1 <= x <= Nx && 1 <= y <= Ny
-                radonim[iφ, is] += im[x, y]*dLz
-            end
-        end
-    end
-    radonim
-end
-
-function backproject(radonim::AbstractArray{T, 2}) where T<:AbstractFloat
-    Nφ, Ns = size(radonim)
-    im = similar(radonim, Ns, Ns)
-    backproject!(im, radonim)
-    im
-end
-
-function backproject!(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2}) where T<:AbstractFloat
-    Nφ, Ns = size(radonim)
-    Nx, Ny = size(im)
-    Ldiag = hypot(Nx, Ny)
-    K = round(Int, Ldiag)
-
-    dLz = Ldiag/(K-1)
-    dLs = Ldiag/(Ns-1)
-    dφ = π/Nφ
-
-    im .= 0.
-    @inbounds for is = 1:Ns, iφ = 1:Nφ
-        φ = (iφ-1)*dφ
-        s = -Ldiag/2 + (is-1)*dLs
-        tx = cos(φ)
-        ty = sin(φ)
-
-        for k=1:K
-            z = -Ldiag/2 + (k-1)*dLz
-            x = round(Int, Nx/2 + z*ty + s*tx)
-            y = round(Int, Ny/2 - z*tx + s*ty)
-
-            if 1 <= x <= Nx && 1 <= y <= Ny
-                im[x, y] += radonim[iφ, is]*dLz
-            end
-        end
-    end
-    im
-end
-
-end # module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/PredictPDPS.jl
--- a/src/PredictPDPS.jl	Thu Apr 25 11:14:41 2024 -0500
+++ b/src/PredictPDPS.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -10,20 +10,8 @@
 # Load external modules
 ########################
 
-using Printf
-using FileIO
-#using JLD2
-using Setfield
-using ImageQualityIndexes: assess_psnr, assess_ssim
-using DelimitedFiles
-import GR
-
+using ImageTools.ImFilter: gaussian
 using AlgTools.Util
-using AlgTools.StructTools
-using AlgTools.LinkedLists
-using AlgTools.Comms
-using ImageTools.Visualise: secs_ns, grayimg, do_visualise
-using ImageTools.ImFilter: gaussian
 
 #####################
 # Load local modules
@@ -32,50 +20,29 @@
 include("OpticalFlow.jl")
 include("Radon.jl")
 include("ImGenerate.jl")
-include("Algorithm.jl")
-include("AlgorithmBoth.jl")
-include("AlgorithmBothGreedyV.jl")
-include("AlgorithmBothCumul.jl")
+include("Run.jl")
+include("AlgorithmProximal.jl")
 include("AlgorithmBothMulti.jl")
-include("AlgorithmBothNL.jl")
 include("AlgorithmFB.jl")
 include("AlgorithmFBDual.jl")
+include("AlgorithmNew.jl")
 include("Stats.jl")
 #include("PlotResults.jl")
-
-
-# Additional
-include("AlgorithmProximal.jl")
-include("AlgorithmGreedy.jl")
-include("AlgorithmRotation.jl")
-include("AlgorithmNoPrediction.jl")
-include("AlgorithmPrimalOnly.jl")
-include("AlgorithmDualScaling.jl")
-include("AlgorithmZeroDual.jl")
 include("PET/PET.jl")
 
 
-import .Algorithm,
-       .AlgorithmBoth,
-       .AlgorithmBothGreedyV,
-       .AlgorithmBothCumul,
-       .AlgorithmBothMulti,
-       .AlgorithmBothNL,
+import .AlgorithmBothMulti,
        .AlgorithmFB,
        .AlgorithmFBDual,
        .AlgorithmProximal,
-       .AlgorithmGreedy,
-       .AlgorithmRotation,
-       .AlgorithmNoPrediction,
-       .AlgorithmPrimalOnly,
-       .AlgorithmDualScaling,
-       .AlgorithmZeroDual
+       .AlgorithmNew
 
 using .ImGenerate
-using .OpticalFlow: DisplacementFull, DisplacementConstant
+using .OpticalFlow
 using .Stats
 #using .PlotResults
 using .PET
+using .Run
 
 ##############
 # Our exports
@@ -101,24 +68,6 @@
 # Parameterisation and experiments
 ###################################
 
-struct Experiment
-    mod :: Module
-    DisplacementT :: Type
-    imgen :: ImGen
-    params :: NamedTuple
-end
-
-function Base.show(io::IO, e::Experiment)
-    displacementname(::Type{DisplacementFull}) = "DisplacementFull"
-    displacementname(::Type{DisplacementConstant}) = "DisplacementConstant"
-    print(io, "
-    mod: $(e.mod)
-    DisplacementT: $(displacementname(e.DisplacementT))
-    imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2])
-    params: $(e.params ⬿ (kernel = "(not shown)",))
-    ")
-end
-
 const default_save_prefix="img/"
 
 const default_params = (
@@ -133,20 +82,17 @@
                              1000, 2000, 2500, 3000, 4000, 5000,
                              6000, 7000, 7500, 8000, 9000, 10000, 8700]),
     pixelwise_displacement=false,
-    dual_flow = true,
-    prox_predict = true,
+    dual_flow = true, # For AlgorithmProximalfrom 2019 paper
     handle_interrupt = true,
     init = :zero,
     plot_movement = false,
     stable_interval = Set(0),
-    ds_exponent = 10,
-    ds_threshold = 10e-12
 )
 
 const square = imgen_square((200, 300))
 const lighthouse = imgen_shake("lighthouse", (200, 300))
 
-const p_known₀_denoising = (
+const p_known₀_denoising = default_params ⬿ (
     noise_level = 0.5,
     shake_noise_level = 0.025,
     shake = 2.0,
@@ -160,7 +106,7 @@
     stable_interval = union(Set(2500:5000),Set(8700:10000)),
 )
 
-const p_known₀ = (
+const p_known₀ = default_params ⬿ (
     noise_level = 0.5,
     shake_noise_level = 0.05,
     shake = 2,
@@ -172,7 +118,7 @@
     τ₀ = 0.01,
 )
 
-const p_unknown₀ = (
+const p_unknown₀ = default_params ⬿ (
     noise_level = 0.3,
     shake_noise_level = 0.05,
     shake = 2,
@@ -190,12 +136,14 @@
 )
 
 
+# Experiments for 2019 paper
+
 const experiments_pdps_known = (
-    Experiment(Algorithm, DisplacementConstant, lighthouse,
+    Experiment(AlgorithmProximal, DisplacementConstant, lighthouse,
                p_known₀_denoising ⬿ (phantom_ρ = 0,)),
-    Experiment(Algorithm, DisplacementConstant, lighthouse,
+    Experiment(AlgorithmProximal, DisplacementConstant, lighthouse,
                p_known₀ ⬿ (phantom_ρ = 100,)),
-    Experiment(Algorithm, DisplacementConstant, square,
+    Experiment(AlgorithmProximal, DisplacementConstant, square,
                p_known₀ ⬿ (phantom_ρ = 0,))
 )
 
@@ -219,182 +167,29 @@
     experiments_fb_known
 ))
 
+# Image stabilisation experiments for 2024 paper. PET experiments are in PET/PET.jl
+
 const denoising_experiments_pdps_known = (
-    Experiment(AlgorithmDualScaling, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
-    Experiment(AlgorithmGreedy, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
-    Experiment(AlgorithmNoPrediction, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
-    Experiment(AlgorithmPrimalOnly, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=DualScaling(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=Greedy(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=nothing,),),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=PrimalOnly(),)),
     Experiment(AlgorithmProximal, DisplacementConstant, lighthouse,
                p_known₀_denoising ⬿ (phantom_ρ = 100,)),
-    Experiment(AlgorithmRotation, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
-    Experiment(AlgorithmZeroDual, DisplacementConstant, lighthouse,
-               p_known₀_denoising),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=Rotation(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=ZeroDual(),)),
 )
 
 const denoising_experiments_all = Iterators.flatten((
     denoising_experiments_pdps_known,
 ))
 
-################
-# Log
-################
-
-struct LogEntry <: IterableStruct
-    iter :: Int
-    time :: Float64
-    function_value :: Float64
-    #v_cumul_true_y :: Float64
-    #v_cumul_true_x :: Float64
-    #v_cumul_est_y :: Float64
-    #v_cumul_est_x :: Float64
-    psnr :: Float64
-    ssim :: Float64
-    #psnr_data :: Float64
-    #ssim_data :: Float64
-end
-
-struct LogEntryHiFi <: IterableStruct
-    iter :: Int
-    v_cumul_true_y :: Float64
-    v_cumul_true_x :: Float64
-end
-
-###############
-# Main routine
-###############
-
-struct State
-    vis :: Union{Channel,Bool,Nothing}
-    start_time :: Union{Real,Nothing}
-    wasted_time :: Real
-    log :: LinkedList{LogEntry}
-    log_hifi :: LinkedList{LogEntryHiFi}
-    aborted :: Bool
-end
-
-function name(e::Experiment, p)
-    ig = e.imgen
-    # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))"
-    return "$(ig.name)_$(e.mod.identifier)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))"
-end
-
-function write_tex(texfile, e_params)
-    open(texfile, "w") do io
-        wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}")
-        wf = (n, s) -> if isdefined(e_params, s)
-                            wp(n, getfield(e_params, s))
-                        end
-        wf("alpha", :α)
-        wf("sigmazero", :σ₀)
-        wf("tauzero", :τ₀)
-        wf("tildetauzero", :τ̃₀)
-        wf("delta", :δ)
-        wf("lambda", :λ)
-        wf("theta", :θ)
-        wf("maxiter", :maxiter)
-        wf("noiselevel", :noise_level)
-        wf("shakenoiselevel", :shake_noise_level)
-        wf("shake", :shake)
-        wf("timestep", :timestep)
-        wf("displacementcount", :displacementcount)
-        wf("phantomrho", :phantom_ρ)
-        if isdefined(e_params, :σ₀)
-            wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}")
-        end
-    end
-end
-
-function run_experiments(;visualise=true,
-                          recalculate=true,
-                          experiments,
-                          save_prefix=default_save_prefix,
-                          fullscreen=false,
-                          kwargs...)
-
-    # Create visualisation
-    if visualise
-        rc = Channel(1)
-        visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen)
-        bind(rc, visproc)
-        vis = rc
-    else
-        vis = false
-    end
-
-    # Run all experiments
-    for e ∈ experiments
-
-        # Parameters for this experiment
-        e_params = default_params ⬿ e.params ⬿ kwargs
-        ename = name(e, e_params)
-        e_params = e_params ⬿ (save_prefix = save_prefix * ename,
-                                dynrange = e.imgen.dynrange,
-                                Λ = e.imgen.Λ)
-
-        if recalculate || !isfile(e_params.save_prefix * ".txt")
-            println("Running experiment \"$(ename)\"")
-
-            # Start data generation task
-            datachannel = Channel{OnlineData{e.DisplacementT}}(2)
-            gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params)
-            bind(datachannel, gentask)
-
-            # Run algorithm
-            iterate = curry(iterate_visualise, datachannel,
-                            State(vis, nothing, 0.0, nothing, nothing, false))
-                            
-            x, y, st = e.mod.solve(e.DisplacementT;
-                                   dim=e.imgen.dim,
-                                   iterate=iterate,
-                                   params=e_params)
-
-            # Clear non-saveable things
-            st = @set st.vis = nothing
-
-            println("Wasted_time: $(st.wasted_time)s")
-
-            if e_params.save_results
-                println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)")
-
-                perffile = e_params.save_prefix * ".txt"
-                hififile = e_params.save_prefix * "_hifi.txt"
-                texfile = e_params.save_prefix * "_params.tex"
-                # datafile = e_params.save_prefix * ".jld2"
-
-                write_log(perffile, st.log, "# params = $(e_params)\n")
-                #write_log(hififile, st.log_hifi, "# params = $(e_params)\n")
-                #write_tex(texfile, e_params)
-                # @save datafile x y st params
-            end
-
-            close(datachannel)
-            wait(gentask)
-
-            if st.aborted
-                break
-            end
-        else
-            println("Skipping already computed experiment \"$(ename)\"")
-            # texfile = e_params.save_prefix * "_params.tex"
-            # write_tex(texfile, e_params)
-        end
-    end
-
-    if visualise
-        # Tell subprocess to finish, and wait
-        put!(rc, nothing)
-        close(rc)
-        wait(visproc)
-    end
-
-    return
-end
-
 #######################
 # Demos and batch runs
 #######################
@@ -403,6 +198,7 @@
     run_experiments(;experiments=(experiment,),
                     save_results=false,
                     save_images=false,
+                    save_prefix=default_save_prefix,
                     visualise=true,
                     recalculate=true,
                     verbose_iter=50,
@@ -428,6 +224,7 @@
 
 function batchrun_article(kwargs...)
     run_experiments(;experiments=experiments_all,
+                    save_prefix=default_save_prefix,
                     save_results=true,
                     save_images=true,
                     visualise=false,
@@ -437,6 +234,7 @@
 
 function batchrun_denoising(;kwargs...)
     run_experiments(;experiments=denoising_experiments_all,
+                    save_prefix=default_save_prefix,
                     save_results=true,
                     save_images=true,
                     visualise=false,
@@ -450,150 +248,6 @@
     batchrun_pet(;kwargs...)
 end
 
-######################################################
-# Iterator that does visualisation and log collection
-######################################################
-
-function iterate_visualise(datachannel::Channel{OnlineData{DisplacementT}},
-                           st :: State,
-                           step :: Function,
-                           params :: NamedTuple) where DisplacementT
-    try
-        sc = nothing
-
-        d = take!(datachannel)
-
-        for iter=1:params.maxiter
-            dnext = take!(datachannel)
-            st = step(d.b_noisy, d.v, dnext.b_noisy) do calc_objective
-                stn = st
-
-                if isnothing(stn.start_time)
-                    # The Julia precompiler is a miserable joke, apparently not crossing module
-                    # boundaries, so only start timing after the first iteration.
-                    stn = @set stn.start_time=secs_ns()
-                end
-
-                verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0
-
-                # Normalise movement to image dimensions so
-                # our TikZ plotting code doesn't need to know
-                # the image pixel size.
-                sc = 1.0./maximum(size(d.b_true))
-                    
-                if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0)
-                    verb_start = secs_ns()
-                    tm = verb_start - stn.start_time - stn.wasted_time
-                    value, x, v, vhist = calc_objective()
-
-                    entry = LogEntry(iter, tm, value,
-                                     #sc*d.v_cumul_true[1],
-                                     #sc*d.v_cumul_true[2],
-                                     #sc*v[1], sc*v[2],
-                                     assess_psnr(x, d.b_true),
-                                     assess_ssim(x, d.b_true),
-                                     #assess_psnr(d.b_noisy, d.b_true),
-                                     #assess_ssim(d.b_noisy, d.b_true)
-                                     )
-
-                    # (**) Collect a singly-linked list of log to avoid array resizing
-                    # while iterating
-                    stn = @set stn.log=LinkedListEntry(entry, stn.log)
-                    
-                    if !isnothing(vhist)
-                        vhist=vhist.*sc
-                    end
-
-                    if verb
-                        @printf("%d/%d J=%f, PSNR=%f, SSIM=%f, avg. FPS=%f\n",
-                                iter, params.maxiter, value, entry.psnr,
-                                entry.ssim, entry.iter/entry.time)
-                        if isa(stn.vis, Channel)
-                            put_onlylatest!(stn.vis, ((d.b_noisy, x),
-                                                      params.plot_movement,
-                                                      stn.log, vhist))
-
-                        end
-                    end
-
-                    if params.save_images && (!haskey(params, :save_images_iters)
-                                              || iter ∈ params.save_images_iters)
-                        fn = (t, ext) -> "$(params.save_prefix)_$(t)_frame$(iter).$(ext)"
-                        # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true))
-                        # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy))
-                        save(File(format"PNG", fn("reco", "png")), grayimg(x))
-                        if !isnothing(vhist)
-                            open(fn("movement", "txt"), "w") do io
-                                writedlm(io, ["est_y" "est_x"])
-                                writedlm(io, vhist)
-                            end
-                        end
-                    end
-
-                    stn = @set stn.wasted_time += (secs_ns() - verb_start)
-
-                    return stn
-                end
-
-                hifientry = LogEntryHiFi(iter, sc*d.v_cumul_true[1], sc*d.v_cumul_true[2])
-                st = @set st.log_hifi=LinkedListEntry(hifientry, st.log_hifi)
-
-                return st
-            end
-            d=dnext
-        end
-    catch ex
-        if params.handle_interrupt && isa(ex, InterruptException)
-            # If SIGINT is received (user pressed ^C), terminate computations,
-            # returning current status. Effectively, we do not call `step()` again,
-            # ending the iterations, but letting the algorithm finish up.
-            # Assuming (**) above occurs atomically, `st.log` should be valid, but
-            # any results returned by the algorithm itself may be partial, as for
-            # reasons of efficiency we do *not* store results of an iteration until
-            # the next iteration is finished.
-            printstyled("\rUser interrupt—finishing up.\n", bold=true, color=202)
-            st = @set st.aborted = true
-        else
-            rethrow(ex)
-        end
-    end
-    
-    return st
-end
-
-function bg_visualise_enhanced(rc; fullscreen=false)
-    process_channel(rc) do d
-        imgs, plot_movement, log, vhist = d
-        do_visualise(imgs, refresh=false, fullscreen=fullscreen)
-        # Overlay movement
-        GR.settextcolorind(5)
-        GR.setcharheight(0.015)
-        GR.settextpath(GR.TEXT_PATH_RIGHT)
-        tx, ty = GR.wctondc(0, 1)
-        GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr)
-        if plot_movement
-            sc=1.0
-            p=unfold_linked_list(log)
-            x=map(e -> 1.5+sc*e.v_cumul_true_x, p)
-            y=map(e -> 0.5+sc*e.v_cumul_true_y, p)
-            GR.setlinewidth(2)
-            GR.setlinecolorind(2)
-            GR.polyline(x, y)
-            x=map(e -> 1.5+sc*e.v_cumul_est_x, p)
-            y=map(e -> 0.5+sc*e.v_cumul_est_y, p)
-            GR.setlinecolorind(3)
-            GR.polyline(x, y)
-            if vhist != nothing
-                GR.setlinecolorind(4)
-                x=map(v -> 1.5+sc*v, vhist[:,2])
-                y=map(v -> 0.5+sc*v, vhist[:,1])
-                GR.polyline(x, y)
-            end
-        end
-        GR.updatews()
-    end
-end
-
 #########################
 # Plotting SSIM and PSNR
 #########################
@@ -604,13 +258,4 @@
 #    fv_plot("lighthouse")
 #end
 
-
-###############
-# Precompiling
-###############
-
-# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}})
-# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool})
-# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}})
-
 end # Module
diff -r 74b1a9f0c35e -r e4a8f662a1ac src/Run.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Run.jl	Thu Apr 25 13:05:40 2024 -0500
@@ -0,0 +1,360 @@
+##################
+# Experiment running and interactive visualisation
+##################
+
+module Run
+
+import GR
+using Setfield
+using Printf
+using ImageQualityIndexes: assess_psnr, assess_ssim
+using FileIO
+#using JLD2
+using DelimitedFiles
+
+using AlgTools.Util
+using AlgTools.StructTools
+using AlgTools.LinkedLists
+using AlgTools.Comms
+
+using ImageTools.Visualise: secs_ns, grayimg, do_visualise
+
+using ..ImGenerate
+using ..OpticalFlow: identifier, DisplacementFull, DisplacementConstant
+
+export run_experiments,
+       Experiment,
+       State,
+       LogEntry,
+       LogEntryHiFi
+
+################
+# Experiment
+################
+
+struct Experiment
+    mod :: Module
+    DisplacementT :: Type
+    imgen :: ImGen
+    params :: NamedTuple
+end
+
+function Base.show(io::IO, e::Experiment)
+    displacementname(::Type{DisplacementFull}) = "DisplacementFull"
+    displacementname(::Type{DisplacementConstant}) = "DisplacementConstant"
+    print(io, "
+    mod: $(e.mod)
+    DisplacementT: $(displacementname(e.DisplacementT))
+    imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2])
+    params: $(e.params ⬿ (kernel = "(not shown)",))
+    ")
+end
+
+################
+# Log
+################
+
+struct LogEntry <: IterableStruct
+    iter :: Int
+    time :: Float64
+    function_value :: Float64
+    #v_cumul_true_y :: Float64
+    #v_cumul_true_x :: Float64
+    #v_cumul_est_y :: Float64
+    #v_cumul_est_x :: Float64
+    psnr :: Float64
+    ssim :: Float64
+    #psnr_data :: Float64
+    #ssim_data :: Float64
+end
+
+struct LogEntryHiFi <: IterableStruct
+    iter :: Int
+    v_cumul_true_y :: Float64
+    v_cumul_true_x :: Float64
+end
+
+###############
+# Main routine
+###############
+
+struct State
+    vis :: Union{Channel,Bool,Nothing}
+    start_time :: Union{Real,Nothing}
+    wasted_time :: Real
+    log :: LinkedList{LogEntry}
+    log_hifi :: LinkedList{LogEntryHiFi}
+    aborted :: Bool
+end
+
+function name(e::Experiment, p)
+    ig = e.imgen
+    id = if haskey(p, :predictor) && ~isnothing(p.predictor)
+        identifier(p.predictor)
+    else
+        e.mod.identifier
+    end
+    # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))"
+    return "$(ig.name)_$(id)_$(Int64(100*p.α))_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))"
+end
+
+function write_tex(texfile, e_params)
+    open(texfile, "w") do io
+        wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}")
+        wf = (n, s) -> if isdefined(e_params, s)
+                            wp(n, getfield(e_params, s))
+                        end
+        wf("alpha", :α)
+        wf("sigmazero", :σ₀)
+        wf("tauzero", :τ₀)
+        wf("tildetauzero", :τ̃₀)
+        wf("delta", :δ)
+        wf("lambda", :λ)
+        wf("theta", :θ)
+        wf("maxiter", :maxiter)
+        wf("noiselevel", :noise_level)
+        wf("shakenoiselevel", :shake_noise_level)
+        wf("shake", :shake)
+        wf("timestep", :timestep)
+        wf("displacementcount", :displacementcount)
+        wf("phantomrho", :phantom_ρ)
+        if isdefined(e_params, :σ₀)
+            wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}")
+        end
+    end
+end
+
+function run_experiments(;visualise=true,
+                          visfn = iterate_visualise,
+                          datatype = OnlineData,
+                          recalculate=true,
+                          experiments,
+                          save_prefix="",
+                          fullscreen=false,
+                          kwargs...)
+
+    # Create visualisation
+    if visualise
+        rc = Channel(1)
+        visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen)
+        bind(rc, visproc)
+        vis = rc
+    else
+        vis = false
+    end
+
+    # Run all experiments
+    for e ∈ experiments
+
+        # Parameters for this experiment
+        e_params = e.params ⬿ kwargs
+        ename = name(e, e_params)
+        e_params = e_params ⬿ (save_prefix = save_prefix * ename,
+                                dynrange = e.imgen.dynrange,
+                                Λ = e.imgen.Λ)
+
+        if recalculate || !isfile(e_params.save_prefix * ".txt")
+            println("Running experiment \"$(ename)\"")
+
+            # Start data generation task
+            datachannel = Channel{datatype{e.DisplacementT}}(2)
+            gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params)
+            bind(datachannel, gentask)
+
+            # Run algorithm
+            iterate = curry(visfn, datachannel,
+                            State(vis, nothing, 0.0, nothing, nothing, false))
+                            
+            x, y, st = e.mod.solve(e.DisplacementT;
+                                   dim=e.imgen.dim,
+                                   iterate=iterate,
+                                   params=e_params)
+
+            # Clear non-saveable things
+            st = @set st.vis = nothing
+
+            println("Wasted_time: $(st.wasted_time)s")
+
+            if e_params.save_results
+                println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)")
+
+                perffile = e_params.save_prefix * ".txt"
+                hififile = e_params.save_prefix * "_hifi.txt"
+                texfile = e_params.save_prefix * "_params.tex"
+                # datafile = e_params.save_prefix * ".jld2"
+
+                write_log(perffile, st.log, "# params = $(e_params)\n")
+                #write_log(hififile, st.log_hifi, "# params = $(e_params)\n")
+                #write_tex(texfile, e_params)
+                # @save datafile x y st params
+            end
+
+            close(datachannel)
+            wait(gentask)
+
+            if st.aborted
+                break
+            end
+        else
+            println("Skipping already computed experiment \"$(ename)\"")
+            # texfile = e_params.save_prefix * "_params.tex"
+            # write_tex(texfile, e_params)
+        end
+    end
+
+    if visualise
+        # Tell subprocess to finish, and wait
+        put!(rc, nothing)
+        close(rc)
+        wait(visproc)
+    end
+
+    return
+end
+
+
+######################################################
+# Iterator that does visualisation and log collection
+######################################################
+
+function iterate_visualise(datachannel::Channel{OnlineData{DisplacementT}},
+                           st :: State,
+                           step :: Function,
+                           params :: NamedTuple) where DisplacementT
+    try
+        sc = nothing
+
+        d = take!(datachannel)
+
+        for iter=1:params.maxiter
+            dnext = take!(datachannel)
+            st = step(d.b_noisy, d.v, dnext.b_noisy) do calc_objective
+                stn = st
+
+                if isnothing(stn.start_time)
+                    # The Julia precompiler is a miserable joke, apparently not crossing module
+                    # boundaries, so only start timing after the first iteration.
+                    stn = @set stn.start_time=secs_ns()
+                end
+
+                verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0
+
+                # Normalise movement to image dimensions so
+                # our TikZ plotting code doesn't need to know
+                # the image pixel size.
+                sc = 1.0./maximum(size(d.b_true))
+                    
+                if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0)
+                    verb_start = secs_ns()
+                    tm = verb_start - stn.start_time - stn.wasted_time
+                    value, x, v, vhist = calc_objective()
+
+                    entry = LogEntry(iter, tm, value,
+                                     #sc*d.v_cumul_true[1],
+                                     #sc*d.v_cumul_true[2],
+                                     #sc*v[1], sc*v[2],
+                                     assess_psnr(x, d.b_true),
+                                     assess_ssim(x, d.b_true),
+                                     #assess_psnr(d.b_noisy, d.b_true),
+                                     #assess_ssim(d.b_noisy, d.b_true)
+                                     )
+
+                    # (**) Collect a singly-linked list of log to avoid array resizing
+                    # while iterating
+                    stn = @set stn.log=LinkedListEntry(entry, stn.log)
+                    
+                    if !isnothing(vhist)
+                        vhist=vhist.*sc
+                    end
+
+                    if verb
+                        @printf("%d/%d J=%f, PSNR=%f, SSIM=%f, avg. FPS=%f\n",
+                                iter, params.maxiter, value, entry.psnr,
+                                entry.ssim, entry.iter/entry.time)
+                        if isa(stn.vis, Channel)
+                            put_onlylatest!(stn.vis, ((d.b_noisy, x),
+                                                      params.plot_movement,
+                                                      stn.log, vhist))
+
+                        end
+                    end
+
+                    if params.save_images && (!haskey(params, :save_images_iters)
+                                              || iter ∈ params.save_images_iters)
+                        fn = (t, ext) -> "$(params.save_prefix)_$(t)_frame$(iter).$(ext)"
+                        # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true))
+                        # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy))
+                        save(File(format"PNG", fn("reco", "png")), grayimg(x))
+                        if !isnothing(vhist)
+                            open(fn("movement", "txt"), "w") do io
+                                writedlm(io, ["est_y" "est_x"])
+                                writedlm(io, vhist)
+                            end
+                        end
+                    end
+
+                    stn = @set stn.wasted_time += (secs_ns() - verb_start)
+
+                    return stn
+                end
+
+                hifientry = LogEntryHiFi(iter, sc*d.v_cumul_true[1], sc*d.v_cumul_true[2])
+                st = @set st.log_hifi=LinkedListEntry(hifientry, st.log_hifi)
+
+                return st
+            end
+            d=dnext
+        end
+    catch ex
+        if params.handle_interrupt && isa(ex, InterruptException)
+            # If SIGINT is received (user pressed ^C), terminate computations,
+            # returning current status. Effectively, we do not call `step()` again,
+            # ending the iterations, but letting the algorithm finish up.
+            # Assuming (**) above occurs atomically, `st.log` should be valid, but
+            # any results returned by the algorithm itself may be partial, as for
+            # reasons of efficiency we do *not* store results of an iteration until
+            # the next iteration is finished.
+            printstyled("\rUser interrupt—finishing up.\n", bold=true, color=202)
+            st = @set st.aborted = true
+        else
+            rethrow(ex)
+        end
+    end
+    
+    return st
+end
+
+function bg_visualise_enhanced(rc; fullscreen=false)
+    process_channel(rc) do d
+        imgs, plot_movement, log, vhist = d
+        do_visualise(imgs, refresh=false, fullscreen=fullscreen)
+        # Overlay movement
+        GR.settextcolorind(5)
+        GR.setcharheight(0.015)
+        GR.settextpath(GR.TEXT_PATH_RIGHT)
+        tx, ty = GR.wctondc(0, 1)
+        GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr)
+        if plot_movement
+            sc=1.0
+            p=unfold_linked_list(log)
+            x=map(e -> 1.5+sc*e.v_cumul_true_x, p)
+            y=map(e -> 0.5+sc*e.v_cumul_true_y, p)
+            GR.setlinewidth(2)
+            GR.setlinecolorind(2)
+            GR.polyline(x, y)
+            x=map(e -> 1.5+sc*e.v_cumul_est_x, p)
+            y=map(e -> 0.5+sc*e.v_cumul_est_y, p)
+            GR.setlinecolorind(3)
+            GR.polyline(x, y)
+            if vhist != nothing
+                GR.setlinecolorind(4)
+                x=map(v -> 1.5+sc*v, vhist[:,2])
+                y=map(v -> 0.5+sc*v, vhist[:,1])
+                GR.polyline(x, y)
+            end
+        end
+        GR.updatews()
+    end
+end
+
+end # module