Initial package setup

.gitignore

@@ -1,4 +1,5 @@
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
-/Manifest.toml
+Manifest.toml
+.vscode

.vscode/settings.json

@@ -0,0 +1 @@
+{}

Project.toml

@@ -3,11 +3,29 @@ uuid = "eacce8a8-16e1-4f0a-a125-bfe7a196bd66"
 authors = ["Laura Beghini"]
 version = "0.1.0"
 
+[deps]
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+MRIBase = "f7771a9a-6e57-4e71-863b-6e4b6a2f17df"
+MRICoilSensitivities = "c57eb701-aafc-44a2-a53c-128049758959"
+MRIFiles = "5a6f062f-bf45-497d-b654-ad17aae2a530"
+MRIReco = "bdf86e05-2d2b-5731-a332-f3fe1f9e047f"
+NIfTI = "a3a9e032-41b5-5fc4-967a-a6b7a19844d3"
+PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
+Scratch = "6c6a2e73-6563-6170-7368-637461726353"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
-julia = "1.1"
+julia = "1.3"
 
 [extras]
+LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["LazyArtifacts", "Scratch", "Test"]

README.md

@@ -2,3 +2,19 @@
 
 [![Build Status](https://github.com/Laura2305/MRINavigator.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Laura2305/MRINavigator.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/Laura2305/MRINavigator.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/Laura2305/MRINavigator.jl)
+
+MRINavigator.jl provides multiple navigator-based correction pipelines for magnetic resonance data. These aim at demodulating time dependent field variations. The package was developed with a focus on spinal cord imaging, however it can be used for multiple imaging applications. 
+The corrections are to be applied on the raw data before the image reconstruction. MRIReco.jl can be used to reconstruct the images.
+
+More details can be found in the Online Documentation (soon to be published).
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+## How to give credit
+
+soon to be poublished
+
+
+## Community Standards
+
+This project is part of the Julia community and follows the [Julia community standards](https://julialang.org/community/standards/). 

src/AdjustData.jl

@@ -0,0 +1,290 @@
+export OrderSlices!, ExtractNoiseData!, ReverseBipolar!, RemoveRef!, CopyTE!, AdjustSubsampleIndices!, ExtractNavigator, ExtractFlags, selectEcho!, selectSlice!
+
+"""
+    OrderSlices!(rawData::RawAcquisitionData)
+
+Spatially order the slices in the MRIReco.jl raw data structure.
+The slices are ordered basing on the position coordinates saved in each profile.
+If these are not present the slices can not be ordered.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function OrderSlices!(rawData::RawAcquisitionData)
+
+    total_num = length(rawData.profiles)
+    slices = zeros(typeof(rawData.profiles[1].head.position[3]), total_num)
+
+    for ii = 1:total_num
+        slices[ii] = rawData.profiles[ii].head.position[3]
+    end
+
+    unique!(slices)
+    slices_indx = sortperm(sortperm(slices))
+
+    for ii = 1:total_num
+        index = rawData.profiles[ii].head.position[3] .== slices
+        rawData.profiles[ii].head.idx.slice = slices_indx[index][1]-1
+    end
+
+end
+
+
+"""
+    flags = ExtractFlags(rawData::RawAcquisitionData) 
+
+Extract the acquisition flags from the MRIReco.jl raw data profiles.
+Return a 31 elements vector for each profile.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function ExtractFlags(rawData::RawAcquisitionData)
+
+    total_num = length(rawData.profiles)
+    flags = zeros(Int64, total_num, 31)
+
+    for ii=1:total_num
+        flags[ii,:] = digits(rawData.profiles[ii].head.flags, base = 2, pad=31)
+    end
+
+    return flags
+
+end
+
+
+"""
+    noisemat = ExtractNoiseData!(rawData::RawAcquisitionData, flags::Array{Int64})
+
+Extract and return the noise acquisition from the MRIReco.jl raw data.
+The noise acquisition is usually the first profile with slice = 0, contrast = 0, repetition = 0.
+The noise profile should have the 19th flag element qual to 1. Check with ExtractFlags if errors occur.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function ExtractNoiseData!(rawData::RawAcquisitionData)
+
+    flags = ExtractFlags(rawData)
+    total_num = length(rawData.profiles)
+    if total_num != size(flags, 1)
+        @error "size of flags and number of profiles in rawData do not match"
+    end
+    noisemat = Matrix{typeof(rawData.profiles[1].data)}
+
+    for ii=1:total_num
+
+        if flags[ii,19] == true
+            noisemat = rawData.profiles[ii].data
+            deleteat!(rawData.profiles, ii)
+            break
+        end
+
+    end
+
+    return noisemat
+
+end
+
+
+"""
+    ReverseBipolar!(rawData::RawAcquisitionData)
+
+Reflect the MRIReco.jl raw data profiles for bipolar acquisition.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function ReverseBipolar!(rawData::RawAcquisitionData)
+
+    flags = ExtractFlags(rawData)
+    total_num = length(rawData.profiles)
+    if total_num != size(flags, 1)
+        @error "size of flags and number of profiles in rawData do not match"
+    end
+
+    for ii=1:total_num
+
+        if flags[ii,22] == true
+            reverse!(rawData.profiles[ii].data, dims=1)
+            rawData.profiles[ii].head.flags=rawData.profiles[ii].head.flags-(2^21)
+        end
+
+    end
+
+end
+
+
+"""
+    RemoveRef!(rawData::RawAcquisitionData, slices::Union{Vector{Int64}, Nothing}, echoes::Union{Vector{Int64}, Nothing})
+
+Remove reference data that are not useful for the navigator-based crrection from acquisitions with phase stabilization on Siemens scanners.
+Not solid to recalls.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function RemoveRef!(rawData::RawAcquisitionData)
+
+    numSlices = rawData.params["enc_lim_slice"].maximum + 1
+    numEchoes = size(rawData.params["TE"],1) + 1
+
+    #Apply this only if using phase stabilizaion
+    removeIndx = numSlices*(numEchoes)
+    deleteat!(rawData.profiles, 1:removeIndx)
+
+end
+
+
+"""
+    CopyTE!(rawData::RawAcquisitionData, acqData::AcquisitionData)
+
+Copy the TE values from the MRIReco.jl raw data structure to the acquisition data structure.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+* `acqData::AcquisitionData` - acquisition data structure obtained converting raw data with MRIReco.jl
+"""
+function CopyTE!(rawData::RawAcquisitionData, acqData::AcquisitionData)
+
+    for ii=1:size(acqData.kdata)[1]
+        acqData.traj[ii].TE = rawData.params["TE"][ii]
+    end
+
+end
+
+
+"""
+    AdjustSubsampleIndices!(acqData::AcquisitionData)
+
+Add subsamples indices in the MRIReco.jl acquisition data structure.
+Needed when conveting data not acquired in the first repetition.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `acqData::AcquisitionData` - acquisition data structure obtained converting raw data with MRIReco.jl
+"""
+function AdjustSubsampleIndices!(acqData::AcquisitionData)
+
+    if isempty(acqData.subsampleIndices[1])
+        for ii = 1:size(acqData.subsampleIndices)[1]
+            acqData.subsampleIndices[ii]=1:size(acqData.kdata[1,1,1])[1]
+        end
+    end
+
+end
+
+
+"""
+    (nav, nav_time) = ExtractNavigator(rawData::RawAcquisitionData, slices::Union{Vector{Int64}, Nothing})
+
+Extract the navigator profiles from the MRIReco.jl raw data structure.
+These are registered with the same indices (contract, slice, encoding step) as the image data for the first echo time.
+Return a navigator array and a navigator time array. The navigator array has four dimensions in order: k-space samples, coils, k-space lines, slices.
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+
+# Arguments
+* `rawData::RawAcquisitionData` - raw data structure obtained loading raw data with MRIReco.jl
+"""
+function ExtractNavigator(rawData::RawAcquisitionData)
+
+    total_num = length(rawData.profiles)
+    numberslices = rawData.params["enc_lim_slice"].maximum +1
+    contrasts = zeros(Int64, total_num)
+    slices = zeros(Int64, total_num)
+    lines = zeros(Int64, total_num)
+    for ii = 1:total_num
+        contrasts[ii] = rawData.profiles[ii].head.idx.contrast
+        slices[ii] = rawData.profiles[ii].head.idx.slice
+        lines[ii] = rawData.profiles[ii].head.idx.kspace_encode_step_1
+    end
+    # keep only the indexes of data saved in the first echo (this includes navigator)
+    contrastsIndx = findall(x->x==0, contrasts)
+    slices = slices[contrastsIndx]
+    lines = lines[contrastsIndx]
+
+    nav = zeros(ComplexF32, size(rawData.profiles[1].data)[1], size(rawData.profiles[1].data)[2],
+        rawData.params["reconSize"][2], numberslices)
+
+    nav_time = zeros(Float64,
+        rawData.params["reconSize"][2], numberslices)
+    #Odd indexes are data first echo, Even indexes are navigator data
+    for ii = 2:2:length(slices)
+        nav[:,:,lines[ii]+1,slices[ii]+1] = rawData.profiles[contrastsIndx[ii]].data
+        nav_time[lines[ii]+1,slices[ii]+1] = rawData.profiles[contrastsIndx[ii]].head.acquisition_time_stamp
+    end
+    #Remove the rows filled with zeroes
+    lines = unique(lines) .+1
+    nav = nav[:,:,lines,:]
+    nav_time = nav_time[lines,:]
+
+    return nav, nav_time
+    #navigator[k-space samples, coils, k-space lines, slices]
+
+end
+
+"""
+    SelectEcho!(acqd, idx_echo)
+
+Extract one or more echoes from the acquisition data structure
+
+# Arguments
+* `acqd::AcquisitionData` - acquisition data structure obtained converting raw data with MRIReco.jl
+* `idx_echo::Vector{Int64}` - vector containing the indexes of the echoes to be selected (starting from 0)
+
+MRIReco reference: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.28792
+"""
+function selectEcho!(acqd::AcquisitionData, idx_echo::Vector{Int64})
+
+    if !isempty(idx_echo)
+        contrasts = size(acqd.kdata)[1]
+        indices = Vector{Int64}(undef, contrasts)
+        for ii=1:contrasts
+            indices[ii] = ii
+        end
+        deleteat!(indices, idx_echo)
+        deleteat!(acqd.subsampleIndices, indices)
+        deleteat!(acqd.traj, indices)
+        acqd.kdata = acqd.kdata[idx_echo,:,:];
+    end
+
+end
+
+"""
+    SelectSlice!(acqd, nav, nav_time, idx_slice)
+
+Extract one or more echoes from the acquisition data structure
+
+# Arguments
+* `acqd::AcquisitionData` - acquisition data structure obtained converting raw data with MRIReco.jl
+* `idx_slice::Vector{Int64}` - vector containing the indexes of the slices to be selected (starting from 0, downer slice)
+
+# Optional arguments with default value = nothing
+* `nav::Union{Array{Complex{T}, 4}, Nothin} = nothing` - navigator profiles obtained with the ExtractNavigator function
+* `nav_time::Union{Array{Complex{Float32}, 2}, Nothing}` - time stamps for the navigator data obtained with ExtractNavigator (in ms from the beginning of the day)
+
+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}
+
+    acqd.kdata = acqd.kdata[:,idx_slice,:]
+    if !isnothing(nav) && !isnothing(nav_time)
+        nav = nav[:,:,:,idx_slice]
+        nav_time = nav_time[:,idx_slice]
+    end
+
+end