Skip to content

Commit d843794

Browse files
authored
Merge pull request #2107 from effigies/enh/hpf_compcor
ENH: Add cosine-basis high-pass-filter to CompCor
2 parents 24fae18 + e03cb13 commit d843794

File tree

5 files changed

+216
-39
lines changed

5 files changed

+216
-39
lines changed

nipype/algorithms/confounds.py

Lines changed: 196 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -323,21 +323,34 @@ class CompCorInputSpec(BaseInterfaceInputSpec):
323323
desc=('Position of mask in `mask_files` to use - '
324324
'first is the default.'))
325325
components_file = traits.Str('components_file.txt', usedefault=True,
326-
desc='Filename to store physiological components')
326+
desc='Filename to store physiological components')
327327
num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
328-
use_regress_poly = traits.Bool(True, usedefault=True,
328+
pre_filter = traits.Enum('polynomial', 'cosine', False, usedefault=True,
329+
desc='Detrend time series prior to component '
330+
'extraction')
331+
use_regress_poly = traits.Bool(True,
332+
deprecated='0.15.0', new_name='pre_filter',
329333
desc=('use polynomial regression '
330334
'pre-component extraction'))
331335
regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
332336
desc='the degree polynomial to use')
333337
header_prefix = traits.Str(desc=('the desired header for the output tsv '
334338
'file (one column). If undefined, will '
335339
'default to "CompCor"'))
340+
high_pass_cutoff = traits.Float(
341+
128, usedefault=True,
342+
desc='Cutoff (in seconds) for "cosine" pre-filter')
343+
repetition_time = traits.Float(
344+
desc='Repetition time (TR) of series - derived from image header if '
345+
'unspecified')
346+
save_pre_filter = traits.Either(
347+
traits.Bool, File, desc='Save pre-filter basis as text file')
336348

337349

338350
class CompCorOutputSpec(TraitedSpec):
339351
components_file = File(exists=True,
340352
desc='text file containing the noise components')
353+
pre_filter_file = File(desc='text file containing high-pass filter basis')
341354

342355

343356
class CompCor(BaseInterface):
@@ -351,7 +364,7 @@ class CompCor(BaseInterface):
351364
>>> ccinterface.inputs.realigned_file = 'functional.nii'
352365
>>> ccinterface.inputs.mask_files = 'mask.nii'
353366
>>> ccinterface.inputs.num_components = 1
354-
>>> ccinterface.inputs.use_regress_poly = True
367+
>>> ccinterface.inputs.pre_filter = 'polynomial'
355368
>>> ccinterface.inputs.regress_poly_degree = 2
356369
357370
"""
@@ -383,17 +396,20 @@ def _run_interface(self, runtime):
383396
self.inputs.merge_method,
384397
self.inputs.mask_index)
385398

399+
if self.inputs.use_regress_poly:
400+
self.inputs.pre_filter = 'polynomial'
401+
402+
# Degree 0 == remove mean; see compute_noise_components
386403
degree = (self.inputs.regress_poly_degree if
387-
self.inputs.use_regress_poly else 0)
404+
self.inputs.pre_filter == 'polynomial' else 0)
388405

389-
imgseries = nb.load(self.inputs.realigned_file,
390-
mmap=NUMPY_MMAP)
406+
imgseries = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP)
391407

392408
if len(imgseries.shape) != 4:
393-
raise ValueError('tCompCor expected a 4-D nifti file. Input {} has '
394-
'{} dimensions (shape {})'.format(
395-
self.inputs.realigned_file, len(imgseries.shape),
396-
imgseries.shape))
409+
raise ValueError('{} expected a 4-D nifti file. Input {} has '
410+
'{} dimensions (shape {})'.format(
411+
self._header, self.inputs.realigned_file,
412+
len(imgseries.shape), imgseries.shape))
397413

398414
if len(mask_images) == 0:
399415
img = nb.Nifti1Image(np.ones(imgseries.shape[:3], dtype=np.bool),
@@ -403,13 +419,41 @@ def _run_interface(self, runtime):
403419

404420
mask_images = self._process_masks(mask_images, imgseries.get_data())
405421

406-
components = compute_noise_components(imgseries.get_data(),
407-
mask_images, degree,
408-
self.inputs.num_components)
422+
TR = 0
423+
if self.inputs.pre_filter == 'cosine':
424+
if isdefined(self.inputs.repetition_time):
425+
TR = self.inputs.repetition_time
426+
else:
427+
# Derive TR from NIfTI header, if possible
428+
try:
429+
TR = imgseries.header.get_zooms()[3]
430+
if imgseries.get_xyzt_units()[1] == 'msec':
431+
TR /= 1000
432+
except (AttributeError, IndexError):
433+
TR = 0
434+
435+
if TR == 0:
436+
raise ValueError(
437+
'{} cannot detect repetition time from image - '
438+
'Set the repetition_time input'.format(self._header))
439+
440+
components, filter_basis = compute_noise_components(
441+
imgseries.get_data(), mask_images, self.inputs.num_components,
442+
self.inputs.pre_filter, degree, self.inputs.high_pass_cutoff, TR)
409443

410444
components_file = os.path.join(os.getcwd(), self.inputs.components_file)
411445
np.savetxt(components_file, components, fmt=b"%.10f", delimiter='\t',
412446
header=self._make_headers(components.shape[1]), comments='')
447+
448+
if self.inputs.pre_filter and self.inputs.save_pre_filter:
449+
pre_filter_file = self._list_outputs()['pre_filter_file']
450+
ftype = {'polynomial': 'poly',
451+
'cosine': 'cos'}[self.inputs.pre_filter]
452+
ncols = filter_basis.shape[1] if filter_basis.size > 0 else 0
453+
header = ['{}{:02d}'.format(ftype, i) for i in range(ncols)]
454+
np.savetxt(pre_filter_file, filter_basis, fmt=b'%.10f',
455+
delimiter='\t', header='\t'.join(header), comments='')
456+
413457
return runtime
414458

415459
def _process_masks(self, mask_images, timeseries=None):
@@ -418,14 +462,19 @@ def _process_masks(self, mask_images, timeseries=None):
418462
def _list_outputs(self):
419463
outputs = self._outputs().get()
420464
outputs['components_file'] = os.path.abspath(self.inputs.components_file)
465+
466+
save_pre_filter = self.inputs.save_pre_filter
467+
if save_pre_filter:
468+
if isinstance(save_pre_filter, bool):
469+
save_pre_filter = os.path.abspath('pre_filter.tsv')
470+
outputs['pre_filter_file'] = save_pre_filter
471+
421472
return outputs
422473

423474
def _make_headers(self, num_col):
424-
headers = []
425475
header = self.inputs.header_prefix if \
426476
isdefined(self.inputs.header_prefix) else self._header
427-
for i in range(num_col):
428-
headers.append(header + '{:02d}'.format(i))
477+
headers = ['{}{:02d}'.format(header, i) for i in range(num_col)]
429478
return '\t'.join(headers)
430479

431480

@@ -473,7 +522,7 @@ class TCompCor(CompCor):
473522
>>> ccinterface.inputs.realigned_file = 'functional.nii'
474523
>>> ccinterface.inputs.mask_files = 'mask.nii'
475524
>>> ccinterface.inputs.num_components = 1
476-
>>> ccinterface.inputs.use_regress_poly = True
525+
>>> ccinterface.inputs.pre_filter = 'polynomial'
477526
>>> ccinterface.inputs.regress_poly_degree = 2
478527
>>> ccinterface.inputs.percentile_threshold = .03
479528
@@ -494,7 +543,7 @@ def _process_masks(self, mask_images, timeseries=None):
494543
for i, img in enumerate(mask_images):
495544
mask = img.get_data().astype(np.bool)
496545
imgseries = timeseries[mask, :]
497-
imgseries = regress_poly(2, imgseries)
546+
imgseries = regress_poly(2, imgseries)[0]
498547
tSTD = _compute_tSTD(imgseries, 0, axis=-1)
499548
threshold_std = np.percentile(tSTD, np.round(100. *
500549
(1. - self.inputs.percentile_threshold)).astype(int))
@@ -569,7 +618,7 @@ def _run_interface(self, runtime):
569618
data = data.astype(np.float32)
570619

571620
if isdefined(self.inputs.regress_poly):
572-
data = regress_poly(self.inputs.regress_poly, data, remove_mean=False)
621+
data = regress_poly(self.inputs.regress_poly, data, remove_mean=False)[0]
573622
img = nb.Nifti1Image(data, img.affine, header)
574623
nb.save(img, op.abspath(self.inputs.detrended_file))
575624

@@ -618,7 +667,7 @@ def _run_interface(self, runtime):
618667
global_signal = in_nii.get_data()[:,:,:,:50].mean(axis=0).mean(axis=0).mean(axis=0)
619668

620669
self._results = {
621-
'n_volumes_to_discard': _is_outlier(global_signal)
670+
'n_volumes_to_discard': is_outlier(global_signal)
622671
}
623672

624673
return runtime
@@ -685,9 +734,10 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False,
685734
func_sd = func_sd[func_sd != 0]
686735

687736
# Compute (non-robust) estimate of lag-1 autocorrelation
688-
ar1 = np.apply_along_axis(AR_est_YW, 1,
689-
regress_poly(0, mfunc, remove_mean=True).astype(
690-
np.float32), 1)[:, 0]
737+
ar1 = np.apply_along_axis(
738+
AR_est_YW, 1,
739+
regress_poly(0, mfunc, remove_mean=True)[0].astype(np.float32),
740+
1)[:, 0]
691741

692742
# Compute (predicted) standard deviation of temporal difference time series
693743
diff_sdhat = np.squeeze(np.sqrt(((1 - ar1) * 2).tolist())) * func_sd
@@ -794,6 +844,27 @@ def is_outlier(points, thresh=3.5):
794844
return timepoints_to_discard
795845

796846

847+
def cosine_filter(data, timestep, period_cut, remove_mean=True, axis=-1):
848+
datashape = data.shape
849+
timepoints = datashape[axis]
850+
851+
data = data.reshape((-1, timepoints))
852+
853+
frametimes = timestep * np.arange(timepoints)
854+
X = _full_rank(_cosine_drift(period_cut, frametimes))[0]
855+
non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
856+
857+
betas = np.linalg.lstsq(X, data.T)[0]
858+
859+
if not remove_mean:
860+
X = X[:, :-1]
861+
betas = betas[:-1]
862+
863+
residuals = data - X.dot(betas).T
864+
865+
return residuals.reshape(datashape), non_constant_regressors
866+
867+
797868
def regress_poly(degree, data, remove_mean=True, axis=-1):
798869
"""
799870
Returns data with degree polynomial regressed out.
@@ -817,6 +888,8 @@ def regress_poly(degree, data, remove_mean=True, axis=-1):
817888
value_array = np.linspace(-1, 1, timepoints)
818889
X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))
819890

891+
non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
892+
820893
# Calculate coefficients
821894
betas = np.linalg.pinv(X).dot(data.T)
822895

@@ -828,7 +901,7 @@ def regress_poly(degree, data, remove_mean=True, axis=-1):
828901
regressed_data = data - datahat
829902

830903
# Back to original shape
831-
return regressed_data.reshape(datashape)
904+
return regressed_data.reshape(datashape), non_constant_regressors
832905

833906

834907
def combine_mask_files(mask_files, mask_method=None, mask_index=None):
@@ -886,37 +959,57 @@ def combine_mask_files(mask_files, mask_method=None, mask_index=None):
886959
return [img]
887960

888961

889-
def compute_noise_components(imgseries, mask_images, degree, num_components):
962+
def compute_noise_components(imgseries, mask_images, num_components,
963+
filter_type, degree, period_cut,
964+
repetition_time):
890965
"""Compute the noise components from the imgseries for each mask
891966
892967
imgseries: a nibabel img
893968
mask_images: a list of nibabel images
894-
degree: order of polynomial used to remove trends from the timeseries
895969
num_components: number of noise components to return
970+
filter_type: type off filter to apply to time series before computing
971+
noise components.
972+
'polynomial' - Legendre polynomial basis
973+
'cosine' - Discrete cosine (DCT) basis
974+
False - None (mean-removal only)
975+
976+
Filter options:
977+
978+
degree: order of polynomial used to remove trends from the timeseries
979+
period_cut: minimum period (in sec) for DCT high-pass filter
980+
repetition_time: time (in sec) between volume acquisitions
896981
897982
returns:
898983
899984
components: a numpy array
985+
basis: a numpy array containing the (non-constant) filter regressors
900986
901987
"""
902988
components = None
989+
basis = np.array([])
903990
for img in mask_images:
904991
mask = img.get_data().astype(np.bool)
905992
if imgseries.shape[:3] != mask.shape:
906-
raise ValueError('Inputs for CompCor, timeseries and mask, '
907-
'do not have matching spatial dimensions '
908-
'({} and {}, respectively)'.format(
909-
imgseries.shape[:3], mask.shape))
993+
raise ValueError(
994+
'Inputs for CompCor, timeseries and mask, do not have '
995+
'matching spatial dimensions ({} and {}, respectively)'.format(
996+
imgseries.shape[:3], mask.shape))
910997

911998
voxel_timecourses = imgseries[mask, :]
912999

9131000
# Zero-out any bad values
9141001
voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0
9151002

916-
# from paper:
917-
# "The constant and linear trends of the columns in the matrix M were
918-
# removed [prior to ...]"
919-
voxel_timecourses = regress_poly(degree, voxel_timecourses)
1003+
# Currently support Legendre-polynomial or cosine or detrending
1004+
# With no filter, the mean is nonetheless removed (poly w/ degree 0)
1005+
if filter_type == 'cosine':
1006+
voxel_timecourses, basis = cosine_filter(
1007+
voxel_timecourses, repetition_time, period_cut)
1008+
elif filter_type in ('polynomial', False):
1009+
# from paper:
1010+
# "The constant and linear trends of the columns in the matrix M were
1011+
# removed [prior to ...]"
1012+
voxel_timecourses, basis = regress_poly(degree, voxel_timecourses)
9201013

9211014
# "Voxel time series from the noise ROI (either anatomical or tSTD) were
9221015
# placed in a matrix M of size Nxm, with time along the row dimension
@@ -936,7 +1029,7 @@ def compute_noise_components(imgseries, mask_images, degree, num_components):
9361029
u[:, :num_components]))
9371030
if components is None and num_components > 0:
9381031
raise ValueError('No components found')
939-
return components
1032+
return components, basis
9401033

9411034

9421035
def _compute_tSTD(M, x, axis=0):
@@ -945,3 +1038,71 @@ def _compute_tSTD(M, x, axis=0):
9451038
stdM[stdM == 0] = x
9461039
stdM[np.isnan(stdM)] = x
9471040
return stdM
1041+
1042+
1043+
# _cosine_drift and _full_rank copied from nipy/modalities/fmri/design_matrix
1044+
#
1045+
# Nipy release: 0.4.1
1046+
# Modified for smooth integration in CompCor classes
1047+
1048+
def _cosine_drift(period_cut, frametimes):
1049+
"""Create a cosine drift matrix with periods greater or equals to period_cut
1050+
1051+
Parameters
1052+
----------
1053+
period_cut: float
1054+
Cut period of the low-pass filter (in sec)
1055+
frametimes: array of shape(nscans)
1056+
The sampling times (in sec)
1057+
1058+
Returns
1059+
-------
1060+
cdrift: array of shape(n_scans, n_drifts)
1061+
cosin drifts plus a constant regressor at cdrift[:,0]
1062+
1063+
Ref: http://en.wikipedia.org/wiki/Discrete_cosine_transform DCT-II
1064+
"""
1065+
len_tim = len(frametimes)
1066+
n_times = np.arange(len_tim)
1067+
hfcut = 1. / period_cut # input parameter is the period
1068+
1069+
# frametimes.max() should be (len_tim-1)*dt
1070+
dt = frametimes[1] - frametimes[0]
1071+
# hfcut = 1/(2*dt) yields len_time
1072+
# If series is too short, return constant regressor
1073+
order = max(int(np.floor(2*len_tim*hfcut*dt)), 1)
1074+
cdrift = np.zeros((len_tim, order))
1075+
nfct = np.sqrt(2.0/len_tim)
1076+
1077+
for k in range(1, order):
1078+
cdrift[:, k-1] = nfct * np.cos((np.pi / len_tim) * (n_times + .5) * k)
1079+
1080+
cdrift[:, order-1] = 1. # or 1./sqrt(len_tim) to normalize
1081+
return cdrift
1082+
1083+
1084+
def _full_rank(X, cmax=1e15):
1085+
"""
1086+
This function possibly adds a scalar matrix to X
1087+
to guarantee that the condition number is smaller than a given threshold.
1088+
1089+
Parameters
1090+
----------
1091+
X: array of shape(nrows, ncols)
1092+
cmax=1.e-15, float tolerance for condition number
1093+
1094+
Returns
1095+
-------
1096+
X: array of shape(nrows, ncols) after regularization
1097+
cmax=1.e-15, float tolerance for condition number
1098+
"""
1099+
U, s, V = np.linalg.svd(X, 0)
1100+
smax, smin = s.max(), s.min()
1101+
c = smax / smin
1102+
if c < cmax:
1103+
return X, c
1104+
IFLOG.warn('Matrix is singular at working precision, regularizing...')
1105+
lda = (smax - cmax * smin) / (cmax - 1)
1106+
s = s + lda
1107+
X = np.dot(U, np.dot(np.diag(s), V))
1108+
return X, cmax

0 commit comments

Comments
 (0)