merge to just one sim_brain_phantom script

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

tests/IVIMpreproc/unit_tests/wip_denoise.py

@@ -7,7 +7,7 @@
 import os
 
 import numpy as np
-from phantoms.brain.sim_brain_phantom_preproc import simulate_brain_phantom
+from phantoms.brain.sim_brain_phantom import simulate_brain_phantom
 import nibabel as nib
 import matplotlib.pyplot as plt
 from src.original.preprocessing.EP_GU.brain_pipeline import denoise_wrap
@@ -19,13 +19,13 @@ def test_denoise_wrap():
     simulate_brain_phantom(snr = snr)
 
     # Denoising
-    im_file = 'phantoms/brain/data/diffusive_snr{}_relax.nii.gz'.format(snr)
+    im_file = 'phantoms/brain/data/diffusive_snr{}.nii.gz'.format(snr)
     im_file_denoised = denoise_wrap(im_file)
     S_denoised = nib.load(im_file_denoised).get_fdata()
 
     # Compute sum of squared residuals before denoising
-    mask = nib.load(os.path.join('phantoms','brain','data','diffusive_snr{}_relax_mask.nii.gz'.format(snr))).get_fdata()
-    ref_file = 'phantoms/brain/data/diffusive_reference_relax.nii.gz'
+    mask = nib.load(os.path.join('phantoms','brain','data','diffusive_snr{}_mask.nii.gz'.format(snr))).get_fdata()
+    ref_file = 'phantoms/brain/data/diffusive_reference.nii.gz'
     S = nib.load(im_file).get_fdata()
     S_ref = nib.load(ref_file).get_fdata()
     ssr_before = np.sum((S[mask!=0,:]-S_ref[mask!=0,:])**2)