ActiveAx analytic formulas

amico/analytic_formulas.py

@@ -0,0 +1,105 @@
+import scipy.special
+import numpy as np
+
+_GAMMA_ = 2.6751525E8
+_am_ = scipy.special.jnp_zeros(1, 60)
+
+def compute_PGSE_Sum(lambda_, beta, delta, DELTA, diff):
+    """Computes PGSE waveform term in eq 8 from [1].
+
+    References
+    ----------
+    .. [1] Ianus et al. (2013) Gaussian phase distribution approximations for oscillating gradient spin echo diffusion MRI. JMR, 227: 25-34
+    """
+    that_sum = 0
+    for n in range(len(lambda_)):
+        term1 = beta[n]/diff
+        term2 = lambda_[n]**2 * diff
+        that_sum += term1/term2 * (lambda_[n]*diff*delta - 1 + np.exp(-lambda_[n]*diff*delta) + np.exp(-lambda_[n]*diff*DELTA) - 0.5*np.exp(-lambda_[n]*diff*(DELTA-delta)) - 0.5*np.exp(-lambda_[n]*diff*(DELTA+delta)))
+    return that_sum
+
+def CylinderGPD(diff,theta,phi,R,scheme,gyro_ratio=_GAMMA_):
+    """Analytic computation of the Cylinder signal for N diffusion gradients of a PGSE protocol [1].
+
+    Attributes
+    ----------
+    diff : float
+        Diffusivity (m^2/s)
+    theta : float
+        Polar angle of the vector defining the fiber orientation (radians)
+    phi : float
+        Azimuth angle of the vector defining the fiber orientation (radians)
+    R : float
+        Cylinder radius (meters)
+    scheme : Scheme object
+        PGSE scheme object encoding diffusion gradient orientations, amplitudes, durations and separations
+    gyro_ratio : float (optional)
+        Gyromagnetic ratio.
+
+    References
+    ----------
+    .. [1] Ianus et al. (2013) Gaussian phase distribution approximations for oscillating gradient spin echo diffusion MRI. JMR, 227: 25-34
+    """
+    n_samples = scheme.nS
+    dir = scheme.raw[:,:3].copy()
+    modG = scheme.raw[:,3].copy()
+    DELTA = scheme.raw[:,4].copy()
+    delta = scheme.raw[:,5].copy()
+    G = dir*np.tile(modG,(3,1)).T
+    modG = np.linalg.norm(G,axis=1)
+    t = DELTA-delta/3.
+    unitGn = np.zeros(n_samples)
+    idx = (modG > 0.0)
+    cyl_vector = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)])[:,None]
+    unitGn[idx] = np.dot(G[idx,:],cyl_vector)[:,0] / (modG[idx] * np.linalg.norm(cyl_vector))
+    alpha = np.arccos(unitGn)
+    lambda_ = (_am_/R)**2
+    beta_factor = 2*(R/_am_)**2 / (_am_**2-1)
+    this_sum = np.zeros(n_samples)
+    for i in range(n_samples):
+        if this_sum[i] == 0:
+            idx = (delta == delta[i]) * (DELTA == DELTA[i])
+            this_sum[idx] = compute_PGSE_Sum(lambda_,beta_factor,delta[i],DELTA[i],diff)
+    Sperp = np.exp(-2 * gyro_ratio**2 * modG**2 * np.sin(alpha)**2 * this_sum)
+    Spar = np.exp(-t  * gyro_ratio**2 * delta**2 * modG**2 * np.cos(alpha)**2 * diff)
+    return Sperp * Spar
+
+def Zeppelin(diffPar,theta,phi,diffPerp,scheme):
+    """Analytic computation of the Zeppelin signal for a given scheme.
+    diffPar : float
+        Parallel diffusivity (m^2/s)
+    theta : float
+        Polar angle of the vector defining the Zeppelin orientation (radians)
+    phi : float
+        Azimuth angle of the vector defining the Zeppelin orientation (radians)
+    diffPerp : float
+        Perpendicular diffusivity (m^2/s)
+    scheme : Scheme object
+        PGSE scheme encoding diffusion gradient orientations and b-values
+    """
+    sinT = np.sin(theta)
+    cosT = np.cos(theta)
+    sinP = np.sin(phi)
+    cosP = np.cos(phi)
+    # Zeppelin orientation
+    n = np.array([cosP * sinT,sinP * sinT,cosT])
+    zepSignal = np.zeros(scheme.nS)
+    for i in range(scheme.nS):
+        dir_norm = np.linalg.norm(scheme.raw[i,:3])
+        if dir_norm > 0:
+            gd = scheme.raw[i,:3]/dir_norm
+        else:
+            gd = scheme.raw[i,:3]
+        #dot product of the cylinder orientation n and the wavenumber q
+        dotqn= np.dot(n,gd)
+        zepSignal[i] = np.exp(-scheme.b[i]*1E6*((diffPar-diffPerp)*dotqn**2 + diffPerp))
+    return zepSignal
+
+def Ball(diff,scheme):
+    """Analytic computation of the Ball signal for a given scheme.
+    diff : float
+        Free diffusivity (m^2/s)
+    scheme : Scheme object
+        PGSE scheme encoding diffusion gradient orientations and b-values
+    """
+    return np.exp(-scheme.b*diff*1E6)

amico/models.py

@@ -3,10 +3,10 @@
 import scipy
 from os.path import exists, join as pjoin
 from os import remove
-import subprocess
 import tempfile
 import amico.lut
 from amico.progressbar import ProgressBar
+from amico.analytic_formulas import CylinderGPD
 from dipy.core.gradients import gradient_table
 from dipy.sims.voxel import single_tensor
 import abc
@@ -329,51 +329,29 @@ def generate( self, out_path, aux, idx_in, idx_out ) :
         scheme_high = amico.lut.create_high_resolution_scheme( self.scheme, b_scale=1E6 )
         filename_scheme = pjoin( out_path, 'scheme.txt' )
         np.savetxt( filename_scheme, scheme_high.raw, fmt='%15.8e', delimiter=' ', header='VERSION: STEJSKALTANNER', comments='' )
-
-        # temporary file where to store "datasynth" output
-        filename_signal = pjoin( tempfile._get_default_tempdir(), next(tempfile._get_candidate_names())+'.Bfloat' )
+
+        gtab = gradient_table( scheme_high.b/1E6, scheme_high.raw[:,0:3] )
 
         nATOMS = len(self.Rs) + len(self.ICVFs) + len(self.d_ISOs)
         progress = ProgressBar( n=nATOMS, prefix="   ", erase=True )
 
         # Cylinder(s)
         for R in self.Rs :
-            CMD = 'datasynth -synthmodel compartment 1 CYLINDERGPD %E 0 0 %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null' % ( self.d_par*1E-6, R, filename_scheme, filename_signal )
-            subprocess.call( CMD, shell=True )
-            if not exists( filename_signal ) :
-                raise RuntimeError( 'Problems generating the signal with "datasynth"' )
-            signal  = np.fromfile( filename_signal, dtype='>f4' )
-            if exists( filename_signal ) :
-                remove( filename_signal )
-
+            signal  = CylinderGPD(self.d_par*1E-6,0.0,0.0,R,scheme_high)
             lm = amico.lut.rotate_kernel( signal, aux, idx_in, idx_out, False )
             np.save( pjoin( out_path, 'A_%03d.npy'%progress.i ), lm )
             progress.update()
 
         # Zeppelin(s)
         for d in [ self.d_par*(1.0-ICVF) for ICVF in self.ICVFs] :
-            CMD = 'datasynth -synthmodel compartment 1 ZEPPELIN %E 0 0 %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null' % ( self.d_par*1E-6, d*1e-6, filename_scheme, filename_signal )
-            subprocess.call( CMD, shell=True )
-            if not exists( filename_signal ) :
-                raise RuntimeError( 'Problems generating the signal with "datasynth"' )
-            signal  = np.fromfile( filename_signal, dtype='>f4' )
-            if exists( filename_signal ) :
-                remove( filename_signal )
-
+            signal = single_tensor( gtab, evals=[d, d, self.d_par] )
             lm = amico.lut.rotate_kernel( signal, aux, idx_in, idx_out, False )
             np.save( pjoin( out_path, 'A_%03d.npy'%progress.i ), lm )
             progress.update()
-
+        
         # Ball(s)
         for d in self.d_ISOs :
-            CMD = 'datasynth -synthmodel compartment 1 BALL %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null' % ( d*1e-6, filename_scheme, filename_signal )
-            subprocess.call( CMD, shell=True )
-            if not exists( filename_signal ) :
-                raise RuntimeError( 'Problems generating the signal with "datasynth"' )
-            signal  = np.fromfile( filename_signal, dtype='>f4' )
-            if exists( filename_signal ) :
-                remove( filename_signal )
-
+            signal = single_tensor( gtab, evals=[d, d, d] )
             lm = amico.lut.rotate_kernel( signal, aux, idx_in, idx_out, True )
             np.save( pjoin( out_path, 'A_%03d.npy'%progress.i ), lm )
             progress.update()