Skip to content

Commit ea62ca6

Browse files
authored
Merge pull request #2198 from effigies/rf/artdetect_norm
RF: Separate mcparam -> affine translation from norm calculation in ArtifactDetect
2 parents 9a988d6 + c74fec0 commit ea62ca6

File tree

1 file changed

+43
-35
lines changed

1 file changed

+43
-35
lines changed

nipype/algorithms/rapidart.py

Lines changed: 43 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -99,6 +99,29 @@ def _calc_norm(mc, use_differences, source, brain_pts=None):
9999
100100
"""
101101

102+
affines = [_get_affine_matrix(mc[i, :], source)
103+
for i in range(mc.shape[0])]
104+
return _calc_norm_affine(affines, use_differences, brain_pts)
105+
106+
107+
def _calc_norm_affine(affines, use_differences, brain_pts=None):
108+
"""Calculates the maximum overall displacement of the midpoints
109+
of the faces of a cube due to translation and rotation.
110+
111+
Parameters
112+
----------
113+
affines : list of [4 x 4] affine matrices
114+
use_differences : boolean
115+
brain_pts : [4 x n_points] of coordinates
116+
117+
Returns
118+
-------
119+
120+
norm : at each time point
121+
displacement : euclidean distance (mm) of displacement at each coordinate
122+
123+
"""
124+
102125
if brain_pts is None:
103126
respos = np.diag([70, 70, 75])
104127
resneg = np.diag([-70, -110, -45])
@@ -107,49 +130,34 @@ def _calc_norm(mc, use_differences, source, brain_pts=None):
107130
else:
108131
all_pts = brain_pts
109132
n_pts = all_pts.size - all_pts.shape[1]
110-
newpos = np.zeros((mc.shape[0], n_pts))
133+
newpos = np.zeros((len(affines), n_pts))
111134
if brain_pts is not None:
112-
displacement = np.zeros((mc.shape[0], int(n_pts / 3)))
113-
for i in range(mc.shape[0]):
114-
affine = _get_affine_matrix(mc[i, :], source)
115-
newpos[i, :] = np.dot(affine,
116-
all_pts)[0:3, :].ravel()
135+
displacement = np.zeros((len(affines), int(n_pts / 3)))
136+
for i, affine in enumerate(affines):
137+
newpos[i, :] = np.dot(affine, all_pts)[0:3, :].ravel()
117138
if brain_pts is not None:
118-
displacement[i, :] = \
119-
np.sqrt(np.sum(np.power(np.reshape(newpos[i, :],
120-
(3, all_pts.shape[1])) -
121-
all_pts[0:3, :],
122-
2),
123-
axis=0))
139+
displacement[i, :] = np.sqrt(np.sum(
140+
np.power(np.reshape(newpos[i, :],
141+
(3, all_pts.shape[1])) - all_pts[0:3, :],
142+
2),
143+
axis=0))
124144
# np.savez('displacement.npz', newpos=newpos, pts=all_pts)
125-
normdata = np.zeros(mc.shape[0])
145+
normdata = np.zeros(len(affines))
126146
if use_differences:
127147
newpos = np.concatenate((np.zeros((1, n_pts)),
128148
np.diff(newpos, n=1, axis=0)), axis=0)
129149
for i in range(newpos.shape[0]):
130150
normdata[i] = \
131-
np.max(np.sqrt(np.sum(np.reshape(np.power(np.abs(newpos[i, :]), 2),
132-
(3, all_pts.shape[1])), axis=0)))
151+
np.max(np.sqrt(np.sum(
152+
np.reshape(np.power(np.abs(newpos[i, :]), 2),
153+
(3, all_pts.shape[1])),
154+
axis=0)))
133155
else:
134156
newpos = np.abs(signal.detrend(newpos, axis=0, type='constant'))
135157
normdata = np.sqrt(np.mean(np.power(newpos, 2), axis=1))
136158
return normdata, displacement
137159

138160

139-
def _nanmean(a, axis=None):
140-
"""Return the mean excluding items that are nan
141-
142-
>>> a = [1, 2, np.nan]
143-
>>> _nanmean(a)
144-
1.5
145-
146-
"""
147-
if axis:
148-
return np.nansum(a, axis) / np.sum(1 - np.isnan(a), axis)
149-
else:
150-
return np.nansum(a) / np.sum(1 - np.isnan(a))
151-
152-
153161
class ArtifactDetectInputSpec(BaseInterfaceInputSpec):
154162
realigned_files = InputMultiPath(File(exists=True),
155163
desc="Names of realigned functional data files",
@@ -376,11 +384,11 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
376384
vol = data[:, :, :, t0]
377385
# Use an SPM like approach
378386
mask_tmp = vol > \
379-
(_nanmean(vol) / self.inputs.global_threshold)
387+
(np.nanmean(vol) / self.inputs.global_threshold)
380388
mask = mask * mask_tmp
381389
for t0 in range(timepoints):
382390
vol = data[:, :, :, t0]
383-
g[t0] = _nanmean(vol[mask])
391+
g[t0] = np.nanmean(vol[mask])
384392
if len(find_indices(mask)) < (np.prod((x, y, z)) / 10):
385393
intersect_mask = False
386394
g = np.zeros((timepoints, 1))
@@ -390,7 +398,7 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
390398
for t0 in range(timepoints):
391399
vol = data[:, :, :, t0]
392400
mask_tmp = vol > \
393-
(_nanmean(vol) / self.inputs.global_threshold)
401+
(np.nanmean(vol) / self.inputs.global_threshold)
394402
mask[:, :, :, t0] = mask_tmp
395403
g[t0] = np.nansum(vol * mask_tmp) / np.nansum(mask_tmp)
396404
elif masktype == 'file': # uses a mask image to determine intensity
@@ -400,15 +408,15 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
400408
mask = mask > 0.5
401409
for t0 in range(timepoints):
402410
vol = data[:, :, :, t0]
403-
g[t0] = _nanmean(vol[mask])
411+
g[t0] = np.nanmean(vol[mask])
404412
elif masktype == 'thresh': # uses a fixed signal threshold
405413
for t0 in range(timepoints):
406414
vol = data[:, :, :, t0]
407415
mask = vol > self.inputs.mask_threshold
408-
g[t0] = _nanmean(vol[mask])
416+
g[t0] = np.nanmean(vol[mask])
409417
else:
410418
mask = np.ones((x, y, z))
411-
g = _nanmean(data[mask > 0, :], 1)
419+
g = np.nanmean(data[mask > 0, :], 1)
412420

413421
# compute normalized intensity values
414422
gz = signal.detrend(g, axis=0) # detrend the signal

0 commit comments

Comments
 (0)