diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index f09fb72..3051326 100644 --- a/.github/workflows/documentation.yml +++ b/.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 diff --git a/docs/make.jl b/docs/make.jl index 24bc915..007bf9e 100644 --- a/docs/make.jl +++ b/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}} diff --git a/docs/src/examples.md b/docs/src/examples.md index 378e303..2687325 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -1,3 +1,110 @@ # [Getting started](@id examples) -aaa \ No newline at end of file +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]) +``` \ No newline at end of file diff --git a/docs/src/functions.md b/docs/src/functions.md index 5c41278..3af75e4 100644 --- a/docs/src/functions.md +++ b/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 diff --git a/docs/src/installation.md b/docs/src/installation.md index 76c72a6..12996a8 100644 --- a/docs/src/installation.md +++ b/docs/src/installation.md @@ -1,4 +1,4 @@ -# Installation instructions +# [Installation instructions](@id install) In the Julia REPL, simply type `]` and ```julia diff --git a/examples/example_basic_usage.jl b/examples/example_basic_usage.jl index 337a43e..82292b6 100644 --- a/examples/example_basic_usage.jl +++ b/examples/example_basic_usage.jl @@ -1,4 +1,4 @@ -using RetrospectiveMotionCorrectionMRI, FastSolversForWeightedTV, UtilitiesForMRI, AbstractProximableFunctions, LinearAlgebra, Test +using RetrospectiveMotionCorrectionMRI, FastSolversForWeightedTV, UtilitiesForMRI, AbstractProximableFunctions, LinearAlgebra, PyPlot # Spatial geometry fov = (1f0, 2f0, 2f0) @@ -6,56 +6,73 @@ 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)/α)) \ No newline at end of file +# 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]) \ No newline at end of file diff --git a/src/image_reconstruction.jl b/src/image_reconstruction.jl index a538761..80d91c4 100644 --- a/src/image_reconstruction.jl +++ b/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,6 +43,8 @@ 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) @@ -50,10 +52,10 @@ Performs image reconstruction by fitting the data `d` through a linear operator 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 \ No newline at end of file diff --git a/src/motion_corrected_reconstruction.jl b/src/motion_corrected_reconstruction.jl index 18790ca..dd99468 100644 --- a/src/motion_corrected_reconstruction.jl +++ b/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