.. 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-25 .. code-block:: Python 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 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 = 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 = 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 = 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.273 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.24.0?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 ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_3d_structure_tensor.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_