X 射线图像处理#

本教程演示了如何使用 NumPy、imageio、Matplotlib 和 SciPy 读取和处理 X 射线图像。您将学习如何加载医学图像、关注特定部位以及使用 高斯拉普拉斯-高斯索贝尔Canny 滤波器进行边缘检测,以直观地比较它们。

例如,当您正在构建一个算法来帮助 检测肺炎 作为 Kaggle 比赛 的一部分时,X 射线图像分析可以成为您数据分析和 机器学习工作流程 的一部分。在医疗保健行业,当图像估计占所有医疗数据的 至少 90% 时,医学图像处理和分析尤为重要。

您将使用来自 ChestX-ray8 数据集的放射学图像,该数据集由 美国国立卫生研究院 (NIH) 提供。ChestX-ray8 包含来自 30,000 多名患者的 100,000 多张去身份识别后的 X 射线图像,格式为 PNG。您可以在 NIH 的公共 Box 存储库/images 文件夹中找到 ChestX-ray8 的文件。(有关更多详细信息,请参阅 2017 年在 CVPR(计算机视觉会议)上发表的 研究论文。)

为了方便起见,一小部分 PNG 图像已保存到本教程的存储库中的 tutorial-x-ray-image-processing/ 下,因为 ChestX-ray8 包含千兆字节的数据,您可能会发现以批次下载它很困难。

A series of 9 x-ray images of the same region of a patient's chest is shown with different types of image processing filters applied to each image. Each x-ray shows different types of biological detail.

先决条件#

读者应该了解一些 Python、NumPy 数组和 Matplotlib 的知识。为了复习,您可以学习 Python 和 Matplotlib PyPlot 教程,以及 NumPy 快速入门

本教程中使用了以下包

  • imageio 用于读取和写入图像数据。医疗保健行业通常使用 DICOM 格式进行医学成像,imageio 非常适合读取这种格式。为了简单起见,在本教程中,您将使用 PNG 文件。

  • Matplotlib 用于数据可视化。

  • SciPy 用于通过 ndimage 进行多维图像处理。

本教程可以在隔离的环境(如 Virtualenvconda)中本地运行。您可以使用 Jupyter Notebook 或 JupyterLab 运行每个笔记本单元格。

目录#

  1. 使用 imageio 检查 X 射线

  2. 将图像组合成多维数组以演示进展

  3. 使用拉普拉斯-高斯、高斯梯度、索贝尔和 Canny 滤波器进行边缘检测

  4. 使用 np.where() 将掩码应用于 X 射线

  5. 比较结果


使用 imageio 检查 X 射线#

让我们从一个简单的例子开始,只使用 ChestX-ray8 数据集中的一个 X 射线图像。

文件 — 00000011_001.png — 已经为您下载并保存在 /tutorial-x-ray-image-processing 文件夹中。

1. 使用 imageio 加载图像

import os
import imageio

DIR = "tutorial-x-ray-image-processing"

xray_image = imageio.v3.imread(os.path.join(DIR, "00000011_001.png"))

2. 检查其形状是否为 1024x1024 像素,并且数组是否由 8 位整数组成

print(xray_image.shape)
print(xray_image.dtype)
(1024, 1024)
uint8

3. 导入 matplotlib 并在灰度颜色图中显示图像

import matplotlib.pyplot as plt

plt.imshow(xray_image, cmap="gray")
plt.axis("off")
plt.show()
../_images/06294911bd707f233c372476eaf58785976e5e34ff65c441b84ecfc4d79c7fe3.png

将图像组合成多维数组以演示进展#

在下一个示例中,您将使用从 ChestX-ray8 数据集中下载并从其中一个数据集文件中提取的 9 张 X 射线 1024x1024 像素图像,而不是 1 张图像。它们的编号从 ...000.png...008.png,我们假设它们属于同一患者。

1. 导入 NumPy,读取每个 X 射线,并创建一个三维数组,其中第一个维度对应于图像编号

import numpy as np
num_imgs = 9

combined_xray_images_1 = np.array(
    [imageio.v3.imread(os.path.join(DIR, f"00000011_00{i}.png")) for i in range(num_imgs)]
)

2. 检查包含 9 个堆叠图像的新 X 射线图像数组的形状

combined_xray_images_1.shape
(9, 1024, 1024)

请注意,第一个维度的形状与 num_imgs 相匹配,因此 combined_xray_images_1 数组可以解释为 2D 图像堆栈。

3. 您现在可以使用 Matplotlib 将每个帧并排绘制来显示“健康进展”

fig, axes = plt.subplots(nrows=1, ncols=num_imgs, figsize=(30, 30))

for img, ax in zip(combined_xray_images_1, axes):
    ax.imshow(img, cmap='gray')
    ax.axis('off')
../_images/4763adcb6df4df7555699dedfab40762077abe0817d3a0d4d502db15b8e998d9.png

4. 此外,将进展显示为动画可能会有所帮助。让我们使用 imageio.mimwrite() 创建一个 GIF 文件,并在笔记本中显示结果

GIF_PATH = os.path.join(DIR, "xray_image.gif")
imageio.mimwrite(GIF_PATH, combined_xray_images_1, format= ".gif", duration=1000)

这将为我们提供:一个动画 gif 反复循环遍历一系列 8 张 X 射线,显示患者胸部的同一视角,只是在不同的时间点。可以直观地比较患者的骨骼和内部器官,从帧到帧。

使用拉普拉斯-高斯、高斯梯度、索贝尔和 Canny 滤波器进行边缘检测#

在处理生物医学数据时,强调 2D “边缘” 以关注图像中的特定特征可能很有用。为此,使用 图像梯度 在检测颜色像素强度变化时可能特别有用。

具有高斯二阶导数的拉普拉斯滤波器#

让我们从一个使用 拉普拉斯 滤波器(“拉普拉斯-高斯”)的 n 维 高斯 二阶导数开始。这种拉普拉斯方法侧重于具有快速强度值变化的像素,并与高斯平滑相结合以 去除噪声。让我们检查它如何在分析 2D X 射线图像时如何有用。

  • 拉普拉斯-高斯滤波器的实现比较简单:1)从 SciPy 导入 ndimage 模块;以及 2)使用 sigma(标量)参数调用 scipy.ndimage.gaussian_laplace(),这会影响高斯滤波器的标准差(您将在下面的示例中使用 1

from scipy import ndimage

xray_image_laplace_gaussian = ndimage.gaussian_laplace(xray_image, sigma=1)

显示原始 X 射线和使用拉普拉斯-高斯滤波器的 X 射线

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))

axes[0].set_title("Original")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Laplacian-Gaussian (edges)")
axes[1].imshow(xray_image_laplace_gaussian, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()
../_images/044727346378eab96e6f59f655cdc24432ae6b01178ada0cd011cf4c1871679c.png

高斯梯度幅度方法#

另一种可能很有用的边缘检测方法是 高斯(梯度)滤波器。它使用高斯导数计算多维梯度幅度,并通过消除 高频 图像成分来提供帮助。

1. 使用 sigma(标量)参数(用于标准差;您将在下面的示例中使用 2)调用 scipy.ndimage.gaussian_gradient_magnitude()

x_ray_image_gaussian_gradient = ndimage.gaussian_gradient_magnitude(xray_image, sigma=2)

2. 显示原始 X 射线和使用高斯梯度滤波器的 X 射线

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))

axes[0].set_title("Original")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Gaussian gradient (edges)")
axes[1].imshow(x_ray_image_gaussian_gradient, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()
../_images/21254b0b3c3b826503acec7f0e3585c0cc6921b40593b110b563895920f13ef2.png

索贝尔-费尔德曼算子(索贝尔滤波器)#

为了在 2D X 光图像的水平和垂直轴上找到高空间频率区域(边缘或边缘图),可以使用 索贝尔-费尔德曼算子(索贝尔滤波器) 技术。索贝尔滤波器通过 卷积 将两个 3x3 核矩阵(每个轴一个)应用于 X 光图像。然后,使用 勾股定理 将这两个点(梯度)组合起来,以产生梯度幅度。

**1.** 使用索贝尔滤波器(scipy.ndimage.sobel())在 X 光图像的 X 轴和 Y 轴上进行操作。然后,使用 勾股定理 和 NumPy 的 np.hypot() 计算应用了索贝尔滤波器后的 xy 之间的距离,以获得幅度。最后,对重新缩放的图像进行归一化,以使像素值在 0 到 255 之间。

图像归一化 遵循 output_channel = 255.0 * (input_channel - min_value) / (max_value - min_value) 公式。因为你正在使用灰度图像,所以只需要归一化一个通道。

x_sobel = ndimage.sobel(xray_image, axis=0)
y_sobel = ndimage.sobel(xray_image, axis=1)

xray_image_sobel = np.hypot(x_sobel, y_sobel)

xray_image_sobel *= 255.0 / np.max(xray_image_sobel)

**2.** 将新的图像数组数据类型从 float16 更改为 32 位浮点数格式,以 使其与 Matplotlib 兼容

print("The data type - before: ", xray_image_sobel.dtype)

xray_image_sobel = xray_image_sobel.astype("float32")

print("The data type - after: ", xray_image_sobel.dtype)
The data type - before:  float16
The data type - after:  float32

**3.** 显示原始 X 光图像和应用了索贝尔“边缘”滤波器的图像。注意,使用灰度和 CMRmap 色图都有助于突出显示边缘。

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 15))

axes[0].set_title("Original")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Sobel (edges) - grayscale")
axes[1].imshow(xray_image_sobel, cmap="gray")
axes[2].set_title("Sobel (edges) - CMRmap")
axes[2].imshow(xray_image_sobel, cmap="CMRmap")
for i in axes:
    i.axis("off")
plt.show()
../_images/e2a53ba9faedba6c37fe48498469bb17319befb5ed7659a7a0e81887a65f7204.png

Canny 滤波器#

你也可以考虑使用另一个著名的边缘检测滤波器,称为 Canny 滤波器

首先,应用一个 高斯 滤波器来去除图像中的噪声。在本例中,你使用的是 傅里叶 滤波器,它通过 卷积 过程对 X 光图像进行平滑处理。接下来,在图像的两个轴上应用 Prewitt 滤波器 来帮助检测一些边缘 - 这将产生两个梯度值。与索贝尔滤波器类似,Prewitt 算子也通过 卷积 将两个 3x3 核矩阵(每个轴一个)应用于 X 光图像。最后,使用 勾股定理 计算两个梯度之间的幅度,并像以前一样对图像进行 归一化

**1.** 使用 SciPy 的傅里叶滤波器(scipy.ndimage.fourier_gaussian())和一个小的 sigma 值来去除 X 光图像中的一些噪声。然后,使用 scipy.ndimage.prewitt() 计算两个梯度。接下来,使用 NumPy 的 np.hypot() 测量梯度之间的距离。最后,像以前一样对重新缩放的图像进行 归一化

fourier_gaussian = ndimage.fourier_gaussian(xray_image, sigma=0.05)

x_prewitt = ndimage.prewitt(fourier_gaussian, axis=0)
y_prewitt = ndimage.prewitt(fourier_gaussian, axis=1)

xray_image_canny = np.hypot(x_prewitt, y_prewitt)

xray_image_canny *= 255.0 / np.max(xray_image_canny)

print("The data type - ", xray_image_canny.dtype)
The data type -  float64

**2.** 绘制原始 X 光图像和使用 Canny 滤波器技术检测边缘的图像。可以使用 prismnipy_spectralterrain Matplotlib 色图来突出显示边缘。

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(20, 15))

axes[0].set_title("Original")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Canny (edges) - prism")
axes[1].imshow(xray_image_canny, cmap="prism")
axes[2].set_title("Canny (edges) - nipy_spectral")
axes[2].imshow(xray_image_canny, cmap="nipy_spectral")
axes[3].set_title("Canny (edges) - terrain")
axes[3].imshow(xray_image_canny, cmap="terrain")
for i in axes:
    i.axis("off")
plt.show()
../_images/821ba442c270ffaf54f34336f10f42ca87fad9db0a281ec9b6d2b0eaf5277cb6.png

使用 np.where() 对 X 光图像应用掩码#

为了筛选出 X 光图像中的某些像素以帮助检测特定特征,可以使用 NumPy 的 np.where(condition: array_like (bool), x: array_like, y: ndarray),该函数在 True 时返回 x,在 False 时返回 y

识别感兴趣区域 - 图像中某些像素集 - 可能有用,掩码充当与原始图像形状相同的布尔数组。

**1.** 获取你一直在使用的原始 X 光图像中像素值的某些基本统计信息。

print("The data type of the X-ray image is: ", xray_image.dtype)
print("The minimum pixel value is: ", np.min(xray_image))
print("The maximum pixel value is: ", np.max(xray_image))
print("The average pixel value is: ", np.mean(xray_image))
print("The median pixel value is: ", np.median(xray_image))
The data type of the X-ray image is:  uint8
The minimum pixel value is:  0
The maximum pixel value is:  255
The average pixel value is:  172.52233219146729
The median pixel value is:  195.0

**2.** 数组数据类型是 uint8,最小值/最大值结果表明 X 光图像中使用了所有 256 种颜色(从 0255)。让我们使用 ndimage.histogram() 和 Matplotlib 可视化原始原始 X 光图像的 *像素强度分布*。

pixel_intensity_distribution = ndimage.histogram(
    xray_image, min=np.min(xray_image), max=np.max(xray_image), bins=256
)

plt.plot(pixel_intensity_distribution)
plt.title("Pixel intensity distribution")
plt.show()
../_images/cc4dfa333bc7edfe45a0dfc4922dd6f19824ddf4b8a737260f95590c55a5a8eb.png

正如像素强度分布所示,存在许多低值(大约在 0 到 20 之间)和非常高的值(大约在 200 到 240 之间)。

**3.** 你可以使用 NumPy 的 np.where() 创建不同的条件掩码 - 例如,让我们只保留那些像素值超过某个阈值的图像值。

# The threshold is "greater than 150"
# Return the original image if true, `0` otherwise
xray_image_mask_noisy = np.where(xray_image > 150, xray_image, 0)

plt.imshow(xray_image_mask_noisy, cmap="gray")
plt.axis("off")
plt.show()
../_images/d8fbefa6b6dbf425b86f93ef71e4dbcb2a30a0af6490080296061a0414c9cb9b.png
# The threshold is "greater than 150"
# Return `1` if true, `0` otherwise
xray_image_mask_less_noisy = np.where(xray_image > 150, 1, 0)

plt.imshow(xray_image_mask_less_noisy, cmap="gray")
plt.axis("off")
plt.show()
../_images/20ce77e6a25ece90b4987a9611be0a98709b265ea06faf39593b323e8051391f.png

比较结果#

让我们显示一些你到目前为止处理过的 X 光图像的结果。

fig, axes = plt.subplots(nrows=1, ncols=9, figsize=(30, 30))

axes[0].set_title("Original")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Laplace-Gaussian (edges)")
axes[1].imshow(xray_image_laplace_gaussian, cmap="gray")
axes[2].set_title("Gaussian gradient (edges)")
axes[2].imshow(x_ray_image_gaussian_gradient, cmap="gray")
axes[3].set_title("Sobel (edges) - grayscale")
axes[3].imshow(xray_image_sobel, cmap="gray")
axes[4].set_title("Sobel (edges) - hot")
axes[4].imshow(xray_image_sobel, cmap="hot")
axes[5].set_title("Canny (edges) - prism)")
axes[5].imshow(xray_image_canny, cmap="prism")
axes[6].set_title("Canny (edges) - nipy_spectral)")
axes[6].imshow(xray_image_canny, cmap="nipy_spectral")
axes[7].set_title("Mask (> 150, noisy)")
axes[7].imshow(xray_image_mask_noisy, cmap="gray")
axes[8].set_title("Mask (> 150, less noisy)")
axes[8].imshow(xray_image_mask_less_noisy, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()
../_images/ff067bd723259aec5922a6bfe4ca7bad495afdaa91e81ec81a0a798c4977feb6.png

下一步#

如果你想使用自己的样本,可以使用 此图像 或者在 Openi 数据库中搜索其他各种图像。Openi 包含许多生物医学图像,如果你带宽较低或受到可以下载的数据量限制,它特别有用。

要详细了解生物医学图像数据的图像处理或边缘检测,你可能会发现以下资料有用: