X射线图像处理#
本教程演示了如何使用 NumPy、imageio、Matplotlib 和 SciPy 读取和处理 X 光图像。您将学习如何加载医学图像、关注某些部分，并使用 高斯、拉普拉斯-高斯、索贝尔 和 Canny 滤波器进行边缘检测来视觉比较它们。
X光图像分析可以成为你数据分析和机器学习工作流程的一部分，例如，当你构建一个帮助检测肺炎的算法作为Kaggle 竞赛的一部分时。在医疗行业中，医学图像处理和分析尤为重要，当图像估计占至少90%的所有医疗数据时。
你将使用来自 ChestX-ray8 数据集的放射影像，该数据集由 美国国立卫生研究院 (NIH) 提供。ChestX-ray8 包含超过 100,000 张去识别化的 PNG 格式 X 光图像，来自超过 30,000 名患者。你可以在 NIH 的公共 Box 仓库 的 /images 文件夹中找到 ChestX-ray8 的文件。（更多详情，请参阅 2017 年在 CVPR（一个计算机视觉会议）上发表的 论文。）
为了方便起见，少量 PNG 图像已保存到本教程的仓库中，位于 tutorial-x-ray-image-processing/ 下，因为 ChestX-ray8 包含数 GB 的数据，您可能会发现批量下载它具有挑战性。
前提条件#
读者应该具备一些 Python、NumPy 数组和 Matplotlib 的知识。为了复习记忆，你可以参加 Python 和 Matplotlib PyPlot 教程，以及 NumPy 快速入门。
以下包在本教程中使用：
imageio 用于读取和写入图像数据。医疗保健行业通常使用 DICOM 格式进行医学成像，而 imageio 应该非常适合读取该格式。为了简单起见，在本教程中，您将使用 PNG 文件。
Matplotlib 用于数据可视化。
本教程可以在一个隔离的环境中本地运行，例如 Virtualenv 或 conda。您可以使用 Jupyter Notebook 或 JupyterLab 来运行每个笔记本单元。
目录#
使用
imageio检查 X 光片将图像组合成多维数组以展示进展
使用拉普拉斯-高斯、高斯梯度、索贝尔和坎尼滤波器的边缘检测
使用
np.where()对 X 光片应用掩码比较结果
使用 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()
将图像组合成多维数组以演示进展#
在下一个示例中，您将使用从 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 数组可以被解释为一系列二维图像的堆叠。
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')
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)
这给了我们：
使用拉普拉斯-高斯、高斯梯度、索贝尔和坎尼滤波器的边缘检测#
处理生物医学数据时，强调图像的二维”边缘”可能有助于聚焦于图像中的特定特征。为此，在检测颜色像素强度变化时，使用图像梯度可能特别有帮助。
带有高斯二阶导数的拉普拉斯滤波器#
让我们从一个n维的 Laplace 滤波器（“拉普拉斯-高斯”）开始，该滤波器使用 Gaussian 的二阶导数。这种拉普拉斯方法专注于值中快速强度变化的像素，并与高斯平滑结合以 去除噪声。让我们看看它在分析2D X射线图像时如何有用。
Laplacian-Gaussian 滤波器的实现相对简单：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光图像和经过拉普拉斯-高斯滤波器的图像：
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()
高斯梯度幅度方法#
另一种有用的边缘检测方法是 Gaussian (梯度) 滤波器。它通过高斯导数计算多维梯度幅值，并通过去除 高频 图像成分来提供帮助。
1. 使用 sigma（标量）参数（用于标准差；你将在下面的示例中使用 2）调用 scipy.ndimage.gaussian_gradient_magnitude()：
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()
Sobel-Feldman 算子（Sobel 滤波器）#
要找到沿2D X射线图像的水平和垂直轴的高空间频率区域（边缘或边缘图），可以使用 Sobel-Feldman 算子（Sobel 滤波器） 技术。Sobel 滤波器通过 卷积 将两个3x3的核矩阵——每个轴一个——应用于X射线图像。然后，这两个点（梯度）使用 勾股定理 组合以产生梯度幅值。
1. 使用 Sobel 滤波器 — (scipy.ndimage.sobel()) — 在 X 射线的 x 轴和 y 轴上。然后，使用 勾股定理 和 NumPy 的 np.hypot() 计算 x 和 y（应用 Sobel 滤波器后）之间的距离以获得幅度。最后，对重缩放的图像进行归一化，使像素值在 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()
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()
使用 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()
正如像素强度分布所示，有许多低（大约在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()
# 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()
比较结果#
让我们展示一些你迄今为止处理过的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()
下一步#
如果你想使用自己的样本，你可以使用这张图片或在Openi数据库中搜索其他各种图片。Openi包含许多生物医学图像，如果你带宽较低和/或受限于可以下载的数据量，它尤其有用。
要了解更多关于生物医学图像数据中的图像处理或简单的边缘检测，您可能会发现以下材料有用：
使用 Scikit-Image 和 pydicom 在 Python 中进行 DICOM 处理和分割（Radiology Data Quest）
使用 Numpy 和 Scipy 进行图像处理 (Scipy 讲义)
强度值 (演示文稿, DataCamp)
使用树莓派和Python进行物体检测 (Maker Portal)
X射线数据准备和分割 使用深度学习 (一个由Kaggle托管的Jupyter笔记本)
图像过滤 (讲座幻灯片, CS6670: 计算机视觉, 康奈尔大学)
Python中的边缘检测 和 NumPy (Towards Data Science)
边缘检测 使用 Scikit-Image (数据 Carpentry)
图像梯度和梯度滤波 (讲座幻灯片, 16-385 计算机视觉, 卡内基梅隆大学)