Merge pull request #48 from NordicMRspine/smallfixes

Smallfixes

.github/workflows/CI.yml

@@ -18,8 +18,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7'
           - '1.8'
+          - '1.9'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -40,7 +40,7 @@ LazyArtifacts = "1"
 MAT = "0.10"
 MRICoilSensitivities = "0.1"
 MRIFiles = "0.1, 0.3"
-MRIReco = "0.7, 0.8"
+MRIReco = "0.7"
 NIfTI = "0.5.9, 0.6"
 PolygonOps = "0.1"
 PyPlot = "2"

docs/src/GettingStarted.md

@@ -29,6 +29,7 @@ params[:trust_SCT] = false
 params[:use_centerline] = true
 params[:corr_type] = "FFT_unwrap"
 params[:FFT_interval] = 35 # millimetres
+params[:mask_thresh] = 0.13
 params[:root_path] = "/Users/me/my_data/"
 
 params[:label] = params[:corr_type] * "_rep_" * string(params[:rep])

src/AdjustData.jl

@@ -236,9 +236,10 @@ function ExtractNavigator(rawData::RawAcquisitionData)
     slices = slices[contrastsIndx]
     lines = lines[contrastsIndx]
 
-    nav = zeros(typeof(rawData.profiles[1].data[1,1]), size(rawData.profiles[1].data)[1], size(rawData.profiles[1].data)[2], rawData.params["reconSize"][2], numberslices)
+    nav = zeros(typeof(rawData.profiles[1].data[1,1]), size(rawData.profiles[1].data)[1], size(rawData.profiles[1].data)[2],
+        rawData.params["enc_lim_kspace_encoding_step_1"].maximum+1, numberslices)
 
-    nav_time = zeros(Float64, rawData.params["reconSize"][2], numberslices)
+    nav_time = zeros(Float64, rawData.params["enc_lim_kspace_encoding_step_1"].maximum+1, numberslices)
 
     #Odd indexes are data first echo, Even indexes are navigator data
     for ii = 2:2:length(slices)
@@ -269,9 +270,9 @@ Extract one or more echoes from the acquisition data structure
 
 MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
 """
-function selectEcho!(acqd::AcquisitionData, idx_echo::Vector{Int64})
+function selectEcho!(acqd::AcquisitionData, idx_echo::Union{Vector{Int64}, Nothing})
 
-    if !isempty(idx_echo)
+    if !isnothing(idx_echo)
 
         contrasts = size(acqd.kdata)[1]
         indices = Vector{Int64}(undef, contrasts)
@@ -303,16 +304,20 @@ Extract one or more echoes from the acquisition data structure
 
 MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
 """
-function selectSlice!(acqd::AcquisitionData, idx_slice::Vector{Int64}, nav::Union{Array{Complex{T}, 4}, Nothing} = nothing, nav_time::Union{Array{Float64, 2}, Nothing} = nothing) where {T}
+function selectSlice!(acqd::AcquisitionData, idx_slice::Union{Vector{Int64}, Nothing}, nav::Union{Array{Complex{T}, 4}, Nothing} = nothing, nav_time::Union{Array{Float64, 2}, Nothing} = nothing) where {T}
 
     # get kdata from slice
-    acqd.kdata = acqd.kdata[:,idx_slice,:]
+    if !isnothing(idx_slice)
+        acqd.kdata = acqd.kdata[:,idx_slice,:]
 
-    if !isnothing(nav) && !isnothing(nav_time)
+        if !isnothing(nav) && !isnothing(nav_time)
 
-        nav = nav[:,:,:,idx_slice]
-        nav_time = nav_time[:,idx_slice]
+            nav = nav[:,:,:,idx_slice]
+            nav_time = nav_time[:,idx_slice]
 
+        end
     end
 
+    return nav, nav_time
+
 end

src/CoilSensMap.jl

@@ -162,35 +162,54 @@ function ResizeSensit!(sensit::Array{Complex{T},4}, acqMap::AcquisitionData, acq
         @warn "The coil sensitivity maps have already been resized, the function cannot be executed."
 
     elseif freq_enc_FoV[1] < freq_enc_FoV[2] || phase_enc_FoV[1] < phase_enc_FoV[2]
-        @error "The reference data field of view is smaller than the image data field of view."
-
-    else
-        freq_enc_FoV_disc = round(Int64, (freq_enc_FoV[1] - freq_enc_FoV[2]) / (freq_enc_FoV[1]/freq_enc_samples[1]) / 2)
-        phase_enc_FoV_disc = round(Int64, (phase_enc_FoV[1] - phase_enc_FoV[2]) / (phase_enc_FoV[1]/phase_enc_samples[1]) / 2)
-
-        # Remove the sensit data that are not included in the image data FoV
-        sensit = sensit[freq_enc_FoV_disc+1:end-freq_enc_FoV_disc, phase_enc_FoV_disc+1:end-phase_enc_FoV_disc, :, :]
-
-        # Interpolate the sensit data to the image data
-        cartes_index = findall(x -> x!=0, sensit)
-        mask = zeros(Float32, size(sensit)) # compute the mask
-        for ii in cartes_index
-            mask[ii] = 1
-        end
-
-        # Linear interpolation
-        sensit = mapslices(x ->imresize(x, (freq_enc_samples[2], phase_enc_samples[2])), sensit, dims=[1,2])
-        mask = mapslices(x ->imresize(x, (freq_enc_samples[2], phase_enc_samples[2])), mask, dims=[1,2])
-
-        # Remove interpolated outline
-        cartes_index = findall(x -> x!=1, mask)
-        for ii in cartes_index
-            mask[ii] = 0
-        end
-
-        sensit = mask .* sensit
+        @warn "The reference data field of view is smaller than the image data field of view."
 
+        expand_freq = 0
+        expand_phase = 0
+        if freq_enc_FoV[1] < freq_enc_FoV[2]
+            expand_freq = ceil(Int64, (freq_enc_FoV[2] - freq_enc_FoV[1]) / (freq_enc_FoV[1]/freq_enc_samples[1]) / 2)
+        elseif phase_enc_FoV[1] < phase_enc_FoV[2]
+            expand_phase = ceil(Int64, (phase_enc_FoV[2] - phase_enc_FoV[1])/ (phase_enc_FoV[1]/phase_enc_samples[1]) / 2)
+        end        
+
+        expand_phase_vec = zeros(typeof(sensit[1,1,1,1]),  sizeSensit[1], expand_phase, sizeSensit[3], sizeSensit[4])
+        sensit = cat(expand_phase_vec, sensit, expand_phase_vec, dims = 2)
+
+        expand_freq_vec = zeros(typeof(sensit[1,1,1,1]), expand_freq, sizeSensit[2], sizeSensit[3], sizeSensit[4])
+        sensit = cat(expand_freq_vec, sensit, expand_freq_vec, dims = 1)
+
+        freq_enc_FoV[1] = freq_enc_FoV[1] + 2 * expand_freq * freq_enc_FoV[1]/freq_enc_samples[1]
+        freq_enc_samples[1] = freq_enc_samples[1] + 2 * expand_freq
+        phase_enc_FoV[1] = phase_enc_FoV[1] + 2 * expand_phase * phase_enc_FoV[1]/phase_enc_samples[1]
+        phase_enc_samples[1] = phase_enc_samples[1] + 2 * expand_phase
+
+    end
+
+    freq_enc_FoV_disc = round(Int64, (freq_enc_FoV[1] - freq_enc_FoV[2]) / (freq_enc_FoV[1]/freq_enc_samples[1]) / 2)
+    phase_enc_FoV_disc = round(Int64, (phase_enc_FoV[1] - phase_enc_FoV[2]) / (phase_enc_FoV[1]/phase_enc_samples[1]) / 2)
+
+    # Remove the sensit data that are not included in the image data FoV
+    sensit = sensit[freq_enc_FoV_disc+1:end-freq_enc_FoV_disc, phase_enc_FoV_disc+1:end-phase_enc_FoV_disc, :, :]
+
+    # Interpolate the sensit data to the image data
+    cartes_index = findall(x -> x!=0, sensit)
+    mask = zeros(Float32, size(sensit)) # compute the mask
+    for ii in cartes_index
+        mask[ii] = 1
     end
+
+    # Linear interpolation
+    sensit = mapslices(x ->imresize(x, (freq_enc_samples[2], phase_enc_samples[2])), sensit, dims=[1,2])
+    mask = mapslices(x ->imresize(x, (freq_enc_samples[2], phase_enc_samples[2])), mask, dims=[1,2])
+
+    # Remove interpolated outline
+    cartes_index = findall(x -> x!=1, mask)
+    for ii in cartes_index
+        mask[ii] = 0
+    end
+
+    sensit = mask .* sensit
+
 end
 
 

src/NavData.jl

@@ -8,7 +8,7 @@ mutable struct additionalNavInput
     numlines::Int64
     TE_nav::Float64
     dt_nav::Float64
-    freq_enc_FoV::Union{Array{Int64}, Nothing}
+    freq_enc_FoV::Union{Array{Float64}, Nothing}
     freq_enc_samples::Union{Array{Int64}, Nothing}
     phase_enc_samples::Union{Array{Int64}, Nothing}
     nav_time::Union{Array{Float64, 2}, Nothing}
@@ -70,6 +70,19 @@ function additionalNavInput(
 
     if !isnothing(acqMap) && !isnothing(acqData)
         (freq_enc_FoV, freq_enc_samples, phase_enc_FoV, phase_enc_samples) = Find_scaling_sensit(acqMap, acqData)
+        expand_freq = 0
+        expand_phase = 0
+        if freq_enc_FoV[1] < freq_enc_FoV[2]
+            @warn "The reference data field of view is smaller than the image data field of view."
+            expand_freq = ceil(Int64, (freq_enc_FoV[2] - freq_enc_FoV[1]) / (freq_enc_FoV[1]/freq_enc_samples[1]))
+        elseif phase_enc_FoV[1] < phase_enc_FoV[2]
+            @warn "The reference data field of view is smaller than the image data field of view."
+            expand_phase = ceil(Int64, (phase_enc_FoV[2] - phase_enc_FoV[1])/ (phase_enc_FoV[1]/phase_enc_samples[1]))
+        end
+        freq_enc_FoV[1] = freq_enc_FoV[1] + expand_freq * freq_enc_FoV[1]/freq_enc_samples[1]
+        freq_enc_samples[1] = freq_enc_samples[1] + expand_freq
+        #phase_enc_FoV[1] = phase_enc_FoV[1] + expand_phase * phase_enc_FoV[1]/phase_enc_samples[1]
+        phase_enc_samples[1] = phase_enc_samples[1] + expand_phase
     end
 
     return additionalNavInput(numslices, numechoes, numsamples, numlines, TE_nav, dt_nav,

src/NavParameters.jl

@@ -6,15 +6,16 @@ export defaultNavParams
 Define default parameters for data loading, navigator correction and image reconstruction.
 
 # Default parameters options are
-* `slices::Union{Nothing, Vector}`    - number of the slices to be loaded, nothing means all slices
-* `echoes::Union{Nothing, Vector}`    - number of the echoes to be loaded, nothing means all echoes
+* `slices::Union{Nothing, Vector}`    - a vector containing the number of the slices to be loaded, nothing means all slices
+* `echoes::Union{Nothing, Vector}`    - a vector containing the number of the echoes to be loaded, nothing means all echoes
 * `rep::Int`                          - repetition to be loaded, the first repetition is 0. It is mandatory to select one
 * `comp_sensit::Bool`                 - compute the sensitivity maps using the reference scan
 * `comp_centerline::Bool`             - use the Spinal Cord Toolbox (SCT) to find the centerlne position
 * `trust_SCT::Bool`                   - trust SCT or display the resutls and wait for user feedback with the julia REPL
 * `use_centerline::Bool`              - use the spinal cord centerline information in the navigator-based correction
 * `corr_type::String`                 - correction type. Options: "none", "knav", "FFT", "FFT_unwrap"
 * `FFT_interval::String`              - interval in mm to be considered for the FFT based approach
+* `mask_thresh::String`               - masking threshold: increase for reduced mask size, decrease for extended mask size
 
 # Additional required parameters are
 * `path_imgData::String`              - path to the image data file in ISMRMRD format
@@ -26,7 +27,7 @@ Define default parameters for data loading, navigator correction and image recon
 # Additional optional parameters are
 * `path_niftiMap::String`             - path to the file where the reconstructed reference data will be saved in nifti format. The file extension must be .nii
 * `path_centerline::String`           - path to the folder where the Spinal Cord Toolbox (SCT) centerline results will be saved
-* `path_physio::String`              - path to the physiological trace recording in .mat format. The variable should be a two columns vector (1:time [ms], 2:trace).
+* `path_physio::String`              - path to the physiological trace recording in .mat format. The variable should be a two columns vector (1:time [ms], 2:trace) called data.
                                         The time should be expressed in seconds from the beginning of the day and contain time points before and after the image acquisiton (at least 4 s).
 
 ISMRMRD reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.26089
@@ -44,6 +45,7 @@ function defaultNavParams()
     params[:use_centerline] = true
     params[:corr_type] = "FFT"
     params[:FFT_interval] = 35 # [millimiters]
+    params[:mask_thresh] = 0.13
 
     return params
   end

src/SpineCenterline.jl

@@ -2,20 +2,21 @@ export ReconstructSaveMap, ReconstructMap, niftiSaveImg, callSCT, findCenterline
 
 
 """
-    ReconstructSaveMap(path_nifti::String, path_ref::String)
+    ReconstructSaveMap(path_nifti::String, path_ref::String, thresh::Float64)
 
 Reconstruct the coil sensitivity map using the MRIReco.jl function and save it in nifti format without spatial information.
 
 # Arguments
 * `path_nifti::String` - path of the nifti file. The file must have .nii extension
 * `path_rep::String` - path of reference data in ISMRMRD format
+* `thresh::Float64` - masking threshold: increase for reduced mask size, decrease for extended mask size
 
 MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
 ISMRMRD reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.26089
 """
-function ReconstructSaveMap(path_nifti::String, path_ref::String)
+function ReconstructSaveMap(path_nifti::String, path_ref::String, thresh = 0.13)
 
-    (img, acq) = ReconstructMap(path_ref)
+    (img, acq) = ReconstructMap(path_ref, thresh)
     start_voxel = div(acq.encodingSize[1] -  acq.encodingSize[2], 2)
     img = img[start_voxel+1 : start_voxel+acq.encodingSize[2], :, :]
 
@@ -31,17 +32,18 @@ Reconstruct the coil sensitivity map using the MRIReco.jl function.
 
 # Arguments
 * `path_rep::String` - path of reference data in ISMRMRD format
+* `thresh::Float64` - masking threshold: increase for reduced mask size, decrease for extended mask size
 
 MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
 ISMRMRD reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.26089
 """
-function ReconstructMap(path_ref::String)
+function ReconstructMap(path_ref::String, thresh = 0.13)
 
     raw = RawAcquisitionData(ISMRMRDFile(path_ref))
     OrderSlices!(raw)
 
     acq = AcquisitionData(raw, estimateProfileCenter=true)
-    sensit = CompSensit(acq)
+    sensit = CompSensit(acq, thresh)
     img = Reconstruct(acq, sensit)
 
     return img, acq
@@ -142,13 +144,13 @@ SCT reference: https://spinalcordtoolbox.com
 """
 function findCenterline(params::Dict{Symbol, Any})
 
-    @info "Reco and Save"
-    # reconstruct and save in nifti the refence data
-    ReconstructSaveMap(params[:path_niftiMap], params[:path_refData])
-
-    @info "Find SC Centerline"
-    # find the spinal cord centerline on the reconstructed reference data
     if params[:comp_centerline] == true
+        @info "Reco reference scan and Save"
+        # reconstruct and save in nifti format the refence data
+        ReconstructSaveMap(params[:path_niftiMap], params[:path_refData], params[:mask_thresh])
+
+        @info "Find SC Centerline"
+        # find the spinal cord centerline on the reconstructed reference data
         callSCT(params)
     end
 

src/main.jl

@@ -27,7 +27,7 @@ function runNavPipeline(params::Dict{Symbol, Any})
 
     # slice and echo selection on acquisition data
     selectEcho!(acqData, params[:echoes])
-    selectSlice!(acqData, params[:slices], nav, nav_time)
+    (nav, nav_time) = selectSlice!(acqData, params[:slices], nav, nav_time)
 
     @info "read ref data"
     # read reference data
@@ -38,13 +38,15 @@ function runNavPipeline(params::Dict{Symbol, Any})
     @info "sensemaps"
     ## compute or load the coil sensitivity map
     if params[:comp_sensit]
-        CompResizeSaveSensit(acqMap, acqData, params[:path_sensit])
+        CompResizeSaveSensit(acqMap, acqData, params[:path_sensit], params[:mask_thresh])
     end
 
     #Load coil sensitivity
     sensit = FileIO.load(params[:path_sensit], "sensit")
-    sensit = reshape(sensit[:,:,params[:slices],:],(size(sensit,1), size(sensit,2),
-        size(params[:slices],1), size(sensit,4)))
+    if !isnothing(params[:slices])
+        sensit = reshape(sensit[:,:,params[:slices],:],(size(sensit,1), size(sensit,2),
+            size(params[:slices],1), size(sensit,4)))
+    end
 
     # Load centerline (ON LINUX: file is centerline.csv, ON WINDOWS AND MAC: is centerline.nii.csv)
     centerline = nothing
@@ -59,13 +61,15 @@ function runNavPipeline(params::Dict{Symbol, Any})
             end
         end
         centerline = centerline.Column1
-        centerline = centerline[params[:slices]]
+        if !isnothing(params[:slices])
+            centerline = centerline[params[:slices]]
+        end
     end
 
     #Load trace
     trace = nothing
     if params[:corr_type] == "FFT_unwrap"
-        trace = read(matopen(params[:path_physio] * string(params[:rep]+1) * ".mat"), "data")
+        trace = read(matopen(params[:path_physio]), "data")
     end
 
     @info "nav corr"
@@ -101,8 +105,7 @@ function saveNoise(path_imgData::String, path_noise::String)
     @info "Load first rep, save noise acquisition"
     # load the first repetition, slice and echo and save the noise acquisition for optimal results
     # the noise acquisition is saved in the first repetition only
-    rawData = RawAcquisitionData(ISMRMRDFile(path_imgData),
-            slice = 0, contrast = 0, repetition = 0)
+    rawData = RawAcquisitionData(ISMRMRDFile(path_imgData), repetition = 0)
     noisemat = ExtractNoiseData!(rawData)
     FileIO.save(path_noise,"noisemat",noisemat)