Update to AbstractImageReconstruction v0.6

Project.toml

@@ -1,15 +1,15 @@
 name = "DaggerImageReconstruction"
 uuid = "b99085f1-2f43-45a7-bd38-d511f4bff3b1"
 authors = ["nHackel <nHackel@users.noreply.github.com> and contributors"]
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 AbstractImageReconstruction = "a4b4fdbf-6459-4ec9-990d-77e1fa24a91b"
 Dagger = "d58978e5-989f-55fb-8d15-ea34adc7bf54"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [compat]
-AbstractImageReconstruction = "0.4, 0.5"
+AbstractImageReconstruction = "0.6"
 Dagger = "0.18, 0.19"
 Distributed = "1"
 julia = "1.10"

docs/make.jl

@@ -31,6 +31,7 @@ makedocs(
             "Introduction" => "example_intro.md",
             "Algoritm Interface" => "generated/example/algorithm.md",
             "RecoPlan Interface" => "generated/example/daggerplan.md",
+            "Parameter Interface" => "generated/example/utility.md"
         ],
         "API Reference" => "API/api.md",
 

docs/src/API/api.md

@@ -6,6 +6,7 @@ REPL one can access this documentation by entering the help mode with `?`
 ```@docs
 DaggerImageReconstruction.DaggerReconstructionAlgorithm
 DaggerImageReconstruction.DaggerReconstructionParameter
+DaggerImageReconstruction.DaggerReconstructionUtility
 ```
 
 ## DaggerRecoPlan

docs/src/literate/example/0_radon_data.jl

@@ -3,7 +3,7 @@
 # In this example we will set up our radon data using RadonKA.jl, ImagePhantoms.jl and ImageGeoms.jl. We will start with simple 2D images and their sinograms and continue with a time series of 3D images and sinograms.
 
 # ## Background
-# The Radon transform is an integral transform that projects the values of a function(or a phantom) along straight lines onto a detector.
+# The Radon transform is an integral transform that projects the values of a function (or a phantom) along straight lines onto a detector.
 # These projections are recorded for a number of different angles and form the so-called sinogram. The Radon transform and its adjoint form the mathematical basis
 # for computed tomography (CT) and other imaging modalities such as single photon emission computed tomography (SPECT) and positron emission tomography (PET).
 
@@ -18,7 +18,7 @@ size(image)
 
 # This produces a 256x256 image of a Shepp-Logan phantom. Next, we will generate the Radon data using `radon` from RadonKA.jl.
 # The arguments of this function are the image or phantom to be transformed, the angles at which the projections are taken, and the used geometry of the system. For this example we will use the default parallel circle geometry. 
-# For more details, we refer to the RadonKA.jl documentation. We will use 256 angles for the projections, between 0 and π.
+# For more details, we refer to the RadonKA.jl documentation. We will use 256 angles for the projections between 0 and π.
 using RadonKA
 angles = collect(range(0, π, 256))
 sinogram = Array(RadonKA.radon(image, angles))
@@ -39,8 +39,8 @@ plot_image(fig[1,2], sinogram, title = "Sinogram")
 resize_to_layout!(fig)
 fig
 
-# ## 3D Pnantom
-# RadonKA.jl also supports 3D Radon transforms. The first two dimensions are interpreted as the XY plane where the transform applied and the last dimensions is the rotational axis z of the projections.
+# ## 3D Phantom
+# RadonKA.jl also supports 3D Radon transforms. The first two dimensions are interpreted as the XY plane where the transform is applied and the last dimension is the rotational axis z of the projections.
 # For that we need to create a 3D Shepp-Logan phantom. First we retrieve the parameters of the ellipsoids of the Shepp-Logan phantom:
 shape = (64, 64, 64)
 params = map(collect, ellipsoid_parameters(; fovs = shape));
@@ -74,7 +74,7 @@ fig
 
 
 # ## Time Series of 3D Phantoms
-# Lastly, we want to add a time dimension to our 3D phantom. For our example we will increase the intensity of the third ellipsoid every time step or frame.
+# Lastly, we want to add a time dimension to our 3D phantom. For our example, we will increase the intensity of the third ellipsoid every time step or frame.
 images = similar(image, size(image)..., 5)
 sinograms = similar(sinogram, size(sinogram)..., 5)
 for (i, intensity) in enumerate(range(params[3][end], 0.3, 5))

docs/src/literate/example/1_interface.jl

@@ -26,24 +26,38 @@ abstract type AbstractDirectRadonAlgorithm <: AbstractRadonAlgorithm end
 abstract type AbstractIterativeRadonAlgorithm <: AbstractRadonAlgorithm end
 
 # ## Internal Interface
-# Reconstruction algorithms in AbstractImageReconstruction.jl are expected to be implemented in the form of distinct processing steps, implemented in their own `process` methods.
-# The `process` function takes an algorithm, parameters, and inputs and returns the result of the processing step.
-# If no function is defined for an instance of an algorithm, the default implementation is called. This method tries to call the function `process` with the type of the algorithm:
+# In AbstractImageReconstruction.jl, reconstruction algorithms are driven by *parameters*. Parameters are callable objects that implement individual processing steps.
+# A parameter type `MyParams` is expected to implement one of:
 # ```julia
-# process(algo::AbstractImageReconstructionAlgorithm, param::AbstractImageReconstructionParameters, inputs...) = process(typeof(algo), param, inputs...)
+#   (param::MyParams)(::Type{<:MyAlgorithm}, inputs...)
+#   (param::MyParams)(algo::MyAlgorithm,      inputs...)
 # ```
-# The implementation of reconstruction algorithms is therefore expected to either implement the `process` function for the algorithm type or for the instance. Dispatch on instances allow an instance to change its state, while dispatch on types allows for pure helper functions.
+# The type-based variant is preferred for pure functions; the instance-based variant allows mutation of the algorithm state. The default implementation of
+# ```julia
+#   (param::AbstractImageReconstructionParameters)(algo::AbstractImageReconstructionAlgorithm, inputs...) = param(algo, inputs...) → param(typeof(algo), inputs...)
+# ```
+# simply forwards to the type-based method.
+
+# A reconstruction algorithm typically stores a *main* parameter. Multiple processing steps can be encoded by
+# composing parameter calls; there is no requirement to implement a strict linear pipeline.
 
-# A `process` itself can invoke other `process` functions to enable multiple processing steps and generally have arbitarry control flow. It is not required to implement a straight-foward pipeline. We will see this later when we implementd our algorithms.
+# To extend an existing algorithm with new behavior, it's enough to implement new parameters or potentially add an algorithm.
+# Later on, we will see more infrastructure of the package which focuses on parameters and their Configuration.
 
-# Let's define a preprocessing step that we can share between our algorithms. We want to allow the user to select certain frames from a time series and average them.
-# We will use the `@kwdef` macro to provide constructor with keyword arguments and default values
+# Let's define a preprocessing step that we can share between our algorithms. We want to
+# allow the user to select certain frames from a time series and average them. We will use
+# the `@parameter` macro. This is similar to `Base.@kwdef` and allows us to provide a constructor with keyword arguments and default values.
+# It also allows us to validate the values of our parameters:
 using Statistics
-Base.@kwdef struct RadonPreprocessingParameters <: AbstractRadonPreprocessingParameters
+@parameter struct RadonPreprocessingParameters <: AbstractRadonPreprocessingParameters
   frames::Vector{Int64} = []
   numAverages::Int64 = 1
+
+  @validate begin
+    @assert numAverages >= 0 "Averages must be a positive integer"
+  end
 end
-function AbstractImageReconstruction.process(::Type{<:AbstractRadonAlgorithm}, params::RadonPreprocessingParameters, data::AbstractArray{T, 4}) where {T}
+function (params::RadonPreprocessingParameters)(::Type{<:AbstractRadonAlgorithm}, data::AbstractArray{T, 4}) where {T}
   frames = isempty(params.frames) ? (1:size(data, 4)) : params.frames
   data = data[:, :, :, frames]
 
@@ -57,4 +71,5 @@ end
 
 # ## User Interface
 # A user of our package should be able to reconstruct images by calling the `reconstruct` function. This function takes an algorithm and an input and returns the reconstructed image.
-# Internally, the `reconstruct` function calls the `put!` and `take!` functions of the algorithm to pass the input and retrieve the output. Algorithms must implement these functions and are expected to have FIFO behavior.
+# Internally, the `reconstruct` function calls the `put!` and `take!` functions of the algorithm to pass the input and retrieve the output. Algorithms must implement these functions and are expected to have FIFO behavior.
+# However, much of this boilerplate can be created via macros, as we will see in this example.

docs/src/literate/example/2_direct.jl

@@ -6,67 +6,51 @@ export AbstractDirectRadonReconstructionParameters, RadonFilteredBackprojectionP
 # To implement our direct reconstruction algorithms we need to define a few more methods and types. We will start by defining the parameters for the backprojection and for the filtered backprojection. Afterwards we can implement the algorithm itself.
 
 # ## Parameters and Processing
-# For convenience we first introduce a new abstract type for the direct reconstruction paramters:
+# For convenience we first introduce a new abstract type for the direct reconstruction parameters:
 abstract type AbstractDirectRadonReconstructionParameters <: AbstractRadonReconstructionParameters end
 # The backprojection parameters are simple and only contain the number of angles:
-Base.@kwdef struct RadonBackprojectionParameters <: AbstractDirectRadonReconstructionParameters
+@parameter struct RadonBackprojectionParameters <: AbstractDirectRadonReconstructionParameters
   angles::Vector{Float64}
 end
 
 # The filtered backprojection parameters are more complex and contain the number of angles and optionally the filter which should be used:
-Base.@kwdef struct RadonFilteredBackprojectionParameters <: AbstractDirectRadonReconstructionParameters
+@parameter struct RadonFilteredBackprojectionParameters <: AbstractDirectRadonReconstructionParameters
   angles::Vector{Float64}
   filter::Union{Nothing, Vector{Float64}} = nothing
 end
 # Since we have defined no default values for the angles, they are required to be set by the user. A more advanced implementation would also allow for the geometry to be set.
 
 # Next we will implement the process steps for both of our backprojection variants. Since RadonKA.jl expects 2D or 3D arrays we have to transform our time series accordingly.
-function AbstractImageReconstruction.process(algoT::Type{<:AbstractDirectRadonAlgorithm}, params::AbstractDirectRadonReconstructionParameters, data::AbstractArray{T, 4}) where {T}
+function (params::AbstractDirectRadonReconstructionParameters)(algoT::Type{<:AbstractDirectRadonAlgorithm}, data::AbstractArray{T, 4}) where {T}
   result = []
   for i = 1:size(data, 4)
-    push!(result, process(algoT, params, view(data, :, :, :, i)))
+    push!(result, params(algoT, view(data, :, :, :, i)))
   end
   return cat(result..., dims = 4)
 end
-AbstractImageReconstruction.process(::Type{<:AbstractDirectRadonAlgorithm}, params::RadonBackprojectionParameters, data::AbstractArray{T, 3}) where {T} = RadonKA.backproject(data, params.angles)
-AbstractImageReconstruction.process(::Type{<:AbstractDirectRadonAlgorithm}, params::RadonFilteredBackprojectionParameters, data::AbstractArray{T, 3}) where {T} = RadonKA.backproject_filtered(data, params.angles; filter = params.filter)
+(params::RadonBackprojectionParameters)(::Type{<:AbstractDirectRadonAlgorithm}, data::AbstractArray{T, 3}) where {T} = RadonKA.backproject(data, params.angles)
+(params::RadonFilteredBackprojectionParameters)(::Type{<:AbstractDirectRadonAlgorithm}, data::AbstractArray{T, 3}) where {T} = RadonKA.backproject_filtered(data, params.angles; filter = params.filter)
 
 # ## Algorithm
 # The direct reconstruction algorithm has essentially no state to store between reconstructions and thus only needs its parameters as fields. We want our algorithm to accept any combination of our preprocessing and direct reconstruction parameters.
 # This we encode in a new type:
-Base.@kwdef struct DirectRadonParameters{P <: AbstractRadonPreprocessingParameters, R <: AbstractDirectRadonReconstructionParameters} <: AbstractRadonParameters
+@chain struct DirectRadonParameters{P <: AbstractRadonPreprocessingParameters, R <: AbstractDirectRadonReconstructionParameters} <: AbstractRadonParameters
   pre::P
   reco::R
 end
-# And the according processing step:
-function AbstractImageReconstruction.process(algoT::Type{<:AbstractDirectRadonAlgorithm}, params::DirectRadonParameters{P, R}, data::AbstractArray{T, 4}) where {T, P<:AbstractRadonPreprocessingParameters, R<:AbstractDirectRadonReconstructionParameters}
-  data = process(algoT, params.pre, data)
-  return process(algoT, params.reco, data)
-end
+# This composite type simply chains the result of our preprocessing step into our reco step. See the API for @chain for more information.
 
-# Now we can define the algorithm type itself. Algorithms are usually constructed with one argument passing in the user parameters:
-mutable struct DirectRadonAlgorithm{D <: DirectRadonParameters} <: AbstractDirectRadonAlgorithm
-  parameter::D
-  output::Channel{Any}
-  DirectRadonAlgorithm(parameter::D) where D = new{D}(parameter, Channel{Any}(Inf))
-end
-# And they implement a method to retrieve the used parameters:
-AbstractImageReconstruction.parameter(algo::DirectRadonAlgorithm) = algo.parameter
+# Now we can define the algorithm type itself using the `@reconstruction` macro. The macro automatically generates:
 
-# Algorithms are assumed to be stateful. To ensure thread safety, we need to implement the `lock` and `unlock` functions. We will use the `output` channel as a lock:
-Base.lock(algo::DirectRadonAlgorithm) = lock(algo.output)
-Base.unlock(algo::DirectRadonAlgorithm) = unlock(algo.output)
+# - A mutable struct with the parameter field and internal infrastructure
+# - A constructor accepting the parameter
+# - Interface methods: `Base.put!`, `Base.take!`, `Base.wait`, `Base.isready`, `Base.lock`, `Base.unlock`
+# - A parameter accessor method
 
-# And implement the `put!` and `take!` functions, mimicking the behavior of a FIFO channel for reconstructions:
-Base.take!(algo::DirectRadonAlgorithm) = Base.take!(algo.output)
-function Base.put!(algo::DirectRadonAlgorithm, data::AbstractArray{T, 4}) where {T} 
-  lock(algo) do
-    put!(algo.output, process(algo, algo.parameter, data))
-  end
-end
+# It's also possible to implement these functions manually, see the corresponding docstrings for information
 
-# The way the behaviour is implemented here, the algorithm does not buffer any inputs and instead blocks until the currenct reconstruction is done. Outputs are stored until they are retrieved.
+@reconstruction struct DirectRadonAlgorithm{D <: DirectRadonParameters} <: AbstractDirectRadonAlgorithm
+  @parameter parameter::D
+end
 
-# With `wait` and `isready` we can check if the algorithm is currently processing data or if it is ready to accept new inputs:
-Base.wait(algo::DirectRadonAlgorithm) = wait(algo.output)
-Base.isready(algo::DirectRadonAlgorithm) = isready(algo.output)
+# The way the behaviour is implemented here, the algorithm does not buffer any inputs and instead blocks until the current reconstruction is done. Outputs are stored until they are retrieved.