X射线图像处理#
本教程演示如何使用 NumPy、imageio、Matplotlib 和 SciPy 读取和处理 X 射线图像。您将学习如何加载医学图像、关注某些部分，并使用 Gaussian 、 Laplacian -Gaussian、 Sobel和 Canny滤波器进行边缘检测来直观地比较它们。
例如，当您 在Kaggle竞赛 中构建一种有助于检测肺炎的算法时， X 射线图像分析可以成为数据分析和 机器学习工作流程的一部分。在医疗保健行业，当图像估计 至少占所有医疗数据的 90% 时，医疗图像处理和分析就显得尤为重要。
您将使用来自
美国国立卫生研究院 (NIH)
提供的ChestX-ray8数据集的放射学图像。 ChestX-ray8 包含来自 30,000 多名患者的 100,000 多张 PNG 格式的go识别化 X 射线图像。您可以在 NIH 的公共 Box存储库的文件夹中找到 ChestX-ray8 的文件
。 （更多详情请参考
2017年CVPR（计算机视觉会议）发表的研究论文。）/images
为了您的方便，少量 PNG 图像已保存到本教程的存储库中tutorial-x-ray-image-processing/，因为 ChestX-ray8 包含千兆字节的数据，您可能会发现批量下载它具有挑战性。
先决条件#
读者应该具备一些 Python、NumPy 数组和 Matplotlib 知识。要刷新记忆，您可以学习 Python和 Matplotlib PyPlot教程以及 NumPy快速入门。
本教程使用以下包：
imageio用于读取和写入图像数据。医疗保健行业通常使用 DICOM格式进行医学成像，而 imageio应该非常适合读取该格式。为简单起见，在本教程中您将使用 PNG 文件。
Matplotlib用于数据可视化。
本教程可以在隔离环境中本地运行，例如 Virtualenv或 conda。您可以使用Jupyter Notebook 或 JupyterLab来运行每个笔记本单元。
目录＃
检查 X 射线
imageio将图像组合成多维数组以展示进展
使用拉普拉斯-高斯、高斯梯度、Sobel 和 Canny 滤波器进行边缘检测
将面罩应用于 X 射线
np.where()比较结果
用#检查 X 射线imageio
让我们从一个简单的示例开始，仅使用 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 张 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')
4.此外，以动画形式显示进度也会很有帮助。让我们创建一个 GIF 文件并imageio.mimwrite()在笔记本中显示结果：
GIF_PATH = os.path.join(DIR, "xray_image.gif")
imageio.mimwrite(GIF_PATH, combined_xray_images_1, format= ".gif", duration=1000)
这给了我们：
使用拉普拉斯-高斯、高斯梯度、Sobel 和 Canny 滤波器进行边缘检测#
处理生物医学数据时，强调 2D “边缘”以关注图像中的特定特征可能很有用 。为此，在检测颜色像素强度的变化时，使用 图像梯度特别有用。
具有高斯二阶导数的拉普拉斯滤波器#
让我们从 使用高斯二阶导数的 n维拉普拉斯滤波器（“拉普拉斯高斯”） 开始。这种拉普拉斯方法侧重于值强度快速变化的像素，并与高斯平滑相结合以 消除噪声。让我们来看看它如何在分析 2D X 射线图像时发挥作用。
拉普拉斯-高斯滤波器的实现相对简单：1）
ndimage从SciPy导入模块； 2)scipy.ndimage.gaussian_laplace()使用 sigma（标量）参数进行调用，该参数会影响高斯滤波器的标准差（您将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()
高斯梯度幅值法#
另一种有用的边缘检测方法是 高斯（梯度）滤波器。它使用高斯导数计算多维梯度幅值，并通过go除 高频 图像分量来提供帮助。
1.scipy.ndimage.gaussian_gradient_magnitude()
使用 sigma（标量）参数进行调用（用于标准差；您将2在下面的示例中使用）：
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()
Sobel-Feldman 算子（Sobel 滤波器）#
要沿 2D X 射线图像的水平轴和垂直轴查找高空间频率区域（边缘或边缘图），可以使用 Sobel-Feldman 算子（Sobel 滤波器） 技术。 Sobel 滤波器通过卷积将两个 3x3 核矩阵（每个轴一个）应用到 X 射线上。然后，使用毕达哥拉斯定理组合这两点（梯度） 以产生梯度幅值。
1.scipy.ndimage.sobel()在 X 射线的 x 轴和 y 轴上使用 Sobel 滤波器 — ( )。然后，使用毕达哥拉斯定理和 NumPy计算x和
之间的距离y（对它们应用 Sobel 滤波器）
以获得幅度。最后，将重新缩放后的图像的像素值标准化为 0 到 255 之间。np.hypot()
图像归一化
遵循以下公式。由于您使用的是灰度图像，因此您只需标准化一个通道。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.将新的图像数组数据类型更改为 32 位浮点格式，以float16使其
与 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“边缘”滤镜的 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()
Canny 过滤器#
您还可以考虑使用另一种著名的边缘检测滤波器，称为Canny 滤波器。
首先，应用高斯 滤波器go除图像中的噪声。在此示例中，您使用的是 傅立叶滤波器，它通过卷积过程平滑 X 射线 。接下来，您在图像的 2 个轴中的每一个上应用Prewitt 滤波器 ，以帮助检测一些边缘 - 这将产生 2 个梯度值。与 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和Matplotlib 颜色图强调边缘
terrain。
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()
使用#将掩模应用于 X 射线np.where()
要仅筛选出 X 射线图像中的某些像素以帮助检测特定特征，您可以使用 NumPy 应用遮罩，
该遮罩会返回when和when 。np.where(condition: array_like (bool), x: array_like, y: ndarray)xTrueyFalse
识别感兴趣的区域（图像中的某些像素集）可能很有用，并且掩模充当与原始图像形状相同的布尔数组。
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让我们使用 Matplotlib 可视化原始 X 射线图像的像素强度分布ndimage.histogram()：
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 处理和分割（放射学数据探索）
使用 Numpy 和 Scipy 进行图像操作和处理（Scipy 讲座笔记）
强度值（演示、DataCamp）
使用 Raspberry Pi 和 Python 进行物体检测（Maker Portal）
通过深度学习进行 X 射线数据准备和分割（Kaggle 托管的 Jupyter 笔记本）
图像过滤（讲座幻灯片，CS6670：计算机视觉，康奈尔大学）
Python 和 NumPy 中的边缘检测（走向数据科学）
使用 Scikit-Image进行边缘检测(Data Carpentry)
图像梯度和梯度过滤（讲座幻灯片，16-385 计算机视觉，卡内基梅隆大学）