fix: calculation of frame timings in the B-Spline scale

src/eddymotion/data/pet.py

@@ -48,8 +48,11 @@ class PET:
     """Best affine for RAS-to-voxel conversion of coordinates (NIfTI header)."""
     brainmask = attr.ib(default=None, repr=_data_repr)
     """A boolean ndarray object containing a corresponding brainmask."""
-    timepoints = attr.ib(default=None, repr=_data_repr)
-    """A 1D numpy array with the timing of each sample."""
+    frame_time = attr.ib(default=None, repr=_data_repr)
+    """A 1D numpy array with the midpoint timing of each sample."""
+    total_duration = attr.ib(default=None, repr=_data_repr)
+    """A float number represaenting the total duration of acquisition."""
+
     em_affines = attr.ib(default=None)
     """
     List of :obj:`nitransforms.linear.Affine` objects that bring
@@ -138,7 +141,7 @@ def from_filename(cls, filename):
 def load(
     filename,
     brainmask_file=None,
-    frame_times=None,
+    frame_time=None,
     frame_duration=None,
 ):
     """Load PET data."""
@@ -152,20 +155,22 @@ def load(
         affine=img.affine,
     )
 
-    if frame_times is not None:
-        retval.timepoints = np.array(frame_times, dtype="float32")
-    elif frame_duration:
-        retval.timepoints = np.array([
-            np.sum(frame_duration[:i])
-            for i in range(1, len(frame_duration) + 1)
-        ])
-    else:
-        raise RuntimeError("Volume timings are necessary")
+    if frame_time is None:
+        raise RuntimeError(
+            "Start time of frames is mandatory (see https://bids-specification.readthedocs.io/"
+            "en/stable/glossary.html#objects.metadata.FrameTimesStart)"
+        )
+
+    frame_time = np.array(frame_time, dtype="float32") - frame_time[0]
+    if frame_duration is None:
+        frame_duration = np.diff(frame_time)
+        if len(frame_duration) == (retval.dataobj.shape[-1] - 1):
+            frame_duration = np.append(frame_duration, frame_duration[-1])
 
-    assert len(retval.timepoints) == retval.dataobj.shape[-1]
+    retval.total_duration = frame_time[-1] + frame_duration[-1]
+    retval.frame_time = frame_time + 0.5 * np.array(frame_duration, dtype="float32")
 
-    # Base at t=0 sec.
-    retval.timepoints = retval.timepoints - retval.timepoints[0]
+    assert len(retval.frame_time) == retval.dataobj.shape[-1]
 
     if brainmask_file:
         mask = nb.load(brainmask_file)

src/eddymotion/estimator.py

@@ -87,8 +87,11 @@ def fit(
             if hasattr(dwdata, "gradients"):
                 kwargs["gtab"] = dwdata.gradients
 
-            if hasattr(dwdata, "timepoints"):
-                kwargs["timepoints"] = dwdata.timepoints
+            if hasattr(dwdata, "frame_time"):
+                kwargs["timepoints"] = dwdata.frame_time
+
+            if hasattr(dwdata, "total_duration"):
+                kwargs["xlim"] = dwdata.total_duration
 
             index_order = np.arange(len(dwdata))
             np.random.shuffle(index_order)

src/eddymotion/model.py

@@ -273,32 +273,38 @@ def predict(self, gradient, **kwargs):
 class PETModel:
     """A PET imaging realignment model based on B-Spline approximation."""
 
-    __slots__ = ("_t", "_x", "_order", "_coeff", "_mask", "_shape", "_n_ctrl")
+    __slots__ = ("_t", "_x", "_xlim", "_order", "_coeff", "_mask", "_shape", "_n_ctrl")
 
-    def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
+    def __init__(self, timepoints=None, xlim=None, n_ctrl=None, mask=None, order=3, **kwargs):
         """
         Create the B-Spline interpolating matrix.
 
         Parameters:
         -----------
         timepoints : :obj:`list`
             The timing (in sec) of each PET volume.
-            E.g., ``[20.,   40.,   60.,  120.,  180.,  240.,  360.,  480.,  600.,
-            900., 1200., 1800., 2400., 3000.]``
+            E.g., ``[15.,   45.,   75.,  105.,  135.,  165.,  210.,  270.,  330.,
+            420.,  540.,  750., 1050., 1350., 1650., 1950., 2250., 2550.]``
 
         n_ctrl : :obj:`int`
             Number of B-Spline control points. If `None`, then one control point every
             six timepoints will be used. The less control points, the smoother is the
             model.
 
         """
-        if timepoints is None:
+        if timepoints is None or xlim is None:
             raise TypeError("timepoints must be provided in initialization")
 
         self._order = order
         self._mask = mask
 
         self._x = np.array(timepoints, dtype="float32")
+        self._xlim = xlim
+
+        if self._x[0] < 1e-2:
+            raise ValueError("First frame midpoint should not be zero or negative")
+        if self._x[-1] > (self._xlim - 1e-2):
+            raise ValueError("Last frame midpoint should not be equal or greater than duration")
 
         # Calculate index coordinates in the B-Spline grid
         self._n_ctrl = n_ctrl or (len(timepoints) // 4) + 1
@@ -314,9 +320,7 @@ def fit(self, data, *args, **kwargs):
         n_jobs = kwargs.pop("n_jobs", None) or 1
 
         timepoints = kwargs.get("timepoints", None) or self._x
-        x = (
-            (np.array(timepoints, dtype="float32") - self._x[0]) / self._x[-1]
-        ) * self._n_ctrl
+        x = (np.array(timepoints, dtype="float32") / self._xlim) * self._n_ctrl
 
         self._shape = data.shape[:3]
 
@@ -351,7 +355,8 @@ def predict(self, timepoint, **kwargs):
         """Return the *b=0* map."""
         from scipy.interpolate import BSpline
 
-        x = ((timepoint - self._x[0]) / self._x[-1]) * self._n_ctrl
+        # Project sample timing into B-Spline coordinates
+        x = (timepoint / self._xlim) * self._n_ctrl
         A = BSpline.design_matrix(x, self._t, k=self._order)
 
         # A is 1 (num. timepoints) x C (num. coeff)