.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/applications/plot_3d_structure_tensor.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_applications_plot_3d_structure_tensor.py>`
        to download the full example code. or to run this example in your browser via Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_applications_plot_3d_structure_tensor.py:


============================================
Estimate anisotropy in a 3D microscopy image
============================================

In this tutorial, we compute the structure tensor of a 3D image.
For a general introduction to 3D image processing, please refer to
:ref:`sphx_glr_auto_examples_applications_plot_3d_image_processing.py`.
The data we use here are sampled from an image of kidney tissue obtained by
confocal fluorescence microscopy (more details at [1]_ under
``kidney-tissue-fluorescence.tif``).

.. [1] https://gitlab.com/scikit-image/data/#data

.. GENERATED FROM PYTHON SOURCE LINES 16-25

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    import plotly.express as px
    import plotly.io

    import skimage as ski









.. GENERATED FROM PYTHON SOURCE LINES 26-29

Load image
==========
This biomedical image is available through `scikit-image`'s data registry.

.. GENERATED FROM PYTHON SOURCE LINES 29-32

.. code-block:: Python


    data = ski.data.kidney()








.. GENERATED FROM PYTHON SOURCE LINES 33-34

What exactly are the shape and size of our 3D multichannel image?

.. GENERATED FROM PYTHON SOURCE LINES 34-39

.. code-block:: Python


    print(f'number of dimensions: {data.ndim}')
    print(f'shape: {data.shape}')
    print(f'dtype: {data.dtype}')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    number of dimensions: 4
    shape: (16, 512, 512, 3)
    dtype: uint16




.. GENERATED FROM PYTHON SOURCE LINES 40-43

For the purposes of this tutorial, we shall consider only the second color
channel, which leaves us with a 3D single-channel image. What is the range
of values?

.. GENERATED FROM PYTHON SOURCE LINES 43-48

.. code-block:: Python


    n_plane, n_row, n_col, n_chan = data.shape
    v_min, v_max = data[:, :, :, 1].min(), data[:, :, :, 1].max()
    print(f'range: ({v_min}, {v_max})')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    range: (68, 4095)




.. GENERATED FROM PYTHON SOURCE LINES 49-50

Let us visualize the middle slice of our 3D image.

.. GENERATED FROM PYTHON SOURCE LINES 50-60

.. code-block:: Python


    fig1 = px.imshow(
        data[n_plane // 2, :, :, 1],
        zmin=v_min,
        zmax=v_max,
        labels={'x': 'col', 'y': 'row', 'color': 'intensity'},
    )

    plotly.io.show(fig1)




.. raw:: html
    :file: images/sphx_glr_plot_3d_structure_tensor_001.html





.. GENERATED FROM PYTHON SOURCE LINES 61-64

Let us pick a specific region, which shows relative X-Y isotropy. In
contrast, the gradient is quite different (and, for that matter, weak) along
Z.

.. GENERATED FROM PYTHON SOURCE LINES 64-76

.. code-block:: Python


    sample = data[5:13, 380:410, 370:400, 1]
    step = 3
    cols = sample.shape[0] // step + 1
    _, axes = plt.subplots(nrows=1, ncols=cols, figsize=(16, 8))

    for it, (ax, image) in enumerate(zip(axes.flatten(), sample[::step])):
        ax.imshow(image, cmap='gray', vmin=v_min, vmax=v_max)
        ax.set_title(f'Plane = {5 + it * step}')
        ax.set_xticks([])
        ax.set_yticks([])




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_3d_structure_tensor_002.png
   :alt: Plane = 5, Plane = 8, Plane = 11
   :srcset: /auto_examples/applications/images/sphx_glr_plot_3d_structure_tensor_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 77-97

To view the sample data in 3D, run the following code:

.. code-block:: python

    import plotly.graph_objects as go

    (n_Z, n_Y, n_X) = sample.shape
    Z, Y, X = np.mgrid[:n_Z, :n_Y, :n_X]

    fig = go.Figure(
        data=go.Volume(
            x=X.flatten(),
            y=Y.flatten(),
            z=Z.flatten(),
            value=sample.flatten(),
            opacity=0.5,
            slices_z=dict(show=True, locations=[4])
        )
    )
    fig.show()

.. GENERATED FROM PYTHON SOURCE LINES 99-104

Compute structure tensor
========================
Let us visualize the bottom slice of our sample data and determine the
typical size for strong variations. We shall use this size as the
'width' of the window function.

.. GENERATED FROM PYTHON SOURCE LINES 104-115

.. code-block:: Python


    fig2 = px.imshow(
        sample[0, :, :],
        zmin=v_min,
        zmax=v_max,
        labels={'x': 'col', 'y': 'row', 'color': 'intensity'},
        title='Interactive view of bottom slice of sample data.',
    )

    plotly.io.show(fig2)




.. raw:: html
    :file: images/sphx_glr_plot_3d_structure_tensor_003.html





.. GENERATED FROM PYTHON SOURCE LINES 116-124

About the brightest region (i.e., at row ~ 22 and column ~ 17), we can see
variations (and, hence, strong gradients) over 2 or 3 (resp. 1 or 2) pixels
across columns (resp. rows). We may thus choose, say, ``sigma = 1.5`` for
the window
function. Alternatively, we can pass sigma on a per-axis basis, e.g.,
``sigma = (1, 2, 3)``. Note that size 1 sounds reasonable along the first
(Z, plane) axis, since the latter is of size 8 (13 - 5). Viewing slices in
the X-Z or Y-Z planes confirms it is reasonable.

.. GENERATED FROM PYTHON SOURCE LINES 124-128

.. code-block:: Python


    sigma = (1, 1.5, 2.5)
    A_elems = ski.feature.structure_tensor(sample, sigma=sigma)








.. GENERATED FROM PYTHON SOURCE LINES 129-130

We can then compute the eigenvalues of the structure tensor.

.. GENERATED FROM PYTHON SOURCE LINES 130-134

.. code-block:: Python


    eigen = ski.feature.structure_tensor_eigenvalues(A_elems)
    eigen.shape





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    (3, 8, 30, 30)



.. GENERATED FROM PYTHON SOURCE LINES 135-136

Where is the largest eigenvalue?

.. GENERATED FROM PYTHON SOURCE LINES 136-141

.. code-block:: Python


    coords = np.unravel_index(eigen.argmax(), eigen.shape)
    assert coords[0] == 0  # by definition
    coords





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    (np.int64(0), np.int64(1), np.int64(22), np.int64(16))



.. GENERATED FROM PYTHON SOURCE LINES 142-148

.. note::
   The reader may check how robust this result (coordinates
   ``(plane, row, column) = coords[1:]``) is to varying ``sigma``.

Let us view the spatial distribution of the eigenvalues in the X-Y plane
where the maximum eigenvalue is found (i.e., ``Z = coords[1]``).

.. GENERATED FROM PYTHON SOURCE LINES 148-158

.. code-block:: Python


    fig3 = px.imshow(
        eigen[:, coords[1], :, :],
        facet_col=0,
        labels={'x': 'col', 'y': 'row', 'facet_col': 'rank'},
        title=f'Eigenvalues for plane Z = {coords[1]}.',
    )

    plotly.io.show(fig3)




.. raw:: html
    :file: images/sphx_glr_plot_3d_structure_tensor_004.html





.. GENERATED FROM PYTHON SOURCE LINES 159-161

We are looking at a local property. Let us consider a tiny neighborhood
around the maximum eigenvalue in the above X-Y plane.

.. GENERATED FROM PYTHON SOURCE LINES 161-164

.. code-block:: Python


    eigen[0, coords[1], coords[2] - 2 : coords[2] + 1, coords[3] - 2 : coords[3] + 1]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([[0.05530323, 0.05929082, 0.06043806],
           [0.05922725, 0.06268274, 0.06354238],
           [0.06190861, 0.06685075, 0.06840962]])



.. GENERATED FROM PYTHON SOURCE LINES 165-167

If we examine the second-largest eigenvalues in this neighborhood, we can
see that they have the same order of magnitude as the largest ones.

.. GENERATED FROM PYTHON SOURCE LINES 167-170

.. code-block:: Python


    eigen[1, coords[1], coords[2] - 2 : coords[2] + 1, coords[3] - 2 : coords[3] + 1]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([[0.03577746, 0.03577334, 0.03447714],
           [0.03819524, 0.04172036, 0.04323701],
           [0.03139592, 0.03587025, 0.03913327]])



.. GENERATED FROM PYTHON SOURCE LINES 171-173

In contrast, the third-largest eigenvalues are one order of magnitude
smaller.

.. GENERATED FROM PYTHON SOURCE LINES 173-176

.. code-block:: Python


    eigen[2, coords[1], coords[2] - 2 : coords[2] + 1, coords[3] - 2 : coords[3] + 1]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([[0.00337661, 0.00306529, 0.00276288],
           [0.0041869 , 0.00397519, 0.00375595],
           [0.00479742, 0.00462116, 0.00442455]])



.. GENERATED FROM PYTHON SOURCE LINES 177-179

Let us visualize the slice of sample data in the X-Y plane where the
maximum eigenvalue is found.

.. GENERATED FROM PYTHON SOURCE LINES 179-190

.. code-block:: Python


    fig4 = px.imshow(
        sample[coords[1], :, :],
        zmin=v_min,
        zmax=v_max,
        labels={'x': 'col', 'y': 'row', 'color': 'intensity'},
        title=f'Interactive view of plane Z = {coords[1]}.',
    )

    plotly.io.show(fig4)




.. raw:: html
    :file: images/sphx_glr_plot_3d_structure_tensor_005.html





.. GENERATED FROM PYTHON SOURCE LINES 191-196

Let us visualize the slices of sample data in the X-Z (left) and Y-Z (right)
planes where the maximum eigenvalue is found. The Z axis is the vertical
axis in the subplots below. We can see the expected relative invariance
along the Z axis (corresponding to longitudinal structures in the kidney
tissue), especially in the Y-Z plane (``longitudinal=1``).

.. GENERATED FROM PYTHON SOURCE LINES 196-204

.. code-block:: Python


    subplots = np.dstack((sample[:, coords[2], :], sample[:, :, coords[3]]))
    fig5 = px.imshow(
        subplots, zmin=v_min, zmax=v_max, facet_col=2, labels={'facet_col': 'longitudinal'}
    )

    plotly.io.show(fig5)




.. raw:: html
    :file: images/sphx_glr_plot_3d_structure_tensor_006.html





.. GENERATED FROM PYTHON SOURCE LINES 205-211

As a conclusion, the region about voxel
``(plane, row, column) = coords[1:]`` is
anisotropic in 3D: There is an order of magnitude between the third-largest
eigenvalues on one hand, and the largest and second-largest eigenvalues on
the other hand. We could see this at first glance in figure `Eigenvalues for
plane Z = 1`.

.. GENERATED FROM PYTHON SOURCE LINES 213-219

The neighborhood in question is 'somewhat isotropic' in a plane (which,
here, would be relatively close to the X-Y plane): There is a factor of
less than 2 between the second-largest and largest eigenvalues.
This description is compatible with what we are seeing in the image, i.e., a
stronger gradient across a direction (which, here, would be relatively close
to the row axis) and a weaker gradient perpendicular to it.

.. GENERATED FROM PYTHON SOURCE LINES 221-225

In an ellipsoidal representation of the 3D structure tensor [2]_,
we would get the pancake situation.

.. [2] https://en.wikipedia.org/wiki/Structure_tensor#Interpretation_2


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.170 seconds)


.. _sphx_glr_download_auto_examples_applications_plot_3d_structure_tensor.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/scikit-image/scikit-image/v0.25.2?filepath=notebooks/auto_examples/applications/plot_3d_structure_tensor.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_3d_structure_tensor.ipynb <plot_3d_structure_tensor.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_3d_structure_tensor.py <plot_3d_structure_tensor.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_3d_structure_tensor.zip <plot_3d_structure_tensor.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_