Merge pull request #80 from MagneticResonanceImaging/bugfix_cmaps

Fix incorrect trj indexing for low-res recon in calculate_coil_maps

Author: IvoMaatman

ext/MRISubspaceReconCUDAExt/BackProjection.jl

@@ -87,6 +87,9 @@ end
 function MRISubspaceRecon.calculate_backprojection(data::CuArray{Tc}, trj::CuArray{<:Integer,3}, cmaps::AbstractVector{<:CuArray}; U=cu(I(size(trj)[end])), sample_mask=CUDA.ones(Bool, size(trj)[2:end])) where {Tc <: Complex}
     Ncoeff = size(U, 2)
     img_shape = size(cmaps[1])
+    trj_cpu = Array(trj)
+    mask_cpu = Array(sample_mask)
+    @assert all(d -> all((@view(trj_cpu[d, :, :])[mask_cpu]) .>= 1) && all((@view(trj_cpu[d, :, :])[mask_cpu]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     img_idx = CartesianIndices(img_shape)
     T = real(eltype(data))
 
@@ -119,6 +122,9 @@ function MRISubspaceRecon.calculate_backprojection(data::CuArray{Tc}, trj::CuArr
 end
 
 function MRISubspaceRecon.calculate_backprojection(data::CuArray{Tc}, trj::CuArray{<:Integer,3}, img_shape; U=cu(I(size(trj)[end])), sample_mask=CUDA.ones(Bool, size(trj)[2:end])) where {Tc <: Complex}
+    trj_cpu = Array(trj)
+    mask_cpu = Array(sample_mask)
+    @assert all(d -> all((@view(trj_cpu[d, :, :])[mask_cpu]) .>= 1) && all((@view(trj_cpu[d, :, :])[mask_cpu]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     Ncoeff = size(U, 2)
     Ncoil = size(data, 3)
     img_idx = CartesianIndices(img_shape)

ext/MRISubspaceReconCUDAExt/CoilMaps.jl

@@ -71,7 +71,9 @@ function MRISubspaceRecon.calculate_coil_maps(
     mask_calib .&= Array(sample_mask)
     mask_calib = cu(mask_calib)
 
-    x = MRISubspaceRecon.reconstruct_coilwise(data, trj, calib_size; U, sample_mask=mask_calib, Niter_cg)
+    # Shift trajectory to 1:calib_size range for reconstruction on the calib grid
+    trj_calib = CuArray{Int16}(trj .- lower_bound)
+    x = MRISubspaceRecon.reconstruct_coilwise(data, trj_calib, calib_size; U, sample_mask=mask_calib, Niter_cg)
 
     imdims = ntuple(i -> i, length(img_shape))
     kbp = fftshift(x, imdims)

ext/MRISubspaceReconCUDAExt/FFTNormalOp.jl

@@ -59,6 +59,9 @@ end
 #############################################################################
 
 function calculate_kernel_cartesian(img_shape, trj::CuArray{<:Integer,3}, U; sample_mask=CUDA.ones(Bool, size(trj)[2:end]), verbose=false)
+    trj_cpu = Array(trj)
+    mask_cpu = Array(sample_mask)
+    @assert all(d -> all((@view(trj_cpu[d, :, :])[mask_cpu]) .>= 1) && all((@view(trj_cpu[d, :, :])[mask_cpu]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     Ncoeff = size(U, 2)
     Nrep = size(U, 3) # number of repetitions (defaults to 1)
 

ext/MRISubspaceReconCUDAExt/NFFTNormalOp.jl

@@ -112,7 +112,6 @@ function calculate_kernel_noncartesian(img_shape_os, trj::CuArray{T,3}, U::CuArr
    # Configure threads and blocks for each kernel within the coefficient loop
     threads_multiply, blocks_multiply, threads_store, blocks_store = launch_config_kernel(nsamp_t, kmask_indcs)
 
-
     verbose && println("calculating non-Cartesian kernel...")
     t = @elapsed CUDA.@sync for ic2 ∈ 1:Ncoeff, ic1 ∈ 1:Ncoeff
         if ic2 >= ic1 # eval. only upper triangular matrix

src/BackProjection.jl

@@ -106,6 +106,7 @@ end
 function calculate_backprojection(data::AbstractArray{Tc}, trj::AbstractArray{<:Integer,3}, cmaps::AbstractVector{<:AbstractArray}; U=I(size(trj)[end]), sample_mask=trues(size(trj)[2:end])) where {Tc <: Complex}
     Ncoeff = size(U, 2)
     img_shape = size(cmaps[1])
+    @assert all(d -> all((@view(trj[d, :, :])[sample_mask]) .>= 1) && all((@view(trj[d, :, :])[sample_mask]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     img_idx = CartesianIndices(img_shape)
 
     dataU = similar(data, img_shape..., Ncoeff)
@@ -131,6 +132,7 @@ function calculate_backprojection(data::AbstractArray{Tc}, trj::AbstractArray{<:
 end
 
 function calculate_backprojection(data::AbstractArray{Tc}, trj::AbstractArray{<:Integer,3}, img_shape; U=I(size(trj)[end]), sample_mask=trues(size(trj)[2:end])) where {Tc <: Complex}
+    @assert all(d -> all((@view(trj[d, :, :])[sample_mask]) .>= 1) && all((@view(trj[d, :, :])[sample_mask]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     Ncoeff = size(U, 2)
     Ncoil = size(data, 3)
     img_idx = CartesianIndices(img_shape)

src/CoilMaps.jl

@@ -80,7 +80,9 @@ function calculate_coil_maps(
     mask_calib = dropdims(all(trj .> lower_bound; dims=1) .& all(trj .< upper_bound; dims=1); dims=1)
     mask_calib .&= sample_mask
 
-    x = reconstruct_coilwise(data, trj, calib_size; U, sample_mask=mask_calib, Niter_cg)
+    # Shift trajectory to 1:calib_size range for reconstruction on the calib grid
+    trj_calib = Int16.(trj .- lower_bound)
+    x = reconstruct_coilwise(data, trj_calib, calib_size; U, sample_mask=mask_calib, Niter_cg)
 
     imdims = ntuple(i -> i, length(img_shape))
     kbp = fftshift(x, imdims)

src/FFTNormalOp.jl

@@ -88,6 +88,7 @@ struct _FFTNormalOp{S,E,F,G,H,I,J,K,L,M}
 end
 
 function calculate_kernel_cartesian(img_shape, trj::AbstractArray{<:Integer,3}, U; sample_mask=trues(size(trj)[2:end]), verbose=false)
+    @assert all(d -> all((@view(trj[d, :, :])[sample_mask]) .>= 1) && all((@view(trj[d, :, :])[sample_mask]) .<= img_shape[d]), 1:size(trj, 1)) "Cartesian trajectory values must be in the range 1:img_shape[d] for each dimension d."
     Ncoeff = size(U, 2)
     Λ = zeros(eltype(U), Ncoeff, Ncoeff, img_shape...)