From aa8491a18d54563d1ad7fed63e1f6501ad7b03d8 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Sat, 10 Dec 2022 16:27:58 +0100
Subject: [PATCH 1/9] enh: initial draft of the PET uptake model

Resolves: #66.
---
 src/eddymotion/data/pet.py   | 170 +++++++++++++++++++++++++++++++++++
 src/eddymotion/model/base.py |  76 ++++++++++++++++
 2 files changed, 246 insertions(+)
 create mode 100644 src/eddymotion/data/pet.py

diff --git a/src/eddymotion/data/pet.py b/src/eddymotion/data/pet.py
new file mode 100644
index 00000000..22328d73
--- /dev/null
+++ b/src/eddymotion/data/pet.py
@@ -0,0 +1,170 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2022 The NiPreps Developers <nipreps@gmail.com>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""PET data representation."""
+from collections import namedtuple
+from pathlib import Path
+from tempfile import mkdtemp
+
+import attr
+import h5py
+import nibabel as nb
+import numpy as np
+from nitransforms.linear import Affine
+
+
+def _data_repr(value):
+    if value is None:
+        return "None"
+    return f"<{'x'.join(str(v) for v in value.shape)} ({value.dtype})>"
+
+
+@attr.s(slots=True)
+class PET:
+    """Data representation structure for dMRI data."""
+
+    dataobj = attr.ib(default=None, repr=_data_repr)
+    """A numpy ndarray object for the data array, without *b=0* volumes."""
+    affine = attr.ib(default=None, repr=_data_repr)
+    """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."""
+    em_affines = attr.ib(default=None)
+    """
+    List of :obj:`nitransforms.linear.Affine` objects that bring
+    PET timepoints into alignment.
+    """
+    _filepath = attr.ib(
+        factory=lambda: Path(mkdtemp()) / "em_cache.h5",
+        repr=False,
+    )
+    """A path to an HDF5 file to store the whole dataset."""
+
+    def __len__(self):
+        """Obtain the number of high-*b* orientations."""
+        return self.dataobj.shape[-1]
+
+    def set_transform(self, index, affine, order=3):
+        """Set an affine, and update data object and gradients."""
+        reference = namedtuple("ImageGrid", ("shape", "affine"))(
+            shape=self.dataobj.shape[:3], affine=self.affine
+        )
+        xform = Affine(matrix=affine, reference=reference)
+
+        if not Path(self._filepath).exists():
+            self.to_filename(self._filepath)
+
+        # read original PET
+        with h5py.File(self._filepath, "r") as in_file:
+            root = in_file["/0"]
+            dframe = np.asanyarray(root["dataobj"][..., index])
+
+        dmoving = nb.Nifti1Image(dframe, self.affine, None)
+
+        # resample and update orientation at index
+        self.dataobj[..., index] = np.asanyarray(
+            xform.apply(dmoving, order=order).dataobj,
+            dtype=self.dataobj.dtype,
+        )
+
+        # update transform
+        if self.em_affines is None:
+            self.em_affines = [None] * len(self)
+
+        self.em_affines[index] = xform
+
+    def to_filename(self, filename, compression=None, compression_opts=None):
+        """Write an HDF5 file to disk."""
+        filename = Path(filename)
+        if not filename.name.endswith(".h5"):
+            filename = filename.parent / f"{filename.name}.h5"
+
+        with h5py.File(filename, "w") as out_file:
+            out_file.attrs["Format"] = "EMC/PET"
+            out_file.attrs["Version"] = np.uint16(1)
+            root = out_file.create_group("/0")
+            root.attrs["Type"] = "pet"
+            for f in attr.fields(self.__class__):
+                if f.name.startswith("_"):
+                    continue
+
+                value = getattr(self, f.name)
+                if value is not None:
+                    root.create_dataset(
+                        f.name,
+                        data=value,
+                        compression=compression,
+                        compression_opts=compression_opts,
+                    )
+
+    def to_nifti(self, filename, insert_b0=False):
+        """Write a NIfTI 1.0 file to disk."""
+        nii = nb.Nifti1Image(self.dataobj, self.affine, None)
+        nii.header.set_xyzt_units("mm")
+        nii.to_filename(filename)
+
+    @classmethod
+    def from_filename(cls, filename):
+        """Read an HDF5 file from disk."""
+        with h5py.File(filename, "r") as in_file:
+            root = in_file["/0"]
+            data = {
+                k: np.asanyarray(v) for k, v in root.items() if not k.startswith("_")
+            }
+        return cls(**data)
+
+
+def load(
+    filename,
+    brainmask_file=None,
+    volume_timings=None,
+    volume_spacings=None,
+):
+    """Load PET data."""
+    filename = Path(filename)
+    if filename.name.endswith(".h5"):
+        return PET.from_filename(filename)
+
+    img = nb.load(filename)
+    retval = PET(
+        dataobj=img.get_fdata(dtype="float32"),
+        affine=img.affine,
+    )
+
+    if volume_timings is not None:
+        retval.timepoints = np.array(volume_timings)
+    elif volume_spacings:
+        x = np.array([
+            np.sum(volume_spacings[:i])
+            for i in range(1, len(volume_spacings) + 1)
+        ])
+        retval.timepoints = x - x[0]
+    else:
+        raise RuntimeError("Volume timings are necessary")
+
+    if brainmask_file:
+        mask = nb.load(brainmask_file)
+        retval.brainmask = np.asanyarray(mask.dataobj)
+
+    return retval
diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py
index 65b76dd6..24da4dce 100644
--- a/src/eddymotion/model/base.py
+++ b/src/eddymotion/model/base.py
@@ -293,6 +293,82 @@ def predict(self, gradient, **kwargs):
         return self._data
 
 
+class PETModel:
+    """A PET imaging realignment model based on B-Spline approximation."""
+
+    __slots__ = ("_t", "_x0", "_x1", "_order", "_coeff", "_mask", "_shape")
+
+    def __init__(self, timepoints, 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.]``
+
+        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.
+
+        """
+        self._order = order
+        self._mask = mask
+
+        x = np.array(timepoints)
+        self._x0 = x[0]
+        x -= self._x0
+
+        self._x1 = x[-1]
+        x /= self._x1
+
+        # Calculate index coordinates in the B-Spline grid
+        n_ctrl = n_ctrl or (len(timepoints) // 6) + 1
+        x *= n_ctrl
+
+        # B-Spline knots
+        self._t = np.linspace(-2.0, float(n_ctrl) + 2.0, n_ctrl + 4)
+
+    def fit(self, data, timepoints, *args, **kwargs):
+        """Do nothing."""
+        from scipy.interpolate import BSpline
+        from scipy.sparse.linalg import cg
+
+        self._shape = data.shape[:3]
+
+        # Convert data into V (voxels) x T (timepoints)
+        data = (
+            data.reshape((-1, data.shape[-1]))
+            if self._mask is None else data[self._mask]
+        )
+
+        x = (np.array(timepoints) - self._x0) / self._x1
+        A = BSpline.design_matrix(x, self._t, k=self._order)
+
+        self._coeff, _ = cg(A.T @ A, A.T @ data)
+
+    def predict(self, timepoint, **kwargs):
+        """Return the *b=0* map."""
+        from scipy.interpolate import BSpline
+
+        x = (timepoint - self._x0) / self._x1
+        A = BSpline.design_matrix(x, self._t, k=self._order)
+
+        # A is 1 (num. timepoints) x C (num. coeff)
+        # self._coeff is C (num. coeff) x V (num. voxels)
+        predicted = (A @ self._coeff).T
+
+        if self._mask is None:
+            return predicted.reshape(self._shape)
+
+        retval = np.zeros(self._shape, dtype="float32")
+        retval[self._mask] = predicted
+        return retval
+
+
 class DTIModel(BaseModel):
     """A wrapper of :obj:`dipy.reconst.dti.TensorModel`."""
 

From 5dad300843c51c22c1283a3075a354b154908c48 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Sat, 10 Dec 2022 20:56:41 +0100
Subject: [PATCH 2/9] enh: revise loading function and object internals

---
 src/eddymotion/data/pet.py | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/src/eddymotion/data/pet.py b/src/eddymotion/data/pet.py
index 22328d73..3bac542d 100644
--- a/src/eddymotion/data/pet.py
+++ b/src/eddymotion/data/pet.py
@@ -138,8 +138,8 @@ def from_filename(cls, filename):
 def load(
     filename,
     brainmask_file=None,
-    volume_timings=None,
-    volume_spacings=None,
+    frame_times=None,
+    frame_duration=None,
 ):
     """Load PET data."""
     filename = Path(filename)
@@ -152,17 +152,21 @@ def load(
         affine=img.affine,
     )
 
-    if volume_timings is not None:
-        retval.timepoints = np.array(volume_timings)
-    elif volume_spacings:
-        x = np.array([
-            np.sum(volume_spacings[:i])
-            for i in range(1, len(volume_spacings) + 1)
+    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)
         ])
-        retval.timepoints = x - x[0]
     else:
         raise RuntimeError("Volume timings are necessary")
 
+    assert len(retval.timepoints) == retval.dataobj.shape[-1]
+
+    # Base at t=0 sec.
+    retval.timepoints = retval.timepoints - retval.timepoints[0]
+
     if brainmask_file:
         mask = nb.load(brainmask_file)
         retval.brainmask = np.asanyarray(mask.dataobj)

From 28ece64c1a988b54d0c3cca127c317f4b1471770 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Sat, 10 Dec 2022 20:57:07 +0100
Subject: [PATCH 3/9] fix: adapt estimator to better process PET models

---
 src/eddymotion/estimator.py  | 33 ++++++++++-------
 src/eddymotion/model/base.py | 70 ++++++++++++++++++++++++------------
 2 files changed, 68 insertions(+), 35 deletions(-)

diff --git a/src/eddymotion/estimator.py b/src/eddymotion/estimator.py
index 6c7a4970..41294af4 100644
--- a/src/eddymotion/estimator.py
+++ b/src/eddymotion/estimator.py
@@ -83,16 +83,6 @@ def fit(
         if seed or seed == 0:
             np.random.seed(20210324 if seed is True else seed)
 
-        bmask_img = None
-        if dwdata.brainmask is not None:
-            _, bmask_img = mkstemp(suffix="_bmask.nii.gz")
-            nb.Nifti1Image(
-                dwdata.brainmask.astype("uint8"), dwdata.affine, None
-            ).to_filename(bmask_img)
-            kwargs["mask"] = dwdata.brainmask
-
-        kwargs["S0"] = _advanced_clip(dwdata.bzero)
-
         if "num_threads" not in align_kwargs and omp_nthreads is not None:
             align_kwargs["num_threads"] = omp_nthreads
 
@@ -103,6 +93,25 @@ def fit(
                 if model.lower() not in ("b0", "s0", "avg", "average", "mean")
                 else "b0"
             )
+
+            # When downsampling these need to be set per-level
+            bmask_img = None
+            if dwdata.brainmask is not None:
+                _, bmask_img = mkstemp(suffix="_bmask.nii.gz")
+                nb.Nifti1Image(
+                    dwdata.brainmask.astype("uint8"), dwdata.affine, None
+                ).to_filename(bmask_img)
+                kwargs["mask"] = dwdata.brainmask
+
+            if hasattr(dwdata, "bzero") and dwdata.bzero is not None:
+                kwargs["S0"] = _advanced_clip(dwdata.bzero)
+
+            if hasattr(dwdata, "gradients"):
+                kwargs["gtab"] = dwdata.gradients
+
+            if hasattr(dwdata, "timepoints"):
+                kwargs["timepoints"] = dwdata.timepoints
+
             index_order = np.arange(len(dwdata))
             np.random.shuffle(index_order)
 
@@ -118,7 +127,6 @@ def fit(
 
                 # Factory creates the appropriate model and pipes arguments
                 dwmodel = ModelFactory.init(
-                    gtab=dwdata.gradients,
                     model=model,
                     **kwargs,
                 )
@@ -137,9 +145,10 @@ def fit(
                         pbar.set_description_str(f"[{grad_str}], {n_jobs} jobs")
 
                         if not single_model:  # A true LOGO estimator
+                            if hasattr(dwdata, "gradients"):
+                                kwargs["gtab"] = data_train[1]
                             # Factory creates the appropriate model and pipes arguments
                             dwmodel = ModelFactory.init(
-                                gtab=data_train[1],
                                 model=model,
                                 n_jobs=n_jobs,
                                 **kwargs,
diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py
index 24da4dce..ff9cc3d7 100644
--- a/src/eddymotion/model/base.py
+++ b/src/eddymotion/model/base.py
@@ -41,7 +41,7 @@ class ModelFactory:
     """A factory for instantiating diffusion models."""
 
     @staticmethod
-    def init(gtab, model="DTI", **kwargs):
+    def init(model="DTI", **kwargs):
         """
         Instatiate a diffusion model.
 
@@ -60,15 +60,14 @@ def init(gtab, model="DTI", **kwargs):
 
         """
         if model.lower() in ("s0", "b0"):
-            return TrivialB0Model(gtab=gtab, S0=kwargs.pop("S0"))
+            return TrivialB0Model(S0=kwargs.pop("S0"))
 
         if model.lower() in ("avg", "average", "mean"):
-            return AverageDWModel(gtab=gtab, **kwargs)
+            return AverageDWModel(**kwargs)
 
-        # Generate a GradientTable object for DIPY
-        if model.lower() in ("dti", "dki"):
+        if model.lower() in ("dti", "dki", "pet"):
             Model = globals()[f"{model.upper()}Model"]
-            return Model(gtab, **kwargs)
+            return Model(**kwargs)
 
         raise NotImplementedError(f"Unsupported model <{model}>.")
 
@@ -217,7 +216,7 @@ class TrivialB0Model:
 
     __slots__ = ("_S0",)
 
-    def __init__(self, gtab, S0=None, **kwargs):
+    def __init__(self, S0=None, **kwargs):
         """Implement object initialization."""
         if S0 is None:
             raise ValueError("S0 must be provided")
@@ -237,7 +236,7 @@ class AverageDWModel:
 
     __slots__ = ("_data", "_th_low", "_th_high", "_bias", "_stat")
 
-    def __init__(self, gtab, **kwargs):
+    def __init__(self, **kwargs):
         r"""
         Implement object initialization.
 
@@ -296,9 +295,9 @@ def predict(self, gradient, **kwargs):
 class PETModel:
     """A PET imaging realignment model based on B-Spline approximation."""
 
-    __slots__ = ("_t", "_x0", "_x1", "_order", "_coeff", "_mask", "_shape")
+    __slots__ = ("_t", "_x", "_order", "_coeff", "_mask", "_shape", "_n_ctrl")
 
-    def __init__(self, timepoints, n_ctrl=None, mask=None, order=3, **kwargs):
+    def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
         """
         Create the B-Spline interpolating matrix.
 
@@ -315,28 +314,36 @@ def __init__(self, timepoints, n_ctrl=None, mask=None, order=3, **kwargs):
             model.
 
         """
+        if timepoints is None:
+            raise TypeError("timepoints must be provided in initialization")
+
         self._order = order
         self._mask = mask
 
-        x = np.array(timepoints)
-        self._x0 = x[0]
-        x -= self._x0
-
-        self._x1 = x[-1]
-        x /= self._x1
+        self._x = np.array(timepoints, dtype="float32")
+        self._x -= self._x[0]
+        self._x /= self._x[-1]
 
         # Calculate index coordinates in the B-Spline grid
-        n_ctrl = n_ctrl or (len(timepoints) // 6) + 1
-        x *= n_ctrl
+        self._n_ctrl = n_ctrl or (len(timepoints) // 8) + 1
 
         # B-Spline knots
-        self._t = np.linspace(-2.0, float(n_ctrl) + 2.0, n_ctrl + 4)
+        self._t = np.arange(-3, float(self._n_ctrl) + 4, dtype="float32")
 
-    def fit(self, data, timepoints, *args, **kwargs):
+    def fit(self, data, *args, **kwargs):
         """Do nothing."""
         from scipy.interpolate import BSpline
         from scipy.sparse.linalg import cg
 
+        n_jobs = kwargs.pop("n_jobs", None) or 1
+
+        timepoints = kwargs.get("timepoints", None)
+        x = (
+            (np.array(timepoints, dtype="float32") - self._x[0]) / self._x[-1]
+            if timepoints is not None
+            else self._x
+        ) * self._n_ctrl
+
         self._shape = data.shape[:3]
 
         # Convert data into V (voxels) x T (timepoints)
@@ -345,10 +352,27 @@ def fit(self, data, timepoints, *args, **kwargs):
             if self._mask is None else data[self._mask]
         )
 
-        x = (np.array(timepoints) - self._x0) / self._x1
+        # A.shape = (T, K - 4); T= n. timepoints, K= n. knots (with padding)
         A = BSpline.design_matrix(x, self._t, k=self._order)
+        AT = A.T
+        ATdotA = AT @ A
+
+        # One single CPU - linear execution (full model)
+        if n_jobs == 1:
+            self._coeff = np.array([
+                cg(ATdotA, AT @ v)[0]
+                for v in data
+            ])
+            return
+
+        # Parallelize process with joblib
+        with Parallel(n_jobs=n_jobs) as executor:
+            results = executor(
+                delayed(cg)(ATdotA, AT @ v)
+                for v in data
+            )
 
-        self._coeff, _ = cg(A.T @ A, A.T @ data)
+        self._coeff = np.array(results)
 
     def predict(self, timepoint, **kwargs):
         """Return the *b=0* map."""
@@ -358,7 +382,7 @@ def predict(self, timepoint, **kwargs):
         A = BSpline.design_matrix(x, self._t, k=self._order)
 
         # A is 1 (num. timepoints) x C (num. coeff)
-        # self._coeff is C (num. coeff) x V (num. voxels)
+        # self._coeff is V (num. voxels) x K - 4
         predicted = (A @ self._coeff).T
 
         if self._mask is None:

From a70fe98fbe65f1a96a6ef04b6d2f1c46291d38a8 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Sat, 10 Dec 2022 22:32:33 +0100
Subject: [PATCH 4/9] fix: final refinements to model

---
 src/eddymotion/model/base.py | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py
index ff9cc3d7..400a1e7a 100644
--- a/src/eddymotion/model/base.py
+++ b/src/eddymotion/model/base.py
@@ -321,11 +321,9 @@ def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
         self._mask = mask
 
         self._x = np.array(timepoints, dtype="float32")
-        self._x -= self._x[0]
-        self._x /= self._x[-1]
 
         # Calculate index coordinates in the B-Spline grid
-        self._n_ctrl = n_ctrl or (len(timepoints) // 8) + 1
+        self._n_ctrl = n_ctrl or (len(timepoints) // 4) + 1
 
         # B-Spline knots
         self._t = np.arange(-3, float(self._n_ctrl) + 4, dtype="float32")
@@ -337,11 +335,9 @@ def fit(self, data, *args, **kwargs):
 
         n_jobs = kwargs.pop("n_jobs", None) or 1
 
-        timepoints = kwargs.get("timepoints", None)
+        timepoints = kwargs.get("timepoints", None) or self._x
         x = (
             (np.array(timepoints, dtype="float32") - self._x[0]) / self._x[-1]
-            if timepoints is not None
-            else self._x
         ) * self._n_ctrl
 
         self._shape = data.shape[:3]
@@ -360,8 +356,7 @@ def fit(self, data, *args, **kwargs):
         # One single CPU - linear execution (full model)
         if n_jobs == 1:
             self._coeff = np.array([
-                cg(ATdotA, AT @ v)[0]
-                for v in data
+                cg(ATdotA, AT @ v)[0] for v in data
             ])
             return
 
@@ -372,18 +367,18 @@ def fit(self, data, *args, **kwargs):
                 for v in data
             )
 
-        self._coeff = np.array(results)
+        self._coeff = np.array([r[0] for r in results])
 
     def predict(self, timepoint, **kwargs):
         """Return the *b=0* map."""
         from scipy.interpolate import BSpline
 
-        x = (timepoint - self._x0) / self._x1
+        x = ((timepoint - self._x[0]) / self._x[-1]) * self._n_ctrl
         A = BSpline.design_matrix(x, self._t, k=self._order)
 
         # A is 1 (num. timepoints) x C (num. coeff)
         # self._coeff is V (num. voxels) x K - 4
-        predicted = (A @ self._coeff).T
+        predicted = np.squeeze(A @ self._coeff.T)
 
         if self._mask is None:
             return predicted.reshape(self._shape)

From 4fdeecf4166ed633e136dbe1acc5e8833177154d Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Sat, 10 Dec 2022 22:37:56 +0100
Subject: [PATCH 5/9] fix: failing test

---
 test/test_model.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/test_model.py b/test/test_model.py
index 22c9b75f..4b3187d1 100644
--- a/test/test_model.py
+++ b/test/test_model.py
@@ -98,7 +98,6 @@ def test_two_initialisations(datadir):
 
     # Direct initialisation
     model1 = model.AverageDWModel(
-        dmri_dataset.gradients,
         S0=dmri_dataset.bzero,
         th_low=100,
         th_high=1000,

From ef5c7f087594eb16c7836a7b206a5deb3595ebe1 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Mon, 12 Dec 2022 15:29:24 +0100
Subject: [PATCH 6/9] Apply suggestions from code review

Co-authored-by: Ariel Rokem <arokem@gmail.com>
---
 src/eddymotion/data/pet.py   | 2 +-
 src/eddymotion/model/base.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/eddymotion/data/pet.py b/src/eddymotion/data/pet.py
index 3bac542d..7de22d6c 100644
--- a/src/eddymotion/data/pet.py
+++ b/src/eddymotion/data/pet.py
@@ -40,7 +40,7 @@ def _data_repr(value):
 
 @attr.s(slots=True)
 class PET:
-    """Data representation structure for dMRI data."""
+    """Data representation structure for PET data."""
 
     dataobj = attr.ib(default=None, repr=_data_repr)
     """A numpy ndarray object for the data array, without *b=0* volumes."""
diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py
index 400a1e7a..977253f7 100644
--- a/src/eddymotion/model/base.py
+++ b/src/eddymotion/model/base.py
@@ -329,7 +329,7 @@ def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
         self._t = np.arange(-3, float(self._n_ctrl) + 4, dtype="float32")
 
     def fit(self, data, *args, **kwargs):
-        """Do nothing."""
+        """Fit the model."""
         from scipy.interpolate import BSpline
         from scipy.sparse.linalg import cg
 

From b41d1057585aa960e69e332c5a1ba010a5f13731 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Mon, 12 Dec 2022 17:54:15 +0100
Subject: [PATCH 7/9] fix: calculation of frame timings in the B-Spline scale

---
 src/eddymotion/data/pet.py   | 35 ++++++++++++++++++++---------------
 src/eddymotion/estimator.py  |  7 +++++--
 src/eddymotion/model/base.py | 23 ++++++++++++++---------
 3 files changed, 39 insertions(+), 26 deletions(-)

diff --git a/src/eddymotion/data/pet.py b/src/eddymotion/data/pet.py
index 7de22d6c..415b477c 100644
--- a/src/eddymotion/data/pet.py
+++ b/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)
diff --git a/src/eddymotion/estimator.py b/src/eddymotion/estimator.py
index 41294af4..9624c1e3 100644
--- a/src/eddymotion/estimator.py
+++ b/src/eddymotion/estimator.py
@@ -109,8 +109,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)
diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py
index 977253f7..59b7a863 100644
--- a/src/eddymotion/model/base.py
+++ b/src/eddymotion/model/base.py
@@ -295,9 +295,9 @@ 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.
 
@@ -305,8 +305,8 @@ def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
         -----------
         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
@@ -314,13 +314,19 @@ def __init__(self, timepoints=None, n_ctrl=None, mask=None, order=3, **kwargs):
             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
@@ -336,9 +342,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]
 
@@ -373,7 +377,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)

From 2cee6a912529088284f99821cc444b933450fef6 Mon Sep 17 00:00:00 2001
From: Martin Norgaard <martin.noergaard@nru.dk>
Date: Thu, 15 Dec 2022 12:02:33 -0800
Subject: [PATCH 8/9] FIX: Update __init__.py to inclyde PETModel

---
 src/eddymotion/model/__init__.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/eddymotion/model/__init__.py b/src/eddymotion/model/__init__.py
index 6fa40bd8..e2eb4565 100644
--- a/src/eddymotion/model/__init__.py
+++ b/src/eddymotion/model/__init__.py
@@ -27,6 +27,7 @@
     DKIModel,
     DTIModel,
     TrivialB0Model,
+    PETModel,
 )
 
 __all__ = (
@@ -35,4 +36,5 @@
     "DKIModel",
     "DTIModel",
     "TrivialB0Model",
+    "PETModel",
 )

From c046fbadb7e04bd126fec93a28905d9a2673c395 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Fri, 16 Dec 2022 09:24:04 +0100
Subject: [PATCH 9/9] fix: pin scipy >= 1.8.0

Co-authored-by: Martin Norgaard <martin.noergaard@nru.dk>
---
 docs/conf.py | 1 +
 setup.cfg    | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/conf.py b/docs/conf.py
index 65645bbd..e3c3419b 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -55,6 +55,7 @@
     "pandas",
     "seaborn",
     "skimage",
+    "sklearn",
     "svgutils",
     "tqdm",
     "transforms3d",
diff --git a/setup.cfg b/setup.cfg
index 327ef964..f093ade8 100755
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,9 +30,10 @@ install_requires =
     joblib
     nipype>= 1.5.1, < 2.0
     nitransforms>=21.0.0
+    numpy>=1.17.3
     nest-asyncio>=1.5.1
     scikit-image>=0.14.2
-    scikit-learn>=1.0.1
+    scipy>=1.8.0
 test_requires =
     codecov
     coverage