Merge pull request #19 from OSIPI/brain_phantom

Brain phantom

Author: oliverchampion

phantoms/brain/README.md

@@ -0,0 +1,7 @@
+# Code for generating ivim phantoms of the brain
+
+Two sets of ground truth data are used for simulation:
+- Simulation in the diffusive regime (IVIM parameters D, f and D*): reference values from Ryghög et al. 2014 and b-values from Federau et al. 2012
+- Simulation in the ballistic regime (IVIM parameters D, f and vd): reference values and sequence parameters from Ahlgren et al. 2016
+
+The segmentation used by the simulations is the ICBM 2009a nonlinear symmetric 3T atlas (https://nist.mni.mcgill.ca/icbm-152-nonlinear-atlases-2009/), the same as in e.g. ASLDRO (https://asldro.readthedocs.io/en/stable/).

phantoms/brain/ground_truth/ballistic.bval

@@ -0,0 +1 @@
+0 10 20 30 40 50 60 70 80 90 100 120 140 160 180 200 10 20 30 40 50 60 70 80 90 100 120 140 160 180 200

phantoms/brain/ground_truth/ballistic.cval

@@ -0,0 +1 @@
+0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.471 0.666 0.816 0.942 1.054 1.154 1.247 1.333 1.414 1.490 1.632 1.763 1.885 1.999 2.107 

phantoms/brain/ground_truth/ballistic_calc_cval.py

@@ -0,0 +1,15 @@
+import os
+import numpy as np
+
+# Data from Ahlgren et al 2016 NMRinBiomed
+Delta = 7.5e-3 # 7.5 ms
+delta = 7.3e-3 # 7.3 ms
+
+# For DDE sequence we have
+# b = 2y^2G^2d^2(D-d/3), c = 2yGdD => c = sqrt(b 2D^2/(D-d/3))
+bval_file = os.path.join(os.path.dirname(__file__),'ballistic.bval')
+b = np.loadtxt(bval_file)
+c = np.sqrt(b*2*Delta**2/(Delta-delta/3))
+c[1:(c.size-1)//2+1] = 0 # flow compensated => c = 0
+cval_file = bval_file.replace('bval','cval')
+np.savetxt(cval_file,c,fmt='%.3f',newline=' ')

phantoms/brain/ground_truth/ballistic_groundtruth.json

@@ -0,0 +1,6 @@
+{
+    "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
+}

phantoms/brain/ground_truth/diffusive.bval

@@ -0,0 +1 @@
+0 10 20 30 40 80 110 140 170 200 300 400 500 600 700 800 900

phantoms/brain/ground_truth/diffusive_groundtruth.json

@@ -0,0 +1,5 @@
+{
+    "D":[0.81e-3,0.86e-3,3e-3],
+    "f":[0.044,0.033,0],
+    "Dstar":[84e-3,76e-3,0]
+}

phantoms/brain/sim_brain_phantom.py

@@ -0,0 +1,61 @@
+import os
+import shutil
+import json
+import numpy as np
+import nibabel as nib
+from scipy.ndimage import zoom
+
+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(folder,'ground_truth','hrgt_icbm_2009a_nls_3t.nii.gz'))
+segmentation = np.squeeze(nii.get_fdata()[...,-1])
+
+with open(os.path.join(folder,'ground_truth',regime+'_groundtruth.json'), 'r') as f:
+       ivim_pars = json.load(f)
+S0 = 1
+
+# Sequence parameters
+bval_file = os.path.join(folder,'ground_truth',regime+'.bval')
+b = np.loadtxt(bval_file)
+if regime == BALLISTIC_REGIME:
+    cval_file = bval_file.replace('bval','cval')
+    c = np.loadtxt(cval_file)
+
+# Calculate signal
+S = np.zeros(list(np.shape(segmentation))+[b.size])
+
+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"])):
+        S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Db-c**2*vd**2))
+else:
+    for i,(D,f,Dstar) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"])):
+        S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Dstar))
+
+# 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
+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 
+nii_out = nib.Nifti1Image(im_noise,np.eye(4))
+base_name = os.path.join(folder,'data','{}_snr{}'.format(regime,snr))
+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')