diff --git a/bids/analysis/analysis.py b/bids/analysis/analysis.py index fd3c28660..d5fcf5666 100644 --- a/bids/analysis/analysis.py +++ b/bids/analysis/analysis.py @@ -386,20 +386,32 @@ def get_design_matrix(self, names=None, format='long', mode='both', kwargs['timing'] = True kwargs['sparse'] = False + acquisition_time = None if sampling_rate == 'TR': trs = {var.run_info[0].tr for var in self.collection.variables.values()} + tas = {var.run_info[0].ta for var in self.collection.variables.values()} if not trs: raise ValueError("Repetition time unavailable; specify sampling_rate " "explicitly") elif len(trs) > 1: raise ValueError("Non-unique Repetition times found ({!r}); specify " - "sampling_rate explicitly") - sampling_rate = 1. / trs.pop() + "sampling_rate explicitly".format(trs)) + TR = trs.pop() + if not tas: + warnings.warn("Acquisition time unavailable; assuming TA = TR") + tas = {TR} + elif len(tas) > 1: + raise ValueError("Non-unique acquisition times found ({!r})".format(tas)) + + sampling_rate = 1. / TR + acquisition_time = tas.pop() elif sampling_rate == 'highest': sampling_rate = None dense_df = coll.to_df(names, format='wide', include_sparse=include_sparse, - sampling_rate=sampling_rate, **kwargs) + sampling_rate=sampling_rate, + integration_window=acquisition_time, + **kwargs) if dense_df is not None: dense_df = dense_df.drop(['onset', 'duration'], axis=1) diff --git a/bids/analysis/tests/test_transformations.py b/bids/analysis/tests/test_transformations.py index 21610289a..e8166d334 100644 --- a/bids/analysis/tests/test_transformations.py +++ b/bids/analysis/tests/test_transformations.py @@ -28,7 +28,7 @@ def sparse_run_variable_with_missing_values(): 'duration': [1.2, 1.6, 0.8, 2], 'amplitude': [1, 1, np.nan, 1] }) - run_info = [RunInfo({'subject': '01'}, 20, 2, 'dummy.nii.gz')] + run_info = [RunInfo({'subject': '01'}, 20, 2, 2, 'dummy.nii.gz')] var = SparseRunVariable(name='var', data=data, run_info=run_info, source='events') return BIDSRunVariableCollection([var]) diff --git a/bids/analysis/transformations/compute.py b/bids/analysis/transformations/compute.py index 153bb846d..3254740d3 100644 --- a/bids/analysis/transformations/compute.py +++ b/bids/analysis/transformations/compute.py @@ -9,6 +9,14 @@ from bids.analysis import hrf from bids.variables import SparseRunVariable, DenseRunVariable +try: + from math import gcd +except ImportError: + def gcd(a, b): + if b == 0: + return a + return gcd(b, a % b) + class Convolve(Transformation): """Convolve the input variable with an HRF. @@ -35,6 +43,25 @@ def _transform(self, var, model='spm', derivative=False, dispersion=False, if isinstance(var, SparseRunVariable): sr = self.collection.sampling_rate + # Resolve diferences in TR and TA to milliseconds + trs = {ri.tr for ri in var.run_info} + tas = {ri.ta for ri in var.run_info} + assert len(trs) == 1 + assert len(tas) == 1 + TR = int(np.round(1000. * trs.pop())) + TA = int(np.round(1000. * tas.pop())) + SR = int(np.round(1000. / sr)) + if TA < TR: + # Use a unit that fits an whole number of times into both + # the interscan interval (TR) and the integration window (TA) + dt = gcd(TR, TA) + if dt > SR: + # Find the nearest whole factor of dt >= SR + # This compromises on the target sampling rate and + # one that neatly bins TR and TA + dt = dt // (dt // SR) + sr = 1000. / dt + var = var.to_dense(sr) df = var.to_df(entities=False) diff --git a/bids/variables/entities.py b/bids/variables/entities.py index a1060f6f3..edde25be5 100644 --- a/bids/variables/entities.py +++ b/bids/variables/entities.py @@ -31,28 +31,28 @@ class RunNode(Node): ''' Represents a single Run in a BIDS project. Args: - id (int): The index of the run. entities (dict): Dictionary of entities for this Node. image_file (str): The full path to the corresponding nifti image. duration (float): Duration of the run, in seconds. repetition_time (float): TR for the run. - task (str): The task name for this run. + acquisition_time (float): TA for the run. ''' - def __init__(self, entities, image_file, duration, repetition_time): + def __init__(self, entities, image_file, duration, repetition_time, acquisition_time=None): self.image_file = image_file self.duration = duration self.repetition_time = repetition_time + self.acquisition_time = acquisition_time or repetition_time super(RunNode, self).__init__('run', entities) def get_info(self): return RunInfo(self.entities, self.duration, self.repetition_time, - self.image_file) + self.acquisition_time, self.image_file) # Stores key information for each Run. -RunInfo = namedtuple('RunInfo', ['entities', 'duration', 'tr', 'image']) +RunInfo = namedtuple('RunInfo', ['entities', 'duration', 'tr', 'ta', 'image']) class NodeIndex(Node): diff --git a/bids/variables/io.py b/bids/variables/io.py index b73cc68fa..e3c265be9 100644 --- a/bids/variables/io.py +++ b/bids/variables/io.py @@ -85,6 +85,27 @@ def load_variables(layout, types=None, levels=None, skip_empty=True, return dataset +def _get_timing_info(img_md): + if 'RepetitionTime' in img_md: + tr = img_md['RepetitionTime'] + if 'DelayTime' in img_md: + ta = tr - img_md['DelayTime'] + elif 'SliceTiming' in img_md: + slicetimes = sorted(img_md['SliceTiming']) + # a, b ... z + # z = final slice onset, b - a = slice duration + ta = np.round(slicetimes[-1] + slicetimes[1] - slicetimes[0], 3) + # If the "silence" is <1ms, consider acquisition continuous + if np.abs(tr - ta) < 1e-3: + ta = tr + else: + ta = tr + elif 'VolumeTiming' in img_md: + return NotImplemented + + return tr, ta + + def _load_time_variables(layout, dataset=None, columns=None, scan_length=None, drop_na=True, events=True, physio=True, stim=True, regressors=True, skip_empty=True, scope='all', @@ -146,8 +167,9 @@ def _load_time_variables(layout, dataset=None, columns=None, scan_length=None, if 'run' in entities: entities['run'] = int(entities['run']) - tr = layout.get_metadata(img_f, suffix='bold', scope=scope, - full_search=True)['RepetitionTime'] + img_md = layout.get_metadata(img_f, suffix='bold', scope=scope, + full_search=True) + tr, ta = _get_timing_info(img_md) # Get duration of run: first try to get it directly from the image # header; if that fails, try to get NumberOfVolumes from the @@ -167,7 +189,8 @@ def _load_time_variables(layout, dataset=None, columns=None, scan_length=None, raise ValueError(msg) run = dataset.get_or_create_node('run', entities, image_file=img_f, - duration=duration, repetition_time=tr) + duration=duration, repetition_time=tr, + acquisition_time=ta) run_info = run.get_info() # Process event files diff --git a/bids/variables/kollekshuns.py b/bids/variables/kollekshuns.py index 7b6687b83..411f46570 100644 --- a/bids/variables/kollekshuns.py +++ b/bids/variables/kollekshuns.py @@ -243,7 +243,7 @@ def _all_dense(self): for v in self.variables.values()]) def resample(self, sampling_rate=None, variables=None, force_dense=False, - in_place=False, kind='linear'): + in_place=False, kind='linear', **kwargs): ''' Resample all dense variables (and optionally, sparse ones) to the specified sampling rate. @@ -276,7 +276,8 @@ def resample(self, sampling_rate=None, variables=None, force_dense=False, # None if in_place; no update needed _var = var.resample(sampling_rate, inplace=in_place, - kind=kind) + kind=kind, + **kwargs) if not in_place: _variables[name] = _var @@ -289,6 +290,7 @@ def resample(self, sampling_rate=None, variables=None, force_dense=False, def to_df(self, variables=None, format='wide', sparse=True, sampling_rate=None, include_sparse=True, include_dense=True, + integration_window=None, **kwargs): ''' Merge columns into a single pandas DataFrame. @@ -343,9 +345,10 @@ def to_df(self, variables=None, format='wide', sparse=True, sampling_rate = sampling_rate or self.sampling_rate # Make sure all variables have the same sampling rate - variables = list(self.resample(sampling_rate, variables, - force_dense=True, - in_place=False).values()) + variables = list( + self.resample(sampling_rate, variables, + force_dense=True, in_place=False, + integration_window=integration_window).values()) return super(BIDSRunVariableCollection, self).to_df(variables, format, **kwargs) diff --git a/bids/variables/tests/test_variables.py b/bids/variables/tests/test_variables.py index 84d220bfb..4cc3fd2b8 100644 --- a/bids/variables/tests/test_variables.py +++ b/bids/variables/tests/test_variables.py @@ -19,7 +19,7 @@ def generate_DEV(name='test', sr=20, duration=480): ent_names = ['task', 'run', 'session', 'subject'] entities = {e: uuid.uuid4().hex for e in ent_names} image = uuid.uuid4().hex + '.nii.gz' - run_info = RunInfo(entities, duration, 2, image) + run_info = RunInfo(entities, duration, 2, 2, image) return DenseRunVariable(name='test', values=values, run_info=run_info, source='dummy', sampling_rate=sr) diff --git a/bids/variables/variables.py b/bids/variables/variables.py index 42bb43392..66764b3cd 100644 --- a/bids/variables/variables.py +++ b/bids/variables/variables.py @@ -425,7 +425,7 @@ def _build_entity_index(self, run_info, sampling_rate): self.timestamps = pd.concat(_timestamps, axis=0, sort=True) return pd.concat(index, axis=0, sort=True).reset_index(drop=True) - def resample(self, sampling_rate, inplace=False, kind='linear'): + def resample(self, sampling_rate, integration_window=None, inplace=False, kind='linear'): '''Resample the Variable to the specified sampling rate. Parameters @@ -442,7 +442,7 @@ def resample(self, sampling_rate, inplace=False, kind='linear'): ''' if not inplace: var = self.clone() - var.resample(sampling_rate, True, kind) + var.resample(sampling_rate, integration_window, True, kind) return var if sampling_rate == self.sampling_rate: @@ -453,18 +453,36 @@ def resample(self, sampling_rate, inplace=False, kind='linear'): self.index = self._build_entity_index(self.run_info, sampling_rate) - x = np.arange(n) num = len(self.index) - from scipy.interpolate import interp1d - f = interp1d(x, self.values.values.ravel(), kind=kind) - x_new = np.linspace(0, n - 1, num=num) - self.values = pd.DataFrame(f(x_new)) + if integration_window is not None: + from scipy.sparse import lil_matrix + old_times = np.arange(n) / old_sr + new_times = np.arange(num) / sampling_rate + integrator = lil_matrix((num, n), dtype=np.uint8) + count = None + for i, new_time in enumerate(new_times): + cols = (old_times >= new_time) & (old_times < new_time + integration_window) + # This could be determined analytically, but dodging off-by-one errors + if count is None: + count = np.sum(cols) + integrator[i, cols] = 1 + + old_vals = self.values.values + self.values = pd.DataFrame(integrator.tocsr().dot(old_vals) / count) + else: + from scipy.interpolate import interp1d + x = np.arange(n) + f = interp1d(x, self.values.values.ravel(), kind=kind) + x_new = np.linspace(0, n - 1, num=num) + self.values = pd.DataFrame(f(x_new)) + assert len(self.values) == len(self.index) self.sampling_rate = sampling_rate - def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None): + def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None, + integration_window=None): '''Convert to a DataFrame, with columns for name and entities. Parameters @@ -480,7 +498,8 @@ def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None): sampled uniformly). If False, omits them. ''' if sampling_rate not in (None, self.sampling_rate): - return self.resample(sampling_rate).to_df(condition, entities) + return self.resample(sampling_rate, + integration_window=integration_window).to_df(condition, entities) df = super(DenseRunVariable, self).to_df(condition, entities)