.. _getting-started:

===============
Getting started
===============

The ``talon`` package, at its core, provides a way to transform a tractogram
into a linear operator, or more precisely a matrix.
This matrix can be used in two ways: to generate data and to explain data.
In both cases, the type of the data is arbitrary and is specified by the user,
not by ``talon``. To quickly get you started, the following examples illustrate
both use cases.

If you haven't already, start by installing ``talon``.
See the :ref:`installation` section.

This short introduction is separated into 4 parts:

* :ref:`creating-test-tractogram`
* :ref:`building-linear-operator`
* :ref:`generating-data`
* :ref:`explaining-data`


.. _creating-test-tractogram:

Creating a test tractogram
==========================

To generate data using ``talon``, we need a tractogram.
Here we will generate streamlines organised into a + sign.

.. code:: python

    import numpy as np
    from scipy.interpolate import interp1d

    # The number of voxels in each dimension of the output image.
    image_size = 25

    center = image_size // 2
    t = np.linspace(0, 1, int(image_size / 0.1))

    # Generate the horizontal and vertical streamlines.
    horizontal_points = np.array([[0, center, center], [image_size - 1, center, center]])
    horizontal_streamline = interp1d([0, 1], horizontal_points, axis=0)(t)

    vertical_points = np.array([[center, 0, center], [center, image_size - 1, center]])
    vertical_streamline = interp1d([0, 1], vertical_points, axis=0)(t)

    # A tractogram is just a collection of streamlines.
    tractogram = [horizontal_streamline, vertical_streamline]

To visualize the geometry of the streamlines, you can display them using ``matplotlib``.

.. code:: python

    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    fig = plt.figure(figsize=(5, 5), dpi=150)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(tractogram[0][:,0], tractogram[0][:,1], tractogram[0][:,2], 'k')
    ax.plot(tractogram[1][:,0], tractogram[1][:,1], tractogram[1][:,2], 'k')
    ax.view_init(90,90)
    ax.set_zticks([])
    plt.show()

You should see the following image:

.. image:: img/getting-started-streamlines.png
    :align: center
    :width: 50em
    :height: 50em


.. _building-linear-operator:

Building the linear operator
============================

Now that we have a tractogram, we can start using ``talon``.
First, we will *voxelize* the tractogram by separating each streamline into
voxel elements. If you are familiar with tractography, streamlines are generated
by following peaks of an image. Voxelizing a tractogram is the opposite i.e.
creating peaks from streamlines. In order to voxelize the tractogram, we first
need to provide a list of *directions* of the possible orientations of the
streamlines represented as an array of unit vectors.

.. code:: python

    import talon

    directions = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=np.float)
    image_shape = (image_size,) * 3
    indices, lengths = talon.voxelize(tractogram, directions, image_shape)

Next we define how each streamline direction is projected onto the data.

.. code:: python

    generators = np.ones((len(directions), 1))

Finally, we build the linear operator :math:`A`.

.. code:: python

    A = talon.operator(generators, indices, lengths)

Note that generators can be multidimensional.
One way to illustrate this is to use the directions as generators.

.. code:: python

    G = talon.operator(directions, indices, lengths)


.. _generating-data:

Generating data with a linear operator
======================================

To generate data simply multiply (using the ``@`` operator) the linear operator
by a weight vector.

.. code:: python

    # Using a vector off all ones gives all streamlines equal weight.
    x = np.ones(A.shape[1])
    b = A @ x

    # We can do the same thing with the multidimensional operator.
    m = G @ x

The data vector ``b`` can be reshaped into an image and visualized.

.. code:: python

    image = b.reshape(image_shape)

    plt.figure(figsize=(5, 5), dpi=150)
    plt.imshow(image[:, :, center])
    plt.colorbar(shrink=0.8)
    plt.show()

An we obtain the following image which corresponds to the streamline density.

.. image:: img/getting-started-density.png
    :align: center
    :width: 50em
    :height: 40em

The second data vector can also be visualized, but requires a bit more
manipulation.

.. code:: python

    rgb_image = m.reshape(image_shape + (3,))

    plt.figure(figsize=(5, 5), dpi=150)
    plt.imshow(rgb_image[:, :, center])
    plt.show()

.. image:: img/getting-started-rgb.png
    :align: center
    :width: 40em
    :height: 40em


.. _explaining-data:

Explaining data with a linear operator
======================================

Considering the case where an error in the tractography algorithm generates a
spurious streamline in our tractogram. In the case of our example, we simply
add a diagonal streamline to *tractogram*.

.. code:: python

    diagonal_points = np.array([[0, center, center], [center, image_size - 1, center]])
    diagonal_streamline = interp1d([0, 1], diagonal_points, axis=0)(t)

    tractogram.append(diagonal_streamline)

    # Visualize the new tractogram.
    fig = plt.figure(figsize=(5, 5), dpi=150)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(tractogram[0][:,0], tractogram[0][:,1], tractogram[0][:,2], 'k')
    ax.plot(tractogram[1][:,0], tractogram[1][:,1], tractogram[1][:,2], 'k')
    ax.plot(tractogram[2][:,0], tractogram[2][:,1], tractogram[1][:,2], 'k')
    ax.view_init(90,90)
    ax.set_zticks([])
    plt.show()

.. image:: img/getting-started-spurious.png
    :align: center
    :width: 50em
    :height: 50em

Given *b*, the data generated using by the original tractogram, we can use
``talon`` to calculate the contribution of each streamline to the data. In order
to do so, we first have to generate a *linear operator* using the new
tractogram. In this case, we use also use a set of 1000 equally spaced unit
vectors as *directions*.

.. code:: python

    directions = talon.utils.directions(1000)
    generators = np.ones((len(directions), 1))
    indices, lengths = talon.voxelize(tractogram, directions, image_shape)
    Z = talon.operator(generators, indices, lengths)

What we want to find are the streamline contributions *x* which minimize

.. math::

    \frac{1}{2} || Z  x - b||^2 + \Omega(x)

In this example it does not make sense to have streamlines with a negative
contribution, therefore, :math:`\Omega(x)` will be set as a positivity
constraint. In ``talon``, we can force positivity constraint  using the
``talon.regularization`` function.

.. code:: python

    positivity_constraint = talon.regularization(non_negativity=True)

The resulting regularization term is then given to the ``talon.solve`` function
in order to obtain the streamlines contributions.

.. code:: python

    solution = talon.solve(Z, b, reg_term=positivity_constraint)
    print('solution.x = [%.2f, %.2f, %.2f]' % tuple(solution.x))

.. code::

    solution.x = [1.00, 1.00, 0.00]

As it is possible to see, the two original streamlines contribute equally to
the data while the third streamline does not contribute.

We can use the ``talon`` solution to filter the tractogram and visualize only
the streamlines presenting a non-zero contribution.

.. code:: python

    # New filtered tractogram.
    filtered_tractogram = []

    fig = plt.figure(figsize=(5, 5), dpi=150)
    ax = fig.add_subplot(111, projection='3d')

    for i,s in enumerate(tractogram):

        # If the current streamline contributes to the data.
        if solution.x[i] > 0.0:

            # Add streamline to filtered tractogram.
            filtered_tractogram.append(s)

            # Visualize the streamline.
            ax.plot(s[:,0], s[:,1], s[:,2], 'k')

    ax.view_init(90,90)
    ax.set_zticks([])
    plt.show()

.. image:: img/getting-started-filtered.png
    :align: center
    :width: 50em
    :height: 50em