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