.. 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 ` 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-27 .. code-block:: default import matplotlib.pyplot as plt import numpy as np import plotly.express as px import plotly.io from skimage import ( data, feature ) .. GENERATED FROM PYTHON SOURCE LINES 28-31 Load image ========== This biomedical image is available through `scikit-image`'s data registry. .. GENERATED FROM PYTHON SOURCE LINES 31-34 .. code-block:: default data = data.kidney() .. GENERATED FROM PYTHON SOURCE LINES 35-36 What exactly are the shape and size of our 3D multichannel image? .. GENERATED FROM PYTHON SOURCE LINES 36-41 .. code-block:: default 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 42-45 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 45-50 .. code-block:: default 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 51-52 Let us visualize the middle slice of our 3D image. .. GENERATED FROM PYTHON SOURCE LINES 52-62 .. code-block:: default 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 63-66 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 66-78 .. code-block:: default 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 79-99 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 101-106 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 106-117 .. code-block:: default 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 118-126 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 126-130 .. code-block:: default sigma = (1, 1.5, 2.5) A_elems = feature.structure_tensor(sample, sigma=sigma) .. GENERATED FROM PYTHON SOURCE LINES 131-132 We can then compute the eigenvalues of the structure tensor. .. GENERATED FROM PYTHON SOURCE LINES 132-136 .. code-block:: default eigen = 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 137-138 Where is the largest eigenvalue? .. GENERATED FROM PYTHON SOURCE LINES 138-143 .. code-block:: default coords = np.unravel_index(eigen.argmax(), eigen.shape) assert coords[0] == 0 # by definition coords .. rst-class:: sphx-glr-script-out .. code-block:: none (0, 1, 22, 16) .. GENERATED FROM PYTHON SOURCE LINES 144-150 .. 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 150-160 .. code-block:: default 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 161-163 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 163-166 .. code-block:: default 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 167-169 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 169-172 .. code-block:: default 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 173-175 In contrast, the third-largest eigenvalues are one order of magnitude smaller. .. GENERATED FROM PYTHON SOURCE LINES 175-178 .. code-block:: default 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 179-181 Let us visualize the slice of sample data in the X-Y plane where the maximum eigenvalue is found. .. GENERATED FROM PYTHON SOURCE LINES 181-192 .. code-block:: default 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 193-198 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 198-210 .. code-block:: default 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 211-217 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 219-225 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 227-231 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 2.021 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.21.x?filepath=notebooks/auto_examples/applications/plot_3d_structure_tensor.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_3d_structure_tensor.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_3d_structure_tensor.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_