Merge pull request #14 from mjt320/develop

handle more negative inputs

README.md

@@ -42,5 +42,6 @@ Most functionality is demonstrated in Jupyter notebook format in ./demo
 ---
 
 ### Updates
-Release 1.0.2 - Changed AIF interpolation method to linear to further reduce oscillations.
-Release 1.0.1 - Changed AIF interpolation method to avoid oscillations. Added demo notebook on interpolation. 
+Release 1.0.3 - Add exception handling for some zero/negative inputs.  
+Release 1.0.2 - Changed AIF interpolation method to linear to further reduce oscillations.  
+Release 1.0.1 - Changed AIF interpolation method to avoid oscillations. Added demo notebook on interpolation.  

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "sepal"
-version = "1.0.2"
+version = "1.0.3"
 description = "Quantitative MRI processing"
 readme = "README.md"
 authors = [{ name = "Michael Thrippleton", email = "mjt320@googlemail.com" }]

src/sepal/dce_fit.py

@@ -112,9 +112,9 @@ def proc(self, enh, t10, k_fa=1):
         Returns:
             ndarray: 1D array of tissue concentrations (mM)
         """
-        if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa):
+        if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa) or k_fa <= 0 or t10 <= 0:
             raise ValueError(
-                f'Unable to calculate concentration: nan arguments received.')
+                f'Unable to calculate concentration: invalid t10/k_fa/enh values received.')
         e_samples = conc_to_enh(self.C_samples, t10, k_fa, self.c_to_r_model,
                                 self.signal_model)
         C_st = self.C_samples[np.concatenate((argrelextrema(e_samples,
@@ -168,9 +168,9 @@ def proc(self, enh, t10, k_fa=1):
         Returns:
             ndarray: 1D array of tissue concentrations (mM)
         """
-        if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa):
+        if any(np.isnan(enh)) or np.isnan(t10) or np.isnan(k_fa) or t10 <= 0 or k_fa <= 0:
             raise ValueError(
-                f'Unable to calculate concentration: nan arguments received.')
+                f'Unable to calculate concentration: invalid enh/t10/k_fa values received.')
         cos_fa_true = np.cos(k_fa * self.fa)
         exp_r10_tr = np.exp(self.tr/t10)
         C_t = -np.log((exp_r10_tr * (enh-100*cos_fa_true-enh*exp_r10_tr+100)) /
@@ -514,7 +514,7 @@ def pkp_to_enh(pk_pars, hct, k_fa, t10_tissue, t10_blood, pk_model,
         Model to predict one or more exponential relaxation components given
         the relaxation rates for each compartment and water exchange behaviour.
     signal_model : SignalModel
-        Model descriibing the relaxation-signal relationship.
+        Model describing the relaxation-signal relationship.
 
     Returns
     -------

src/sepal/fitting.py

@@ -82,7 +82,7 @@ def proc_image(self, input_images, arg_images=None, mask=None,
                 contain 1 or 0 only. 1 indicates voxels to be processed.
                 Defaults to None (process all voxels).
             threshold (float): Voxel is processed if max input value in
-                series (e.g. flip angle or time series) is >= threshold.
+                series (e.g. flip angle series or time series) is >= threshold.
                 Defaults to -np.inf
             dir (str): Directory for output images. If None, no output
                 images are written. Defaults to None.

src/sepal/t1_fit.py

@@ -58,9 +58,9 @@ def proc(self, s, k_fa=1):
                 t1 (float): T1 (s)
 
         """
-        if any(np.isnan(s)):
+        if any(np.isnan(s)) or np.isnan(k_fa) or k_fa <= 0:
             raise ValueError(
-                f'Unable to calculate T1: nan signal values received.')
+                f'Unable to calculate T1: nan signal or nan/<=0 k_fa values received.')
         with np.errstate(divide='ignore', invalid='ignore'):
             fa_true = k_fa * self.fa_rad
             sr = s[0] / s[1]
@@ -112,9 +112,9 @@ def proc(self, s, k_fa=1):
                 t1 (float): T1 (s)
 
         """
-        if any(np.isnan(s)) or np.isnan(k_fa):
+        if any(np.isnan(s)) or np.isnan(k_fa) or k_fa <= 0:
             raise ArithmeticError(
-                f'Unable to calculate T1: nan signal or k_fa values received.')
+                f'Unable to calculate T1: nan signal or nan/<=0 k_fa values received.')
         fa_true = k_fa * self.fa_rad
         y = s / np.sin(fa_true)
         x = s / np.tan(fa_true)
@@ -165,9 +165,9 @@ def proc(self, s, k_fa=1):
                 t1 (float): T1 (s)
 
         """
-        if any(np.isnan(s)) or np.isnan(k_fa):
+        if any(np.isnan(s)) or np.isnan(k_fa) or k_fa <= 0:
             raise ValueError(
-                f'Unable to calculate T1: nan signal or k_fa values received.')
+                f'Unable to calculate T1: nan signal or nan/<=0 k_fa values received.')
         # use linear fit to obtain initial guess, otherwise start with T1=1
         try:
             x0 = np.array(self.linear_fitter.proc(s, k_fa=k_fa))