Added motion parameters utilities

README.md

@@ -1,3 +1,3 @@
 # UtilitiesForMRI
 
-Utilities for MRI
+Collection of utilities for 3D MRI, first and foremost comprising a NFFT linear operator which may be perturbed by phase-encoding dependent rigid-body motion. The perturbed operator is differentiable with respect to rigid-body motion parameters (translation and rotations).

src/UtilitiesForMRI.jl

@@ -8,5 +8,6 @@ include("./kspace_sampling.jl")
 include("./translations.jl")
 include("./rotations.jl")
 include("./nfft.jl")
+include("./motion_parameter_utils.jl")
 
 end

src/motion_parameter_utils.jl

@@ -0,0 +1,65 @@
+# Motion parameter utilities
+
+export derivative1d_linop, derivative1d_motionpars_linop, interpolation1d_linop, interpolation1d_motionpars_linop
+
+function derivative1d_linop(t::AbstractVector{T}, order::Integer) where {T<:Real}
+    (order != 1) && (order != 2) && error("Order not supported (only 1 or 2)")
+    nt = length(t)
+    Δt = t[2:end]-t[1:end-1]
+    if order == 1
+        d_0 = [-T(1)./Δt; T(0)]
+        d_p1 = T(1)./Δt;
+        return spdiagm(nt, nt, 0 => d_0, 1 => d_p1)
+    elseif order == 2
+        Δt_m1 = Δt[1:end-1]
+        Δt_p1 = Δt[2:end]
+        Δt_mean = (Δt_m1+Δt_p1)/T(2)
+        d_m1 = T(1)./(Δt_m1.*Δt_mean)
+        d_p1 = [T(0); T(1)./(Δt_p1.*Δt_mean)]
+        d_0  = [T(0); -d_m1-d_p1[2:end]]
+        return spdiagm(nt, nt, -1 => d_m1, 0 => d_0, 1 => d_p1)
+    end
+end
+
+function derivative1d_motionpars_linop(t::AbstractVector{T}, order::Integer; pars::NTuple{6,Bool}=(true, true, true, true, true, true)) where {T<:Real}
+    D = derivative1d_linop(t, order)
+    Id = sparse(I, length(t), length(t))
+    A = []
+    for i = 1:6
+        pars[i] ? push!(A, D) : push!(A, Id)
+    end
+    return cat(A...; dims=(1,2))
+end
+
+function interpolation1d_linop(t::AbstractVector{T}, ti::AbstractVector{T}; interp::Symbol=:linear) where {T<:Real}
+    (interp != :linear) && error("Only linear interpolation supported")
+    nt = length(t)
+    nti = length(ti)
+    I = repeat(reshape(1:nti, :, 1); outer=(1,2))
+    J = Array{Int64,2}(undef, nti, 2)
+    V = Array{T,2}(undef, nti, 2)
+    for i = 1:nti
+        if t[1] == ti[i]
+            J[i,:] = [1 2]
+            V[i,:] = [T(1) T(0)]
+        else
+            idx = maximum(findall(t.<ti[i]))
+            J[i,:] = [idx idx+1]
+            Δt = t[idx+1]-t[idx]
+            V[i,:] = [(t[idx+1]-ti[i])/Δt (ti[i]-t[idx])/Δt]
+        end
+    end
+    return sparse(vec(I),vec(J),vec(V),nti,nt)
+end
+
+function interpolation1d_motionpars_linop(t::NTuple{6,AbstractVector{T}}, ti::NTuple{6,AbstractVector{T}}; interp::Symbol=:linear) where {T<:Real}
+    Ip = []
+    for i=1:6
+        push!(Ip, interpolation1d_linop(t[i], ti[i]; interp=interp))
+    end
+    return cat(Ip...; dims=(1,2))
+end
+
+interpolation1d_motionpars_linop(t::AbstractVector{T}, ti::NTuple{6,AbstractVector{T}}; interp::Symbol=:linear) where {T<:Real} = interpolation1d_motionpars_linop((t,t,t,t,t,t), ti; interp=interp)
+
+interpolation1d_motionpars_linop(t::NTuple{6,AbstractVector{T}}, ti::AbstractVector{T}; interp::Symbol=:linear) where {T<:Real} = interpolation1d_motionpars_linop(t, (ti,ti,ti,ti,ti,ti); interp=interp)

src/nfft.jl

@@ -1,7 +1,7 @@
 # NFFT utilities
 
 export NFFTLinOp, NFFTParametericLinOp, NFFTParametericDelayedEval, JacobianNFFTEvaluated
-export nfft_linop, nfft, Jacobia, ∂
+export nfft_linop, nfft, Jacobia, ∂, sparse_matrix_GaussNewton
 
 
 ## NFFT linear operator
@@ -105,4 +105,22 @@ AbstractLinearOperators.domain_size(∂Fu::JacobianNFFTEvaluated) = (size(∂Fu.
 AbstractLinearOperators.range_size(∂Fu::JacobianNFFTEvaluated) = size(∂Fu.J)[1:2]
 AbstractLinearOperators.matvecprod(∂Fu::JacobianNFFTEvaluated{T}, Δθ::AbstractArray{Complex{T},2}) where {T<:Real} = sum(∂Fu.J.*reshape(Δθ, :, 1, 6); dims=3)[:,:,1]
 Base.:*(∂Fu::JacobianNFFTEvaluated{T}, Δθ::AbstractArray{T,2}) where {T<:Real} = ∂Fu*complex(Δθ)
-AbstractLinearOperators.matvecprod_adj(∂Fu::JacobianNFFTEvaluated{T}, Δd::AbstractArray{Complex{T},2}) where {T<:Real} = sum(conj(∂Fu.J).*reshape(Δd, size(Δd)..., 1); dims=2)[:,1,:]
+AbstractLinearOperators.matvecprod_adj(∂Fu::JacobianNFFTEvaluated{T}, Δd::AbstractArray{Complex{T},2}) where {T<:Real} = sum(conj(∂Fu.J).*reshape(Δd, size(Δd)..., 1); dims=2)[:,1,:]
+
+
+## Other utilities
+
+function sparse_matrix_GaussNewton(J::JacobianNFFTEvaluated{T}) where {T<:Real}
+    J = J.J
+    GN = similar(J, T, size(J, 1), 6, 6)
+    @inbounds for i = 1:6, j = 1:6
+        GN[:,i,j] = real(sum(conj(J[:,:,i]).*J[:,:,j]; dims=2)[:,1])
+    end
+    h(i,j) = spdiagm(0 => GN[:,i,j])
+    return [h(1,1) h(1,2) h(1,3) h(1,4) h(1,5) h(1,6);
+            h(2,1) h(2,2) h(2,3) h(2,4) h(2,5) h(2,6);
+            h(3,1) h(3,2) h(3,3) h(3,4) h(3,5) h(3,6);
+            h(4,1) h(4,2) h(4,3) h(4,4) h(4,5) h(4,6);
+            h(5,1) h(5,2) h(5,3) h(5,4) h(5,5) h(5,6);
+            h(6,1) h(6,2) h(6,3) h(6,4) h(6,5) h(6,6)]
+end

test/runtests.jl

@@ -5,4 +5,5 @@ using UtilitiesForMRI, Test
     include("./test_translations.jl")
     include("./test_rotations.jl")
     include("./test_nfft.jl")
+    include("./test_motionpars_utils.jl")
 end

test/test_motionpars_utils.jl

@@ -0,0 +1,14 @@
+using UtilitiesForMRI, LinearAlgebra, Test, PyPlot
+
+nt = 256
+t = Float64.(1:nt)
+D = derivative1d_motionpars_linop(t, 2; pars=(true,false,true,false,true,true))
+θ = zeros(Float64, nt, 6); θ[129-30:129+30,:] .= 1; θ[1:10,1] .= 1; θ[end-9:end,1] .= 1
+Dθ = reshape(D*vec(θ), :, 6)
+
+nt = 256+1
+t = Float64.(0:nt-1)
+ti = Float64.(0:2:nt-1); nti = length(ti)
+Itp = interpolation1d_motionpars_linop((ti,ti,t,ti,ti,t), t)
+θ = randn(Float64, 4*nti+2*nt)
+Iθ = reshape(Itp*θ, nt, 6)

test/test_nfft.jl

@@ -37,4 +37,9 @@ d, _, ∂Fθu = ∂(F()*u, θ)
 t = 1e-6
 Fθu_p1 = F(θ+0.5*t*Δθ)*u
 Fθu_m1 = F(θ-0.5*t*Δθ)*u
-@test (Fθu_p1-Fθu_m1)/t ≈ ∂Fθu*Δθ rtol=1e-3
+@test (Fθu_p1-Fθu_m1)/t ≈ ∂Fθu*Δθ rtol=1e-3
+
+# Sparse matrix Gauss-Newton
+H = sparse_matrix_GaussNewton(∂Fθu) # ∂Fθu'*∂Fθu
+Δθ = randn(Float64, nt, 6); Δθ *= norm(θ)/norm(Δθ)
+@test H*vec(Δθ) ≈ vec(real(∂Fθu'*(∂Fθu*Δθ))) rtol=1e-6