From 7f3172ca28a04aeb3145a40d223ed46109387052 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 12 Dec 2022 17:54:15 +0100 Subject: [PATCH] fix: calculation of frame timings in the B-Spline scale --- src/eddymotion/data/pet.py | 35 ++++++++++++++++++++--------------- src/eddymotion/estimator.py | 7 +++++-- src/eddymotion/model.py | 19 ++++++++++++------- 3 files changed, 37 insertions(+), 24 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 16f756d6..455a7204 100644 --- a/src/eddymotion/estimator.py +++ b/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) diff --git a/src/eddymotion/model.py b/src/eddymotion/model.py index 15830514..f59e4b90 100644 --- a/src/eddymotion/model.py +++ b/src/eddymotion/model.py @@ -273,9 +273,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. @@ -292,13 +292,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 @@ -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)