initial commit

.github/workflows/CI.yml

@@ -24,10 +24,11 @@ jobs:
       matrix:
         version:
           - '1.10'
-          - '1.9'
           - 'pre'
         os:
           - ubuntu-latest
+          - macos-latest
+          - windows-latest
         arch:
           - x64
     steps:

Project.toml

@@ -3,8 +3,17 @@ uuid = "4ba9dd0f-9419-4edc-b963-b8b02e27d02f"
 authors = ["aTrotier <[email protected]> and contributors"]
 version = "0.0.1"
 
+[deps]
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+MRICoilSensitivities = "c57eb701-aafc-44a2-a53c-128049758959"
+MRIFiles = "5a6f062f-bf45-497d-b654-ad17aae2a530"
+MRIReco = "bdf86e05-2d2b-5731-a332-f3fe1f9e047f"
+NIfTI = "a3a9e032-41b5-5fc4-967a-a6b7a19844d3"
+QuantitativeMRI = "659767e3-31a7-4dc1-8563-6e03f484b231"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/BIDS.jl

@@ -0,0 +1,90 @@
+export write_bids_MP2RAGE
+
+"""
+      write_bids_MP2RAGE(d::Dict,subname::AbstractString,folder="")
+
+This function writes data from a dictionary (`d`) in BIDS (Brain Imaging Data Structure) format for MP2RAGE acquisitions.
+
+**Arguments:**
+
+* `d` (Dict): A dictionary containing the data to be written. Expected keys:
+    * `im_reco` (Array): 5D array containing the reconstructed images.
+    * `MP2RAGE` (Array): T1 map image.
+    * `T1maps` (Array): Additional T1 map data (optional).
+    * `params_prot` (Dict): Dictionary containing acquisition parameters.
+    * `params_MP2RAGE` (Dict): Dictionary containing MP2RAGE specific parameters.
+* `subname` (AbstractString): The name of the subject.
+* `folder` (AbstractString, optional): The folder where the BIDS data will be written. Defaults to the current directory.
+
+**Functionality:**
+
+1. Creates a directory structure for the anatomical data under `folder/subname/anat`.
+2. Defines a list of file paths for different image types associated with MP2RAGE acquisitions.
+3. Extracts relevant data from the dictionary `d` for each image type.
+4. Creates NIfTI volumes (`NIVolume`) with the extracted data and specified voxel size from `d["params_prot"]`.
+5. Writes each NIfTI volume to a compressed file (`.nii.gz`) in the corresponding directory.
+6. Extracts acquisition parameters from `d`.
+7. Creates a dictionary (`JSON_dict`) containing these parameters in BIDS format:
+    * `InversionTime`: List of inversion times (TI1, TI2) in seconds.
+    * `RepetitionTimeExcitation`: Repetition time (TR) in seconds.
+    * `RepetitionTimePreparation`: MP2RAGE specific repetition time (MP2RAGE_TR) in seconds.
+    * `NumberShots`: Echo train length (ETL).
+    * `FlipAngle`: List of flip angles (alpha1, alpha2) in degrees.
+    * `MagneticFieldStrength`: Magnetic field strength in Tesla.
+    * `Units`: Units for the data (set to "arbitrary" in this case).
+8. Writes the JSON dictionary to a file named `MP2RAGE.json` in the `folder/subname` directory.
+
+**Note:** This function assumes the dictionary `d` contains the necessary data in the specified format. 
+"""
+function write_bids_MP2RAGE(d::Dict,subname::AbstractString,folder="")
+
+  path_anat = joinpath(folder,subname,"anat")
+  mkpath(path_anat)
+
+
+  path_type = ["_inv-1-mag_MP2RAGE",
+              "_inv-1-phase_MP2RAGE",
+              "_inv-1-complex_MP2RAGE",
+              "_inv-2-mag_MP2RAGE",
+              "_inv-2-phase_MP2RAGE",
+              "_inv-2-complex_MP2RAGE",
+              "_UNIT1",
+              "_T1map"]
+
+  data_ = [ abs.(d["im_reco"][:,:,:,1,1,:]),
+            angle.(d["im_reco"][:,:,:,1,1,:]),
+            d["im_reco"][:,:,:,1,1,:],
+            abs.(d["im_reco"][:,:,:,2,1,:]),
+            angle.(d["im_reco"][:,:,:,2,1,:]),
+            d["im_reco"][:,:,:,2,1,:],
+            d["MP2RAGE"],
+            d["T1maps"]]
+
+  voxel_size = tuple(parse.(Float64,d["params_prot"]["PVM_SpatResol"])...) #mm
+  for (name,data) in zip(path_type, data_)
+    ni = NIVolume(data,voxel_size=voxel_size)
+    niwrite(joinpath(path_anat,subname*name*".nii.gz"),ni)
+  end
+
+
+  # pass parameters
+  b["ACQ_operator"] # required to read ACQ
+  MagneticField = parse(Float64,b["BF1"]) / 42.576
+  p_MP2 = d["params_MP2RAGE"]
+
+  # define JSON dict
+  JSON_dict = Dict{Any,Any}()
+  JSON_dict["InversonTime"]= [p_MP2.TI₁,p_MP2.TI₂] #s
+  JSON_dict["RepetitionTimeExcitation"]= p_MP2.TR
+  JSON_dict["RepetitionTimePreparation"]= p_MP2.MP2RAGE_TR
+  JSON_dict["NumberShots"]= p_MP2.ETL
+  JSON_dict["FlipAngle"]= [p_MP2.α₁,p_MP2.α₂]
+  JSON_dict["MagneticFieldStrength"] = MagneticField
+  JSON_dict["Units"] = "arbitrary"
+
+  # Write the dictionary to a JSON file
+  open(joinpath(folder,subname,"MP2RAGE.json"), "w") do f
+    JSON.print(f, JSON_dict, 4)  # Indent 4 spaces for readability
+  end
+
+end

src/SEQ_BRUKER_a_MP2RAGE_CS_360.jl

@@ -1,5 +1,16 @@
 module SEQ_BRUKER_a_MP2RAGE_CS_360
 
-# Write your package code here.
+# required module
+using MRIReco
+using MRIFiles
+using MRICoilSensitivities
+using QuantitativeMRI
+using Statistics
+using JSON
+using NIfTI
 
+# Write your package code here.
+include("bruker_sequence.jl")
+include("reconstruction.jl")
+include("BIDS.jl")
 end

src/bruker_sequence.jl

@@ -0,0 +1,112 @@
+export RawAcquisitionData_MP2RAGE
+
+"""
+    RawAcquisitionData_MP2RAGE(b::BrukerFile)
+
+Convert a Bruker dataset acquired with the a_MP2RAGE_CS_360 sequence into a 
+`RawAcquisitionData` object compatible with the MRIReco functions.
+
+Input : 
+    - b::BrukerFile
+
+Output :
+    - raw::RawAcquisitionData
+"""
+function RawAcquisitionData_MP2RAGE(b::BrukerFile)
+  T = Complex{MRIFiles.acqWordSize(b)}
+
+  filename = joinpath(b.path, "rawdata.job0")
+
+  N = MRIFiles.pvmMatrix(b)
+  factor_AntiAlias = MRIFiles.pvmAntiAlias(b)
+
+  N = round.(Int,N .* factor_AntiAlias)
+
+  if length(N) < 3
+    N_ = ones(Int,3)
+    N_[1:length(N)] .= N
+    N = N_
+  end
+
+  numChannel = parse.(Int,b["PVM_EncNReceivers"])
+  numAvailableChannel = MRIFiles.pvmEncAvailReceivers(b)
+  numSlices = MRIFiles.acqNumSlices(b)
+  numEchos = parse(Int,b["PVM_NEchoImages"])
+  ETL = parse(Int,b["MP2_EchoTrainLength"]) 
+  numRep = MRIFiles.acqNumRepetitions(b)
+  numEncSteps = MRIFiles.acqSpatialSize1(b) 
+
+  profileLength = N[1]*numChannel#MRIFiles.acqSize(b)[1] #Int((ceil(N[1]*numChannel*sizeof(dtype)/1024))*1024/sizeof(dtype)) # number of points + zeros
+
+  I = open(filename,"r") do fd
+    read!(fd,Array{T,6}(undef, profileLength,
+                                  ETL,
+                                   numEchos,
+                                   numSlices,
+                                   div(numEncSteps, ETL),
+                                   numRep))
+  end
+
+  encSteps1 = parse.(Int,b["PVM_EncGenSteps1"]).+round(Int,N[2]/2)
+  encSteps2 = parse.(Int,b["PVM_EncGenSteps2"]).+round(Int,N[3]/2)
+
+  objOrd = MRIFiles.acqObjOrder(b)
+  objOrd = objOrd.-minimum(objOrd)
+
+  gradMatrix = MRIFiles.acqGradMatrix(b)
+
+  offset1 = MRIFiles.acqReadOffset(b)
+  offset2 = MRIFiles.acqPhase1Offset(b)
+
+  offset3 = ndims(b) == 2 ? MRIFiles.acqSliceOffset(b) : MRIFiles.pvmEffPhase2Offset(b)
+
+  profiles = Profile[]
+  for nR = 1:numRep
+    for nEnc = 1:div(numEncSteps, ETL)
+      for nSl = 1:numSlices
+        for nEcho=1:numEchos
+          for nETL = 1:ETL
+              counter = EncodingCounters(kspace_encode_step_1=encSteps1[nETL+ETL*(nEnc-1)],
+                                          kspace_encode_step_2=encSteps2[nETL+ETL*(nEnc-1)],
+                                          average=0,
+                                          slice=objOrd[nSl],
+                                          contrast=nEcho-1,
+                                          phase=0,
+                                          repetition=nR-1,
+                                          set=0,
+                                          segment=0 )
+
+              G = gradMatrix[:,:,nSl]
+              read_dir = (G[1,1],G[2,1],G[3,1])
+              phase_dir = (G[1,2],G[2,2],G[3,2])
+              slice_dir = (G[1,3],G[2,3],G[3,3])
+
+              # Not sure if the following is correct...
+              pos = offset1[nSl]*G[:,1] +
+                    offset2[nSl]*G[:,2] +
+                    offset3[nSl]*G[:,3]
+
+              position = (pos[1], pos[2], pos[3])
+
+              head = AcquisitionHeader(number_of_samples=N[1], idx=counter,
+                                        read_dir=read_dir, phase_dir=phase_dir,
+                                        slice_dir=slice_dir, position=position,
+                                        center_sample=div(N[1],2),
+                                        available_channels = numChannel, #numAvailableChannel ?
+                                        active_channels = numChannel)
+              traj = Matrix{Float32}(undef,0,0)
+              dat = map(T, reshape(I[:,nETL,nEcho,nSl,nEnc,nR],:,numChannel))
+              push!(profiles, Profile(head,traj,dat) )
+          end
+        end
+      end
+    end
+  end
+
+  params = MRIFiles.brukerParams(b)
+  params["trajectory"] = "cartesian"
+
+  params["encodedSize"] = N
+
+  return RawAcquisitionData(params, profiles)
+end

src/reconstruction.jl

@@ -0,0 +1,57 @@
+export reconstruction_MP2RAGE
+
+function reconstruction_MP2RAGE(path_bruker;mean_NR::Bool = false)
+  b = BrukerFile(path_bruker)
+  raw = RawAcquisitionData_MP2RAGE(b)
+  acq = AcquisitionData(raw,OffsetBruker=true)
+
+  ncalib = minimum(parse.(Int,b["PVM_EncCSNumRef"]))
+
+  if ncalib > 24
+    ncalib = 24
+  end
+  sens_fully = espirit(acq,eigThresh_2 = 0,(6,6,6),ncalib);
+
+  params = Dict{Symbol, Any}()
+  params[:reco] = "multiCoil"
+  params[:reconSize] = acq.encodingSize
+  params[:senseMaps] = sens_fully
+  params[:iterations] = 1
+
+  x_approx = reconstruction(acq, params).data
+  x_approx = permutedims(x_approx,(1,2,3,5,6,4))
+  if mean_NR
+    x_approx = mean(x_approx,dims=6)
+  end
+
+  ## process data to extract T1 maps
+  MP2 = mp2rage_comb(x_approx[:,:,:,:,1,:])
+
+  p = params_from_seq(b)
+
+  T1,range_T1 = mp2rage_T1maps(MP2,p)
+
+  d = Dict{Any,Any}()
+  d["im_reco"] = x_approx
+  d["MP2RAGE"] = MP2
+  d["T1maps"] = T1
+  d["range_T1"] = range_T1
+  d["params_reco"] = params
+  d["params_MP2RAGE"] = p
+  d["params_prot"] = b
+
+  return d
+end
+
+function params_from_seq(b::BrukerFile)
+  return ParamsMP2RAGE(
+      parse(Float64,b["EffectiveTI"][1]),
+      parse(Float64,b["EffectiveTI"][2]),
+      parse(Float64,b["PVM_RepetitionTime"]),
+      parse(Float64,b["MP2_RecoveryTime"]),
+      parse(Int,b["MP2_EchoTrainLength"]),
+      parse.(Float64,split(b["ExcPulse1"],", ")[3]),
+      parse.(Float64,split(b["ExcPulse2"],", ")[3]) # for now only one pulse is used
+  )
+end
+