Mercurial > repos > PredictPDPS / changeset

--- a/src/ImGenerate.jl	Sun Apr 21 19:05:44 2024 +0300
+++ b/src/ImGenerate.jl	Sun Apr 21 19:26:50 2024 +0300
@@ -25,10 +25,11 @@
 # Added for PET
 import PoissonRandom: pois_rand
 import Random: shuffle
-import Images: center, warp
+import Images: center, warp, imresize
 import CoordinateTransformations: recenter
 import Rotations: RotMatrix
 import Interpolations: Flat
+import MAT: matread
 
 
 ##############
@@ -40,7 +41,8 @@
        imgen_square,
        imgen_shake,
        PetOnlineData,
-       imgen_shepplogan_radon
+       imgen_shepplogan_radon,
+       imgen_brainphantom_radon
 
 ##################
 # Data structures
@@ -108,15 +110,15 @@
 
 
 function rotatebytheta(params)
-    r = params.rotation_factor*(2*rand(rng)-1)
+    r = params.rotation_factor*randn(rng)
     return r
 end
 
 function generate_radonmask(params)
     imdim = params.radondims
     sino_sparse = params.sino_sparsity
-    numone = Int64(round(sino_sparse*imdim[1]*imdim[2]))
-    numzero = imdim[1]*imdim[2]-numone
+    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
@@ -180,11 +182,17 @@
                               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,9182737465)   
 
     nextv = shake(params)
-    v_true = nextv()
+    v_true = zero_factor.*nextv()
     v_cumul = copy(v_true)
 
     while true
@@ -199,8 +207,11 @@
             return
         end
         # Next step shake
-        v_true = nextv()
+        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
 
@@ -217,20 +228,25 @@
 # PETscan
 ########################################################################
 function generate_radon(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,9182737465)  
+    Random.seed!(rng,9182737465)
 
     nextv = shake(params)
-    v_true = nextv()
+    v_true = zero_factor.*nextv()
     v_cumul = copy(v_true)
 
     S_true = generate_radonmask(params)
 
-    theta_true = rotatebytheta(params)
+    theta_true = zero_factor*rotatebytheta(params)
     theta_cumul = copy(theta_true)
 
     while true
@@ -261,22 +277,31 @@
         end
 
         # Next step shake
-        v_true = nextv()
+        v_true = zero_factor.*nextv()
         v_cumul .+= v_true
         # Next theta
-        theta_true = rotatebytheta(params)
+        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(origsize,sz)
-    im = convert(Array{Float64},TestImages.shepp_logan(origsize, highContrast=true))
+function imgen_shepplogan_radon(sz)
+    im = convert(Array{Float64},TestImages.shepp_logan(sz[1], highContrast=true))
     dynrange = maximum(im)
     return ImGen(curry(generate_radon, 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_radon, im, sz), sz, 1, dynrange, "brainphantom$(sz[1])x$(sz[2])")
+end
 
+normalise = (data) -> data./maximum(data)
 end # Module

--- a/src/PET/ImGenerate.jl	Sun Apr 21 19:05:44 2024 +0300
+++ b/src/PET/ImGenerate.jl	Sun Apr 21 19:26:50 2024 +0300
@@ -182,11 +182,17 @@
                               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)   
+    Random.seed!(rng,9182737465)   
 
     nextv = shake(params)
-    v_true = nextv()
+    v_true = zero_factor.*nextv()
     v_cumul = copy(v_true)
 
     while true
@@ -201,8 +207,11 @@
             return
         end
         # Next step shake
-        v_true = nextv()
+        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
 

--- a/src/PET/PET.jl	Sun Apr 21 19:05:44 2024 +0300
+++ b/src/PET/PET.jl	Sun Apr 21 19:26:50 2024 +0300
@@ -111,7 +111,7 @@
     shake_noise_level = 0.1,
     shake = 1.0, 
     rotation_factor = 0.075,                  
-    rotation_noise_level = 0.00075,         
+    rotation_noise_level = 0.0075,         
     α = 1.0,
     ρ̃₀ = 1.0,
     σ̃₀ = 1.0,

--- a/src/PredictPDPS.jl	Sun Apr 21 19:05:44 2024 +0300
+++ b/src/PredictPDPS.jl	Sun Apr 21 19:26:50 2024 +0300
@@ -138,6 +138,8 @@
     handle_interrupt = true,
     init = :zero,
     plot_movement = false,
+    # stable_interval = Set(),
+    stable_interval = union(Set(1000:2000),Set(3000:4000),Set(5000:7000),Set(9000:10000)),
 )
 
 const square = imgen_square((200, 300))