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

\[\begin{split}V = \begin{bmatrix} A \\ B \end{bmatrix} \qquad H = \begin{bmatrix} A & B \end{bmatrix}.\end{split}\]

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,)