Note
Go to the end to download the full example code. or to run this example in your browser via Binder
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
Explore 3D images (of cells).
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).
import matplotlib.pyplot as plt
import numpy as np
import plotly.express as px
import plotly.io
from skimage import data, feature
Load image#
This biomedical image is available through scikit-image’s data registry.
data = data.kidney()
What exactly are the shape and size of our 3D multichannel image?
print(f'number of dimensions: {data.ndim}')
print(f'shape: {data.shape}')
print(f'dtype: {data.dtype}')
number of dimensions: 4
shape: (16, 512, 512, 3)
dtype: uint16
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?
range: (68, 4095)
Let us visualize the middle slice of our 3D image.
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.
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([])
To view the sample data in 3D, run the following code:
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()
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.
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.
We can then compute the eigenvalues of the structure tensor.
(3, 8, 30, 30)
Where is the largest eigenvalue?
coords = np.unravel_index(eigen.argmax(), eigen.shape)
assert coords[0] == 0 # by definition
coords
(np.int64(0), np.int64(1), np.int64(22), np.int64(16))
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]).