# HG changeset patch # User Neil Dizon # Date 1713716810 -10800 # Node ID a30a8745a83bde5941caa8229a15ac09d1437c94 # Parent 97737e4e71973cff53567ff92897aa907f31039e added stable interval for denoising diff -r 97737e4e7197 -r a30a8745a83b src/ImGenerate.jl --- 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 diff -r 97737e4e7197 -r a30a8745a83b src/PET/ImGenerate.jl --- 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 diff -r 97737e4e7197 -r a30a8745a83b src/PET/PET.jl --- 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, diff -r 97737e4e7197 -r a30a8745a83b src/PredictPDPS.jl --- 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))