ENH: Add a Bandpass filter interface under algorithms.filters
#2915
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,109 @@ | ||
| # -*- coding: utf-8 -*- | ||
| # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- | ||
| # vi: set ft=python sts=4 ts=4 sw=4 et: | ||
| """ | ||
| Signal processing tools | ||
| """ | ||
| from __future__ import (print_function, division, unicode_literals, | ||
| absolute_import) | ||
|
|
||
| from ..utils.filemanip import fname_presuffix | ||
| from ..interfaces.base import (traits, TraitedSpec, SimpleInterface, | ||
| BaseInterfaceInputSpec, File) | ||
|
|
||
|
|
||
| class BandpassInputSpec(BaseInterfaceInputSpec): | ||
| in_file = File(exists=True, mandatory=True, desc='functional data') | ||
| freq_low = traits.Range(0.0, min=0.0, usedefault=True, | ||
| desc='low frequency cutoff (in Hz); ' | ||
| 'the default of 0 sets low pass cutoff to Nyquist') | ||
| freq_hi = traits.Range(0.0, min=0.0, usedefault=True, | ||
| desc='high frequency cutoff (in Hz); ' | ||
| 'the default of 0 sets high pass cutoff to 0') | ||
| repetition_time = traits.Either(None, traits.Range( | ||
| 0.0, min=0.0, exclude_low=True), desc='repetition_time') | ||
|
|
||
|
|
||
| class BandpassOutputSpec(TraitedSpec): | ||
| out_file = File(exists=True, desc='bandpass filtered functional data') | ||
|
|
||
|
|
||
| class Bandpass(SimpleInterface): | ||
| """ | ||
| Bandpass filtering for functional MRI timeseries | ||
| """ | ||
| input_spec = BandpassInputSpec | ||
| output_spec = BandpassOutputSpec | ||
|
|
||
| def _run_interface(self, runtime): | ||
| self._results['out_file'] = _bandpass_filter( | ||
| in_file=self.inputs.in_file, | ||
| tr=self.inputs.repetition_time, | ||
| freq_low=self.inputs.freq_low, | ||
| freq_hi=self.inputs.freq_hi, | ||
| out_file=fname_presuffix( | ||
| self.inputs.in_file, | ||
| suffix='_filtered', newpath=runtime.cwd), | ||
| ) | ||
| return runtime | ||
|
|
||
|
|
||
| def _bandpass_filter(in_file, tr=None, freq_low=0, freq_hi=0, out_file=None): | ||
| """ | ||
| Bandpass filter the input files | ||
|
|
||
| Parameters | ||
| ---------- | ||
| files : str | ||
| 4D NIfTI file | ||
| freq_low : float | ||
| cutoff frequency for the low pass filter (in Hz) | ||
| the default of 0 sets low pass cutoff to Nyquist | ||
| freq_hi : float | ||
| cutoff frequency for the high pass filter (in Hz) | ||
| the default of 0 sets high pass cutoff to 0 | ||
|
|
||
| tr : float | ||
| repetition time (in seconds) | ||
| out_file : str | ||
| output file name | ||
|
|
||
| """ | ||
| import numpy as np | ||
| import nibabel as nb | ||
|
|
||
| if freq_hi > 0 and freq_low >= freq_hi: | ||
| raise ValueError("Low-cutoff frequency can't be greater than the high-cutoff") | ||
|
|
||
| img = nb.load(in_file, mmap=NUMPY_MMAP) | ||
| timepoints = img.shape[-1] | ||
| F = np.zeros((timepoints)) | ||
|
|
||
| if tr is None: # If TR is not set, find in the image file header | ||
| tr = img.header.get_zooms()[3] | ||
|
|
||
| sampling_rate = 1. / tr | ||
|
|
||
| lowidx = timepoints // 2 + 1 # "/" replaced by "//" | ||
| if freq_low > 0: | ||
| # "np.round(..." replaced by "int(np.round(..." | ||
| lowidx = int(np.round(freq_low / sampling_rate * timepoints)) | ||
|
|
||
| highidx = 0 | ||
| if freq_hi > 0: | ||
| highidx = int(np.round(freq_hi / sampling_rate * timepoints)) # same | ||
|
|
||
| F[highidx:lowidx] = 1 | ||
| F = ((F + F[::-1]) > 0).astype(int) | ||
| try: | ||
| data = img.get_fdata() | ||
| except AttributeError: | ||
| data = img.get_data() | ||
|
|
||
| if np.all(F): | ||
| filtered_data = data | ||
| return in_file | ||
|
|
||
| filtered_data = np.real(np.fft.ifftn(np.fft.fftn(data) * F)) | ||
| img_out = nb.Nifti1Image(filtered_data, img.affine, img.header) | ||
| img_out.to_filename(out_file) | ||
| return out_file | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,30 @@ | ||
| # AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT | ||
| from __future__ import unicode_literals | ||
| from ..filters import Bandpass | ||
|
|
||
|
|
||
| def test_Bandpass_inputs(): | ||
| input_map = dict( | ||
| freq_hi=dict( | ||
| min=0.0, | ||
| usedefault=True, | ||
| ), | ||
| freq_low=dict( | ||
| min=0.0, | ||
| usedefault=True, | ||
| ), | ||
| in_file=dict(mandatory=True, ), | ||
| repetition_time=dict(), | ||
| ) | ||
| inputs = Bandpass.input_spec() | ||
|
|
||
| for key, metadata in list(input_map.items()): | ||
| for metakey, value in list(metadata.items()): | ||
| assert getattr(inputs.traits()[key], metakey) == value | ||
| def test_Bandpass_outputs(): | ||
| output_map = dict(out_file=dict(), ) | ||
| outputs = Bandpass.output_spec() | ||
|
|
||
| for key, metadata in list(output_map.items()): | ||
| for metakey, value in list(metadata.items()): | ||
| assert getattr(outputs.traits()[key], metakey) == value |
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Necessarily NIfTI?