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)