NumPy 应用:X射线图像处理

本文详细介绍了如何使用Python的NumPy、imageio、Matplotlib和SciPy库读取、处理和分析X射线图像,包括加载图像、边缘检测(拉普拉斯-高斯、高斯梯度、Sobel和Canny滤波器)以及使用掩码进行数据筛选。通过ChestX-ray8数据集展示了这些技术在医学图像分析中的应用。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >


本教程演示了如何使用 NumPy、imageio、Matplotlib 和 SciPy 读取和处理 X 射线图像。您将学习如何加载医学图像,关注特定部分,并使用高斯、拉普拉斯-高斯、Sobel 和 Canny 过滤器进行边缘检测,通过视觉比较它们。

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

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

为了方便起见,本教程的存储库中保存了一小部分 PNG 图像,位于 tutorial-x-ray-image-processing/ 文件夹中,因为 ChestX-ray8 包含大量数据,您可能会发现批量下载它很具有挑战性。

一系列显示同一患者胸部相同区域的 9 张 X 射线图像,每张图像都应用了不同类型的图像处理滤波器。每张 X 射线图像显示了不同类型的生物细节。

先决条件

读者应具备一些 Python、NumPy 数组和 Matplotlib 的知识。为了恢复记忆,您可以参考 Python 和 Matplotlib 的教程,以及 NumPy 的快速入门。

本教程使用了以下软件包:

  • imageio 用于读取和写入图像数据。医疗保健行业通常使用 DICOM 格式进行医学成像,而 imageio 应该非常适合读取该格式。为了简单起见,在本教程中,您将使用 PNG 文件进行操作。
  • Matplotlib 用于数据可视化。
  • SciPy 用于通过 ndimage 进行多维图像处理。

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

目录

  1. 使用 imageio 检查 X 射线图像
  2. 将图像组合成多维数组以演示进展
  3. 使用拉普拉斯-高斯、高斯梯度、Sobel 和 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/8823375d747db0de213a6103b7b33d27ca2ac725bf6eab6b09f5b22a69ed406c.png

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

在下一个示例中,您将使用 ChestX-ray8 数据集中的 9 张 1024x1024 像素的 X 射线图像,这些图像已从数据集文件中下载并提取出来。它们的编号从 ...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/68b0dec3e74b8e1f78074af65ff2677edff488d44473a9610c40e51d24898563.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 射线图像,显示了患者胸部在不同时间点的相同视角。可以从帧到帧地直观比较患者的骨骼和内部器官。

使用拉普拉斯-高斯、高斯梯度、Sobel 和 Canny 过滤器进行边缘检测

在处理生物医学数据时,强调 2D “边缘” 可以有助于关注图像中的特定特征。为了做到这一点,使用 图像梯度 在检测颜色像素强度的变化时特别有帮助。

使用拉普拉斯-高斯、高斯梯度、Sobel 和 Canny 过滤器

让我们从使用 拉普拉斯 过滤器(“拉普拉斯-高斯”)开始,该过滤器使用 高斯 二阶导数。这种拉普拉斯方法侧重于值快速变化的像素,并与高斯平滑相结合以去除噪声。让我们看看它在分析 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 射线图像和应用拉普拉斯-高斯过滤器后的图像:

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

axes[0].set_title("原始图像")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("拉普拉斯-高斯(边缘)")
axes[1].imshow(xray_image_laplace_gaussian, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()

../_images/bd95089bfc08b32657c867ab41aba09e1874081ed0f713849ce2b85c3356e009.png

使用高斯梯度幅值法

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

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

x_ray_image_gaussian_gradient = ndimage.gaussian_gradient_magnitude(xray_image, sigma=2)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))

axes[0].set_title("原始图像")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("高斯梯度(边缘)")
axes[1].imshow(x_ray_image_gaussian_gradient, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()

../_images/975915cba39c6c45a7a1aefdfb47b82ab33f75f2ccc197684f7fbaf9d96bf382.png

Sobel-Feldman 算子(Sobel 滤波器)

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

1. 使用 Sobel 滤波器(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("更改前的数据类型:", xray_image_sobel.dtype)

xray_image_sobel = xray_image_sobel.astype("float32")

print("更改后的数据类型:", xray_image_sobel.dtype)
更改前的数据类型: float16
更改后的数据类型: float32

3. 显示原始 X 射线图像和应用 Sobel “边缘”滤波器的图像。注意,这里使用了灰度和 CMRmap 颜色映射来帮助强调边缘:

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

axes[0].set_title("原始图像")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Sobel(边缘)- 灰度")
axes[1].imshow(xray_image_sobel, cmap="gray")
axes[2].set_title("Sobel(边缘)- CMRmap")
axes[2].imshow(xray_image_sobel, cmap="CMRmap")
for i in axes:
    i.axis("off")
plt.show()

../_images/60943aa49a2b8fb72b2aaa7e2badee9422d24a7e80e54a8ef098b5be694bdaba.png

Canny 滤波器

您还可以考虑使用另一种用于边缘检测的著名滤波器,称为 Canny 滤波器

首先,应用 高斯 滤波器来消除图像中的噪声。在本例中,您使用的是 傅里叶 滤波器,它通过 卷积 过程平滑 X 射线图像。接下来,对图像的每个轴应用 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("数据类型 - ", xray_image_canny.dtype)
数据类型 -  float64

2. 绘制原始 X 射线图像以及使用 Canny 滤波器技术检测到的边缘图像。可以使用 prismnipy_spectralterrain Matplotlib 颜色映射来强调边缘。

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

axes[0].set_title("原始图像")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Canny(边缘)- prism")
axes[1].imshow(xray_image_canny, cmap="prism")
axes[2].set_title("Canny(边缘)- nipy_spectral")
axes[2].imshow(xray_image_canny, cmap="nipy_spectral")
axes[3].set_title("Canny(边缘)- terrain")
axes[3].imshow(xray_image_canny, cmap="terrain")
for i in axes:
    i.axis("off")
plt.show()

../_images/16aa720c547c2df3db628bda5e6892cbe3b383beff4f9662f849c73266f6f419.png

使用 np.where() 对 X 射线图像应用掩码

为了仅筛选出 X 射线图像中的特定像素以帮助检测特定特征,可以使用 NumPy 的 np.where(condition: array_like (bool), x:array_like, y: ndarray) 应用掩码,当条件为 True 时返回 x,否则返回 y

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

1. 获取有关您正在处理的原始 X 射线图像中像素值的一些基本统计信息:

print("X 射线图像的数据类型是:", xray_image.dtype)
print("最小像素值是:", np.min(xray_image))
print("最大像素值是:", np.max(xray_image))
print("平均像素值是:", np.mean(xray_image))
print("中位数像素值是:", np.median(xray_image))
X 射线图像的数据类型是: uint8
最小像素值是: 0
最大像素值是: 255
平均像素值是: 172.52233219146729
中位数像素值是: 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("像素强度分布")
plt.show()

../_images/4f458fe6c8214c2a5a6e9b5bb7941e68854dda5fbde3196d2e533bf3b4649f55.png

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

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

# 阈值为“大于 150”
# 如果为真,则返回原始图像,否则返回 `0`
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/28fbefa54339c06d139ecc19469207fa11163ad2ea8e2abd975f37eeb516994e.png

# 阈值为“大于 150”
# 如果为真,则返回 `1`,否则返回 `0`
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/e640c76062bbbe8b32546c848a1d0e5ad82d1eb08662f7fb4470233ef4cd7866.png

比较结果

让我们显示一些您迄今为止处理过的处理后的 X 射线图像的结果:

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

axes[0].set_title("原始图像")
axes[0].imshow(xray_image, cmap="gray")
axes[1].set_title("Laplace-Gaussian(边缘)")
axes[1].imshow(xray_image_laplace_gaussian, cmap="gray")
axes[2].set_title("高斯梯度(边缘)")
axes[2].imshow(x_ray_image_gaussian_gradient, cmap="gray")
axes[3].set_title("Sobel(边缘)- 灰度")
axes[3].imshow(xray_image_sobel, cmap="gray")
axes[4].set_title("Sobel(边缘)- 热图")
axes[4].imshow(xray_image_sobel, cmap="hot")
axes[5].set_title("Canny(边缘)- prism")
axes[5].imshow(xray_image_canny, cmap="prism")
axes[6].set_title("Canny(边缘)- nipy_spectral")
axes[6].imshow(xray_image_canny, cmap="nipy_spectral")
axes[7].set_title("掩码(> 150,有噪声)")
axes[7].imshow(xray_image_mask_noisy, cmap="gray")
axes[8].set_title("掩码(> 150,较少噪声)")
axes[8].imshow(xray_image_mask_less_noisy, cmap="gray")
for i in axes:
    i.axis("off")
plt.show()

下一步

如果你想使用自己的样本,你可以使用[这张图片](https://openi.nlm.nih.gov/detailedresult?img=CXR3666_IM-1824-1001&query=chest infection&it=xg&req=4&npos=32)或在Openi数据库中搜索其他各种图片。Openi包含许多生物医学图像,如果你的带宽较低和/或受到可下载数据量限制,它可能特别有帮助。

要了解更多关于生物医学图像数据处理或简单边缘检测的内容,你可能会发现以下材料有用:

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

数智笔记

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值