Skip to content

ENH: Speed up tractogram.apply_affine #1000

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Feb 10, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 17 additions & 3 deletions nibabel/affines.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ class AffineError(ValueError):
pass


def apply_affine(aff, pts):
def apply_affine(aff, pts, inplace=False):
""" Apply affine matrix `aff` to points `pts`

Returns result of application of `aff` to the *right* of `pts`. The
Expand All @@ -39,6 +39,10 @@ def apply_affine(aff, pts):
pts : (..., N-1) array-like
Points, where the last dimension contains the coordinates of each
point. For 3D, the last dimension will be length 3.
inplace : bool, optional
If True, attempt to apply the affine directly to ``pts``.
If False, or in-place application fails, a freshly allocated
array will be returned.

Returns
-------
Expand Down Expand Up @@ -80,8 +84,18 @@ def apply_affine(aff, pts):
# rzs == rotations, zooms, shears
rzs = aff[:-1, :-1]
trans = aff[:-1, -1]
res = np.dot(pts, rzs.T) + trans[None, :]
return res.reshape(shape)

if inplace:
try:
np.dot(pts, rzs.T, out=pts)
except ValueError:
inplace = False
else:
pts += trans[None, :]
if not inplace:
pts = pts @ rzs.T + trans[None, :]

return pts.reshape(shape)


def to_matvec(transform):
Expand Down
4 changes: 4 additions & 0 deletions nibabel/streamlines/array_sequence.py
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,10 @@ def __init__(self, iterable=None, buffer_size=4):

self.extend(iterable)

@property
def is_sliced_view(self):
return self._lengths.sum() != self._data.shape[0]

@property
def is_array_sequence(self):
return True
Expand Down
8 changes: 6 additions & 2 deletions nibabel/streamlines/tractogram.py
Original file line number Diff line number Diff line change
Expand Up @@ -429,8 +429,12 @@ def apply_affine(self, affine, lazy=False):
if np.all(affine == np.eye(4)):
return self # No transformation.

for i in range(len(self.streamlines)):
self.streamlines[i] = apply_affine(affine, self.streamlines[i])
if self.streamlines.is_sliced_view:
# Apply affine only on the selected streamlines.
for i in range(len(self.streamlines)):
self.streamlines[i] = apply_affine(affine, self.streamlines[i])
else:
self.streamlines._data = apply_affine(affine, self.streamlines._data, inplace=True)

if self.affine_to_rasmm is not None:
# Update the affine that brings back the streamlines to RASmm.
Expand Down
8 changes: 7 additions & 1 deletion nibabel/tests/test_affines.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,12 @@ def test_apply_affine():
pts = pts.reshape((2, 2, 3))
exp_res = exp_res.reshape((2, 2, 3))
assert_array_equal(apply_affine(aff, pts), exp_res)

# Check inplace modification.
res = apply_affine(aff, pts, inplace=True)
assert_array_equal(res, exp_res)
assert np.shares_memory(res, pts)

# That ND also works
for N in range(2, 6):
aff = np.eye(N)
Expand Down Expand Up @@ -203,7 +209,7 @@ def test_rescale_affine():
orig_zooms = voxel_sizes(orig_aff)
orig_axcodes = aff2axcodes(orig_aff)
orig_centroid = apply_affine(orig_aff, (orig_shape - 1) // 2)

for new_shape in (None, tuple(orig_shape), (256, 256, 256), (64, 64, 40)):
for new_zooms in ((1, 1, 1), (2, 2, 3), (0.5, 0.5, 0.5)):
new_aff = rescale_affine(orig_aff, orig_shape, new_zooms, new_shape)
Expand Down