Merge pull request #11 from mjt320/develop

modify AIF interpolation method

README.md

@@ -8,6 +8,7 @@ Created 28 September 2020
 @authors: Michael Thrippleton  
 @email: [email protected]  
 @institution: University of Edinburgh, UK
+---
 
 ### Installation:
 pip install sepal
@@ -38,3 +39,7 @@ Most functionality is demonstrated in Jupyter notebook format in ./demo
 - R2/R2* effects not included in fitting of enhancement curves (but is included for enhancement-to-concentration conversion)
 - Compartment-specific relaxivity parameters/models
 - Fitting water exchange parameters
+---
+
+### Updates
+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.0"
+version = "1.0.1"
 description = "Quantitative MRI processing"
 readme = "README.md"
 authors = [{ name = "Michael Thrippleton", email = "[email protected]" }]

src/sepal/aifs.py

@@ -17,7 +17,7 @@
 from abc import ABC, abstractmethod
 
 import numpy as np
-from scipy.interpolate import interp1d
+from scipy.interpolate import interp1d, Akima1DInterpolator
 
 
 class AIF(ABC):
@@ -71,14 +71,19 @@ def __init__(self, t_data, c_ap_data):
         """
         self.t_data = t_data
         self.c_ap_data = c_ap_data
-        self.c_ap_func = interp1d(t_data, c_ap_data,
-                                  kind='quadratic', bounds_error=False,
-                                  fill_value=(0, c_ap_data[-1]))
+        self.c_ap_func = Akima1DInterpolator(t_data, c_ap_data)
 
     def c_ap(self, t):
         """Get AIF plasma concentration(t). Overrides superclass method."""
         #  calculate concentration(t) using interpolation function
         c_ap = self.c_ap_func(t)
+        #  replace values outside time range with first/last values
+        idx_before = t < min(self.t_data)
+        idx_after = t > max(self.t_data)
+        c_ap_before = self.c_ap_data[0]
+        c_ap_after = self.c_ap_data[-1]
+        c_ap[idx_before] = c_ap_before
+        c_ap[idx_after] = c_ap_after
         return c_ap