From 6786921973a013011cbed2816a9655146b88f052 Mon Sep 17 00:00:00 2001 From: David Romascano Date: Wed, 2 Aug 2017 21:20:43 +0200 Subject: [PATCH 1/4] NF: added analytical formula to compute Cylinder signal in LUT --- amico/analytic_formulas.py | 104 +++++++++++++++++++++++++++++++++++++ amico/models.py | 36 +++---------- 2 files changed, 111 insertions(+), 29 deletions(-) create mode 100644 amico/analytic_formulas.py diff --git a/amico/analytic_formulas.py b/amico/analytic_formulas.py new file mode 100644 index 0000000..f65e12b --- /dev/null +++ b/amico/analytic_formulas.py @@ -0,0 +1,104 @@ +import scipy.special +import numpy as np + +_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): + """Analytic computation of the Cylinder signal for N diffusion gradients [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 + + 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() + for k in range(scheme.nS): + dir[k,:] /= np.linalg.norm(dir[k,:]) + G = dir*np.tile(scheme.raw[:,3],(3,1)).T + S = np.zeros(n_samples) + GAMMA = 2.6751525E8 + t = scheme.raw[:,4]-scheme.raw[:,5]/3. + unitGn = np.zeros(n_samples) + idx = (scheme.raw[:,3] > 0) + cyl_vector = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)]) + unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T) / (scheme.raw[:,3][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 = (scheme.raw[:,5] == scheme.raw[i,5]) * (scheme.raw[:,4] == scheme.raw[i,4]) + this_sum[idx] = compute_PGSE_Sum(lambda_,beta_factor,scheme.raw[i,5],scheme.raw[i,4],diff) + Sperp = np.exp(-2 * GAMMA**2 * scheme.raw[:,3]**2 * np.sin(alpha)**2 * this_sum) + # the restricted parallel signal + Spar = np.exp(-t * GAMMA**2 * scheme.raw[i,5]**2 * scheme.raw[:,3]**2 * np.cos(alpha)**2 * diff) + # the restricted signal + 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*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) \ No newline at end of file diff --git a/amico/models.py b/amico/models.py index 4c4056e..371eb84 100644 --- a/amico/models.py +++ b/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() From 40074e2fc90d66bb0f5d5c47bb40b7f2dc8c19c4 Mon Sep 17 00:00:00 2001 From: David Romascano Date: Wed, 2 Aug 2017 23:25:05 +0200 Subject: [PATCH 2/4] Fixed Cylinders so signal matches output of Camino; Zeppelin output from dipy doesn't match Camino though... --- amico/analytic_formulas.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/amico/analytic_formulas.py b/amico/analytic_formulas.py index f65e12b..103a6f8 100644 --- a/amico/analytic_formulas.py +++ b/amico/analytic_formulas.py @@ -39,28 +39,28 @@ def CylinderGPD(diff,theta,phi,R,scheme): """ n_samples = scheme.nS dir = scheme.raw[:,:3].copy() - for k in range(scheme.nS): + for k in range(n_samples): dir[k,:] /= np.linalg.norm(dir[k,:]) - G = dir*np.tile(scheme.raw[:,3],(3,1)).T - S = np.zeros(n_samples) + modG = scheme.raw[:,3].copy() + DELTA = scheme.raw[:,4].copy() + delta = scheme.raw[:,5].copy() + G = dir*np.tile(modG,(3,1)).T GAMMA = 2.6751525E8 - t = scheme.raw[:,4]-scheme.raw[:,5]/3. + t = DELTA-delta/3. unitGn = np.zeros(n_samples) - idx = (scheme.raw[:,3] > 0) + idx = (modG > 0) cyl_vector = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)]) - unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T) / (scheme.raw[:,3][idx] * np.linalg.norm(cyl_vector)) + unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T) / (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 = (scheme.raw[:,5] == scheme.raw[i,5]) * (scheme.raw[:,4] == scheme.raw[i,4]) - this_sum[idx] = compute_PGSE_Sum(lambda_,beta_factor,scheme.raw[i,5],scheme.raw[i,4],diff) - Sperp = np.exp(-2 * GAMMA**2 * scheme.raw[:,3]**2 * np.sin(alpha)**2 * this_sum) - # the restricted parallel signal - Spar = np.exp(-t * GAMMA**2 * scheme.raw[i,5]**2 * scheme.raw[:,3]**2 * np.cos(alpha)**2 * diff) - # the restricted signal + 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 * GAMMA**2 * modG**2 * np.sin(alpha)**2 * this_sum) + Spar = np.exp(-t * GAMMA**2 * delta**2 * modG**2 * np.cos(alpha)**2 * diff) return Sperp * Spar def Zeppelin(diffPar,theta,phi,diffPerp,scheme): From 80479cb340f54803b75bde8eb8a33525178512cf Mon Sep 17 00:00:00 2001 From: David Romascano Date: Thu, 31 May 2018 10:10:11 +0200 Subject: [PATCH 3/4] Fixed bug for new version of Numpy that complained about division of (1,N) by (N,) vectors --- amico/analytic_formulas.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/amico/analytic_formulas.py b/amico/analytic_formulas.py index 103a6f8..8c77815 100644 --- a/amico/analytic_formulas.py +++ b/amico/analytic_formulas.py @@ -50,7 +50,7 @@ def CylinderGPD(diff,theta,phi,R,scheme): unitGn = np.zeros(n_samples) idx = (modG > 0) cyl_vector = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)]) - unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T) / (modG[idx] * np.linalg.norm(cyl_vector)) + unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T)[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) @@ -91,7 +91,7 @@ def Zeppelin(diffPar,theta,phi,diffPerp,scheme): 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*1E6*((diffPar-diffPerp)*dotqn**2 + diffPerp)) + zepSignal[i] = np.exp(-scheme.b[i]*1E6*((diffPar-diffPerp)*dotqn**2 + diffPerp)) return zepSignal def Ball(diff,scheme): @@ -101,4 +101,4 @@ def Ball(diff,scheme): scheme : Scheme object PGSE scheme encoding diffusion gradient orientations and b-values """ - return np.exp(-scheme.b*diff*1E6) \ No newline at end of file + return np.exp(-scheme.b*diff*1E6) From ed1927784b4d0d0f01c79cf9d2d443f9ed08906b Mon Sep 17 00:00:00 2001 From: David Romascano Date: Wed, 11 Jul 2018 10:45:44 +0200 Subject: [PATCH 4/4] Modified CylinderGPD to match output of Camino for R = [0.01E-6,0.1E-6,0.5E-6,1.0E-6,5.0E-6,10.0E-6] --- amico/analytic_formulas.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/amico/analytic_formulas.py b/amico/analytic_formulas.py index 8c77815..a22e16d 100644 --- a/amico/analytic_formulas.py +++ b/amico/analytic_formulas.py @@ -1,11 +1,12 @@ 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 @@ -17,9 +18,9 @@ def compute_PGSE_Sum(lambda_, beta, delta, DELTA, 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): - """Analytic computation of the Cylinder signal for N diffusion gradients [1]. - +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 @@ -32,25 +33,25 @@ def CylinderGPD(diff,theta,phi,R,scheme): 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() - for k in range(n_samples): - dir[k,:] /= np.linalg.norm(dir[k,:]) modG = scheme.raw[:,3].copy() DELTA = scheme.raw[:,4].copy() delta = scheme.raw[:,5].copy() G = dir*np.tile(modG,(3,1)).T - GAMMA = 2.6751525E8 + modG = np.linalg.norm(G,axis=1) t = DELTA-delta/3. unitGn = np.zeros(n_samples) - idx = (modG > 0) - cyl_vector = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)]) - unitGn[idx] = np.dot(cyl_vector[None,:],G[idx].T)[0,:] / (modG[idx] * np.linalg.norm(cyl_vector)) + 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) @@ -59,8 +60,8 @@ def CylinderGPD(diff,theta,phi,R,scheme): 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 * GAMMA**2 * modG**2 * np.sin(alpha)**2 * this_sum) - Spar = np.exp(-t * GAMMA**2 * delta**2 * modG**2 * np.cos(alpha)**2 * 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):