Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on Sep 7th, 2012  |  syntax: None  |  size: 3.60 KB  |  hits: 13  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. import numpy as np
  2.  
  3. import tensor as ozt
  4.  
  5.  
  6. # Global constants for this module:
  7. AD = 1.5
  8. RD = 0.5
  9.  
  10. # This converts b values from the data, so that it matches the units of ADC we
  11. # use in the Stejskal/Tanner equation:
  12. SCALE_FACTOR = 1000
  13.  
  14.  
  15.  
  16. class ODF(object):
  17.     """
  18.     This class represents a distribution of fiber orientations within a voxel
  19.    
  20.     """
  21.     def __init__(self,
  22.                  bvecs,
  23.                  weights):
  24.  
  25.         self.bvecs = bvecs
  26.  
  27.         # Enforce that the weights are a pdf:
  28.         if np.sum(weights)!=1:
  29.             e_s = "The ODF is a distribution function, so"
  30.             e_s += "the weights must sum to 1"
  31.             raise ValueError(e_s)
  32.            
  33.         self.weights = weights
  34.        
  35.        
  36. class Voxel(object):
  37.     """
  38.     The representation of the measurement in a single voxel with some
  39.     configuration of fibers crossing through it, represented by an ODF
  40.     """
  41.     def __init__(self,
  42.                  bvecs,
  43.                  bvals,
  44.                  odf,
  45.                  scaling_factor=SCALE_FACTOR,
  46.                  axial_diffusivity=AD,
  47.                  radial_diffusivity=RD):
  48.         """
  49.         The signal in a single voxel is computed from the convolution of a
  50.         response function with the ODF
  51.        
  52.         Parameters
  53.         ----------
  54.  
  55.         bvecs: 3 by n array
  56.              unit vectors on the sphere.
  57.        
  58.         bvals: 1 by n array
  59.              The measurement parameter defining where on the curve we measure
  60.              the exponential decay of the signal
  61.  
  62.         odf: an ODF class instance
  63.             The representation of the orientation distribution function. Note
  64.             that this class also has bvecs, but these bvecs don't have to be
  65.             the same bvecs as the ones of this class. bvecs In this class refer
  66.             to measurement directions, while the bvecs in the ODF class refer
  67.             to directions of fibers within the voxel (not necessarily the
  68.             same...)
  69.  
  70.         S0: float (optional)
  71.             The signal measured in the non diffusion-weighted scans. Default: 1
  72.  
  73.         scaling_factor: float (optional)
  74.             To get the right units on the ADC, sometimes the b value needs to
  75.             be scaled. Typically, divided by 1000 (default).
  76.  
  77.         axial_diffusivity, radial_diffusivity: float
  78.             These parameters
  79.            
  80.         """
  81.  
  82.         self.bvecs = bvecs
  83.         self.bvals = bvals/scaling_factor
  84.         self.odf = odf
  85.        
  86.         # We assume that the response function is a cigar shaped tensor:
  87.         self.Q = np.array([[axial_diffusivity, 0, 0],
  88.                            [0, radial_diffusivity, 0],
  89.                            [0, 0, radial_diffusivity]])
  90.         self.response_function = ozt.Tensor(self.Q, self.bvecs, self.bvals)
  91.  
  92.     def signal(self,S0=1):
  93.         """
  94.         The signal generated by the ODF in each bvec direction
  95.         """
  96.         # Initialize to zero:
  97.         signal = np.zeros(self.bvecs.shape[-1])
  98.  
  99.         # Decompose the response function to its eigen-decomposition:
  100.         evals, evecs = self.response_function.decompose
  101.         # We generate the signal for each one of the fibers in the ODF:
  102.         for odf_idx, odf_vec in enumerate(self.odf.bvecs):
  103.             this_response = ozt.rotate_to_vector(odf_vec,
  104.                                                  evals,
  105.                                                  evecs,
  106.                                                  self.bvecs,
  107.                                                  self.bvals)
  108.  
  109.             signal += (self.odf.weights[odf_idx] *
  110.                        this_response.predicted_signal(S0))
  111.  
  112.         return signal
create a new version of this paste RAW Paste Data