Develop

asltk/asldata.py

@@ -82,13 +82,21 @@ def __init__(
         if kwargs.get('dw_values'):
             self._parameters['dw'] = kwargs.get('dw_values')
 
-    def set_image(self, full_path: str, spec: str):
+    def set_image(self, image, spec: str):
         """Insert a image necessary to define de ASL data processing.
 
         The `spec` parameters specifies what is the type of image to be used in
         ASL processing step. Choose one of the options: `m0` for the M0 volume,
         `pcasl` for the pCASL data.
 
+        Note:
+            The image can be a full path to the image file or a numpy array.
+            In case the image parameter is a path, then the method will load
+            the image file directly and associate it with the `spec` parameter.
+            However, if the image is a numpy array, then the method will
+            pass it to the ASLData object image data regarding the `spec`
+            parameter as well.
+
         Examples:
             >>> data = ASLData()
             >>> path_m0 = './tests/files/m0.nii.gz' # M0 file with shape (5,35,35)
@@ -97,13 +105,19 @@ def set_image(self, full_path: str, spec: str):
             (5, 35, 35)
 
         Args:
-            full_path (str): The full path with filename to be loaded
+            image (str): The image to be used.
             spec (str): The type of image being used in the ASL processing.
         """
-        if spec == 'm0':
-            self._m0_image = load_image(full_path)
-        elif spec == 'pcasl':
-            self._asl_image = load_image(full_path)
+        if isinstance(image, str) and os.path.exists(image):
+            if spec == 'm0':
+                self._m0_image = load_image(image)
+            elif spec == 'pcasl':
+                self._asl_image = load_image(image)
+        elif isinstance(image, np.ndarray):
+            if spec == 'm0':
+                self._m0_image = image
+            elif spec == 'pcasl':
+                self._asl_image = image
 
     def get_ld(self):
         """Obtain the LD array values"""

asltk/models/signal_dynamic.py

@@ -0,0 +1,336 @@
+import math
+
+import numpy as np
+
+
+def asl_model_buxton(
+    tau: list,
+    w: list,
+    m0: float,
+    cbf: float,
+    att: float,
+    lambda_value: float = 0.98,
+    t1b: float = 1650.0,
+    alpha: float = 0.85,
+):
+    """Buxton model to calculate the ASL magnetization values.
+
+    It is assumed that the LD and PLD values are coherent with the ASl Buxton
+    model, i.e. the both has the same array size.
+
+    The calculations is given assuming a voxel value. Hence, all the `tau`,
+    `w`, `cbf` and `att` values must representas a voxel in the image.
+
+    Note:
+        The CBF value is the original scale, without assuming the normalized
+        CBF value. See more details at the CBFMapping class documentation.
+
+    Args:
+        tau (list): LD values
+        w (list): PLD values
+        m0 (float): The M0 magnetization value
+        cbf (float): The CBF value, not been assumed as normalized.
+        att (float): The ATT value
+        lambda_value (float, optional): The blood-brain partition coefficient (0 to 1.0). Defaults to 0.98.
+        t1b (float, optional): The T1 relaxation value of the blood. Defaults to 1650.0.
+        alpha (float, optional): The labeling efficiency. Defaults to 0.85.
+
+    Returns:
+        (numpy.ndarray): A numpy array with the magnetization values calculated
+    """
+    tau = tau.tolist() if isinstance(tau, np.ndarray) else tau
+    w = w.tolist() if isinstance(w, np.ndarray) else w
+
+    if not (isinstance(tau, list) ^ isinstance(tau, tuple)):
+        raise ValueError('tau parameter must be a list or tuple of values.')
+
+    if not isinstance(w, list) ^ isinstance(w, tuple):
+        raise ValueError('w parameter must be a list or tuple of values.')
+
+    for v in tau:
+        if not isinstance(v, float) ^ isinstance(v, int):
+            raise ValueError('tau list must contain float or int values')
+
+    for v in w:
+        if not isinstance(v, float) ^ isinstance(v, int):
+            raise ValueError('w list must contain float or int values')
+
+    # if len(tau) != len(w):
+    #     raise SyntaxError("tau and w parameters must be at the same size.")
+
+    t = np.add(tau, w).tolist()
+
+    t1bp = 1 / ((1 / t1b) + (cbf / lambda_value))
+    m_values = np.zeros(len(tau))
+
+    for i in range(0, len(tau)):
+        try:
+            if t[i] < att:
+                m_values[i] = 0.0
+            elif (att <= t[i]) and (t[i] < tau[i] + att):
+                q = 1 - math.exp(-(t[i] - att) / t1bp)
+                m_values[i] = (
+                    2.0 * m0 * cbf * t1bp * alpha * q * math.exp(-att / t1b)
+                )
+            else:
+                q = 1 - math.exp(-tau[i] / t1bp)
+                m_values[i] = (
+                    2.0
+                    * m0
+                    * cbf
+                    * t1bp
+                    * alpha
+                    * q
+                    * math.exp(-att / t1b)
+                    * math.exp(-(t[i] - tau[i] - att) / t1bp)
+                )
+        except OverflowError:   # pragma: no cover
+            m_values[i] = 0.0
+
+    return m_values
+
+
+def asl_model_multi_te(
+    tau: list,
+    w: list,
+    te: list,
+    m0: float,
+    cbf: float,
+    att: float,
+    t2b: float = 165.0,
+    t2csf: float = 75.0,
+    tblcsf: float = 1400.0,
+    alpha: float = 0.85,
+    t1b: float = 1650.0,
+    t1csf: float = 1400.0,
+):
+    """Multi Time of Echos (TE) ASL model to calculate the T1 relaxation time for
+    blood and Grey Matter exchange.
+
+    This model is directly used on the MultiTE_ASLMapping class.
+
+    Reference: Ultra-long-TE arterial spin labeling reveals rapid and
+    brain-wide blood-to-CSF water transport in humans, NeuroImage,
+    doi: 10.1016/j.neuroimage.2021.118755
+
+    Args:
+        tau (list): The LD values
+        w (list): The PLD values
+        te (list): The TE values
+        m0 (float): The M0 voxel value
+        cbf (float): The CBF voxel value
+        att (float): The ATT voxel value
+        t2b (float, optional): The T2 relaxation value for blood. Defaults to 165.0.
+        t2csf (float, optional): The T2 relaxation value for CSF. Defaults to 75.0.
+        tblcsf (float, optional): The T1 relaxation value between blood and CSF. Defaults to 1400.0.
+        alpha (float, optional): The pulse labeling efficiency. Defaults to 0.85.
+        t1b (float, optional): The T1 relaxation value for blood. Defaults to 1650.0.
+        t1csf (float, optional): The T1 relaxation value for CSF. Defaults to 1400.0.
+
+    Returns:
+        (nd.ndarray): The magnetization values for T1-Blood-GM
+    """
+    t1bp = 1 / ((1 / t1b) + (1 / tblcsf))
+    t1csfp = 1 / ((1 / t1csf) + (1 / tblcsf))
+
+    t2bp = 1 / ((1 / t2b) + (1 / tblcsf))
+    t2csfp = 1 / ((1 / t2csf) + (1 / tblcsf))
+
+    t = np.add(tau, w).tolist()
+
+    mag_total = np.zeros(len(tau))
+
+    for i in range(0, len(tau)):
+        try:
+            if t[i] < att:
+                S1b = 0.0
+                S1csf = 0.0
+                if te[i] < (att - t[i]):
+                    Sb = 0
+                    Scsf = 0
+                elif (att - t[i]) <= te[i] and te[i] < (att + tau[i] - t[i]):
+                    Sb = (
+                        2
+                        * alpha
+                        * m0
+                        * cbf
+                        * t2bp
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * (1 - math.exp(-(te[i] - att + t[i]) / t2bp))
+                    )   #% measured signal = S2
+                    Scsf = (
+                        2
+                        * alpha
+                        * m0
+                        * cbf
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * (
+                            t2csf
+                            * (1 - math.exp(-(te[i] - att + t[i]) / t2csf))
+                            - t2csfp
+                            * (1 - math.exp(-(te[i] - att + t[i]) / t2csfp))
+                        )
+                    )
+                else:   #% att + tau - t <= te
+                    Sb = (
+                        2
+                        * alpha
+                        * m0
+                        * cbf
+                        * t2bp
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * math.exp(-(te[i] - att + t[i]) / t2bp)
+                        * (math.exp(tau[i] / t2bp) - 1)
+                    )
+                    Scsf = (
+                        2
+                        * alpha
+                        * m0
+                        * cbf
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * (
+                            t2csf
+                            * math.exp(-(te[i] - att + t[i]) / t2csf)
+                            * (math.exp(tau[i] / t2csf) - 1)
+                            - t2csfp
+                            * math.exp(-(te[i] - att + t[i]) / t2csfp)
+                            * (math.exp(tau[i] / t2csfp) - 1)
+                        )
+                    )
+            elif (att <= t[i]) and (t[i] < (att + tau[i])):
+                S1b = (
+                    2
+                    * alpha
+                    * m0
+                    * cbf
+                    * t1bp
+                    * math.exp(-att / t1b)
+                    * (1 - math.exp(-(t[i] - att) / t1bp))
+                )
+                S1csf = (
+                    2
+                    * alpha
+                    * m0
+                    * cbf
+                    * math.exp(-att / t1b)
+                    * (
+                        t1csf * (1 - math.exp(-(t[i] - att) / t1csf))
+                        - t1csfp * (1 - math.exp(-(t[i] - att) / t1csfp))
+                    )
+                )
+                if te[i] < (att + tau[i] - t[i]):
+                    Sb = S1b * math.exp(
+                        -te[i] / t2bp
+                    ) + 2 * alpha * m0 * cbf * t2bp * math.exp(
+                        -att / t1b
+                    ) * math.exp(
+                        -te[i] / t2b
+                    ) * (
+                        1 - math.exp(-te[i] / t2bp)
+                    )
+                    Scsf = (
+                        S1b
+                        * (1 - math.exp(-te[i] / tblcsf))
+                        * math.exp(-te[i] / t2csf)
+                        + S1csf * math.exp(-te[i] / t2csf)
+                        + 2
+                        * alpha
+                        * m0
+                        * cbf
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * (
+                            t2csf * (1 - math.exp(-te[i] / t2csf))
+                            - t2csfp * (1 - math.exp(-te[i] / t2csfp))
+                        )
+                    )
+                else:   # att + tau - t <= te
+                    Sb = S1b * math.exp(
+                        -te[i] / t2bp
+                    ) + 2 * alpha * m0 * cbf * t2bp * math.exp(
+                        -att / t1b
+                    ) * math.exp(
+                        -te[i] / t2b
+                    ) * math.exp(
+                        -te[i] / t2bp
+                    ) * (
+                        math.exp((att + tau[i] - t[i]) / t2bp) - 1
+                    )
+                    Scsf = (
+                        S1b
+                        * (1 - math.exp(-te[i] / tblcsf))
+                        * math.exp(-te[i] / t2csf)
+                        + S1csf * math.exp(-te[i] / t2csf)
+                        + 2
+                        * alpha
+                        * m0
+                        * cbf
+                        * math.exp(-att / t1b)
+                        * math.exp(-te[i] / t2b)
+                        * (
+                            t2csf
+                            * math.exp(-te[i] / t2csf)
+                            * (math.exp((att + tau[i] - t[i]) / t2csf) - 1)
+                            - t2csfp
+                            * math.exp(-te[i] / t2csfp)
+                            * (math.exp((att + tau[i] - t[i]) / t2csfp) - 1)
+                        )
+                    )
+            else:   # att+tau < t
+                S1b = (
+                    2
+                    * alpha
+                    * m0
+                    * cbf
+                    * t1bp
+                    * math.exp(-att / t1b)
+                    * math.exp(-(t[i] - att) / t1bp)
+                    * (math.exp(tau[i] / t1bp) - 1)
+                )
+                S1csf = (
+                    2
+                    * alpha
+                    * m0
+                    * cbf
+                    * math.exp(-att / t1b)
+                    * (
+                        t1csf
+                        * math.exp(-(t[i] - att) / t1csf)
+                        * (math.exp(tau[i] / t1csf) - 1)
+                        - t1csfp
+                        * math.exp(-(t[i] - att) / t1csfp)
+                        * (math.exp(tau[i] / t1csfp) - 1)
+                    )
+                )
+
+                Sb = S1b * math.exp(-te[i] / t2bp)
+                Scsf = S1b * (1 - math.exp(-te[i] / tblcsf)) * math.exp(
+                    -te[i] / t2csf
+                ) + S1csf * math.exp(-te[i] / t2csf)
+        except (OverflowError, RuntimeError):   # pragma: no cover
+            Sb = 0.0
+            Scsf = 0.0
+
+        mag_total[i] = Sb + Scsf
+
+    return mag_total
+
+
+def asl_model_multi_dw(
+    b_values: list, A1: list, D1: float, A2: list, D2: float
+):
+    mag_total = np.zeros(len(b_values))
+
+    for i in range(0, len(b_values)):
+        try:
+            mag_total[i] = A1 * math.exp(-b_values[i] * D1) + A2 * math.exp(
+                -b_values[i] * D2
+            )
+        except (OverflowError, RuntimeError):   # pragma: no cover
+            mag_total[i] = 0.0
+
+    return mag_total

asltk/reconstruction/cbf_mapping.py

@@ -8,8 +8,8 @@
 
 from asltk.asldata import ASLData
 from asltk.aux_methods import _check_mask_values
+from asltk.models.signal_dynamic import asl_model_buxton
 from asltk.mri_parameters import MRIParameters
-from asltk.utils import asl_model_buxton
 
 # Global variables to assist multi cpu threading
 cbf_map = None