diff -r 6ce506e8ab02 -r e6a2ebaf467d src/PET/PET.jl
--- a/src/PET/PET.jl	Sun Apr 21 13:44:02 2024 +0300
+++ b/src/PET/PET.jl	Sun Apr 21 18:45:09 2024 +0300
@@ -39,6 +39,8 @@
 include("AlgorithmPrimalOnly.jl")
 include("AlgorithmProximal.jl")
 include("AlgorithmRotation.jl")
+include("AlgorithmZeroDual.jl")
+include("PlotResults.jl")
 
 using .Radon 
 using .ImGenerate
@@ -49,13 +51,20 @@
 using .AlgorithmPrimalOnly
 using .AlgorithmProximal
 using .AlgorithmRotation
+using .AlgorithmZeroDual
+using .PlotResults
 
 
 ##############
 # Our exports
 ##############
 
-export demo_pet1, demo_pet2, demo_pet3, demo_pet4, demo_pet5, demo_pet6, batchrun_pet
+export demo_petS1, demo_petS2, demo_petS3, 
+       demo_petS4, demo_petS5, demo_petS6, demo_petS7,
+       demo_petB1, demo_petB2, demo_petB3, 
+       demo_petB4, demo_petB5, demo_petB6, demo_petB7,
+       batchrun_shepplogan, batchrun_brainphantom, batchrun_pet,
+       plot_pet
 
 ###################################
 # Parameterisation and experiments
@@ -109,8 +118,7 @@
     δ = 0.9,
     σ₀ = 0.5,
     τ₀ = 0.9,
-    λ = 1,
-    origsize = 256,             
+    λ = 1,            
     radondims = [128,64],
     sz = (256,256),
     scale = 1,
@@ -118,28 +126,54 @@
     sino_sparsity = 0.5,
     L = 300.0,
     L_experiment = false,
-    stable_interval = union(Set(1500:200),Set(3500:4000)),
+    stable_interval = union(Set(1500:2000),Set(3500:4000)),
+)
+
+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(AlgorithmProximal, DisplacementConstant, shepplogan,
+               p_known₀_pet ⬿ (phantom_ρ = 100,)),
+    Experiment(AlgorithmRotation, DisplacementConstant, shepplogan,
+               p_known₀_pet),   
+    Experiment(AlgorithmZeroDual, DisplacementConstant, shepplogan,
+               p_known₀_pet),                       
 )
 
-const petscan = imgen_shepplogan_radon(p_known₀_pet.sz)
-
-const pet_experiments_pdps_known = (
-    Experiment(AlgorithmDualScaling, DisplacementConstant, petscan,
+const brainphantom_experiments_pdps_known = (
+    Experiment(AlgorithmDualScaling, DisplacementConstant, brainphantom,
                p_known₀_pet),
-    Experiment(AlgorithmGreedy, DisplacementConstant, petscan,
+    Experiment(AlgorithmGreedy, DisplacementConstant, brainphantom,
                p_known₀_pet),    
-    Experiment(AlgorithmNoPrediction, DisplacementConstant, petscan,
+    Experiment(AlgorithmNoPrediction, DisplacementConstant, brainphantom,
                p_known₀_pet),          
-    Experiment(AlgorithmPrimalOnly, DisplacementConstant, petscan,
+    Experiment(AlgorithmPrimalOnly, DisplacementConstant, brainphantom,
                p_known₀_pet), 
-    Experiment(AlgorithmProximal, DisplacementConstant, petscan,
+    Experiment(AlgorithmProximal, DisplacementConstant, brainphantom,
                p_known₀_pet ⬿ (phantom_ρ = 100,)),
-    Experiment(AlgorithmRotation, DisplacementConstant, petscan,
-               p_known₀_pet),                         
+    Experiment(AlgorithmRotation, DisplacementConstant, brainphantom,
+               p_known₀_pet),   
+    Experiment(AlgorithmZeroDual, DisplacementConstant, brainphantom,
+               p_known₀_pet),                       
 )
 
-const pet_experiments_all = Iterators.flatten((
-    pet_experiments_pdps_known,
+
+const shepplogan_experiments_all = Iterators.flatten((
+    shepplogan_experiments_pdps_known,
+))
+
+const brainphantom_experiments_all = Iterators.flatten((
+    brainphantom_experiments_pdps_known,
 ))
 
 ################
@@ -312,16 +346,25 @@
                     kwargs...)
 end
 
-demo_pet1 = () -> demo(pet_experiments_pdps_known[1]) # Dual scaling
-demo_pet2 = () -> demo(pet_experiments_pdps_known[2]) # Greedy
-demo_pet3 = () -> demo(pet_experiments_pdps_known[3]) # No Prediction
-demo_pet4 = () -> demo(pet_experiments_pdps_known[4]) # Primal only
-demo_pet5 = () -> demo(pet_experiments_pdps_known[5]) # Proximal (old)
-demo_pet6 = () -> demo(pet_experiments_pdps_known[6]) # Rotation
+demo_petS1 = () -> demo(shepplogan_experiments_pdps_known[1]) # Dual scaling
+demo_petS2 = () -> demo(shepplogan_experiments_pdps_known[2]) # Greedy
+demo_petS3 = () -> demo(shepplogan_experiments_pdps_known[3]) # No Prediction
+demo_petS4 = () -> demo(shepplogan_experiments_pdps_known[4]) # Primal only
+demo_petS5 = () -> demo(shepplogan_experiments_pdps_known[5]) # Proximal (old)
+demo_petS6 = () -> demo(shepplogan_experiments_pdps_known[6]) # Rotation
+demo_petS7 = () -> demo(shepplogan_experiments_pdps_known[7]) # Zero dual
+
+demo_petB1 = () -> demo(brainphantom_experiments_pdps_known[1]) # Dual scaling
+demo_petB2 = () -> demo(brainphantom_experiments_pdps_known[2]) # Greedy
+demo_petB3 = () -> demo(brainphantom_experiments_pdps_known[3]) # No Prediction
+demo_petB4 = () -> demo(brainphantom_experiments_pdps_known[4]) # Primal only
+demo_petB5 = () -> demo(brainphantom_experiments_pdps_known[5]) # Proximal (old)
+demo_petB6 = () -> demo(brainphantom_experiments_pdps_known[6]) # Rotation
+demo_petB7 = () -> demo(brainphantom_experiments_pdps_known[7]) # Zero dual
 
 
-function batchrun_pet(kwargs...)
-    run_experiments(;experiments=pet_experiments_all,
+function batchrun_shepplogan(kwargs...)
+    run_experiments(;experiments=shepplogan_experiments_all,
                     save_results=true,
                     save_images=true,
                     visualise=false,
@@ -329,6 +372,20 @@
                     kwargs...)
 end
 
+function batchrun_brainphantom(kwargs...)
+    run_experiments(;experiments=brainphantom_experiments_all,
+                    save_results=true,
+                    save_images=true,
+                    visualise=false,
+                    recalculate=false,
+                    kwargs...)
+end
+
+function batchrun_pet(kwargs...)
+    batchrun_shepplogan()
+    batchrun_brainphantom()
+end
+
 ######################################################
 # Iterator that does visualisation and log collection
 ######################################################
@@ -336,11 +393,8 @@
 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
 
@@ -496,6 +550,19 @@
     return apply.(im)
 end
 
+
+#########################
+# Plotting SSIM and PSNR
+#########################
+
+function plot_pet(kwargs...)
+    ssim_plot("shepplogan")
+    ssim_plot("brainphantom")
+    psnr_plot("shepplogan")
+    psnr_plot("brainphantom")
+end
+
+
 ###############
 # Precompiling
 ###############