估计 3D 显微镜图像中的各向异性

在本教程中，我们计算 3D 图像的结构张量。有关 3D 图像处理的概览，请参考 探索 3D 图像（细胞）。我们在这里使用的数据取样自共聚焦荧光显微镜获得的肾脏组织图像（有关详细信息，请参阅 [1]，位于 kidney-tissue-fluorescence.tif 下）。

[1]

scikit-image/data/#data

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

from skimage import data, feature

加载图像

此生物医学图像可通过 scikit-image 的数据注册表获取。

data = data.kidney()

我们的 3D 多通道图像的形状和大小到底是什么？

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

在本教程中，我们将只考虑第二个颜色通道，这将留下一个 3D 单通道图像。值的范围是多少？

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})')

range: (68, 4095)

让我们可视化 3D 图像的中间切片。

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)

让我们选择一个特定的区域，该区域显示相对 X-Y 各向同性。相反，梯度沿 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([])

要以 3D 方式查看示例数据，请运行以下代码

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()

计算结构张量

让我们可视化示例数据的底部切片并确定强变异的典型大小。我们将使用此大小作为窗口函数的“宽度”。

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)

关于最亮的区域（即，在行 ~ 22 和列 ~ 17 处），我们可以看到跨列（分别为行）2 或 3 个（分别为 1 或 2 个）像素的变异（因此，强梯度）。因此，我们可以选择，例如，sigma = 1.5 用于窗口函数。或者，我们可以按轴传递 sigma，例如，sigma = (1, 2, 3)。请注意，大小 1 沿第一个（Z，平面）轴听起来很合理，因为后者的尺寸为 8（13 - 5）。查看 X-Z 或 Y-Z 平面中的切片证实了这一点。

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

然后我们可以计算结构张量的特征值。

eigen = feature.structure_tensor_eigenvalues(A_elems)
eigen.shape

(3, 8, 30, 30)

最大的特征值在哪里？

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))

注意

读者可以检查此结果（坐标 (plane, row, column) = coords[1:]）对于不同的 sigma 的鲁棒性。

让我们查看特征值在找到最大特征值的 X-Y 平面的空间分布（即，Z = coords[1]）。

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)

我们正在查看一个局部属性。让我们考虑上述 X-Y 平面中最大特征值周围的一个微小邻域。

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

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

如果我们检查此邻域中的第二大特征值，我们可以看到它们的量级与最大特征值相同。

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

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

相反，第三大特征值小了一个数量级。

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

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

让我们可视化示例数据在找到最大特征值的 X-Y 平面的切片。

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)

让我们可视化示例数据在找到最大特征值的 X-Z（左）和 Y-Z（右）平面中的切片。Z 轴是下面子图中的垂直轴。我们可以看到沿 Z 轴（对应于肾脏组织中的纵向结构）的预期相对不变性，特别是在 Y-Z 平面 (longitudinal=1) 中。

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)

总之，关于体素 (plane, row, column) = coords[1:] 的区域在 3D 中是各向异性的：一方面，第三大特征值与最大特征值和第二大特征值之间存在一个数量级。我们可以在图 平面 Z = 1 的特征值 中一目了然地看到这一点。

所讨论的邻域在一个平面上是“有点各向同性”的（这里，将相对接近 X-Y 平面）：第二大特征值与最大特征值之间存在不到 2 的因子。此描述与我们在图像中看到的内容相符，即跨一个方向（这里，将相对接近行轴）的较强梯度和垂直于它的较弱梯度。

在 3D 结构张量的椭球表示 [2] 中，我们将得到薄饼情况。

[2]

https://en.wikipedia.org/wiki/Structure_tensor#Interpretation_2