This tutorial will give you an impression of how to process and manipulate your neurons' morphology. See the :ref:`API reference <api_morph>` for a complete list of available functions.

As you might imagine some manipulations (e.g. smoothing or simplification) will work on all/most neuron types while others will only work on specific types. For example rerooting only makes sense on a ``TreeNeuron``.

The rule of thumb is this: if a function is called e.g. ``simplify_neuron`` it should work with multiple if not all neuron types while specialized functions will be called e.g. ``reroot_skeleton``. So depending on what data you are working with and what you want to get out of it, you might have to explicitly convert between neuron types. See the respective function's docstring for details!

Rerooting
=========

``TreeNeurons`` are hierarchical trees and as such typically have a single "root" node (fragmented neurons will have multiple roots). The root is important because it is used as the reference/origin for a bunch of analyses such as Strahler order. Typically, you want the root to be the soma. Because the root is so important, :class:`TreeNeuron` can be rerooted:

.. code:: ipython3

    import navis
    
    n = navis.example_neurons(1, kind='skeleton')
    print(n.soma)


.. parsed-literal::

    [4177]


``.soma`` returns the node ID of the soma (if there is one) and can be used to reroot

.. code:: ipython3

    navis.reroot_skeleton(n, n.soma)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>type</th>
          <td>navis.TreeNeuron</td>
        </tr>
        <tr>
          <th>name</th>
          <td>DA1_lPN_R</td>
        </tr>
        <tr>
          <th>id</th>
          <td>1734350788</td>
        </tr>
        <tr>
          <th>n_nodes</th>
          <td>4465</td>
        </tr>
        <tr>
          <th>n_connectors</th>
          <td>2705</td>
        </tr>
        <tr>
          <th>n_branches</th>
          <td>598</td>
        </tr>
        <tr>
          <th>n_leafs</th>
          <td>619</td>
        </tr>
        <tr>
          <th>cable_length</th>
          <td>266476.875</td>
        </tr>
        <tr>
          <th>soma</th>
          <td>[4177]</td>
        </tr>
        <tr>
          <th>units</th>
          <td>8 nanometer</td>
        </tr>
        <tr>
          <th>created_at</th>
          <td>2022-07-27 21:16:48.240027</td>
        </tr>
        <tr>
          <th>id</th>
          <td>1734350788</td>
        </tr>
        <tr>
          <th>origin</th>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
      </tbody>
    </table>
    </div>



Note that the root is implicitly also important for :class:`~navis.MeshNeuron` because we're using their skeletons for a couple operations/analyses. 

Simplifying
===========

If you work with large lists of neurons you may want to downsample/simplifiy before e.g. trying to plot them. This is one of the things that - in principle work - with all neuron types. The implementation, however, depends on the neuron type. Lookup the respective function's help (e.g. via the :ref:`API <api>`) for details. 

For ``TreeNeurons`` downsampling means skipping N nodes (here 10):

.. code:: ipython3

    sk = navis.example_neurons(n=1, kind='skeleton')
    print(sk.n_nodes)


.. parsed-literal::

    4465


.. code:: ipython3

    sk_downsampled = navis.downsample_neuron(sk, downsampling_factor=10, inplace=False)
    print(sk_downsampled.n_nodes)


.. parsed-literal::

    1304


For ``MeshNeurons`` downsampling will reduce the number of faces by a factor of N:

.. code:: ipython3

    me = navis.example_neurons(n=1, kind='mesh')
    print(me.n_faces)


.. parsed-literal::

    13054


.. code:: ipython3

    me_downsampled = navis.downsample_neuron(me, downsampling_factor=10, inplace=False)
    print(me_downsampled.n_faces)


.. parsed-literal::

    1474


.. note::
   Under the hood :func:`~navis.downsample_neuron` calls :func:`~navis.simplify_mesh` for ``MeshNeurons``.
   That function then requires one of the supported backends for mesh operations to be installed: Blender
   3D, ``pymeshlab`` or ``open3d``.

``TreeNeurons`` can also be "resampled" (up or down) to a given resolution (i.e. distance between nodes):

.. code:: ipython3

    sk = navis.example_neurons(n=1, kind='skeleton')
    print(sk.sampling_resolution * sk.units)


.. parsed-literal::

    477.4501679731243 nanometer


.. code:: ipython3

    # Note that we can provide a unit ("1 micron") here because our neuron has units set
    sk_resampled = navis.resample_skeleton(sk, resample_to='1 micron', inplace=False)
    print(sk_resampled.sampling_resolution * sk_resampled.units)


.. parsed-literal::

    1072.2341923485653 nanometer


Let's visualize what we did there:

.. code:: ipython3

    import matplotlib.pyplot as plt
    
    nodes_original = sk.nodes[['x', 'y' ,'z']].values 
    nodes_downsampled = sk_downsampled.nodes[['x', 'y' ,'z']].values 
    nodes_resampled = sk_resampled.nodes[['x', 'y' ,'z']].values 
    
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    _ = navis.plot2d(nodes_original, method='2d', view=('x', '-z'), scatter_kws=dict(c='blue'), ax=axes[0])
    _ = navis.plot2d(nodes_resampled, method='2d', view=('x', '-z'), scatter_kws=dict(c='red'), ax=axes[1])
    _ = navis.plot2d(nodes_downsampled, method='2d', view=('x', '-z'), scatter_kws=dict(c='green'), ax=axes[2])
    
    for ax, title in zip(axes, ['Original', 'Resampled to 1um', 'Downsampled 10x']):
        ax.set_title(title, color='lightgrey')
        ax.set_axis_off()
    
    plt.show()



.. image:: morph_processing_files/morph_processing_16_0.png


As you can see the resampling increased the node density in the backbone and decreased it in the finer neurites to bring things on par. Downsampling just thinned out the nodes across the board.

.. note::

   Resampling has a caveat you need to be aware of: nodes are not merely moved around to match the 
   desired resolution - they are regenerated from scratch. As consequence, the original node IDs
   are - with a few exceptions - all gone.

Smoothing
=========

Smoothing is one of those things that work on all neurons but the approaches are so vastly different that there are separate functions: :func:`navis.smooth_skeleton`, :func:`navis.smooth_mesh` and :func:`navis.smooth_voxels`:

.. code:: ipython3

    # smooth_skeleton uses a rolling window along the linear segments
    sk = navis.example_neurons(n=1, kind='skeleton')
    sk_smoothed = navis.smooth_skeleton(sk, window=5, inplace=False)






.. code:: ipython3

    # smooth_mesh uses a iterative rounds of Laplacian smoothing
    me = navis.example_neurons(n=1, kind='mesh')
    me_smoothed = navis.smooth_mesh(me, iterations=5, inplace=False)

Cutting & Pruning
=================

Cutting and pruning work best if there is a sense of topology which implicitly requires a skeleton. Many functions will also work on ``MeshNeurons`` though. That's because the operation is performed on their skeleton and changes are propagated back to the mesh. Fair warning though: this may not be perfect (e.g. the resulting mesh might not be watertight) - should be good enough for a first pass though.

Let's start with something easy: cutting a skeleton in two at a given node.

.. code:: ipython3

    # Load the neuron
    n = navis.example_neurons(1, kind='skeleton')
    
    # Pick a node ID
    cut_node_id = n.nodes.node_id.values[333]
    distal, proximal = navis.cut_skeleton(n, cut_node_id)

Plot the two fragments:

.. code:: ipython3

    # Note that we are using method='2d' here because that makes annotating the plot easier
    fig, ax = distal.plot2d(color='cyan', method='2d', view=('x', '-z'))
    fig, ax = proximal.plot2d(color='green', ax=ax, method='2d', view=('x', '-z'))
    
    # Annotate cut point
    cut_coords = distal.nodes.set_index('node_id').loc[distal.root, ['x', 'z']].values[0]
    ax.annotate('cut point',
                xy=(cut_coords[0], -cut_coords[1]),
                color='lightgrey',
                xytext=(cut_coords[0], -cut_coords[1]-2000), va='center', ha='center',
                arrowprops=dict(shrink=0.1, width=2, color='lightgrey'),
                )
    
    plt.show()



.. image:: morph_processing_files/morph_processing_23_0.png


If instead of a node ID, you have an x/y/z coordinate where you want to cut: use the ``.snap`` method to find the closest node to that location.

.. code:: ipython3

    node_id, dist = n.snap([14000, 16200, 12000])
    print(f'Closest node: {node_id} at distance {dist * n.units:.2f} {n.units.units}')


.. parsed-literal::

    Closest node: 334 at distance 565.69 nanometer nanometer


Instead of cutting a neuron in two, we can also just prune bits off:

.. code:: ipython3

    n_pruned = n.prune_distal_to(cut_node_id, inplace=False)
    
    cut_coords = n.nodes.set_index('node_id').loc[cut_node_id, ['x', 'z']].values
    
    # Plot original neuron in red and with dotted line
    fig, ax = n.plot2d(color='red', method='2d', linestyle=(0, (5, 10)), view=('x', '-z'))
    
    # Plot remaining neurites in red
    fig, ax = n_pruned.plot2d(color='green', method='2d', ax=ax, view=('x', '-z'), lw=1.2)
    
    # Annotate cut point
    ax.annotate('cut point',
                xy=(cut_coords[0], -cut_coords[1]),
                color='lightgrey',
                xytext=(cut_coords[0], -cut_coords[1]-2000), va='center', ha='center',
                arrowprops=dict(shrink=0.1, width=2, color='lightgrey'),
                )
    
    plt.show()



.. image:: morph_processing_files/morph_processing_27_0.png


:func:`~navis.cut_skeleton` also takes multiple cut nodes, in case you want to chop your neuron into multiple pieces.

As an (extreme) example, let's cut a neuron at every single branch point:

.. code:: ipython3

    n = navis.example_neurons(1, kind='skeleton')
    
    branch_points = n.nodes[n.nodes.type=='branch'].node_id.values
    
    cut = navis.cut_skeleton(n, branch_points)
    cut.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>type</th>
          <th>name</th>
          <th>id</th>
          <th>n_nodes</th>
          <th>n_connectors</th>
          <th>n_branches</th>
          <th>n_leafs</th>
          <th>cable_length</th>
          <th>soma</th>
          <th>units</th>
          <th>created_at</th>
          <th>id</th>
          <th>origin</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>navis.TreeNeuron</td>
          <td>DA1_lPN_R</td>
          <td>1734350788</td>
          <td>4</td>
          <td>16</td>
          <td>0</td>
          <td>2</td>
          <td>373.565735</td>
          <td>None</td>
          <td>8 nanometer</td>
          <td>2022-07-27 21:19:57.584201</td>
          <td>1734350788</td>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
        <tr>
          <th>1</th>
          <td>navis.TreeNeuron</td>
          <td>DA1_lPN_R</td>
          <td>1734350788</td>
          <td>5</td>
          <td>10</td>
          <td>0</td>
          <td>2</td>
          <td>431.518494</td>
          <td>None</td>
          <td>8 nanometer</td>
          <td>2022-07-27 21:19:57.584201</td>
          <td>1734350788</td>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
        <tr>
          <th>2</th>
          <td>navis.TreeNeuron</td>
          <td>DA1_lPN_R</td>
          <td>1734350788</td>
          <td>6</td>
          <td>12</td>
          <td>0</td>
          <td>2</td>
          <td>388.590637</td>
          <td>None</td>
          <td>8 nanometer</td>
          <td>2022-07-27 21:19:57.584201</td>
          <td>1734350788</td>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
        <tr>
          <th>3</th>
          <td>navis.TreeNeuron</td>
          <td>DA1_lPN_R</td>
          <td>1734350788</td>
          <td>8</td>
          <td>22</td>
          <td>0</td>
          <td>2</td>
          <td>665.534912</td>
          <td>None</td>
          <td>8 nanometer</td>
          <td>2022-07-27 21:19:57.584201</td>
          <td>1734350788</td>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
        <tr>
          <th>4</th>
          <td>navis.TreeNeuron</td>
          <td>DA1_lPN_R</td>
          <td>1734350788</td>
          <td>28</td>
          <td>8</td>
          <td>0</td>
          <td>2</td>
          <td>1534.385498</td>
          <td>None</td>
          <td>8 nanometer</td>
          <td>2022-07-27 21:19:57.584201</td>
          <td>1734350788</td>
          <td>/Users/philipps/Google Drive/Cloudbox/Github/n...</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # Plot neuron fragments
    fig, ax = navis.plot2d(cut, linewidth=1.5)
    
    # Rotate to front view
    ax.azim, ax.elev = -90, -90
    ax.dist = 6
    
    plt.show()







.. image:: morph_processing_files/morph_processing_30_1.png


Let's try something more sophisticated: pruning by Strahler index.

.. code:: ipython3

    # Load a fresh skeleton
    n = navis.example_neurons(1, kind='skeleton')
    
    # Reroot to soma
    n = n.reroot(n.soma)
    
    # This will prune off terminal branches (the lowest two Strahler indices)
    n_pruned = n.prune_by_strahler(to_prune = [1, 2], inplace=False)
    
    # Plot original neurons in red
    fig, ax = n.plot2d(color='red')
    
    # Plot remaining neurites in green
    fig, ax = n_pruned.plot2d(color='green', ax=ax, linewidth=1)
    
    # Rotate to front view
    ax.azim, ax.elev = -90, -90
    ax.dist = 6
    
    plt.show()



.. image:: morph_processing_files/morph_processing_32_0.png


We can also turn this around and remove only the higher order branches. Let's use this example to show that we can also do this with ``MeshNeurons``:

.. code:: ipython3

    # Load an example mesh neuron
    m = navis.example_neurons(1, kind='mesh')
    
    # This will prune to the just terminal branches
    m_pruned = navis.prune_by_strahler(m, to_prune=range(3, 100), inplace=False)
    
    # Plot original neuron in cyan
    fig, ax = m.plot2d(color='cyan', figsize=(10, 10))
    
    # Plot remaining neurites red
    fig, ax = m_pruned.plot2d(color='red', ax=ax)
    
    # Rotate to front view
    ax.azim, ax.elev = -90, -90
    ax.dist = 6
    
    plt.show()



.. image:: morph_processing_files/morph_processing_34_0.png


Alternatively, we can prune terminal branches based on size:

.. code:: ipython3

    # This will prune all branches smaller than 10 microns
    m_pruned = navis.prune_twigs(m, size='10 microns', inplace=False)
    
    # Plot original neuron in red
    fig, ax = m.plot2d(color='red', figsize=(10, 10))
    
    # Plot remaining neurites in cyan
    fig, ax = m_pruned.plot2d(color='cyan', ax=ax, linewidth=.75, alpha=.5)
    
    # Rotate to front view
    ax.azim, ax.elev = -90, -90
    ax.dist = 6
    
    plt.show()



.. image:: morph_processing_files/morph_processing_36_0.png


Intersecting with Volumes
=========================

We can also intersect neurons with :class:`navis.Volume` (and ``trimesh.Trimesh`` for that matter).

.. code:: ipython3

    # Load an example navis.Volume
    lh = navis.example_volume('LH')
    
    # Prune by volume
    m_lh = navis.in_volume(m, lh, inplace=False)
    m_outside_lh = navis.in_volume(m, lh, mode='OUT', inplace=False)

And plot!

.. code:: ipython3

    # Plot pruned branchs neuron in green
    fig, ax = navis.plot2d([m_lh, m_outside_lh, lh], color=['red', 'green'], figsize=(10, 10))
    
    # Rotate to front view
    ax.azim, ax.elev = -90, -90
    ax.dist = 6
    
    plt.show()







.. image:: morph_processing_files/morph_processing_40_1.png


`Does this work with all neuron types`? There is no simple answer unfortunately. In theory, anything that works on skeletons should also work on meshes, and vice versa. However, :class:`~navis.Dotprops` and :class:`~navis.VoxelNeuron` are so fundamentally different that certain operations just don't make sense. For example we can't cut them but we can subset them to a given volume. Check out the :ref:`API <api_func_matrix>` docs for an overview of what works with which neuron type.

Note that :func:`navis.in_volume` also works with arbitrary spatial data (i.e. `(N, 3)` arrays of x/y/z locations):

.. code:: ipython3

    # Get the connectors for one of our above skeletons
    cn = sk.connectors
    
    # Add a column that tells us which connectors are in the LH volume
    cn['in_lh'] = navis.in_volume(cn[['x', 'y', 'z']].values, lh)
    cn.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>connector_id</th>
          <th>node_id</th>
          <th>type</th>
          <th>x</th>
          <th>y</th>
          <th>z</th>
          <th>roi</th>
          <th>confidence</th>
          <th>in_lh</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0</td>
          <td>1436</td>
          <td>pre</td>
          <td>6444</td>
          <td>21608</td>
          <td>14516</td>
          <td>LH(R)</td>
          <td>0.959</td>
          <td>True</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1</td>
          <td>1436</td>
          <td>pre</td>
          <td>6457</td>
          <td>21634</td>
          <td>14474</td>
          <td>LH(R)</td>
          <td>0.997</td>
          <td>True</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2</td>
          <td>2638</td>
          <td>pre</td>
          <td>4728</td>
          <td>23538</td>
          <td>14179</td>
          <td>LH(R)</td>
          <td>0.886</td>
          <td>True</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3</td>
          <td>1441</td>
          <td>pre</td>
          <td>5296</td>
          <td>22059</td>
          <td>16048</td>
          <td>LH(R)</td>
          <td>0.967</td>
          <td>True</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4</td>
          <td>1872</td>
          <td>pre</td>
          <td>4838</td>
          <td>23114</td>
          <td>15146</td>
          <td>LH(R)</td>
          <td>0.990</td>
          <td>True</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # Count the number of connectors (pre and post) in- and outside the LH:
    cn.groupby(['type', 'in_lh']).size()




.. parsed-literal::

    type  in_lh
    post  False    1978
          True      106
    pre   False     325
          True      296
    dtype: int64



About half the presynapses are in the LH (most of the rest will be in the MB calyx). The large majority of postsynapses are outside the LH in the antennal lobe where this neuron has its dendrites.

That's it for now! Please see the :ref:`NBLAST tutorial <nblast_intro>` for morphological comparisons using NBLAST and the :ref:`API<api_morph>` for a full list of morphology-related functions.