Updated docs

.github/workflows/documentation.yml

@@ -4,7 +4,7 @@ on:
   push:
     branches:
       - main # update to match your development branch (master, main, dev, trunk, ...)
-    tags: '*'
+    tags: '1.0.0'
   pull_request:
 
 jobs:
@@ -22,6 +22,7 @@ jobs:
                                        Pkg.add(url="https://github.com/grizzuti/AbstractLinearOperators.git");
                                        Pkg.add(url="https://github.com/grizzuti/AbstractProximableFunctions.git");
                                        Pkg.add(url="https://github.com/grizzuti/FastSolversForWeightedTV.git");
+                                       Pkg.add(url="https://github.com/grizzuti/UtilitiesForMRI.git");
                                        Pkg.add(PackageSpec(path=pwd()));
                                        Pkg.instantiate()'
       - name: Build and deploy

docs/make.jl

@@ -1,4 +1,4 @@
-using Documenter, AbstractLinearOperators, AbstractProximableFunctions, FastSolversForWeightedTV, RetrospectiveMotionCorrectionMRI
+using Documenter, AbstractLinearOperators, AbstractProximableFunctions, FastSolversForWeightedTV, UtilitiesForMRI, RetrospectiveMotionCorrectionMRI
 
 const RealOrComplex{T<:Real} = Union{T, Complex{T}}
 

docs/src/examples.md

@@ -1,3 +1,110 @@
 # [Getting started](@id examples)
 
-aaa
+In this section, we present a brief step-by-step tutorial on how to perform motion correction for MRI using the tools provided by `RetrospectiveMotionCorrectionMRI`.
+
+Make sure to have installed the package `RetrospectiveMotionCorrectionMRI` and its dependency by following the instructions highlighted [here](@ref install). For plotting, we are also going to make use of `PyPlot`, to install it press `]` in the Julia REPL and type:
+```julia
+(@v1.8) pkg> add PyPlot
+```
+
+Let's start by loading all the needed modules:
+```julia
+using RetrospectiveMotionCorrectionMRI, FastSolversForWeightedTV, UtilitiesForMRI, AbstractProximableFunctions, LinearAlgebra, PyPlot
+```
+
+We can define a simple spatial discretization `X` and ``k``-space trajectory `K` (in this case, corresponding to Cartesian dense sampling). The associated Fourier transform will be `F`:
+```julia
+# Spatial geometry
+fov = (1f0, 2f0, 2f0)
+n = (64, 64, 64)
+o = (0.5f0, 1f0, 1f0)
+X = spatial_geometry(fov, n; origin=o)
+
+# k-space trajectory (Cartesian dense sampling)
+phase_encoding = (1, 2)
+K = kspace_sampling(X, phase_encoding)
+nt, nk = size(K)
+
+# Fourier operator
+F = nfft_linop(X, K)
+```
+
+We are going now to define a simple 3D image and some rigid motion parameter evolution, which corresponds to a step function in time (e.g. the object moves once halfway through the acquisition). The motion-corrupted data is obtained by perturbing the conventional Fourier transform:
+```julia
+# Setting image ground-truth (= not motion corrupted)
+ground_truth = zeros(ComplexF32, n)
+ground_truth[33-10:33+10, 33-10:33+10, 33-10:33+10] .= 1
+
+# Setting simple rigid motion parameters
+nt, _ = size(K)
+θ_true = zeros(Float32, nt, 6)
+θ_true[div(nt,2)+51:end,:] .= reshape([0.0f0, 0.0f0, 0.0f0, Float32(pi)/180*10, 0f0, 0f0], 1, 6)
+
+# Motion-corrupted data
+d = F(θ_true)*ground_truth
+```
+
+To perform motion correction we need to define several options related to the alternating minimization scheme described in [this methodological section](@ref theory). Refer to that section for in-depth discussion of each of these options:
+```julia
+# Optimization options:
+    ## Image reconstruction options
+    h = spacing(X); L = 4f0*sum(1 ./h.^2)
+    opt_inner = FISTA_options(L; Nesterov=true, niter=10)
+    # η = 1f-2*structural_mean(ground_truth)   # reference guide
+    # P = structural_weight(ground_truth; η=η) # reference guide
+    P = nothing                                # no reference
+    g = gradient_norm(2, 1, size(ground_truth), h; complex=true, weight=P, options=opt_inner)
+    ε = 0.8f0*g(ground_truth)
+    opt_imrecon = image_reconstruction_options(; prox=indicator(g ≤ ε), Nesterov=true, niter=5, niter_estimate_Lipschitz=3, verbose=true, fun_history=true)
+
+    ## Parameter estimation options
+    ti = Float32.(range(1, nt; length=16))
+    t = Float32.(1:nt)
+    Ip = interpolation1d_motionpars_linop(ti, t)
+    D = derivative1d_motionpars_linop(t, 2; pars=(true, true, true, true, true, true))/4f0
+    opt_parest = parameter_estimation_options(; niter=5, steplength=1f0, λ=0f0, scaling_diagonal=1f-3, scaling_mean=1f-4, scaling_id=0f0, reg_matrix=D, interp_matrix=Ip, verbose=true, fun_history=true)
+
+    ## Overall options
+    options = motion_correction_options(; image_reconstruction_options=opt_imrecon, parameter_estimation_options=opt_parest, niter=40, verbose=true, fun_history=true)
+```
+Notably, we can opt to perform motion correction with or without a reference guide. This specific behavior is determined by the weight linear operator `P` (comment/uncomment the corresponding line to study the effect on the final results).
+
+The conventional reconstruction is simply performed by:
+```julia
+# Conventional reconstruction
+u_conventional = F'*d
+```
+The results will be highly corrupted.
+
+To run motion correction:
+```julia
+# Motion-corrected reconstruction
+θ0 = zeros(Float32, length(ti), 6)
+u0 = zeros(ComplexF32, n)
+u, θ = motion_corrected_reconstruction(F, d, u0, θ0, options)
+θ = reshape(Ip*vec(θ), length(t), 6)
+```
+
+Finally, we can compare the results by:
+```julia
+# Plotting
+    ## Image
+    figure()
+    x, y, _ = coord(X)
+    extent = (x[1], x[end], y[1], y[end])
+    vmin = -0.1; vmax = 1.1
+    subplot(1, 3, 1)
+    imshow(abs.(u_conventional[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Conventional")
+    subplot(1, 3, 2)
+    imshow(abs.(u[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Corrected")
+    subplot(1, 3, 3)
+    imshow(abs.(ground_truth[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Ground-truth")
+
+    ## Motion parameters
+    figure()
+    plot(θ[:, 4])
+    plot(θ_true[:, 4])
+```

docs/src/functions.md

@@ -13,7 +13,7 @@ motion_correction_options(; image_reconstruction_options::ImageReconstructionOpt
 ## Image reconstruction
 
 ```@docs
-image_reconstruction(F::AbstractLinearOperator{CT,3,2}, d::AbstractArray{CT,2}, initial_estimate::AbstractArray{CT,3}, options::RetrospectiveMotionCorrectionMRI.ImageReconstructionOptionsFISTA) where {CT<:RealOrComplex}
+image_reconstruction(F::AbstractLinearOperator{CT,3,2}, d::AbstractArray{CT,2}, initial_estimate::AbstractArray{CT,3}, options::RetrospectiveMotionCorrectionMRI.ImageReconstructionOptionsFISTA) where {T<:Real,CT<:RealOrComplex{T}}
 ```
 
 ```@docs

docs/src/installation.md

@@ -1,4 +1,4 @@
-# Installation instructions
+# [Installation instructions](@id install)
 
 In the Julia REPL, simply type `]` and
 ```julia

examples/example_basic_usage.jl

@@ -1,61 +1,78 @@
-using RetrospectiveMotionCorrectionMRI, FastSolversForWeightedTV, UtilitiesForMRI, AbstractProximableFunctions, LinearAlgebra, Test
+using RetrospectiveMotionCorrectionMRI, FastSolversForWeightedTV, UtilitiesForMRI, AbstractProximableFunctions, LinearAlgebra, PyPlot
 
 # Spatial geometry
 fov = (1f0, 2f0, 2f0)
 n = (64, 64, 64)
 o = (0.5f0, 1f0, 1f0)
 X = spatial_geometry(fov, n; origin=o)
 
-# Cartesian sampling in k-space
-phase_encoding = (1,2)
+# k-space trajectory (Cartesian dense sampling)
+phase_encoding = (1, 2)
 K = kspace_sampling(X, phase_encoding)
 nt, nk = size(K)
 
 # Fourier operator
 F = nfft_linop(X, K)
 
-# Setting image ground-truth (=not motion corrupted)
-ground_truth = zeros(ComplexF32, n); ground_truth[div(64,2)+1-10:div(64,2)+1+10, div(64,2)+1-10:div(64,2)+1+10, div(64,2)+1-10:div(64,2)+1+10] .= 1
+# Setting image ground-truth (= not motion corrupted)
+ground_truth = zeros(ComplexF32, n)
+ground_truth[33-10:33+10, 33-10:33+10, 33-10:33+10] .= 1
 
-# Setting simple motion corruption
+# Setting simple rigid motion parameters
 nt, _ = size(K)
 θ_true = zeros(Float32, nt, 6)
 θ_true[div(nt,2)+51:end,:] .= reshape([0.0f0, 0.0f0, 0.0f0, Float32(pi)/180*10, 0f0, 0f0], 1, 6)
+
+# Motion-corrupted data
 d = F(θ_true)*ground_truth
 
+# Optimization options:
+    ## Image reconstruction options
+    h = spacing(X); L = 4f0*sum(1 ./h.^2)
+    opt_inner = FISTA_options(L; Nesterov=true, niter=10)
+    # η = 1f-2*structural_mean(ground_truth)   # reference guide
+    # P = structural_weight(ground_truth; η=η) # reference guide
+    P = nothing                                # no reference
+    g = gradient_norm(2, 1, size(ground_truth), h; complex=true, weight=P, options=opt_inner)
+    ε = 0.8f0*g(ground_truth)
+    opt_imrecon = image_reconstruction_options(; prox=indicator(g ≤ ε), Nesterov=true, niter=5, niter_estimate_Lipschitz=3, verbose=true, fun_history=true)
+
+    ## Parameter estimation options
+    ti = Float32.(range(1, nt; length=16))
+    t = Float32.(1:nt)
+    Ip = interpolation1d_motionpars_linop(ti, t)
+    D = derivative1d_motionpars_linop(t, 2; pars=(true, true, true, true, true, true))/4f0
+    opt_parest = parameter_estimation_options(; niter=5, steplength=1f0, λ=0f0, scaling_diagonal=1f-3, scaling_mean=1f-4, scaling_id=0f0, reg_matrix=D, interp_matrix=Ip, verbose=true, fun_history=true)
+
+    ## Overall options
+    options = motion_correction_options(; image_reconstruction_options=opt_imrecon, parameter_estimation_options=opt_parest, niter=40, verbose=true, fun_history=true)
+
 # Conventional reconstruction
 u_conventional = F'*d
 
-# Image reconstruction options
-h = spacing(X); LD = 4f0*sum(1 ./h.^2)
-opt_inner = FISTA_options(LD; Nesterov=true, niter=10)
-g = gradient_norm(2, 1, size(ground_truth), h; complex=true, options=opt_inner)
-ε = 0.8f0*g(ground_truth)
-x = reshape(range(-1f0, 1f0; length=64), :, 1, 1)
-y = reshape(range(-1f0, 1f0; length=64), 1, :, 1)
-z = reshape(range(-1f0, 1f0; length=64), 1, 1, :)
-r = sqrt.(x.^2 .+y.^2 .+z.^2)
-C = zero_set(ComplexF32, r .< 0.1f0)
-h = indicator(set_options(g+indicator(C), opt_inner) ≤ ε)
-opt_imrecon = image_reconstruction_options(; prox=h, Nesterov=true, niter=5, verbose=true, fun_history=true)
-
-# Parameter estimation options
-ti = Float32.(range(1, nt; length=16))
-# ti = Float32.(1:nt)
-t = Float32.(1:nt)
-Ip = interpolation1d_motionpars_linop(ti, t)
-D = derivative1d_motionpars_linop(t, 2; pars=(true, true, true, true, true, true))/4f0
-opt_parest = parameter_estimation_options(; niter=5, steplength=1f0, λ=0f0, scaling_diagonal=1f-5, scaling_id=1f-1, reg_matrix=D, interp_matrix=Ip, verbose=true, fun_history=true)
-
-# Solution
-opt = motion_correction_options(; image_reconstruction_options=opt_imrecon, parameter_estimation_options=opt_parest, niter=40, niter_estimate_Lipschitz=3, verbose=true, fun_history=true)
+# Motion-corrected reconstruction
 θ0 = zeros(Float32, length(ti), 6)
 u0 = zeros(ComplexF32, n)
-u, θ = motion_corrected_reconstruction(F, d, u0, θ0, opt)
+u, θ = motion_corrected_reconstruction(F, d, u0, θ0, options)
 θ = reshape(Ip*vec(θ), length(t), 6)
 
-# Comparison
-u_conventional = F'*d
-α = norm(ground_truth, Inf)
-@info string("Conventional: (psnr,ssim)=(", psnr(abs.(u_conventional)/α, abs.(ground_truth)/α), ",", ssim(abs.(u_conventional)/α, abs.(ground_truth)/α),")")
-@info string("Motion-corrected: (psnr,ssim)=(", psnr(abs.(u)/α, abs.(ground_truth)/α), ",", ssim(abs.(u)/α, abs.(ground_truth)/α))
+# Plotting
+    ## Image
+    figure()
+    x, y, _ = coord(X)
+    extent = (x[1], x[end], y[1], y[end])
+    vmin = -0.1; vmax = 1.1
+    subplot(1, 3, 1)
+    imshow(abs.(u_conventional[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Conventional")
+    subplot(1, 3, 2)
+    imshow(abs.(u[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Corrected")
+    subplot(1, 3, 3)
+    imshow(abs.(ground_truth[:, end:-1:1, 33])'; vmin=vmin, vmax=vmax, extent=extent)
+    title("Ground-truth")
+
+    ## Motion parameters
+    figure()
+    plot(θ[:, 4])
+    plot(θ_true[:, 4])

src/image_reconstruction.jl

@@ -20,7 +20,7 @@ end
                                    fun_history=false)
 
 Returns image reconstruction options for the routine [`image_reconstruction`](@ref):
-- `prox`: proximal operator associated to a chosen regularization function
+- `prox`: regularization function (for which a proximal operator is implemented)
 - `Lipschitz_constant`: Lipschitz constant of a smooth objective
 - `niter_estimate_Lipschitz`: when set to an `Integer`, the iterative power method is invoked in order to estimate the Lipschitz constant with the specified number of iteration
 - `Nesterov`: allows Nesterov acceleration (default)
@@ -43,17 +43,19 @@ end
 
 AbstractProximableFunctions.fun_history(options::ImageReconstructionOptionsFISTA) = fun_history(options.options)
 
+AbstractProximableFunctions.set_Lipschitz_constant(options::ImageReconstructionOptionsFISTA, L::Real) = ImageReconstructionOptionsFISTA(options.prox, set_Lipschitz_constant(options.options, L), options.niter_estimate_Lipschitz)
+
 """
     image_reconstruction(F, d, initial_estimate, options)
 
 Performs image reconstruction by fitting the data `d` through a linear operator `F` (e.g. the Fourier transform). `initial_estimate` sets the initial guess. The minimization `options` are passed via the routine [`image_reconstruction_options`](@ref).
 
 See this [section](@ref imrecon) for more details on the solution algorithm.
 """
-function image_reconstruction(F::AbstractLinearOperator{CT,3,2}, d::AbstractArray{CT,2}, initial_estimate::AbstractArray{CT,3}, options::ImageReconstructionOptionsFISTA) where {CT<:RealOrComplex}
+function image_reconstruction(F::AbstractLinearOperator{CT,3,2}, d::AbstractArray{CT,2}, initial_estimate::AbstractArray{CT,3}, options::ImageReconstructionOptionsFISTA) where {T<:Real,CT<:RealOrComplex{T}}
     if ~isnothing(options.niter_estimate_Lipschitz)
         L = T(1.1)*spectral_radius(F*F'; niter=options.niter_estimate_Lipschitz)
-        set_Lipschitz_constant(options.options, L)
+        opt_FISTA = set_Lipschitz_constant(options.options, L)
     end
-    return leastsquares_solve(F, d, options.prox, initial_estimate, options.options)
+    return leastsquares_solve(F, d, options.prox, initial_estimate, opt_FISTA)
 end

src/motion_corrected_reconstruction.jl

@@ -55,7 +55,7 @@ function motion_corrected_reconstruction(F::StructuredNFFTtype2LinOp{T}, d::Abst
         # Image reconstruction
         options.verbose && (@info string("--- Image reconstruction..."))
         flag_interp ? (θ_ = reshape(Ip*vec(θ), :, 6)) : (θ_ = θ); Fθ = F(θ_)
-        u = image_reconstruction(Fθ, d, u, options_imrecon)
+        u = image_reconstruction(Fθ, d, u, options.options_imrecon)
 
         # Motion-parameter estimation
         (n == options.niter) && break