Simulation

simulation.py

import numpy as np

import tensor as ozt


# Global constants for this module:
AD = 1.5
RD = 0.5

# This converts b values from the data, so that it matches the units of ADC we
# use in the Stejskal/Tanner equation:
SCALE_FACTOR = 1000



class ODF(object):
    """
    This class represents a distribution of fiber orientations within a voxel
    
    """
    def __init__(self,
                 bvecs,
                 weights):

        self.bvecs = bvecs

        # Enforce that the weights are a pdf: 
        if np.sum(weights)!=1:
            e_s = "The ODF is a distribution function, so"
            e_s += "the weights must sum to 1"
            raise ValueError(e_s)
            
        self.weights = weights
        
        
class Voxel(object):
    """
    The representation of the measurement in a single voxel with some
    configuration of fibers crossing through it, represented by an ODF
    """
    def __init__(self,
                 bvecs,
                 bvals,
                 odf,
                 scaling_factor=SCALE_FACTOR,
                 axial_diffusivity=AD,
                 radial_diffusivity=RD):
        """
        The signal in a single voxel is computed from the convolution of a
        response function with the ODF
        
        Parameters
        ----------

        bvecs: 3 by n array
             unit vectors on the sphere.
        
        bvals: 1 by n array
             The measurement parameter defining where on the curve we measure
             the exponential decay of the signal

        odf: an ODF class instance
            The representation of the orientation distribution function. Note
            that this class also has bvecs, but these bvecs don't have to be
            the same bvecs as the ones of this class. bvecs In this class refer
            to measurement directions, while the bvecs in the ODF class refer
            to directions of fibers within the voxel (not necessarily the
            same...) 

        S0: float (optional)
            The signal measured in the non diffusion-weighted scans. Default: 1

        scaling_factor: float (optional)
            To get the right units on the ADC, sometimes the b value needs to
            be scaled. Typically, divided by 1000 (default).

        axial_diffusivity, radial_diffusivity: float
            These parameters
            
        """

        self.bvecs = bvecs
        self.bvals = bvals/scaling_factor
        self.odf = odf
        
        # We assume that the response function is a cigar shaped tensor:
        self.Q = np.array([[axial_diffusivity, 0, 0],
                           [0, radial_diffusivity, 0],
                           [0, 0, radial_diffusivity]])
        self.response_function = ozt.Tensor(self.Q, self.bvecs, self.bvals)

    def signal(self,S0=1):
        """
        The signal generated by the ODF in each bvec direction
        """
        # Initialize to zero:
        signal = np.zeros(self.bvecs.shape[-1])

        # Decompose the response function to its eigen-decomposition:
        evals, evecs = self.response_function.decompose
        # We generate the signal for each one of the fibers in the ODF: 
        for odf_idx, odf_vec in enumerate(self.odf.bvecs):
            this_response = ozt.rotate_to_vector(odf_vec,
                                                 evals,
                                                 evecs,
                                                 self.bvecs,
                                                 self.bvals)

            signal += (self.odf.weights[odf_idx] *
                       this_response.predicted_signal(S0))

        return signal