Add RadonOp based on RadonKA

Author: nHackel

Project.toml

@@ -19,6 +19,7 @@ julia = "1.9"
 GPUArrays = "8, 9, 10"
 NFFT = "0.13"
 LinearOperators = "2.3.3"
+RadonKA = "0.6"
 Wavelets = "0.9, 0.10"
 Reexport = "1.0"
 FFTW = "1.0"
@@ -28,6 +29,7 @@ GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 NFFT = "efe261a4-0d2b-5849-be55-fc731d526b0d"
 Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+RadonKA = "86de8297-835b-47df-b249-c04e8db91db5"
 
 [targets]
 test = ["Test", "FFTW", "Wavelets", "NFFT"]
@@ -36,4 +38,5 @@ test = ["Test", "FFTW", "Wavelets", "NFFT"]
 LinearOperatorNFFTExt = ["NFFT", "FFTW"]
 LinearOperatorFFTWExt = "FFTW"
 LinearOperatorWaveletExt = "Wavelets"
-LinearOperatorGPUArraysExt = "GPUArrays"
+LinearOperatorGPUArraysExt = "GPUArrays"
+LinearOperatorRadonKAExt = "RadonKA"

ext/LinearOperatorRadonKAExt/LinearOperatorRadonKAExt.jl

@@ -0,0 +1,7 @@
+module LinearOperatorRadonKAExt
+
+using LinearOperatorCollection, RadonKA
+
+include("RadonOp.jl")
+
+end

ext/LinearOperatorRadonKAExt/RadonOp.jl

@@ -0,0 +1,46 @@
+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}
+  nrow :: Int
+  ncol :: Int
+  symmetric :: Bool
+  hermitian :: Bool
+  prod! :: Function
+  tprod! :: Nothing
+  ctprod! :: Function
+  nprod :: Int
+  ntprod :: Int
+  nctprod :: Int
+  args5 :: Bool
+  use_prod5! :: Bool
+  allocated5 :: Bool
+  Mv5 :: vecT
+  Mtu5 :: vecT
+  angles :: vecT2
+  geometry :: G
+  μ :: A
+end
+
+LinearOperators.storage_type(op::RadonOpImpl) = typeof(op.Mv5)
+
+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, true, false, true, 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; μ))

src/LinearOperatorCollection.jl

@@ -30,7 +30,7 @@ end
 
 export linearOperatorList, createLinearOperator
 export AbstractLinearOperatorFromCollection, WaveletOp, FFTOp, DCTOp, DSTOp, NFFTOp,
-       SamplingOp, NormalOp, WeightingOp, GradientOp
+       SamplingOp, NormalOp, WeightingOp, GradientOp, RadonOp
 
 abstract type AbstractLinearOperatorFromCollection{T} <: AbstractLinearOperator{T} end
 abstract type WaveletOp{T} <: AbstractLinearOperatorFromCollection{T} end
@@ -41,6 +41,7 @@ abstract type NFFTOp{T} <: AbstractLinearOperatorFromCollection{T} end
 abstract type SamplingOp{T} <: AbstractLinearOperatorFromCollection{T} end
 abstract type NormalOp{T} <: AbstractLinearOperatorFromCollection{T} end
 abstract type GradientOp{T} <: AbstractLinearOperatorFromCollection{T} end
+abstract type RadonOp{T} <: AbstractLinearOperatorFromCollection{T} end
 
 
 """