Preproc-main into main

.gitignore

@@ -40,5 +40,10 @@ coverage.xml
 *.pyc
 phantoms/MR_XCAT_qMRI/*.json
 phantoms/MR_XCAT_qMRI/*.txt
+phantoms/brain/data/*.nii.gz
+phantoms/brain/ground_truth/*.nii.gz
 tests/IVIMmodels/unit_tests/models
 models
+temp
+.vscode
+.env

doc/code_contributions_record.csv

@@ -29,4 +29,5 @@ IVIM,Fitting,LSQ fitting,Part of the SUPER-IVIM-DC python package - trust region
 IVIM,Fitting,LSQ fitting,Part of the SUPER-IVIM-DC python package - segmented algorithm,TCML_TechnionIIT,TF2.4_IVIM-MRI_CodeCollection/src/original/TCML_TechnionIIT/SUPER-IVIM-DC,Noam Korngut/Elad Rotman/Onur Afacan Sila Kurugol/Yael Zaffrani-Reznikov/Shira Nemirovsky-Rotman/Simon Warfield/Moti Freiman,TCML Technion IIT,super_ivim_dc.source.Classsic_ivim_fit.IVIM_fit_sls,https://doi.org/10.1007/978-3-031-16434-7_71,TCML_TechnionIIT_lsqSLS
 IVIM,Fitting,LSQ fitting,Part of the SUPER-IVIM-DC python package - Levenberg Marquart lsq algorithm initialized with segmented fitting,TCML_TechnionIIT,TF2.4_IVIM-MRI_CodeCollection/src/original/TCML_TechnionIIT/SUPER-IVIM-DC,Noam Korngut/Elad Rotman/Onur Afacan Sila Kurugol/Yael Zaffrani-Reznikov/Shira Nemirovsky-Rotman/Simon Warfield/Moti Freiman,TCML Technion IIT,super_ivim_dc.source.Classsic_ivim_fit.IVIM_fit_sls_lm,https://doi.org/10.1007/978-3-031-16434-7_71,TCML_TechnionIIT_lsq_sls_lm
 IVIM,Fitting,LSQ fitting,Part of the SUPER-IVIM-DC python package - trust region lsq algorithm initialized with segmented fitting,TCML_TechnionIIT,TF2.4_IVIM-MRI_CodeCollection/src/original/TCML_TechnionIIT/SUPER-IVIM-DC,Noam Korngut/Elad Rotman/Onur Afacan Sila Kurugol/Yael Zaffrani-Reznikov/Shira Nemirovsky-Rotman/Simon Warfield/Moti Freiman,TCML Technion IIT,super_ivim_dc.source.Classsic_ivim_fit.IVIM_fit_sls_trf,https://doi.org/10.1007/978-3-031-16434-7_71,TCML_TechnionIIT_lsq_sls_trf
-IVIM,Fitting,Segmented linear fitting,Consensus algorithm,TF_reference,TF2.4_IVIM-MRI_CodeCollection/src/original/TF_reference/,Ben Neijndorff/OSIPI taskforce 2.4,OSIPI,segmented_IVIMfit,tba,TF_reference_IVIMfit
+IVIM,Fitting,Segmented linear fitting,Consensus algorithm,TF_reference,TF2.4_IVIM-MRI_CodeCollection/src/original/TF_reference/,Ben Neijndorff/OSIPI taskforce 2.4,OSIPI,segmented_IVIMfit,tba,TF_reference_IVIMfit
+IVIM,Preprocessing,Full pipeline,Brain pipeline,EP_GU,TF2.4_IVIM-MRI_CodeCollection/src/original/EP_GU,Elina Petersson,University of Gothenburg,brain_preproc,tba,

phantoms/brain/ground_truth/ballistic_groundtruth.json

@@ -1,6 +1,9 @@
 {
+    "tissues": ["WM","GM","CSF"],
     "D":[1.2e-3,0.98e-3,3e-3],
     "f":[0.0243,0.0164,0],
     "vd":[1.71,1.73,0],
-    "Db":1.75e-3
+    "Db":1.75e-3,
+    "T1":[1.0, 1.5, 3.4],
+    "T2":[0.046,0.068,0.45]
 }

phantoms/brain/ground_truth/diffusive_groundtruth.json

@@ -1,5 +1,8 @@
 {
+    "tissues": ["WM","GM","CSF"],
     "D":[0.81e-3,0.86e-3,3e-3],
     "f":[0.044,0.033,0],
-    "Dstar":[84e-3,76e-3,0]
-}
+    "Dstar":[84e-3,76e-3,0],
+    "T1":[1.0, 1.5, 3.4],
+    "T2":[0.046,0.068,0.45]
+} 

phantoms/brain/sim_brain_phantom.py

@@ -6,23 +6,28 @@
 from scipy.ndimage import zoom
 from utilities.data_simulation.Download_data import download_data
 
-if __name__ == "__main__":
+DIFFUSIVE_REGIME = 'diffusive'
+BALLISTIC_REGIME = 'ballistic'
+
+def simulate_brain_phantom(regime=DIFFUSIVE_REGIME, snr=100, sim_relaxation=True, TE=60e-3, TR=5, resolution=[3,3,3]):
+    '''
+    Simulation parameters can be set by changing the default values of the function arguments. 
+    The default values are chosen to be suitable for a diffusive regime phantom, but can be changed to simulate a ballistic regime phantom as well. 
+    The simulated image is saved in phantoms/brain/data with the name 'diffusive_sn{snr}.nii.gz' or 'ballistic_sn{snr}.nii.gz' depending on the regime. 
+    The corresponding bvals and cvals (if applicable) are also saved in the same folder.
+
+    regime: 'diffusive' or 'ballistic' to choose the type of phantom to simulate
+    snr: signal-to-noise ratio of the simulated image
+    sim_relaxation: whether to simulate T1 and T2 relaxation effects in the signal
+    TE: echo time in seconds
+    TR: repetition time in seconds
+    resolution: voxel size in mm
+    '''
     download_data()
 
-    DIFFUSIVE_REGIME = 'diffusive'
-    BALLISTIC_REGIME = 'ballistic'
 
     folder = os.path.dirname(__file__)
 
-    ###########################
-    #  Simulation parameters  #
-    regime = DIFFUSIVE_REGIME
-    snr = 200
-    resolution = [3,3,3]
-    #                         #
-    ###########################
-
-
     # Ground truth
     nii = nib.load(os.path.join(os.path.split(os.path.split(folder)[0])[0],'download','Phantoms','brain','ground_truth','hrgt_icbm_2009a_nls_3t.nii.gz'))
     segmentation = np.squeeze(nii.get_fdata()[...,-1])
@@ -40,20 +45,30 @@
 
     # Calculate signal
     S = np.zeros(list(np.shape(segmentation))+[b.size])
-
+    print(ivim_pars)
     if regime == BALLISTIC_REGIME:
         Db = ivim_pars["Db"]
-        for i,(D,f,vd) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["vd"])):
+        for i,(D,f,vd,T1,T2) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["vd"],ivim_pars["T1"],ivim_pars["T2"])):
             S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Db-c**2*vd**2))
+            if sim_relaxation:
+                S[segmentation==i+1,:] *= np.exp(-TE/T2)*(1-np.exp(-TR/T1))
     else:
-        for i,(D,f,Dstar) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"])):
+        for i,(D,f,Dstar,T1,T2) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"],ivim_pars["T1"],ivim_pars["T2"])):
             S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Dstar))
+            if sim_relaxation:
+                S[segmentation==i+1,:] *= np.exp(-TE/T2)*(1-np.exp(-TR/T1))
 
     # Resample to suitable resolution
     im = zoom(S,np.append(np.diag(nii.affine)[:3]/np.array(resolution),1),order=1)
     sz = im.shape
 
-    # Add noise
+    # Save image without noise for reference
+    nii_out = nib.Nifti1Image(im,np.eye(4))
+    base_name = os.path.join(folder,'data','{}_reference'.format(regime,snr))
+    nib.save(nii_out,base_name+'.nii.gz')
+    shutil.copyfile(bval_file,base_name+'.bval')
+
+    # Add Rician noise
     im_noise = np.abs(im + S0/snr*(np.random.randn(sz[0],sz[1],sz[2],sz[3])+1j*np.random.randn(sz[0],sz[1],sz[2],sz[3])))
 
     # Save as image and sequence parameters 
@@ -62,4 +77,14 @@
     nib.save(nii_out,base_name+'.nii.gz')
     shutil.copyfile(bval_file,base_name+'.bval')
     if regime == BALLISTIC_REGIME:
-        shutil.copyfile(cval_file,base_name+'.cval')
+        shutil.copyfile(cval_file,base_name+'.cval')
+
+    # Resample and save segmentation
+    segmentation = zoom(segmentation,np.diag(nii.affine)[:3]/np.array(resolution),order=0)
+    nii_out = nib.Nifti1Image(segmentation,np.eye(4))
+    base_name = os.path.join(folder,'data','{}_snr{}_mask'.format(regime,snr))
+    nib.save(nii_out,base_name+'.nii.gz')
+
+
+if __name__ == "__main__":
+    simulate_brain_phantom()

pytest.ini

@@ -3,6 +3,7 @@ markers =
     slow: marks tests as slow (deselect with '-m "not slow"')
 addopts =  
     -m 'not slow'
+    --ignore-glob=**/wip_*.py
 testpaths = tests
 filterwarnings =
     ignore::Warning

requirements.txt

@@ -18,4 +18,7 @@ sphinx_rtd_theme
 pytest-json-report
 statsmodels
 ivimnet
-nlopt
+nlopt
+nipype
+itk
+#ivim-mri

src/original/preprocessing/EP_GU/brain_pipeline.py

@@ -0,0 +1,184 @@
+'''
+Preprocessing of brain IVIM data
+Requirements: 
+- ivim
+- fsl
+- freesurfer
+- itk
+- itk-elastix
+'''
+#%%
+import numpy as np
+import nibabel as nib
+import subprocess, os
+from nipype.interfaces import fsl
+import itk
+from ivim.io.base import read_bval, write_bval
+from ivim.preproc.base  import extract, combine, average
+from ivim.preproc import signal_drift
+
+
+def denoise_wrap(im_file):
+    ''' Denoising '''
+    im_file_denoised = im_file.replace('.nii.gz','-pca.nii.gz')
+    noise_level = im_file.replace('.nii.gz','-noise_level.nii.gz')
+    print('Denoising...')
+    subprocess.run(['dwidenoise', '-noise', noise_level, im_file, im_file_denoised])
+    subprocess.run(['cp',im_file.replace('.nii.gz','.bval'),im_file_denoised.replace('.nii.gz','.bval')])
+    subprocess.run(['cp',im_file.replace('.nii.gz','.bvec'),im_file_denoised.replace('.nii.gz','.bvec')])
+    return im_file_denoised
+
+def degibbs_wrap(im_file):
+    ''' Gibbs ringing artifact removal '''
+    im_file_degibbs = im_file.replace('.nii.gz','-gib.nii.gz')
+    print('Degibbsing...')
+    subprocess.run(['mrdegibbs', im_file, im_file_degibbs])
+    subprocess.run(['cp',im_file.replace('.nii.gz','.bval'),im_file_degibbs.replace('.nii.gz','.bval')])
+    subprocess.run(['cp',im_file.replace('.nii.gz','.bvec'),im_file_degibbs.replace('.nii.gz','.bvec')])
+    return im_file_degibbs
+
+def extract_volume(in_file, t_index):
+    '''
+    Extracts a 3D volume at time index t_index from a 4D image.
+    Returns the path to the temporary file.
+    '''
+    out_file = in_file.replace('.nii.gz',f'-tmp{t_index}.nii.gz')
+    subprocess.run(['fslroi',in_file,out_file,str(t_index),'1'])
+    return out_file
+
+def register_3D_to_ref(in_file, ref_idx, out_file,method='rigid'):
+    '''
+    Registers a 3D image to a reference image using ITK Elastix.
+    '''
+    PixelType = itk.F
+    Dimension3D = 3
+    Image3DType = itk.Image[PixelType, Dimension3D]
+
+    # Built-in default parameter map seems to work well for brain
+    parameter_object = itk.ParameterObject.New()
+    parameter_object.AddParameterMap(parameter_object.GetDefaultParameterMap(method))
+
+    # Load reference image
+    tmp_ref = extract_volume(in_file, ref_idx)
+    ref_itk = itk.imread(tmp_ref, pixel_type=PixelType)
+
+    # Loop through all volumes and register to reference
+    n_vols = nib.load(in_file).shape[3]
+    registered_vols = []
+    for t in range(n_vols):
+        if t == ref_idx:
+            registered_vols.append(tmp_ref)  # No need to register reference to itself
+            continue
+
+        # Load moving image
+        tmp_moving = extract_volume(in_file, t)
+        moving_itk = itk.imread(tmp_moving, pixel_type=PixelType)
+
+        print(f'Registering volume {t} to reference volume {ref_idx}...')
+        elastix = itk.ElastixRegistrationMethod.New(ref_itk, moving_itk)
+        # elastix.SetFixedImage(ref_itk)
+        # elastix.SetMovingImage(moving_itk)
+        elastix.SetParameterObject(parameter_object)
+        elastix.Update()
+        result = elastix.GetOutput()
+        tmp_out = out_file.replace('.nii.gz',f'-tmp{t}.nii.gz')
+        itk.imwrite(result, tmp_out)
+        os.remove(tmp_moving)  # Clean up temp moving image
+        registered_vols.append(tmp_out)
+
+    out_file = in_file.replace('.nii.gz','-rigid.nii.gz')
+    # Stack back to 4D
+    subprocess.run(['fslmerge','-t',out_file]+registered_vols)
+
+    # Clean up
+    os.remove(tmp_ref)  # Clean up temp reference image
+    for t in range(1,n_vols):
+        os.remove(out_file.replace('-rigid.nii.gz',f'-reg-tmp{t}.nii.gz'))  # Clean up temp registered images
+    return out_file
+
+def topup_wrap(im_file,bval_file, b0rev_file, b0rev_bval_file):
+    ''' Susceptibility distortion correction '''
+    print('Running topup...')
+    im_file_unwarp = im_file.replace('.nii.gz','-unwarp.nii.gz')
+
+    # Prepare b0:s for topup
+    extract(im_file=im_file, bval_file=bval_file, outbase=im_file.replace('.nii.gz','-b0'),b_ex=0) # we dont perform averaging since it seems preferable to not do it
+    combine(dwi_files=[im_file.replace('.nii.gz','-b0.nii.gz'),b0rev_file],
+            bval_files=[im_file.replace('.nii.gz','-b0.bval'),b0rev_bval_file], 
+            outbase=im_file.replace('.nii.gz','-b0-b0rev'))
+    b = read_bval(im_file.replace('.nii.gz','-b0.bval'))
+    brev = read_bval(b0rev_bval_file)
+    interbase = im_file.replace('.nii.gz','-b0-b0rev')
+    acqp_file = interbase + '_acqparams.txt'
+    with open(acqp_file,'w') as f:
+        for _ in range(np.sum(b==0)):
+            f.write('0 1 0 0.050\n')
+        for _ in range(np.sum(brev==0)):
+            f.write('0 -1 0 0.050\n')
+
+    topup_cmd = f'topup --imain={im_file.replace(".nii.gz","-b0-b0rev.nii.gz")} --datain={acqp_file} --config=b02b0.cnf --subsamp=1 --out={interbase} --verbose'
+    os.system(topup_cmd)
+
+    applytopup_cmd = f'applytopup --imain={im_file} --datain={acqp_file} --inindex=1 --topup={interbase} --out={im_file_unwarp} --method=jac --verbose'
+    os.system(applytopup_cmd)   
+    subprocess.run(['cp',im_file.replace('.nii.gz','.bval'),im_file_unwarp.replace('.nii.gz','.bval')])
+
+    return im_file_unwarp
+
+def brain_preproc(im_file: str, bval_file: str, b0rev_file: str, b0rev_bval_file: str, temp_folder: str):
+
+    # Denoising
+    im_file_denoised = denoise_wrap(im_file)
+
+    # Degibbsing
+    im_file_degibbs = degibbs_wrap(im_file_denoised)
+
+    # Registration to b0
+    im_file_reg = register_3D_to_ref(im_file_degibbs, ref_idx=0, out_file=im_file_degibbs.replace('.nii.gz','-reg.nii.gz'))
+
+    # Suseptibility distortion correction 
+    im_file_unwarp = topup_wrap(im_file_reg, bval_file, b0rev_file, b0rev_bval_file)
+
+    # Brain extraction
+    brain_mask = im_file_unwarp.replace('.nii.gz','-brain-mask.nii.gz')
+    extract(im_file_unwarp, im_file_unwarp.replace('.nii.gz','.bval'), outbase=im_file_unwarp.replace('.nii.gz','-b0'), b_ex=0)
+    average(im_file_unwarp.replace('.nii.gz','-b0.nii.gz'), im_file_unwarp.replace('.nii.gz','-b0.bval'), im_file_unwarp.replace('.nii.gz','-b0-avr'))
+    subprocess.run(['mri_synthstrip', '-i', im_file_unwarp.replace('.nii.gz','-b0-avr.nii.gz'), '-o', im_file_unwarp.replace('.nii.gz', '-brain.nii.gz'), '-m', im_file_unwarp.replace('.nii.gz', '-brain-mask.nii.gz')])
+    os.remove(im_file_unwarp.replace('.nii.gz','-b0.nii.gz'))
+    os.remove(im_file_unwarp.replace('.nii.gz','-b0.bval'))
+    os.remove(im_file_unwarp.replace('.nii.gz','-b0-avr.nii.gz'))
+
+    # Signal drift correction
+    signal_drift.spatiotemporal(im_file = im_file_unwarp, 
+                                bval_file = im_file_unwarp.replace('.nii.gz','.bval'), 
+                                outbase = im_file_unwarp.replace('.nii.gz','-sigdri'), 
+                                roi_file = brain_mask)
+    subprocess.run(['cp',im_file_unwarp.replace('.nii.gz','.bval'),im_file_unwarp.replace('.nii.gz','-sigdri_corr.bval')])
+#%%
+# im_file = os.getcwd()+'/temp/diffusive_snr200.nii.gz'
+# bval_file = os.getcwd()+'/temp/diffusive_snr200.bval'
+# b0rev_file = os.getcwd()+'/temp/diffusive_snr200_b0.nii.gz'
+# b0rev_bval_file = os.getcwd()+'/temp/diffusive_snr200_b0.bval'
+
+# im_file_denoised = denoise_wrap(im_file)
+# im_file_degibbs = degibbs_wrap(im_file_denoised)
+# im_file_reg = register_3D_to_ref(im_file_degibbs, ref_idx=0, out_file=im_file_degibbs.replace('.nii.gz','-reg.nii.gz'))
+# im_file_unwarp = topup_wrap(im_file_reg, bval_file, b0rev_file, b0rev_bval_file)
+
+# Brain extraction
+# brain_mask = im_file_unwarp.replace('.nii.gz','-brain-mask.nii.gz')
+# extract(im_file_unwarp, bval_file, outbase=im_file_unwarp.replace('.nii.gz','-b0'), b_ex=0)
+# average(im_file_unwarp.replace('.nii.gz','-b0.nii.gz'), im_file_unwarp.replace('.nii.gz','-b0.bval'), im_file_unwarp.replace('.nii.gz','-b0-avr'))
+# subprocess.run(['mri_synthstrip', '-i', im_file_unwarp.replace('.nii.gz','-b0-avr.nii.gz'), '-o', im_file_unwarp.replace('.nii.gz', '-brain.nii.gz'), '-m', im_file_unwarp.replace('.nii.gz', '-brain-mask.nii.gz')])
+# os.remove(im_file_unwarp.replace('.nii.gz','-b0.nii.gz'))
+# os.remove(im_file_unwarp.replace('.nii.gz','-b0.bval'))
+# os.remove(im_file_unwarp.replace('.nii.gz','-b0-avr.nii.gz'))
+
+# # Signal drift correction
+# signal_drift.spatiotemporal(im_file = im_file_unwarp, 
+#                             bval_file = bval_file, 
+#                             outbase = im_file_unwarp.replace('.nii.gz','-sigdri'), 
+#                             roi_file = brain_mask)
+# subprocess.run(['cp',im_file_unwarp.replace('.nii.gz','.bval'),im_file_unwarp.replace('.nii.gz','-sigdri_corr.bval')])
+
+# %%

src/original/preprocessing/EP_GU/info.json

@@ -0,0 +1,4 @@
+{
+    "developers": "Oscar Jalnefjord",
+    "institutions": ["University of Gothenburg, Gothenburg, Sweden", "Sahlgrenska University Hospital, Gothenburg, Sweden"]
+}