Mercurial > repos > PredictPDPS / file comparison
comparison: src/PET/ImGenerate.jl
src/PET/ImGenerate.jl
- changeset 36
- e4a8f662a1ac
- parent 35
- 74b1a9f0c35e
- child 37
- bba159cf1636
equal
deleted
inserted
replaced
| 35:74b1a9f0c35e | 36:e4a8f662a1ac |
|---|---|
| 1 ################### | |
| 2 # Image generation | |
| 3 ################### | |
| 4 | |
| 5 module ImGenerate | |
| 6 | |
| 7 using ColorTypes: Gray | |
| 8 import TestImages | |
| 9 # We don't really *directly* depend on QuartzImageIO. The import here is | |
| 10 # merely a workaround to suppress warnings when loading TestImages. | |
| 11 # Something is broken in those packages. | |
| 12 import QuartzImageIO | |
| 13 | |
| 14 using AlgTools.Util | |
| 15 using AlgTools.Comms | |
| 16 using ImageTools.Translate | |
| 17 | |
| 18 using ..OpticalFlow: Image, DisplacementConstant, DisplacementFull | |
| 19 using ..Radon | |
| 20 | |
| 21 # Added for reproducibility | |
| 22 import StableRNGs: StableRNG, Random | |
| 23 const rng = StableRNG(314159) | |
| 24 | |
| 25 # Added for PET | |
| 26 import PoissonRandom: pois_rand | |
| 27 import Random: shuffle | |
| 28 import Images: center, warp, imresize | |
| 29 import CoordinateTransformations: recenter | |
| 30 import Rotations: RotMatrix | |
| 31 import Interpolations: Flat | |
| 32 import MAT: matread | |
| 33 | |
| 34 | |
| 35 ############## | |
| 36 # Our exports | |
| 37 ############## | |
| 38 | |
| 39 export ImGen, | |
| 40 OnlineData, | |
| 41 imgen_square, | |
| 42 imgen_shake, | |
| 43 PetOnlineData, | |
| 44 imgen_shepplogan_radon, | |
| 45 imgen_brainphantom_radon | |
| 46 | |
| 47 ################## | |
| 48 # Data structures | |
| 49 ################## | |
| 50 | |
| 51 struct ImGen | |
| 52 f :: Function | |
| 53 dim :: Tuple{Int64,Int64} | |
| 54 Λ :: Float64 | |
| 55 dynrange :: Float64 | |
| 56 name :: String | |
| 57 end | |
| 58 | |
| 59 struct OnlineData{DisplacementT} | |
| 60 b_true :: Image | |
| 61 b_noisy :: Image | |
| 62 v :: DisplacementT | |
| 63 v_true :: DisplacementT | |
| 64 v_cumul_true :: DisplacementT | |
| 65 end | |
| 66 | |
| 67 struct PetOnlineData{DisplacementT} | |
| 68 b_true :: Image | |
| 69 sinogram_true :: Image | |
| 70 sinogram_noisy :: Image | |
| 71 v :: DisplacementT | |
| 72 v_true :: DisplacementT | |
| 73 v_cumul_true :: DisplacementT | |
| 74 theta :: DisplacementT # theta = thetaknown, theta_cumul | |
| 75 S :: Image | |
| 76 end | |
| 77 | |
| 78 ################### | |
| 79 # Shake generation | |
| 80 ################### | |
| 81 | |
| 82 function make_const_v(displ, sz) | |
| 83 v = zeros(2, sz...) | |
| 84 v[1, :, :] .= displ[1] | |
| 85 v[2, :, :] .= displ[2] | |
| 86 return v | |
| 87 end | |
| 88 | |
| 89 function shake(params) | |
| 90 if !haskey(params, :shaketype) || params.shaketype == :gaussian | |
| 91 return () -> params.shake.*randn(rng,2) | |
| 92 elseif params.shaketype == :disk | |
| 93 return () -> begin | |
| 94 θ = 2π*rand(rng,Float64) | |
| 95 r = params.shake*√(rand(rng,Float64)) | |
| 96 return [r*cos(θ), r*sin(θ)] | |
| 97 end | |
| 98 elseif params.shaketype == :circle | |
| 99 return () -> begin | |
| 100 θ = 2π*rand(rng,Float64) | |
| 101 r = params.shake | |
| 102 return [r*cos(θ), r*sin(θ)] | |
| 103 end | |
| 104 else | |
| 105 error("Unknown shaketype $(params.shaketype)") | |
| 106 end | |
| 107 end | |
| 108 | |
| 109 pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) | |
| 110 | |
| 111 | |
| 112 function rotatebytheta(params) | |
| 113 r = params.rotation_factor*randn(rng) | |
| 114 return r | |
| 115 end | |
| 116 | |
| 117 function generate_radonmask(params) | |
| 118 imdim = params.radondims | |
| 119 sino_sparse = params.sino_sparsity | |
| 120 numzero = Int64(round(sino_sparse*imdim[1]*imdim[2])) | |
| 121 numone = imdim[1]*imdim[2]-numzero | |
| 122 A = shuffle(rng,reshape([ones(numone); zeros(numzero)],(imdim[1],imdim[2]))) | |
| 123 return A | |
| 124 end | |
| 125 | |
| 126 ################ | |
| 127 # Moving square | |
| 128 ################ | |
| 129 | |
| 130 function generate_square(sz, | |
| 131 :: Type{DisplacementT}, | |
| 132 datachannel :: Channel{OnlineData{DisplacementT}}, | |
| 133 params) where DisplacementT | |
| 134 | |
| 135 if false | |
| 136 v₀ = make_const_v(0.1.*(-1, 1), sz) | |
| 137 nextv = () -> v₀ | |
| 138 elseif DisplacementT == DisplacementFull | |
| 139 nextv = pixelwise(shake(params), sz) | |
| 140 elseif DisplacementT == DisplacementConstant | |
| 141 nextv = shake(params) | |
| 142 else | |
| 143 @error "Invalid DisplacementT" | |
| 144 end | |
| 145 | |
| 146 # Constant linear displacement everywhere has Jacobian determinant one | |
| 147 # (modulo the boundaries which we ignore here) | |
| 148 m = round(Int, sz[1]/5) | |
| 149 b_orig = zeros(sz...) | |
| 150 b_orig[sz[1].-(2*m:3*m), 2*m:3*m] .= 1 | |
| 151 | |
| 152 v_true = nextv() | |
| 153 v_cumul = copy(v_true) | |
| 154 | |
| 155 while true | |
| 156 # Flow original data and add noise | |
| 157 b_true = zeros(sz...) | |
| 158 translate_image!(b_true, b_orig, v_cumul; threads=true) | |
| 159 b = b_true .+ params.noise_level.*randn(rng,sz...) | |
| 160 v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) | |
| 161 # Pass true data to iteration routine | |
| 162 data = OnlineData{DisplacementT}(b_true, b, v, v_true, v_cumul) | |
| 163 if !put_unless_closed!(datachannel, data) | |
| 164 return | |
| 165 end | |
| 166 # Next step shake | |
| 167 v_true = nextv() | |
| 168 v_cumul .+= v_true | |
| 169 end | |
| 170 end | |
| 171 | |
| 172 function imgen_square(sz) | |
| 173 return ImGen(curry(generate_square, sz), sz, 1, 1, "square$(sz[1])x$(sz[2])") | |
| 174 end | |
| 175 | |
| 176 ################ | |
| 177 # Shake a photo | |
| 178 ################ | |
| 179 | |
| 180 function generate_shake_image(im, sz, | |
| 181 :: Type{DisplacementConstant}, | |
| 182 datachannel :: Channel{OnlineData{DisplacementConstant}}, | |
| 183 params :: NamedTuple) | |
| 184 | |
| 185 | |
| 186 # Set up counter and zero factor for stabilisation | |
| 187 @assert(params.maxiter ≥ maximum(params.stable_interval)) | |
| 188 indx = 1 | |
| 189 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 | |
| 190 | |
| 191 # Restart the seed to enable comparison across predictors | |
| 192 Random.seed!(rng,314159) | |
| 193 | |
| 194 | |
| 195 nextv = shake(params) | |
| 196 v_true = zero_factor.*nextv() | |
| 197 v_cumul = copy(v_true) | |
| 198 | |
| 199 while true | |
| 200 # Extract subwindow of original image and add noise | |
| 201 b_true = zeros(sz...) | |
| 202 extract_subimage!(b_true, im, v_cumul; threads=true) | |
| 203 b = b_true .+ params.noise_level.*randn(rng,sz...) | |
| 204 v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) | |
| 205 # Pass data to iteration routine | |
| 206 data = OnlineData{DisplacementConstant}(b_true, b, v, v_true, v_cumul) | |
| 207 if !put_unless_closed!(datachannel, data) | |
| 208 return | |
| 209 end | |
| 210 # Next step shake | |
| 211 v_true = zero_factor.*nextv() | |
| 212 v_cumul .+= v_true | |
| 213 # Update indx and zero factor | |
| 214 indx += 1 | |
| 215 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 | |
| 216 end | |
| 217 end | |
| 218 | |
| 219 function imgen_shake(imname, sz) | |
| 220 im = Float64.(Gray.(TestImages.testimage(imname))) | |
| 221 dynrange = maximum(im) | |
| 222 return ImGen(curry(generate_shake_image, im, sz), sz, 1, dynrange, | |
| 223 "$(imname)$(sz[1])x$(sz[2])") | |
| 224 end | |
| 225 | |
| 226 | |
| 227 | |
| 228 ######################################################################## | |
| 229 # PETscan | |
| 230 ######################################################################## | |
| 231 function generate_sinogram(im, sz, | |
| 232 :: Type{DisplacementConstant}, | |
| 233 datachannel :: Channel{PetOnlineData{DisplacementConstant}}, | |
| 234 params :: NamedTuple) | |
| 235 | |
| 236 # Set up counter and zero factor for stabilisation | |
| 237 @assert(params.maxiter ≥ maximum(params.stable_interval)) | |
| 238 indx = 1 | |
| 239 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 | |
| 240 | |
| 241 # Restart the seed to enable comparison across predictors | |
| 242 Random.seed!(rng,314159) | |
| 243 | |
| 244 nextv = shake(params) | |
| 245 v_true = zero_factor.*nextv() | |
| 246 v_cumul = copy(v_true) | |
| 247 | |
| 248 S_true = generate_radonmask(params) | |
| 249 | |
| 250 theta_true = zero_factor*rotatebytheta(params) | |
| 251 theta_cumul = copy(theta_true) | |
| 252 | |
| 253 while true | |
| 254 # Define the transformation matrix | |
| 255 center_point = (sz[1]/2 + v_true[1], sz[2]/2 + v_true[2]) | |
| 256 tform = recenter(RotMatrix(theta_cumul), center_point) | |
| 257 | |
| 258 # Apply the transformation to the image using warp | |
| 259 b_true = copy(warp(im, tform, axes(im), fillvalue=Flat())) | |
| 260 | |
| 261 v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) | |
| 262 theta = theta_true*(1.0 + params.rotation_noise_level.*randn(rng)) | |
| 263 | |
| 264 # Generate the true and noisy sinograms | |
| 265 sinogram_true = zeros(params.radondims...) | |
| 266 sinogram_true .*= params.scale | |
| 267 radon!(sinogram_true, b_true) | |
| 268 sinogram_noisy = copy(sinogram_true) | |
| 269 | |
| 270 for i=1:params.radondims[1], j=1:params.radondims[2] | |
| 271 sinogram_noisy[i, j] += pois_rand(rng,params.noise_level) | |
| 272 end | |
| 273 | |
| 274 # Pass data to iteration routine | |
| 275 data = PetOnlineData{DisplacementConstant}(b_true, sinogram_true, sinogram_noisy, v, v_true, v_cumul, [theta, theta_cumul], S_true) | |
| 276 if !put_unless_closed!(datachannel, data) | |
| 277 return | |
| 278 end | |
| 279 | |
| 280 # Next step shake | |
| 281 v_true = zero_factor.*nextv() | |
| 282 v_cumul .+= v_true | |
| 283 # Next theta | |
| 284 theta_true = zero_factor*rotatebytheta(params) | |
| 285 theta_cumul += theta_true | |
| 286 # Next sinogram mask | |
| 287 S_true = generate_radonmask(params) | |
| 288 # Update indx and zero factor | |
| 289 indx += 1 | |
| 290 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 | |
| 291 end | |
| 292 end | |
| 293 | |
| 294 function imgen_shepplogan_radon(sz) | |
| 295 im = convert(Array{Float64},TestImages.shepp_logan(sz[1], highContrast=true)) | |
| 296 dynrange = maximum(im) | |
| 297 return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])") | |
| 298 end | |
| 299 | |
| 300 function imgen_brainphantom_radon(sz) | |
| 301 data = matread("src/PET/phantom_slice.mat") | |
| 302 im = normalise(imresize(convert(Array{Float64},data["square_data"]),sz)) | |
| 303 dynrange = maximum(im) | |
| 304 return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "brainphantom$(sz[1])x$(sz[2])") | |
| 305 end | |
| 306 | |
| 307 normalise = (data) -> data./maximum(data) | |
| 308 end # Module |
