Concatenating linear operators¶
It is possible to concatenate linear operators in a way that imitates the
numpy.concatenate
function. The only concatenations that are allowed are
in the vertical and horizontal directions.
The talon.concatenate
function requires an iterable containing the linear
operators to concatenate and the axis along which they have to be concatenated.
The following code shows the correct syntax to concatenate two linear operators \(A\) and \(B\) vertically and horizontally:
V = talon.concatenate((A, B), axis=0) # vertical (default)
H = talon.concatenate((A, B), axis=1) # horizontal
which correspond to the following
Examples¶
Build a tractogram with two crossing bundles of fibers and the corresponding linear operator.
import numpy as np
import talon
from scipy.interpolate import interp1d
# Set seed for reproducibility
np.random.seed(1992)
# The number of voxels in each dimension of the output image.
image_size = 25
center = image_size // 2
n_points = int(image_size / 0.01)
t = np.linspace(0, 1, n_points)
streamlines_per_bundle = 50
def generate_crossing_tractogram():
tractogram = []
horizontal_points = np.array([[0, center, center],
[image_size - 1, center, center]])
horizontal_streamline = interp1d([0, 1], horizontal_points, axis=0)(t)
for k in range(streamlines_per_bundle):
new_streamline = horizontal_streamline.copy()
new_streamline[:,1] += (np.random.rand(1) - 0.5)
tractogram.append(new_streamline)
vertical_points = np.array([[center, 0, center],
[center, image_size - 1, center]])
vertical_streamline = interp1d([0, 1], vertical_points, axis=0)(t)
for k in range(streamlines_per_bundle):
new_streamline = vertical_streamline.copy()
new_streamline[:,0] += (np.random.rand(1) - 0.5)
tractogram.append(new_streamline)
return tractogram
cross_tractogram = generate_crossing_tractogram()
directions = talon.utils.directions(1000)
generators = np.ones((len(directions), 1))
image_shape = (image_size,) * 3
indices, lengths = talon.voxelize(cross_tractogram, directions, image_shape)
A = talon.operator(generators, indices, lengths)
Vertical concatenation¶
If multiple features for each streamline are encoded in different linear operators we can concatenate different linear operators vertically. If \(A\) encodes the linear operator for the set of streamlines \(\alpha\) and generators \(G_1\) and \(B\) encodes the linear operator for the same streamlines but with generators \(G_2\), instead of rebuilding the linear operator from scratch we can concatenate \(A\) and \(B\) vertically to obtain the same result.
G2 = np.random.rand(len(directions), 5) # New generators
B = talon.operator(G2, indices, lengths)
V = talon.concatenate((A,B), axis=0)
print('Shape of A: {}'.format(A.shape))
print('Shape of B: {}'.format(B.shape))
print('Shape of V: {}'.format(V.shape))
print('Check: {} + {} = {}'.format(A.shape[0], B.shape[0], A.shape[0] + B.shape[0]))
Notice that the axis=0
argument is redundant since it is the default.
The output should be the following:
Shape of A: (15625, 100)
Shape of B: (78125, 100)
Shape of V: (93750, 100)
Check: 15625 + 78125 = 93750
Horizontal concatenation¶
One (but not the only) reason to concatenate two linear operators horizontally is to add a set of streamlines to the system. If \(A\) encodes the linear operator for the set of streamlines \(\alpha\) and \(C\) for set \(\beta\), instead of rebuilding the linear operator from scratch we can concatenate \(A\) and \(C\) horizontally to obtain the same result.
def generate_diagonal_tractogram():
tractogram = []
diagonal_points = np.array([[0, center, center],
[center, image_size - 1, center]])
diagonal_streamline = interp1d([0, 1], diagonal_points, axis=0)(t)
for k in range(streamlines_per_bundle):
new_streamline = diagonal_streamline.copy()
new_streamline[:,0] += (np.random.rand(1) - 0.5)
new_streamline[:,1] += (np.random.rand(1) - 0.5)
tractogram.append(new_streamline)
return tractogram
diag_tractogram = generate_diagonal_tractogram()
indices, lengths = talon.voxelize(diag_tractogram, directions, image_shape)
C = talon.operator(generators, indices, lengths) # diagonal
The concatenation of the two linear operators is performed as follows:
H = talon.concatenate([A, C], axis=1)
print('Shape of A: {}'.format(A.shape))
print('Shape of C: {}'.format(C.shape))
print('Shape of H: {}'.format(H.shape))
The output should be the following:
Shape of A: (15625, 100)
Shape of C: (15625, 50)
Shape of H: (15625, 150)
The matrix multiplication and transposition operations work as usual:
x = H @ np.random.rand(H.shape[1])
y = H.T @ np.random.rand(H.shape[0])
print('Shape of x: {}'.format(x.shape))
print('Shape of y: {}'.format(y.shape))
The output should be the following:
Shape of x: (15625,)
Shape of y: (150,)