diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ab4df39..eff9a2f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -48,3 +48,25 @@ jobs: ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 + + docs: + name: Documentation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: '1' + - run: | + julia --project=docs -e ' + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' + - run: | + julia --project=docs -e ' + using Documenter: DocMeta, doctest + using AbstractImageReconstruction' + - run: julia --project=docs docs/make.jl + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file diff --git a/.gitignore b/.gitignore index b067edd..49c4447 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,10 @@ +*.jl.cov +*.jl.*.cov +*.jl.mem +docs/build/ +docs/site/ +docs/src/generated/ + +Manifest.toml + /Manifest.toml diff --git a/README.md b/README.md index 4c92640..d3cf344 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,9 @@ [![Build Status](https://github.com/JuliaImageRecon/AbstractImageReconstruction.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaImageRecon/AbstractImageReconstruction.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![](https://img.shields.io/badge/docs-latest-blue.svg)](https://JuliaImageRecon.github.io/AbstractImageReconstruction.jl/latest) + + This package contains an interface and type hierarchy for image reconstruction algorithms and their parameters, together with associated utility tools. ## Installation diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..aa5ebf1 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,14 @@ +[deps] +AbstractImageReconstruction = "a4b4fdbf-6459-4ec9-990d-77e1fa24a91b" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ImageGeoms = "9ee76f2b-840d-4475-b6d6-e485c9297852" +ImagePhantoms = "71a99df6-f52c-4da1-bd2a-69d6f37f3252" +LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +RadonKA = "86de8297-835b-47df-b249-c04e8db91db5" +RegularizedLeastSquares = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23" + +[compat] +Documenter = "1" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..a4dcf46 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,51 @@ +using Documenter, Literate, AbstractImageReconstruction + +# Generate examples +OUTPUT_BASE = joinpath(@__DIR__(), "src/generated") +INPUT_BASE = joinpath(@__DIR__(), "src/literate") +for (_, dirs, _) in walkdir(INPUT_BASE) + for dir in dirs + OUTPUT = joinpath(OUTPUT_BASE, dir) + INPUT = joinpath(INPUT_BASE, dir) + for file in filter(f -> endswith(f, ".jl"), readdir(INPUT)) + Literate.markdown(joinpath(INPUT, file), OUTPUT) + end + end +end + +makedocs( + format = Documenter.HTML(; + prettyurls=get(ENV, "CI", "false") == "true", + canonical="https://github.com/JuliaImageRecon/AbstractImageReconstruction.jl", + assets=String[], + collapselevel=1, + ), + repo="https://github.com/JuliaImageRecon/AbstractImageReconstruction.jl/blob/{commit}{path}#{line}", + modules = [AbstractImageReconstruction], + sitename = "AbstractImageReconstruction.jl", + authors = "Niklas Hackelberg, Tobias Knopp", + pages = [ + "Home" => "index.md", + "Example: Radon Reconstruction Package" => Any[ + "Introduction" => "example_intro.md", + "Radon Data" => "generated/example/0_radon_data.md", + "Interface" => "generated/example/1_interface.md", + "Direct Reconstruction" => "generated/example/2_direct.md", + "Direct Reconstruction Result" => "generated/example/3_direct_result.md", + ], + "How to" => Any[ + #"Construct RecoPlan" => "generated/howto/reco_plan.md", + #"Caching" => "generated/howto/caching.md", + #"Listeners" => "generated/howto/listeners.md", + ], + #"API Reference" => Any["Solvers" => "API/solvers.md", + #"Regularization Terms" => "API/regularization.md"], + + ], + pagesonly = true, + checkdocs = :none, + doctest = false, + doctestfilters = [r"(\d*)\.(\d{4})\d+"] + ) + +deploydocs(repo = "github.com/JuliaImageRecon/AbstractImageReconstruction.jl.git", push_preview = true) \ No newline at end of file diff --git a/docs/src/API/algorithm.md b/docs/src/API/algorithm.md new file mode 100644 index 0000000..3887043 --- /dev/null +++ b/docs/src/API/algorithm.md @@ -0,0 +1,52 @@ +# API for Regularizers +This page contains documentation of the public API of the RegularizedLeastSquares. In the Julia +REPL one can access this documentation by entering the help mode with `?` + +```@docs +RegularizedLeastSquares.L1Regularization +RegularizedLeastSquares.L2Regularization +RegularizedLeastSquares.L21Regularization +RegularizedLeastSquares.LLRRegularization +RegularizedLeastSquares.NuclearRegularization +RegularizedLeastSquares.TVRegularization +``` + +## Projection Regularization +```@docs +RegularizedLeastSquares.PositiveRegularization +RegularizedLeastSquares.RealRegularization +``` + +## Nested Regularization +```@docs +RegularizedLeastSquares.innerreg(::AbstractNestedRegularization) +RegularizedLeastSquares.sink(::AbstractNestedRegularization) +RegularizedLeastSquares.sinktype(::AbstractNestedRegularization) +``` + +## Scaled Regularization +```@docs +RegularizedLeastSquares.AbstractScaledRegularization +RegularizedLeastSquares.scalefactor +RegularizedLeastSquares.NormalizedRegularization +RegularizedLeastSquares.NoNormalization +RegularizedLeastSquares.MeasurementBasedNormalization +RegularizedLeastSquares.SystemMatrixBasedNormalization +RegularizedLeastSquares.FixedParameterRegularization +``` + +## Misc. Nested Regularization +```@docs +RegularizedLeastSquares.MaskedRegularization +RegularizedLeastSquares.TransformedRegularization +RegularizedLeastSquares.PlugAndPlayRegularization +``` + +## Miscellaneous Functions +```@docs +RegularizedLeastSquares.prox!(::AbstractParameterizedRegularization, ::AbstractArray) +RegularizedLeastSquares.prox!(::Type{<:AbstractParameterizedRegularization}, ::Any, ::Any) +RegularizedLeastSquares.norm(::AbstractParameterizedRegularization, ::AbstractArray) +RegularizedLeastSquares.λ(::AbstractParameterizedRegularization) +RegularizedLeastSquares.norm(::Type{<:AbstractParameterizedRegularization}, ::Any, ::Any) +``` \ No newline at end of file diff --git a/docs/src/API/plan.md b/docs/src/API/plan.md new file mode 100644 index 0000000..5774048 --- /dev/null +++ b/docs/src/API/plan.md @@ -0,0 +1,60 @@ +# API for Solvers +This page contains documentation of the public API of the RegularizedLeastSquares. In the Julia +REPL one can access this documentation by entering the help mode with `?` + +## solve! +```@docs +RegularizedLeastSquares.solve!(::AbstractLinearSolver, ::Any) +RegularizedLeastSquares.init!(::AbstractLinearSolver, ::Any) +RegularizedLeastSquares.init!(::AbstractLinearSolver, ::AbstractSolverState, ::AbstractMatrix) + +``` + +## ADMM +```@docs +RegularizedLeastSquares.ADMM +``` + +## CGNR +```@docs +RegularizedLeastSquares.CGNR +``` + +## Kaczmarz +```@docs +RegularizedLeastSquares.Kaczmarz +``` + +## FISTA +```@docs +RegularizedLeastSquares.FISTA +``` + +## OptISTA +```@docs +RegularizedLeastSquares.OptISTA +``` + +## POGM +```@docs +RegularizedLeastSquares.POGM +``` + +## SplitBregman +```@docs +RegularizedLeastSquares.SplitBregman +``` + +## Miscellaneous +```@docs +RegularizedLeastSquares.solverstate +RegularizedLeastSquares.solversolution +RegularizedLeastSquares.solverconvergence +RegularizedLeastSquares.StoreSolutionCallback +RegularizedLeastSquares.StoreConvergenceCallback +RegularizedLeastSquares.CompareSolutionCallback +RegularizedLeastSquares.linearSolverList +RegularizedLeastSquares.createLinearSolver +RegularizedLeastSquares.applicableSolverList +RegularizedLeastSquares.isapplicable +``` \ No newline at end of file diff --git a/docs/src/example_intro.md b/docs/src/example_intro.md new file mode 100644 index 0000000..353fcc5 --- /dev/null +++ b/docs/src/example_intro.md @@ -0,0 +1,27 @@ +# Small Reconstruction Package for Radon projections +In this example we will implement a small image reconstruction package with the help of `AbstractImageReconstruction.jl`. Our example reconstruction package aims to provide direct and iterative reconstruction algorithms for Radon projection data with optional GPU support. + +Most of the desired functionality is already implemented in various Julia packages. Our reconstruction packages now needs to properly connect these packages and transform the data into the appropriate formats for each package. + +## Installation +In addition to AbstractImageReconstruction.jl, we will need a few more packages to get started. We can install these packages using the Julia package manager. Open a Julia REPL and run the following command: + +```julia +using Pkg +Pkg.add("AbstractImageReconstruction") +``` +This will download and install AbstractImageReconstruction.jl and its dependencies. To install a different version, please consult the [Pkg documentation](https://pkgdocs.julialang.org/dev/managing-packages/#Adding-packages). + + +[RadonKA.jl](https://github.com/roflmaostc/RadonKA.jl/tree/main) provides us with fast Radon forward and backprojections, which we can use for direct reconstructions and preparing example data for our package. + +[LinearOperatorCollection.jl](https://github.com/JuliaImageRecon/LinearOperatorCollection.jl) wraps the functionality of RadonKA.jl in a matrix-free linear operator, which can be used in iterative solvers. + +[RegularizedLeastSquares.jl](https://github.com/JuliaImageRecon/RegularizedLeastSquares.jl) offers a variety of iterative solver and regularization options. + +[ImagePhantoms.jl](https://github.com/JuliaImageRecon/ImagePhantoms.jl) and [ImageGeoms.jl](https://github.com/JuliaImageRecon/ImageGeoms.jl) allow us to define digital software "phantoms", which we will use to test our reconstruction algorithms. + +Lastly, we will use [CairoMakie.jl](https://docs.makie.org/stable/) to visualize our results. + +## Outline + diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..6b3ddaa --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,47 @@ +# AbstractImageReconstruction.jl + +*Abstract Interface for Medical Image Reconstruction Packages* + +## Introduction + +AbstractImageReconstruction.jl is a Julia package that serves as the core API for medical imaging packages. It provides implementations an interface and type hierarchy with which one can represent and implement image reconstruction algorithms, their parameters and runtime behaviour. In particular, this package serves as the API of the Julia packages [MPIReco.jl](https://github.com/MagneticParticleImaging/MPIReco.jl). + +## Features + +* Variety of optimization algorithms optimized for least squares problems +* Storing, loading and manipulating of reconstruction algorithms with (partially) set parameters +* Attaching callback listeners to parameters + +## Installation + +Within Julia, use the package manager: +```julia +using Pkg +Pkg.add("AbstractImageReconstruction") +``` +AbstractImageReconstruction is not intended to be used alone, but together with an image reconstruction package that implements the provided interface, such as [MPIReco.jl](https://github.com/MagneticParticleImaging/MPIReco.jl). + +## Usage +Concrete construction of reconstruction algorithms depend on the implementation of the reconstruction package. Once an algorithms is constructed with the given paramters, images can be reconstructed as follows: +```julia +using AbstractImageReconstruction, MPIReco + +params = ... # Setup reconstruction paramter +algo = ... # Setup chosen algorithm with params +raw = ... # Setup raw data + +image = reconstruct(algo, raw) +``` +Once an algorithm is constructed it can be transformed into a `RecoPlan`. These are mutable and transparent wrappers around the nested types of the algorithm and its paramters, that can be stored and restored to and from TOML files. + +```julia +plan = toPlan(algo) +savePlan(MPIReco, "Example", plan) +plan = loadPlan(MPIReco, "Example", [MPIReco, RegularizedLeastSquares, MPIFiles]) + +algo2 = build(plan) +algo == algo2 # true +``` +Unlike concrete algorithm instances, a `RecoPlan` may still be missing certain values of its fields and it can encode the structure of an image reconstruction algorithm without concrete parameterization. + +It is also possible to attach `Listeners` to `RecoPlan` fields, that call user-specified functions if they are changed. This allows specific `RecoPlans` to provide smart default paramter choices or embedding a plan into a GUI. \ No newline at end of file diff --git a/docs/src/literate/example/0_radon_data.jl b/docs/src/literate/example/0_radon_data.jl new file mode 100644 index 0000000..5765a1c --- /dev/null +++ b/docs/src/literate/example/0_radon_data.jl @@ -0,0 +1,98 @@ +# # Radon Data + +# In this example we will setup our Radon data using RadonKA.jl, ImagePhantoms.jl and ImageGeoms.jl. We will start with simple 2D images and their sinograms and move up to a time series of 3D images and sinograms. + +# ## Background +# The Radon transform is a integral transform which 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). + +# ## 2D Phantom +# We will use a simple Shepp-Logan phantom to generate our Radon data. The phantom is a standard test image for medical imaging and consists of a number of ellipses with different intensities. +# It can be generated with ImagePhantoms.jl and ImageGeoms.jl. as follows: + +using ImagePhantoms, ImageGeoms +N = 256 +image = shepp_logan(N, SheppLoganToft()) +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 image or phantom under transformation, the anlges 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 π. +using RadonKA +angles = collect(range(0, π, 256)) +sinogram = Array(RadonKA.radon(image, angles)) +size(sinogram) + + +# To visualize our progress so far, we will use CairoMakie.jl and a small helper function: +using CairoMakie +function plot_image(figPos, img; title = "", width = 150, height = 150) + ax = CairoMakie.Axis(figPos[1, 1]; yreversed=true, title, width, height) + hidedecorations!(ax) + hm = heatmap!(ax, img) + Colorbar(figPos[1, 2], hm) +end +fig = Figure() +plot_image(fig[1,1], image, title = "Image") +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. +# 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 = (128, 128, 128) +params = map(collect, ellipsoid_parameters(; fovs = shape)); + +# We then scale the intensities of the ellipsoids to [0.0, ..., 1.0]: +toft_settings = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] +for idx in eachindex(toft_settings) + params[idx][10] = toft_settings[idx] +end + +# Finally, we create the 3D Shepp-Logan phantom by defining and sampling our image geometry: +ob = ellipsoid(map(Tuple, params)) +ig = ImageGeom(;dims = shape) +image = phantom(axes(ig)..., ob) +size(image) + +# Now we can compute the 3D Radon transform of our phantom: +sinogram = Array(RadonKA.radon(image, angles)) +size(sinogram) + +# Let's visualize the 3D Radon data: +fig = Figure() +plot_image(fig[1,1], reverse(image[52,:,:]), title = "Slice YZ at z=48") +plot_image(fig[1,2], image[:,80,:], title = "Slice XZ at y=80") +plot_image(fig[2,1], reverse(image[:, :, 48]), title = "Slice XY at z=48") +plot_image(fig[2,2], reverse(sinogram[:,:,48]), title = "Sinogram at z=48") +plot_image(fig[3,1], reverse(image[:, :, 32]), title = "Slice XY at z=32") +plot_image(fig[3,2], reverse(sinogram[:,:,32]), title = "Sinogram at z=32") +resize_to_layout!(fig) +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. +images = similar(image, size(image)..., 5) +sinograms = similar(sinogram, size(sinogram)..., 5) +for (i, intensity) in enumerate(range(params[3][end], 0.3, 5)) + params[3][end] = intensity + local ob = ellipsoid(map(Tuple, params)) + local ig = ImageGeom(;dims = shape) + images[:, :, :, i] = phantom(axes(ig)..., ob) + sinograms[:, :, :, i] = Array(RadonKA.radon(images[:, :, :, i], angles)) +end +size(sinograms) + +fig = Figure() +for i = 1:5 + plot_image(fig[i,1], reverse(images[:, :, 48, i])) + plot_image(fig[i,2], sinograms[:, :, 48, i]) +end +resize_to_layout!(fig) +fig + +# The goal of our reconstruction package is now to recover an approximation of the time-series 3D phantoms from a given time-series of sinograms. +# In the next section we will introduce our class hierarchies and explore the API of AbstractImageReconstruction.jl. \ No newline at end of file diff --git a/docs/src/literate/example/1_interface.jl b/docs/src/literate/example/1_interface.jl new file mode 100644 index 0000000..dafede0 --- /dev/null +++ b/docs/src/literate/example/1_interface.jl @@ -0,0 +1,60 @@ +# # Interface +# This section introduces the abstract types we need to implement for our reconstruction package and how they relate to interface and types of AbstractImageReconstruction.jl. + +# ## Abstract Types +# We start by defining the abstract types we need to implement for our reconstruction package. AbstractImageReconstruction.jl provides two abstract types: +# ```julia +# abstract type AbstractImageReconstructionAlgorithm end +# abstract type AbstractImageReconstructionParameters end +# ``` +# `AbstractImageReconstructionAlgorithms` represent a given reconstruction algorithm, while `AbstractImageReconstructionParameters` represent the parameters an algorithm was constructed with. +# Once constructed, algorithms can be used to reconstruct images repeatly and idealy without unecessary recomputations. + +# For our package we extend these abstract types with our own abstract subtypes: +using AbstractImageReconstruction +abstract type AbstractRadonAlgorithm <: AbstractImageReconstructionAlgorithm end +abstract type AbstractRadonParameters <: AbstractImageReconstructionParameters end + +# Later on we will have parameters that are shared between different algorithms and parameters for different processings steps of a reconstruction. +# In our case we will have preprocessing parameters and reconstruction parameters. For these we introduce the following abstract types: +abstract type AbstractRadonPreprocessingParameters <: AbstractRadonParameters end +abstract type AbstractRadonReconstructionParameters <: AbstractRadonParameters end + +# Since we want to implement both direct and iterative methods for our reconstruction, we introduce the following abstract types: +abstract type AbstractDirectRadonAlgorithm <: AbstractRadonAlgorithm end +abstract type AbstractIterativeRadonAlgorithm <: AbstractRadonAlgorithm end + +# ## Internal Interface +# Reconstruction algorithms in AbstractImageReconstruction.jl are expected to be implemented in terms 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 invoke the function `process` with the type of the algorithm: +# ```julia +# process(algo::AbstractImageReconstructionAlgorithm, param::AbstractImageReconstructionParameters, inputs...) = process(typeof(algo), param, inputs...) +# ``` +# The implementation of reconstruction algorithms is thus expected to either implement the `process` function for the algorithm type or instance. Dispatch on instances allow an instance to change its state, while dispatch on types allows for pure helper functions. +# 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. + +# Let's define a preprocessing step we can share between our algorithms. We want to allow the user to select certain frames from a time data and average them. +# We will use the `@kwdef` macro to provide constructor with keyword arguments and default values +using Statistics +Base.@kwdef struct RadonPreprocessingParameters <: AbstractRadonPreprocessingParameters + frames::Vector{Int64} = [] + numAverages::Int64 = 1 +end +function AbstractImageReconstruction.process(::Type{<:AbstractRadonAlgorithm}, params::RadonPreprocessingParameters, data::AbstractArray{T, 4}) where {T} + frames = isempty(params.frames) ? (1:size(data, 4)) : params.frames + data = data[:, :, :, frames] + + if params.numAverages > 1 + data = reshape(data, size(data)[1:3]..., params.numAverage, :) + data = dropdims(mean(data, dims = 4), dims = 4) + end + + return data +end + +# ## External 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. + +# Algorithms can internally construct and use further algorithms to which they pass data. \ No newline at end of file diff --git a/docs/src/literate/example/2_direct.jl b/docs/src/literate/example/2_direct.jl new file mode 100644 index 0000000..f3d25f4 --- /dev/null +++ b/docs/src/literate/example/2_direct.jl @@ -0,0 +1,63 @@ +include("../../literate/example/1_interface.jl") #hide +using RadonKA #hide + +# # Direct Reconstruction +# 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: +abstract type AbstractDirectRadonReconstructionParameters <: AbstractRadonReconstructionParameters end +# The backprojection parameters are simple and only contain the number of angles: +Base.@kwdef 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 + 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} + result = [] + for i = 1:size(data, 4) + push!(result, process(algoT, params, 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) + +# ## 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, R} <: AbstractRadonParameters where {P<:AbstractRadonPreprocessingParameters, R<:AbstractDirectRadonReconstructionParameters} + 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 + +# Now we can define the algorithm type itself. Algorithms are usually constructed with one argument passing in the user parameters: +mutable struct DirectRadonAlgorithm{D} <: AbstractDirectRadonAlgorithm where D <: DirectRadonParameters + 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.parameters(algo::DirectRadonAlgorithm) = algo.parameter + +# And implement the `put!` and `take!` functions, mimicking the behavior of a FIFO channel: +Base.take!(algo::DirectRadonAlgorithm) = Base.take!(algo.output) +function Base.put!(algo::DirectRadonAlgorithm, data::AbstractArray{T, 4}) where {T} + lock(algo.output) do + put!(algo.output, process(algo, algo.parameter, data)) + end +end + +# 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. diff --git a/docs/src/literate/example/3_direct_result.jl b/docs/src/literate/example/3_direct_result.jl new file mode 100644 index 0000000..5cd3723 --- /dev/null +++ b/docs/src/literate/example/3_direct_result.jl @@ -0,0 +1,62 @@ +include("../../literate/example/1_interface.jl") #hide +include("../../literate/example/2_direct.jl") #hide +using RadonKA, ImagePhantoms, ImageGeoms, CairoMakie, AbstractImageReconstruction #hide +using CairoMakie #hide +function plot_image(figPos, img; title = "", width = 150, height = 150) #hide + ax = CairoMakie.Axis(figPos[1, 1]; yreversed=true, title, width, height) #hide + hidedecorations!(ax) #hide + hm = heatmap!(ax, img) #hide + Colorbar(figPos[1, 2], hm) #hide +end #hide +angles = collect(range(0, π, 256)) #hide +shape = (128, 128, 128) #hide +params = map(collect, ellipsoid_parameters(; fovs = shape)) #hide +toft_settings = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] #hide +for idx in eachindex(toft_settings) #hide + params[idx][10] = toft_settings[idx] #hide +end #hide +ob = ellipsoid(map(Tuple, params)) #hide +ig = ImageGeom(;dims = shape) #hide +image = phantom(axes(ig)..., ob) #hide +sinogram = Array(RadonKA.radon(image, angles)) #hide +sinograms = similar(sinogram, size(sinogram)..., 5) #hide +images = similar(image, size(image)..., 5) #hide +for (i, intensity) in enumerate(range(params[3][end], 0.3, 5)) #hide + params[3][end] = intensity #hide + local ob = ellipsoid(map(Tuple, params)) #hide + local ig = ImageGeom(;dims = shape) #hide + images[:, :, :, i] = phantom(axes(ig)..., ob) #hide + sinograms[:, :, :, i] = Array(RadonKA.radon(images[:, :, :, i], angles)) #hide +end #hide + + +# # Direct Reconstruction Result +# Now that we have implemented our direct reconstruction algorithm, we can use it to reconstruct for example the first three images of our time series. +# We first prepare our parameters: +pre = RadonPreprocessingParameters(frames = collect(1:3)) +back_reco = RadonBackprojectionParameters(;angles) +filter_reco = RadonFilteredBackprojectionParameters(;angles); + +# Then we can construct the algorithms: +algo_back = DirectRadonAlgorithm(DirectRadonParameters(pre, back_reco)) +algo_filter = DirectRadonAlgorithm(DirectRadonParameters(pre, filter_reco)); + +# And apply them to our sinograms: +imag_back = reconstruct(algo_back, sinograms) +imag_filter = reconstruct(algo_filter, sinograms); + +# Finally we can visualize the results: +fig = Figure() +for i = 1:3 + plot_image(fig[i,1], reverse(images[:, :, 48, i])) + plot_image(fig[i,2], sinograms[:, :, 48, i]) + plot_image(fig[i,3], reverse(imag_back[:, :, 48, i])) + plot_image(fig[i,4], reverse(imag_filter[:, :, 48, i])) +end +resize_to_layout!(fig) +fig + +# To add new functionality to our direct reconstruction algorithm, we can simply add new preprocessing or reconstruction parameters and implement the according processing steps. This way we can easily extend our algorithm to support new features. +# We can also define our own file format to store the results of our reconstruction algorithms or dispatch our preprocessing on files and pass files instead of the sinograms directly. + +# The way we use the algorithms here, requires the user to reconstruct the algorithm for each changed parameter. After we implement the more complex iterative reconstructions, we will take a look at how to make this process more manageable. \ No newline at end of file