X光图像处理

X光图像处理#

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

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

你将使用来自ChestX-ray8数据集的放射影像，该数据集由美国国立卫生研究院 (NIH)提供。ChestX-ray8 包含超过 100,000 张来自 30,000 多名患者的去身份信息 PNG 格式 X 光图像。你可以在 NIH 的公共 Box 代码库的 /images 文件夹中找到 ChestX-ray8 的文件。（更多详细信息，请参考在 2017 年 CVPR（计算机视觉会议）上发表的研究论文。）

为方便起见，由于 ChestX-ray8 包含数 GB 的数据，你可能会发现难以批量下载，因此本教程的代码库中已保存少量 PNG 图像，位于 tutorial-x-ray-image-processing/ 下。

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进行多维图像处理。

本教程可以在隔离的环境中本地运行，例如Virtualenv或conda。你可以使用Jupyter Notebook 或 JupyterLab来运行每个 notebook 单元格。

目录#

  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/0d251fa796d1d9ac5e39be882e59c03b3529373b191d5e5ed89f08306eed901c.png

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

在下一个示例中，你将使用来自 ChestX-ray8 数据集的 9 张 1024x1024 像素的 X 光图像（而不是 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/9f0762d8cedbcecd08e8c8ae83b251f53af20ce16c1e866f4f3a94b1d06175b6.png

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

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/32a128991a775f52a7706843f93fde44e9095eaf2bc29c50d0475ad3b5a487a0.png

高斯梯度幅度法#

另一种有用的边缘检测方法是高斯(梯度)滤波器。它使用高斯导数计算多维梯度幅度，并有助于去除高频图像成分。

1. 使用scipy.ndimage.gaussian_gradient_magnitude()，并设置sigma (标量)参数（表示标准差；在下面的示例中将使用2）

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

2. 显示原始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/c5082edba518a9edae9b9cddb07b12bf16f75bd81b08cc11b3deaa5e819d500f.png

Sobel-Feldman算子（Sobel滤波器）#

为了找到二维X射线图像水平和垂直轴上的高空间频率区域（边缘或边缘图），可以使用Sobel-Feldman算子（Sobel滤波器）技术。Sobel滤波器通过卷积将两个3x3的卷积核矩阵（每个轴一个）应用于X射线图像。然后，使用勾股定理组合这两个点（梯度）以产生梯度幅度。

1. 使用Sobel滤波器 — (scipy.ndimage.sobel()) — 作用于X射线的x轴和y轴。然后，使用勾股定理和NumPy的np.hypot()计算应用了Sobel滤波器的x和y之间的距离以获得幅度。最后，对重新缩放的图像进行归一化，使像素值介于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射线图像和应用了Sobel“边缘”滤波器的图像。请注意，使用了灰度和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/66de24fb5e7e46f9acf9c294891a585fd9343e0443645710327b076270173e39.png

Canny滤波器#

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

首先，应用高斯滤波器去除图像中的噪声。在这个例子中，你使用的是傅里叶滤波器，它通过卷积过程来平滑X射线图像。接下来，在图像的两个轴上应用Prewitt滤波器来帮助检测一些边缘——这将产生两个梯度值。与Sobel滤波器类似，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滤波器技术检测到的边缘的图像。可以使用prism、nipy_spectral和terrain 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/8a1aa4842d6eb57d02a00a198ad83fb6a5347156fb75e17a4405a90cfc108d38.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种颜色（从0到255）。让我们使用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/c8b8c88a02c374237669f5662bb50e0e0a21e8fa8859a30de50320129463468c.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/9b6b15b250a702152d26fa062ab445de875e817eb986cc02ffa5eb5f56d003d7.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/b2bc018b0270827fbc35f69463c51694e970da7eced07958dda8f198758755f4.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/bbc80d80d61309ce0c5d0611fc36d529a1fe75078c69d34caa1bf27f1ab9ac4c.png

后续步骤#

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

要了解有关生物医学图像数据中的图像处理或仅仅是边缘检测的更多信息，你可能会发现以下资料有用：