Merge pull request #161 from OSIPI/preproc-main

Preproc-main into main

Author: oliverchampion

.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,

…original/DK_OGC_AmsterdamUMC/__init__.py phantoms/__init__.pysrc/original/DK_OGC_AmsterdamUMC/__init__.py renamed to phantoms/__init__.py src/original/DK_OGC_AmsterdamUMC/__init__.py renamed to phantoms/__init__.py

@@ -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/ballistic_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/ground_truth/diffusive_groundtruth.json

@@ -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()

phantoms/brain/sim_brain_phantom.py

@@ -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

pytest.ini

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