LinearOperatorCollection.jl/ext/LinearOperatorRadonKAExt/RadonOp.jl at e690a25609fa4b1b1a6f7d234d4a2fb1c7b5c5ec · JuliaImageRecon/LinearOperatorCollection.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
"""
    RadonOp(::Type{T}; shape::NTuple{N, Int}, angles, geometry = RadonParallelCircle(shape[1], -(shape[1]-1)÷2:(shape[1]-1)÷2), μ = nothing, S = Vector{T}) where {T, N}

Generates a `RadonOp` which evaluates the Radon transform operator and its adjoint (backprojection) for a given geometry and projection angles.

# Arguments:
* `T`                       - element type for the operator (e.g., `Float64`, `ComplexF32`)
* `shape::NTuple{N, Int}`   - size of the image
* `angles`                  - array of projection angles
* `geometry`                - Radon geometry descriptor (default: parallel beam circle)
* `μ`                       - optional attenuation map (for attenuated Radon transform)
* `S`                       - storage type for internal vectors (default: `Vector{T}`)
"""
function LinearOperatorCollection.RadonOp(::Type{T}; shape::NTuple{N, Int}, angles,
   geometry = RadonParallelCircle(shape[1], -(shape[1]-1)÷2:(shape[1]-1)÷2), μ = nothing, S = Vector{T}) where {T, N}
  return RadonOpImpl(T; shape, angles, geometry, μ, S)
end

mutable struct RadonOpImpl{T, vecT <: AbstractVector{T}, vecT2, G, A} <: RadonOp{T}
  const nrow :: Int
  const ncol :: Int
  const symmetric :: Bool
  const hermitian :: Bool
  const prod! :: Function
  const tprod! :: Nothing
  const ctprod! :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  Mv :: vecT
  Mtu :: vecT
  const angles :: vecT2
  const geometry :: G
  const μ :: A
end

LinearOperators.storage_type(::RadonOpImpl{T, vecT}) where {T, vecT} = vecT

function RadonOpImpl(T::Type; shape::NTuple{N, Int64}, angles, geometry, μ, S) where N
  N_sinogram = length(geometry.in_height)
  N_angles = length(angles)
  d = length(shape) == 3 ? shape[3] : 1
  nrow = N_sinogram * N_angles * d
  ncol = prod(shape)
  return RadonOpImpl(nrow, ncol, false, false,
  (res, x) -> prod_radon!(res, x, shape, angles, geometry, μ),
  nothing,
  (res, x) -> ctprod_radon!(res, x, (N_sinogram, N_angles, d), angles, geometry, μ),
  0, 0, 0, S(undef, 0), S(undef, 0), angles, geometry, μ)
end

prod_radon!(res::vecT, x::vecT, shape, angles::vecT2, geometry::G, μ::A) where {vecT, vecT2, G, A} = copyto!(res, radon(reshape(x, shape), angles; geometry, μ))
prod_radon!(res::vecT, x::vecT, shape, angles::vecT2, ::Nothing, μ::A) where {vecT, vecT2, A} = copyto!(res, radon(reshape(x, shape), angles; μ))

ctprod_radon!(res::vecT, x::vecT, shape, angles::vecT2, geometry::G, μ::A) where {vecT, vecT2, G, A} = copyto!(res, RadonKA.backproject(reshape(x, shape), angles; geometry, μ))
ctprod_radon!(res::vecT, x::vecT, shape, angles::vecT2, ::Nothing, μ::A) where {vecT, vecT2, A} = copyto!(res, RadonKA.backproject(reshape(x, shape), angles; μ))