# HG changeset patch
# User Neil Dizon <neil.dizon@helsinki.fi>
# Date 1714482025 -10800
# Node ID 70fd92ac9da097f7fd198794eac6b8ec49faa90c
# Parent  c67ba8ef595e2b6156008dbc6152e1b84f24d047
experimental dual scaling

diff -r c67ba8ef595e -r 70fd92ac9da0 src/OpticalFlow.jl
--- a/src/OpticalFlow.jl	Thu Apr 25 15:21:45 2024 -0500
+++ b/src/OpticalFlow.jl	Tue Apr 30 16:00:25 2024 +0300
@@ -42,7 +42,7 @@
        HornSchunckData,
        filter_hs,
        petpdflow!,
-       DualScaling, Greedy, Rotation, ZeroDual, PrimalOnly,
+       DualScaling, Greedy, Rotation, ZeroDual, PrimalOnly, ActivatedDual,
        identifier
 
 ###############################################
@@ -64,13 +64,14 @@
 struct DualScaling
     exponent :: Integer
     threshold :: Real
-    DualScaling(e = 10, t = 1e-12) = new(e, t)
+    DualScaling(e = 50, t = 1e-12) = new(e, t)
 end
 
 struct Greedy end
 struct Rotation end
 struct ZeroDual end
 struct PrimalOnly end
+struct ActivatedDual end
 
 function identifier(::DualScaling)
     "pdps_known_dualscaling"
@@ -92,6 +93,10 @@
     "pdps_known_primalonly"
 end
 
+function identifier(::ActivatedDual)
+    "pdps_known_activateddual"
+end
+
 #################################
 # Displacement field based flow
 #################################
@@ -199,6 +204,29 @@
     y .= C.*y
 end
 
+function activation(x)
+    # return (1/(1 + exp(-(x - 1.5))))                          # Works well with dual norm activation for lighthouse and brainphantom
+    # return 0.2/(1 + exp(-5(x-0.5))) + 0.3                     # Works well with dual norm activation for lighthouse and brainphantom
+    return  1.0 - 0.5/(1+ exp(-100(x-0.7))) 
+end
+
+# Experimental predictor for dual scaling based on activation
+function pdflow!(x, Δx, y, Δy, u, flow :: ActivatedDual; 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 = activation.(1 .- cc./ cm)
+    # C[1,:,:] .= c
+    # C[2,:,:] .= c
+    Δx .= activation.((sqrt.(abs.(y[1,:,:]).^2 + abs.(y[2,:,:]).^2))./0.25) 
+    C[1,:,:] .= Δx
+    C[2,:,:] .= Δx
+    y .= C.*y
+end
+
 function pdflow!(x, Δx, y, Δy, u, :: ZeroDual; threads=:none)
     @assert(size(u)==(2,))
     flow!(x, u; threads=(threads==:inner))
@@ -265,13 +293,32 @@
     @. x = Δx
     C = similar(y)
     cc = abs.(x-oldx)
-    cm = max(1e-12,maximum(cc))
-    c = 1 .* (1 .- cc./ cm) .^(10)
+    cm = max(flow.threshold,maximum(cc))
+    c = 1 .* (1 .- cc./ cm) .^flow.exponent
     C[1,:,:] .= c
     C[2,:,:] .= c
     y .= C.*y
 end
 
+# Experimental predictor for dual scaling based on activation
+function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: ActivatedDual; 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 = activation.(1 .- cc./ cm)
+    # C[1,:,:] .= c
+    # C[2,:,:] .= c
+    Δx .= activation.((sqrt.(abs.(y[1,:,:]).^2 + abs.(y[2,:,:]).^2))./0.25) 
+    C[1,:,:] .= Δx
+    C[2,:,:] .= Δx
+    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
diff -r c67ba8ef595e -r 70fd92ac9da0 src/PET/PET.jl
--- a/src/PET/PET.jl	Thu Apr 25 15:21:45 2024 -0500
+++ b/src/PET/PET.jl	Tue Apr 30 16:00:25 2024 +0300
@@ -122,6 +122,8 @@
                p_known₀_pets ⬿ (predictor=Rotation(),)),
     Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
                p_known₀_pets ⬿ (predictor=ZeroDual(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, shepplogan,
+               p_known₀_pets ⬿ (predictor=ActivatedDual(),)),
 )
 
 const brainphantom_experiments_pdps_known = (
@@ -139,6 +141,8 @@
                p_known₀_petb ⬿ (predictor=Rotation(),)),
     Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
                p_known₀_petb ⬿ (predictor=ZeroDual(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, brainphantom,
+               p_known₀_petb ⬿ (predictor=ActivatedDual(),)),
 )
 
 
diff -r c67ba8ef595e -r 70fd92ac9da0 src/PredictPDPS.jl
--- a/src/PredictPDPS.jl	Thu Apr 25 15:21:45 2024 -0500
+++ b/src/PredictPDPS.jl	Tue Apr 30 16:00:25 2024 +0300
@@ -183,6 +183,8 @@
                p_known₀_denoising ⬿ (predictor=Rotation(),)),
     Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
                p_known₀_denoising ⬿ (predictor=ZeroDual(),)),
+    Experiment(AlgorithmNew, DisplacementConstant, lighthouse,
+               p_known₀_denoising ⬿ (predictor=ActivatedDual(),)),
 )
 
 const denoising_experiments_all = Iterators.flatten((
diff -r c67ba8ef595e -r 70fd92ac9da0 src/Stats.jl
--- a/src/Stats.jl	Thu Apr 25 15:21:45 2024 -0500
+++ b/src/Stats.jl	Tue Apr 30 16:00:25 2024 +0300
@@ -42,7 +42,7 @@
 end
 
 const default_imnames = ("lighthouse200x300", "shepplogan256x256", "brainphantom256x256")
-const default_algnames = ("dualscaling", "greedy","noprediction","primalonly","proximal","rotation","zerodual")
+const default_algnames = ("dualscaling", "greedy","noprediction","primalonly","proximal","rotation","zerodual","activateddual")
 
 function calculate_statistics(;imnames=default_imnames,
                                algnames=default_algnames,
@@ -99,7 +99,7 @@
             end
         end
     end
-    sort!(results, [:experiment, :α])
+    sort!(results, [:experiment, :α, :algorithm])
     if isfile(csv_path)
     	rm(csv_path)
     end