掌握Scipy医学图像处理核心技术从图像预处理到特征提取的完整工作流程详解Python科学计算库在医学影像领域的创新应用

引言

医学图像处理在现代医疗诊断、治疗规划和研究中扮演着至关重要的角色。随着医学成像技术的快速发展，CT、MRI、X光、超声等多种成像方式产生的医学图像数据量呈指数级增长。如何高效、准确地处理这些海量数据，从中提取有价值的信息，成为医学影像领域的关键挑战。Python科学计算库Scipy凭借其强大的数值计算能力和丰富的图像处理功能，为医学图像处理提供了高效、灵活的解决方案。本文将详细介绍使用Scipy进行医学图像处理的完整工作流程，从图像预处理到特征提取，并探讨Python科学计算库在医学影像领域的创新应用。

Scipy库概述

Scipy（Scientific Python）是一个开源的Python算法库和数学工具包，基于NumPy构建，提供了更多的高级数学函数、信号处理工具、优化算法以及图像处理功能。在医学图像处理领域，Scipy的ndimage模块（多维图像处理）尤其重要，它提供了一系列用于图像处理的功能，包括滤波、插值、形态学操作、图像测量等。

Scipy的主要优势在于：

高效性：底层使用C和Fortran实现，计算效率高
易用性：简洁的API设计，易于学习和使用
兼容性：与NumPy、Matplotlib等科学计算库无缝集成
开源免费：遵循BSD许可证，可自由使用和分发

以下是一个简单的Scipy导入示例：

# 导入必要的库 import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt from scipy import misc, signal

医学图像基础

医学图像具有一些独特的特性，了解这些特性对于正确处理和分析医学图像至关重要。

常见医学图像类型

CT（计算机断层扫描）：基于X射线，提供高分辨率的骨骼和软组织图像，通常以Hounsfield单位（HU）表示密度值。
MRI（磁共振成像）：利用磁场和射频脉冲，产生高对比度的软组织图像，有不同的序列（T1、T2、DWI等）。
X光片：传统的二维投影图像，主要用于检查骨折和肺部疾病。
超声：基于声波反射，实时成像，常用于产科和心脏检查。
PET（正电子发射断层扫描）：功能性成像，显示代谢活动，常与CT结合（PET-CT）。

医学图像特点

高动态范围：医学图像通常具有较宽的强度范围，需要特殊处理。
噪声特性：不同成像方式有不同的噪声模型（如CT中的泊松噪声）。
各向异性：体素在三个维度上的分辨率可能不同。
标准化需求：医学图像通常需要标准化处理，以便进行比较和定量分析。

以下代码展示了如何加载和显示医学图像：

import numpy as np import matplotlib.pyplot as plt from scipy import misc import skimage.io as io # 加载示例医学图像（这里使用scipy自带的图像作为示例） # 实际应用中，可以使用专门的医学图像加载库如SimpleITK或nibabel image = misc.face(gray=True) # 显示图像 plt.figure(figsize=(10, 8)) plt.imshow(image, cmap='gray') plt.title('示例医学图像') plt.colorbar() plt.show() # 显示图像直方图 plt.figure(figsize=(10, 6)) plt.hist(image.ravel(), bins=256, range=(0, 256)) plt.title('图像直方图') plt.xlabel('像素强度') plt.ylabel('频率') plt.show()

图像预处理技术

图像预处理是医学图像分析的第一步，目的是改善图像质量，增强相关特征，为后续处理做准备。

图像读取与显示

医学图像通常以特殊格式存储，如DICOM（.dcm）、NIfTI（.nii）等。Scipy本身不直接支持这些格式，但可以与其他库结合使用。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage import nibabel as nib # 用于读取NIfTI格式医学图像 # 读取NIfTI格式的医学图像 # 注意：实际使用时需要替换为真实文件路径 try: img = nib.load('example_brain.nii.gz') image_data = img.get_fdata() print(f"图像维度: {image_data.shape}") # 显示中间切片 mid_slice = image_data[:, :, image_data.shape[2]//2] plt.figure(figsize=(10, 8)) plt.imshow(mid_slice.T, cmap='gray', origin='lower') plt.title('脑部MRI中间切片') plt.colorbar() plt.show() except: # 如果没有真实医学图像，使用模拟数据 print("未找到医学图像文件，使用模拟数据") # 创建模拟的3D医学图像数据 x, y, z = np.ogrid[-10:10:100j, -10:10:100j, -10:10:100j] sphere = x**2 + y**2 + z**2 <= 10**2 image_data = sphere.astype(float) # 添加一些噪声 image_data += np.random.normal(0, 0.1, image_data.shape) # 显示中间切片 mid_slice = image_data[:, :, 50] plt.figure(figsize=(10, 8)) plt.imshow(mid_slice.T, cmap='gray', origin='lower') plt.title('模拟脑部MRI中间切片') plt.colorbar() plt.show()

图像增强

图像增强的目的是改善图像的视觉效果，突出特定特征。

对比度调整

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import exposure # 加载图像（使用scipy自带图像作为示例） image = misc.face(gray=True) # 对比度拉伸 p2, p98 = np.percentile(image, (2, 98)) img_rescale = exposure.rescale_intensity(image, in_range=(p2, p98)) # 直方图均衡化 img_eq = exposure.equalize_hist(image) # 自适应直方图均衡化 img_adapteq = exposure.equalize_adapthist(image, clip_limit=0.03) # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(image, cmap='gray') plt.title('原始图像') plt.subplot(222) plt.imshow(img_rescale, cmap='gray') plt.title('对比度拉伸') plt.subplot(223) plt.imshow(img_eq, cmap='gray') plt.title('直方图均衡化') plt.subplot(224) plt.imshow(img_adapteq, cmap='gray') plt.title('自适应直方图均衡化') plt.tight_layout() plt.show() # 显示直方图 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.hist(image.ravel(), bins=256, range=(0, 256)) plt.title('原始图像直方图') plt.subplot(222) plt.hist(img_rescale.ravel(), bins=256, range=(0, 256)) plt.title('对比度拉伸后直方图') plt.subplot(223) plt.hist(img_eq.ravel(), bins=256, range=(0, 1)) plt.title('直方图均衡化后直方图') plt.subplot(224) plt.hist(img_adapteq.ravel(), bins=256, range=(0, 1)) plt.title('自适应直方图均衡化后直方图') plt.tight_layout() plt.show()

图像滤波

图像滤波用于去除噪声、增强特定特征或提取图像信息。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage, signal from skimage import filters, img_as_float # 加载图像并添加噪声 image = misc.face(gray=True) image = img_as_float(image) noisy = image + 0.1 * np.random.randn(*image.shape) noisy = np.clip(noisy, 0, 1) # 确保值在[0,1]范围内 # 应用不同的滤波器 # 高斯滤波 gaussian_filtered = ndimage.gaussian_filter(noisy, sigma=1.0) # 中值滤波 median_filtered = ndimage.median_filter(noisy, size=3) # 维纳滤波（需要scipy.signal） wiener_filtered = signal.wiener(noisy, (5, 5)) # 双边滤波（需要skimage） bilateral_filtered = filters.denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=1) # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(231) plt.imshow(image, cmap='gray') plt.title('原始图像') plt.subplot(232) plt.imshow(noisy, cmap='gray') plt.title('添加噪声后的图像') plt.subplot(233) plt.imshow(gaussian_filtered, cmap='gray') plt.title('高斯滤波') plt.subplot(234) plt.imshow(median_filtered, cmap='gray') plt.title('中值滤波') plt.subplot(235) plt.imshow(wiener_filtered, cmap='gray') plt.title('维纳滤波') plt.subplot(236) plt.imshow(bilateral_filtered, cmap='gray') plt.title('双边滤波') plt.tight_layout() plt.show()

图像配准

图像配准是将不同时间、不同视角或不同模态的图像对齐的过程。在医学图像处理中，配准对于比较患者前后的扫描、融合不同模态的图像等非常重要。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import transform, data, exposure # 加载参考图像和待配准图像 reference_image = data.camera() moving_image = transform.rotate(reference_image, 15) moving_image = ndimage.zoom(moving_image, 0.9) moving_image = ndimage.shift(moving_image, (20, -30)) # 确保图像大小一致 h, w = reference_image.shape moving_image_resized = transform.resize(moving_image, (h, w)) # 显示原始图像和待配准图像 plt.figure(figsize=(12, 6)) plt.subplot(121) plt.imshow(reference_image, cmap='gray') plt.title('参考图像') plt.subplot(122) plt.imshow(moving_image_resized, cmap='gray') plt.title('待配准图像') plt.show() # 简单的相位相关配准（仅用于平移） def phase_correlation(reference, moving): # 计算傅里叶变换 f1 = np.fft.fft2(reference) f2 = np.fft.fft2(moving) # 计算交叉功率谱 cross_power_spectrum = (f1 * np.conj(f2)) / (np.abs(f1) * np.abs(f2)) # 计算逆傅里叶变换 result = np.fft.ifft2(cross_power_spectrum) # 找到峰值位置 peak = np.unravel_index(np.argmax(np.abs(result)), result.shape) # 计算平移量 shift_y = peak[0] if peak[0] <= h//2 else peak[0] - h shift_x = peak[1] if peak[1] <= w//2 else peak[1] - w return shift_y, shift_x # 计算平移量 shift_y, shift_x = phase_correlation(reference_image, moving_image_resized) print(f"检测到的平移量: (y={shift_y}, x={shift_x})") # 应用平移 registered_image = ndimage.shift(moving_image_resized, (-shift_y, -shift_x)) # 显示配准结果 plt.figure(figsize=(12, 6)) plt.subplot(121) plt.imshow(reference_image, cmap='gray') plt.title('参考图像') plt.subplot(122) plt.imshow(registered_image, cmap='gray') plt.title('配准后图像') plt.show() # 显示差异图像 plt.figure(figsize=(12, 6)) plt.subplot(121) plt.imshow(np.abs(reference_image - moving_image_resized), cmap='hot') plt.title('配准前差异') plt.subplot(122) plt.imshow(np.abs(reference_image - registered_image), cmap='hot') plt.title('配准后差异') plt.show()

图像分割技术

图像分割是将图像划分为多个区域或对象的过程，是医学图像分析中的关键步骤。准确的分割有助于识别和量化解剖结构或病变区域。

阈值分割

阈值分割是最简单的分割方法之一，适用于目标与背景有明显对比度差异的情况。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import filters, data, color # 加载图像 image = data.brain() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image = color.rgb2gray(image) # 应用不同的阈值方法 # 全局阈值 thresh_global = filters.threshold_otsu(image) binary_global = image > thresh_global # 自适应阈值 block_size = 35 binary_adaptive = filters.threshold_local(image, block_size, offset=10) binary_adaptive = image > binary_adaptive # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(image, cmap='gray') plt.title('原始图像') plt.subplot(222) plt.hist(image.ravel(), bins=256) plt.axvline(thresh_global, color='r') plt.title('直方图与Otsu阈值') plt.subplot(223) plt.imshow(binary_global, cmap='gray') plt.title(f'全局阈值分割 (阈值={thresh_global:.2f})') plt.subplot(224) plt.imshow(binary_adaptive, cmap='gray') plt.title('自适应阈值分割') plt.tight_layout() plt.show()

区域生长

区域生长是一种基于相似性的分割方法，从种子点开始，逐步将相邻的相似像素合并到区域中。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import data, color, segmentation # 加载图像 image = data.coins() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image = color.rgb2gray(image) # 定义区域生长算法 def region_growing(image, seed, threshold): """ 简单的区域生长算法实现 参数: image: 输入图像 seed: 种子点坐标 (row, col) threshold: 相似性阈值 返回: 分割后的二值图像 """ # 创建与输入图像大小相同的标记矩阵 rows, cols = image.shape segmented = np.zeros_like(image, dtype=np.uint8) # 获取种子点值 seed_value = image[seed[0], seed[1]] # 创建队列并将种子点加入 queue = [seed] segmented[seed[0], seed[1]] = 1 # 定义8邻域 neighbors = [(-1,-1), (-1,0), (-1,1), (0,-1), (0,1), (1,-1), (1,0), (1,1)] while queue: current_point = queue.pop(0) # 检查所有邻域 for d in neighbors: r = current_point[0] + d[0] c = current_point[1] + d[1] # 检查边界 if r < 0 or r >= rows or c < 0 or c >= cols: continue # 如果像素未被访问且与种子值相似 if segmented[r, c] == 0 and abs(int(image[r, c]) - int(seed_value)) < threshold: segmented[r, c] = 1 queue.append((r, c)) return segmented # 选择种子点（这里手动选择，实际应用中可以交互式选择） seed_points = [(150, 150), (250, 300), (350, 200)] threshold = 30 # 相似性阈值 # 对每个种子点应用区域生长 segmented_image = np.zeros_like(image, dtype=np.uint8) for seed in seed_points: region = region_growing(image, seed, threshold) segmented_image = np.maximum(segmented_image, region) # 显示结果 plt.figure(figsize=(15, 5)) plt.subplot(131) plt.imshow(image, cmap='gray') plt.title('原始图像') for seed in seed_points: plt.plot(seed[1], seed[0], 'ro') plt.subplot(132) plt.imshow(segmented_image, cmap='gray') plt.title('区域生长分割结果') plt.subplot(133) plt.imshow(color.label2rgb(segmented_image, image, bg_label=0)) plt.title('彩色叠加显示') plt.tight_layout() plt.show()

边缘检测

边缘检测是识别图像中亮度变化明显的点，这些点通常对应于物体的边界。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage, signal from skimage import filters, feature, data # 加载图像 image = data.brain() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image = color.rgb2gray(image) # 应用不同的边缘检测方法 # Sobel算子 sobel_h = ndimage.sobel(image, axis=0) sobel_v = ndimage.sobel(image, axis=1) sobel = np.sqrt(sobel_h**2 + sobel_v**2) # Canny边缘检测 edges_canny = feature.canny(image, sigma=1.0) # LoG (Laplacian of Gaussian) log = ndimage.gaussian_laplace(image, sigma=1.0) # 寻找零交叉点 zero_crossing = np.zeros_like(log) for i in range(1, log.shape[0]-1): for j in range(1, log.shape[1]-1): block = log[i-1:i+2, j-1:j+2] if (block[1,1] < 0 and (block[0,1] > 0 or block[2,1] > 0 or block[1,0] > 0 or block[1,2] > 0)) or (block[1,1] > 0 and (block[0,1] < 0 or block[2,1] < 0 or block[1,0] < 0 or block[1,2] < 0)): zero_crossing[i, j] = 1 # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(image, cmap='gray') plt.title('原始图像') plt.subplot(222) plt.imshow(sobel, cmap='gray') plt.title('Sobel边缘检测') plt.subplot(223) plt.imshow(edges_canny, cmap='gray') plt.title('Canny边缘检测') plt.subplot(224) plt.imshow(zero_crossing, cmap='gray') plt.title('LoG零交叉边缘检测') plt.tight_layout() plt.show()

基于机器学习的分割

基于机器学习的分割方法利用训练数据学习图像特征与分割结果之间的关系，可以实现更复杂的分割任务。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import data, color, segmentation, feature, future from sklearn.ensemble import RandomForestClassifier from functools import partial # 加载图像 image = data.brain() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image_gray = color.rgb2gray(image) else: image_gray = image.copy() # 创建训练标签（这里使用简单的阈值创建示例标签） # 在实际应用中，这些标签应该由专家手动标注 labels = np.zeros(image_gray.shape, dtype=np.uint8) labels[50:150, 50:150] = 1 # 前景区域 labels[200:300, 200:300] = 2 # 背景区域 # 提取特征 def _features(image): """提取图像特征""" features = [] # 原始强度 features.append(image) # 高斯滤波后的强度 for sigma in [1, 2, 4]: features.append(ndimage.gaussian_filter(image, sigma=sigma)) # 梯度幅值 features.append(ndimage.sobel(image)) # 拉普拉斯 features.append(ndimage.laplace(image)) return np.stack(features, axis=-1) # 提取特征 features = _features(image_gray) # 训练随机森林分类器 clf = RandomForestClassifier(n_estimators=50, n_jobs=-1, max_depth=10, random_state=0) # 将标签展平为1D数组 training_labels = labels[labels > 0] training_features = features[labels > 0] # 训练分类器 clf.fit(training_features, training_labels) # 预测整个图像 result = clf.predict(features.reshape(-1, features.shape[-1])) result = result.reshape(image_gray.shape) # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(image_gray, cmap='gray') plt.title('原始图像') plt.subplot(222) plt.imshow(labels, cmap='jet') plt.title('训练标签') plt.subplot(223) plt.imshow(result, cmap='jet') plt.title('随机森林分割结果') plt.subplot(224) plt.imshow(color.label2rgb(result, image_gray, bg_label=0)) plt.title('彩色叠加显示') plt.tight_layout() plt.show()

特征提取

特征提取是从分割后的图像区域中提取定量信息的过程，这些特征可以用于分类、检测、诊断等任务。

形态学特征

形态学特征描述了目标区域的形状和大小特性。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import measure, data, color, draw # 创建一个示例二值图像 image = np.zeros((300, 300), dtype=np.uint8) # 添加一些不同形状的对象 # 圆形 rr, cc = draw.circle(100, 100, 50) image[rr, cc] = 1 # 椭圆 rr, cc = draw.ellipse(200, 200, 60, 30) image[rr, cc] = 1 # 矩形 image[50:100, 200:250] = 1 # 不规则形状 rr, cc = draw.polygon([150, 200, 250, 200, 150], [50, 30, 50, 70, 90]) image[rr, cc] = 1 # 标记连通区域 label_image = measure.label(image) regions = measure.regionprops(label_image) # 显示原始图像和标记结果 plt.figure(figsize=(12, 6)) plt.subplot(121) plt.imshow(image, cmap='gray') plt.title('原始二值图像') plt.subplot(122) plt.imshow(label_image, cmap='jet') plt.title('标记后的图像') plt.show() # 提取并显示每个区域的形态学特征 print("形态学特征分析:") print("-" * 80) for i, region in enumerate(regions): print(f"区域 {i+1}:") print(f" 面积: {region.area} 像素") print(f" 周长: {region.perimeter} 像素") print(f" 离心率: {region.eccentricity:.4f}") print(f" 圆度: {(4 * np.pi * region.area / region.perimeter**2):.4f}") print(f" 长轴长度: {region.major_axis_length:.2f} 像素") print(f" 短轴长度: {region.minor_axis_length:.2f} 像素") print(f" 方向: {region.orientation:.4f} 弧度") print(f" 质心: ({region.centroid[0]:.2f}, {region.centroid[1]:.2f})") print(f" 边界框: {region.bbox}") print("-" * 80) # 可视化一些特征 plt.figure(figsize=(15, 10)) for i, region in enumerate(regions): plt.subplot(2, 2, i+1) plt.imshow(image, cmap='gray') # 绘制质心 y0, x0 = region.centroid plt.plot(x0, y0, 'ro') # 绘制边界框 minr, minc, maxr, maxc = region.bbox plt.plot([minc, maxc, maxc, minc, minc], [minr, minr, maxr, maxr, minr], 'r-') # 绘制主轴 y0, x0 = region.centroid orientation = region.orientation length = region.major_axis_length / 2 x1 = x0 + length * np.cos(orientation) y1 = y0 + length * np.sin(orientation) x2 = x0 - length * np.cos(orientation) y2 = y0 - length * np.sin(orientation) plt.plot([x1, x2], [y1, y2], 'g-') plt.title(f'区域 {i+1} 特征可视化') plt.tight_layout() plt.show()

纹理特征

纹理特征描述了图像区域的表面特性，对于区分不同组织类型非常有用。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from skimage import data, color, feature from skimage.filters import gabor # 加载图像 image = data.brick() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image_gray = color.rgb2gray(image) else: image_gray = image.copy() # 提取不同的纹理特征 # 1. 局部二值模式 (LBP) radius = 3 n_points = 8 * radius lbp = feature.local_binary_pattern(image_gray, n_points, radius, method='uniform') # 2. Gabor滤波器 freqs = [0.1, 0.3, 0.5] gabor_responses = [] for freq in freqs: filt_real, filt_imag = gabor(image_gray, frequency=freq) gabor_responses.append(np.sqrt(filt_real**2 + filt_imag**2)) # 3. 灰度共生矩阵 (GLCM) 特征 def glcm_features(image, distances=[1], angles=[0, np.pi/4, np.pi/2, 3*np.pi/4], levels=256): """计算灰度共生矩阵特征""" from skimage.feature import greycomatrix, greycoprops # 计算灰度共生矩阵 glcm = greycomatrix(image, distances=distances, angles=angles, levels=levels, symmetric=True, normed=True) # 计算特征 features = {} properties = ['contrast', 'dissimilarity', 'homogeneity', 'energy'] for prop in properties: features[prop] = greycoprops(glcm, prop).mean() return features # 计算GLCM特征 glcm_feats = glcm_features(image_gray) # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(231) plt.imshow(image_gray, cmap='gray') plt.title('原始图像') plt.subplot(232) plt.imshow(lbp, cmap='gray') plt.title('局部二值模式 (LBP)') for i, (freq, resp) in enumerate(zip(freqs, gabor_responses)): plt.subplot(233 + i) plt.imshow(resp, cmap='gray') plt.title(f'Gabor滤波器 (频率={freq})') plt.tight_layout() plt.show() # 打印GLCM特征 print("灰度共生矩阵 (GLCM) 特征:") print("-" * 40) for feat_name, feat_value in glcm_feats.items(): print(f"{feat_name}: {feat_value:.4f}") print("-" * 40) # 显示LBP直方图 plt.figure(figsize=(10, 6)) plt.hist(lbp.ravel(), bins=n_points + 2, range=(0, n_points + 2)) plt.title('LBP直方图') plt.xlabel('LBP值') plt.ylabel('频率') plt.show()

强度特征

强度特征描述了图像区域的亮度分布特性。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage, stats from skimage import data, color, measure # 加载图像 image = data.brain() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image_gray = color.rgb2gray(image) else: image_gray = image.copy() # 创建一个掩码（这里使用简单的阈值） threshold = filters.threshold_otsu(image_gray) mask = image_gray > threshold # 标记连通区域 label_image = measure.label(mask) regions = measure.regionprops(label_image, intensity_image=image_gray) # 显示原始图像和掩码 plt.figure(figsize=(12, 6)) plt.subplot(121) plt.imshow(image_gray, cmap='gray') plt.title('原始图像') plt.subplot(122) plt.imshow(mask, cmap='gray') plt.title('二值掩码') plt.show() # 提取并显示每个区域的强度特征 print("强度特征分析:") print("-" * 80) for i, region in enumerate(regions[:5]): # 只显示前5个区域 print(f"区域 {i+1}:") print(f" 平均强度: {region.mean_intensity:.4f}") print(f" 最大强度: {region.max_intensity:.4f}") print(f" 最小强度: {region.min_intensity:.4f}") print(f" 强度标准差: {region.intensity_std:.4f}") # 计算强度直方图特征 intensity_values = region.intensity_image[region.image] hist, bin_edges = np.histogram(intensity_values, bins=10) # 计算偏度和峰度 skewness = stats.skew(intensity_values) kurtosis = stats.kurtosis(intensity_values) print(f" 强度偏度: {skewness:.4f}") print(f" 强度峰度: {kurtosis:.4f}") print("-" * 80) # 可视化一些区域的强度分布 plt.figure(figsize=(15, 10)) for i, region in enumerate(regions[:4]): plt.subplot(2, 2, i+1) # 获取该区域的强度值 intensity_values = region.intensity_image[region.image] # 绘制直方图 plt.hist(intensity_values, bins=20, alpha=0.7) plt.axvline(region.mean_intensity, color='r', linestyle='dashed', linewidth=1, label='平均值') plt.axvline(region.mean_intensity + region.intensity_std, color='g', linestyle='dashed', linewidth=1, label='标准差') plt.axvline(region.mean_intensity - region.intensity_std, color='g', linestyle='dashed', linewidth=1) plt.title(f'区域 {i+1} 强度分布') plt.xlabel('强度值') plt.ylabel('频率') plt.legend() plt.tight_layout() plt.show()

频域特征

频域特征描述了图像在频率域中的特性，对于分析周期性模式和纹理非常有用。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage, fftpack from skimage import data, color # 加载图像 image = data.camera() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image_gray = color.rgb2gray(image) else: image_gray = image.copy() # 计算傅里叶变换 f_transform = fftpack.fft2(image_gray) f_shift = fftpack.fftshift(f_transform) magnitude_spectrum = 20 * np.log(np.abs(f_shift)) # 创建不同的滤波器 rows, cols = image_gray.shape crow, ccol = rows // 2, cols // 2 # 低通滤波器 low_pass_filter = np.zeros((rows, cols), np.uint8) r = 30 # 半径 y, x = np.ogrid[:rows, :cols] mask_area = (x - ccol) ** 2 + (y - crow) ** 2 <= r*r low_pass_filter[mask_area] = 1 # 高通滤波器 high_pass_filter = np.ones((rows, cols), np.uint8) high_pass_filter[mask_area] = 0 # 带通滤波器 band_pass_filter = np.zeros((rows, cols), np.uint8) r_inner, r_outer = 10, 50 mask_area_inner = (x - ccol) ** 2 + (y - crow) ** 2 >= r_inner*r_inner mask_area_outer = (x - ccol) ** 2 + (y - crow) ** 2 <= r_outer*r_outer band_pass_filter[mask_area_inner & mask_area_outer] = 1 # 应用滤波器 f_shift_low = f_shift * low_pass_filter f_shift_high = f_shift * high_pass_filter f_shift_band = f_shift * band_pass_filter # 逆变换 img_back_low = np.abs(fftpack.ifft2(fftpack.ifftshift(f_shift_low))) img_back_high = np.abs(fftpack.ifft2(fftpack.ifftshift(f_shift_high))) img_back_band = np.abs(fftpack.ifft2(fftpack.ifftshift(f_shift_band))) # 显示结果 plt.figure(figsize=(15, 12)) plt.subplot(331) plt.imshow(image_gray, cmap='gray') plt.title('原始图像') plt.subplot(332) plt.imshow(magnitude_spectrum, cmap='gray') plt.title('频谱') plt.subplot(333) plt.imshow(low_pass_filter, cmap='gray') plt.title('低通滤波器') plt.subplot(334) plt.imshow(img_back_low, cmap='gray') plt.title('低通滤波结果') plt.subplot(335) plt.imshow(high_pass_filter, cmap='gray') plt.title('高通滤波器') plt.subplot(336) plt.imshow(img_back_high, cmap='gray') plt.title('高通滤波结果') plt.subplot(337) plt.imshow(band_pass_filter, cmap='gray') plt.title('带通滤波器') plt.subplot(338) plt.imshow(img_back_band, cmap='gray') plt.title('带通滤波结果') plt.tight_layout() plt.show() # 提取频域特征 def extract_frequency_features(image): """提取频域特征""" # 计算傅里叶变换 f_transform = fftpack.fft2(image) f_shift = fftpack.fftshift(f_transform) magnitude_spectrum = np.abs(f_shift) # 计算径向分布 rows, cols = image.shape crow, ccol = rows // 2, cols // 2 # 创建径向距离矩阵 y, x = np.ogrid[:rows, :cols] r = np.sqrt((x - ccol) ** 2 + (y - crow) ** 2) r = r.astype(int) # 计算径向平均 tbin = np.bincount(r.ravel(), magnitude_spectrum.ravel()) nr = np.bincount(r.ravel()) radial_profile = tbin / nr # 计算特征 features = { 'dc_component': magnitude_spectrum[crow, ccol], 'spectral_centroid': np.sum(np.arange(len(radial_profile)) * radial_profile) / np.sum(radial_profile), 'spectral_spread': np.sqrt(np.sum(((np.arange(len(radial_profile)) - features['spectral_centroid']) ** 2) * radial_profile) / np.sum(radial_profile)), 'spectral_energy': np.sum(radial_profile ** 2), 'spectral_entropy': -np.sum(radial_profile * np.log2(radial_profile + 1e-10)) } return features, radial_profile # 提取频域特征 freq_features, radial_profile = extract_frequency_features(image_gray) # 打印频域特征 print("频域特征分析:") print("-" * 40) for feat_name, feat_value in freq_features.items(): print(f"{feat_name}: {feat_value:.4f}") print("-" * 40) # 显示径向分布 plt.figure(figsize=(10, 6)) plt.plot(radial_profile) plt.title('频谱径向分布') plt.xlabel('频率半径') plt.ylabel('平均幅度') plt.grid(True) plt.show()

实际案例分析

让我们通过一个完整的案例来展示如何使用Scipy处理实际的医学图像数据。在这个案例中，我们将处理脑部MRI图像，从预处理到特征提取的完整流程。

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage, signal, fftpack from skimage import filters, feature, measure, segmentation, color, exposure import nibabel as nib # 用于读取NIfTI格式医学图像 # 注意：由于我们可能没有真实的医学图像文件，这里将创建模拟数据 # 在实际应用中，应该使用真实的医学图像 def create_simulated_brain_mri(): """创建模拟的脑部MRI数据""" # 创建3D体积 size = 128 x, y, z = np.ogrid[-size//2:size//2, -size//2:size//2, -size//2:size//2] # 创建脑部形状（椭球） brain = (x**2 / (size/2.5)**2 + y**2 / (size/2.5)**2 + z**2 / (size/3)**2) <= 1 # 创建脑室（较小的椭球） ventricle = (x**2 / (size/8)**2 + y**2 / (size/8)**2 + z**2 / (size/6)**2) <= 1 # 创建肿瘤（不规则形状） tumor = ((x-size/4)**2 / (size/10)**2 + (y-size/4)**2 / (size/10)**2 + (z-size/4)**2 / (size/10)**2) <= 1 # 创建灰质和白质 gray_matter = brain & ~ventricle & ~tumor white_matter = (x**2 / (size/3)**2 + y**2 / (size/3)**2 + z**2 / (size/3.5)**2) <= 1 white_matter = white_matter & ~ventricle & ~tumor # 创建MRI数据 mri_data = np.zeros((size, size, size)) # 设置不同组织的强度值 mri_data[brain] = 30 # 脑脊液 mri_data[gray_matter] = 70 # 灰质 mri_data[white_matter] = 110 # 白质 mri_data[ventricle] = 20 # 脑室 mri_data[tumor] = 150 # 肿瘤 # 添加噪声 mri_data += np.random.normal(0, 5, mri_data.shape) # 应用平滑 mri_data = ndimage.gaussian_filter(mri_data, sigma=1.0) return mri_data # 创建模拟脑部MRI数据 brain_mri = create_simulated_brain_mri() # 显示中间切片 mid_slice = brain_mri[:, :, brain_mri.shape[2]//2] plt.figure(figsize=(10, 8)) plt.imshow(mid_slice.T, cmap='gray', origin='lower') plt.title('模拟脑部MRI中间切片') plt.colorbar() plt.show() # 1. 图像预处理 # 1.1 强度标准化 def intensity_normalization(image): """强度标准化""" min_val = np.min(image) max_val = np.max(image) return (image - min_val) / (max_val - min_val) normalized_mri = intensity_normalization(brain_mri) # 1.2 噪声去除 denoised_mri = ndimage.median_filter(normalized_mri, size=3) # 1.3 对比度增强 enhanced_mri = exposure.equalize_adapthist(denoised_mri, clip_limit=0.03) # 显示预处理结果 mid_slice_original = normalized_mri[:, :, normalized_mri.shape[2]//2] mid_slice_denoised = denoised_mri[:, :, denoised_mri.shape[2]//2] mid_slice_enhanced = enhanced_mri[:, :, enhanced_mri.shape[2]//2] plt.figure(figsize=(15, 5)) plt.subplot(131) plt.imshow(mid_slice_original.T, cmap='gray', origin='lower') plt.title('原始图像') plt.subplot(132) plt.imshow(mid_slice_denoised.T, cmap='gray', origin='lower') plt.title('去噪后') plt.subplot(133) plt.imshow(mid_slice_enhanced.T, cmap='gray', origin='lower') plt.title('对比度增强后') plt.tight_layout() plt.show() # 2. 图像分割 # 2.1 使用多阈值分割 def multi_threshold_segmentation(image, thresholds): """多阈值分割""" segmented = np.zeros_like(image, dtype=np.uint8) for i, thresh in enumerate(thresholds): if i == 0: segmented[image <= thresh] = i + 1 else: segmented[(image > thresholds[i-1]) & (image <= thresh)] = i + 1 segmented[image > thresholds[-1]] = len(thresholds) + 1 return segmented # 计算阈值 thresholds = filters.threshold_multiotsu(mid_slice_enhanced, classes=4) segmented_slice = multi_threshold_segmentation(mid_slice_enhanced, thresholds) # 2.2 区域生长分割 def region_growing_3d(image, seed, threshold): """3D区域生长算法""" # 创建与输入图像大小相同的标记矩阵 dims = image.shape segmented = np.zeros(dims, dtype=np.uint8) # 获取种子点值 seed_value = image[seed[0], seed[1], seed[2]] # 创建队列并将种子点加入 queue = [seed] segmented[seed[0], seed[1], seed[2]] = 1 # 定义6邻域（3D） neighbors = [ (-1,0,0), (1,0,0), (0,-1,0), (0,1,0), (0,0,-1), (0,0,1) ] while queue: current_point = queue.pop(0) # 检查所有邻域 for d in neighbors: r = current_point[0] + d[0] c = current_point[1] + d[1] z = current_point[2] + d[2] # 检查边界 if (r < 0 or r >= dims[0] or c < 0 or c >= dims[1] or z < 0 or z >= dims[2]): continue # 如果像素未被访问且与种子值相似 if (segmented[r, c, z] == 0 and abs(image[r, c, z] - seed_value) < threshold): segmented[r, c, z] = 1 queue.append((r, c, z)) return segmented # 选择种子点（肿瘤区域） seed_point = (brain_mri.shape[0]//2 + brain_mri.shape[0]//4, brain_mri.shape[1]//2 + brain_mri.shape[1]//4, brain_mri.shape[2]//2 + brain_mri.shape[2]//4) # 应用区域生长 tumor_segmented = region_growing_3d(enhanced_mri, seed_point, 0.1) # 显示分割结果 plt.figure(figsize=(15, 5)) plt.subplot(131) plt.imshow(mid_slice_enhanced.T, cmap='gray', origin='lower') plt.title('预处理后图像') plt.subplot(132) plt.imshow(segmented_slice.T, cmap='jet', origin='lower') plt.title('多阈值分割') plt.subplot(133) plt.imshow(tumor_segmented[:, :, brain_mri.shape[2]//2].T, cmap='gray', origin='lower') plt.title('肿瘤区域生长分割') plt.tight_layout() plt.show() # 3. 特征提取 # 3.1 形态学特征 label_image = measure.label(tumor_segmented) regions = measure.regionprops(label_image) if regions: region = regions[0] print("肿瘤形态学特征:") print(f" 体积: {region.area} 体素") print(f" 表面积: {region.perimeter} 像素") print(f" 等效直径: {region.equivalent_diameter:.2f} 体素") print(f" 离心率: {region.eccentricity:.4f}") print(f" 长轴长度: {region.major_axis_length:.2f} 体素") print(f" 短轴长度: {region.minor_axis_length:.2f} 体素") print(f" 方向: {region.orientation:.4f} 弧度") print(f" 质心: ({region.centroid[0]:.2f}, {region.centroid[1]:.2f}, {region.centroid[2]:.2f})") else: print("未检测到肿瘤区域") # 3.2 强度特征 if regions: region = regions[0] intensity_values = enhanced_mri[tumor_segmented > 0] print("n肿瘤强度特征:") print(f" 平均强度: {region.mean_intensity:.4f}") print(f" 最大强度: {region.max_intensity:.4f}") print(f" 最小强度: {region.min_intensity:.4f}") print(f" 强度标准差: {region.intensity_std:.4f}") print(f" 强度偏度: {stats.skew(intensity_values):.4f}") print(f" 强度峰度: {stats.kurtosis(intensity_values):.4f}") # 3.3 纹理特征 def glcm_features_3d(image, mask, distances=[1], angles=[0, np.pi/4, np.pi/2, 3*np.pi/4], levels=64): """计算3D图像的灰度共生矩阵特征""" # 提取ROI区域 roi = image[mask > 0] # 量化强度值 min_val = np.min(roi) max_val = np.max(roi) roi_quantized = ((roi - min_val) / (max_val - min_val) * (levels - 1)).astype(np.uint8) # 创建2D切片用于GLCM计算 # 这里简化处理，只计算中间切片 mid_slice = image.shape[2] // 2 slice_image = image[:, :, mid_slice] slice_mask = mask[:, :, mid_slice] if np.sum(slice_mask) == 0: return {} # 量化切片 slice_roi = slice_image[slice_mask > 0] slice_min = np.min(slice_roi) slice_max = np.max(slice_roi) slice_quantized = ((slice_image - slice_min) / (slice_max - slice_min) * (levels - 1)).astype(np.uint8) # 计算GLCM glcm = feature.greycomatrix(slice_quantized, distances=distances, angles=angles, levels=levels, symmetric=True, normed=True) # 计算特征 features = {} properties = ['contrast', 'dissimilarity', 'homogeneity', 'energy', 'correlation'] for prop in properties: features[prop] = feature.greycoprops(glcm, prop).mean() return features # 计算纹理特征 texture_features = glcm_features_3d(enhanced_mri, tumor_segmented) if texture_features: print("n肿瘤纹理特征:") for feat_name, feat_value in texture_features.items(): print(f" {feat_name}: {feat_value:.4f}") # 3.4 频域特征 def frequency_features_3d(image, mask): """计算3D图像的频域特征""" # 提取ROI区域 roi = image[mask > 0] if len(roi) == 0: return {} # 计算ROI的傅里叶变换 # 这里简化处理，只计算中间切片 mid_slice = image.shape[2] // 2 slice_image = image[:, :, mid_slice] slice_mask = mask[:, :, mid_slice] if np.sum(slice_mask) == 0: return {} # 提取ROI区域 roi_slice = slice_image * slice_mask # 计算傅里叶变换 f_transform = fftpack.fft2(roi_slice) f_shift = fftpack.fftshift(f_transform) magnitude_spectrum = np.abs(f_shift) # 计算径向分布 rows, cols = roi_slice.shape crow, ccol = rows // 2, cols // 2 # 创建径向距离矩阵 y, x = np.ogrid[:rows, :cols] r = np.sqrt((x - ccol) ** 2 + (y - crow) ** 2) r = r.astype(int) # 计算径向平均 tbin = np.bincount(r.ravel(), magnitude_spectrum.ravel()) nr = np.bincount(r.ravel()) radial_profile = tbin / (nr + 1e-10) # 避免除以零 # 计算特征 features = { 'dc_component': magnitude_spectrum[crow, ccol], 'spectral_energy': np.sum(radial_profile ** 2), 'spectral_entropy': -np.sum(radial_profile * np.log2(radial_profile + 1e-10)) } return features, radial_profile # 计算频域特征 freq_features, radial_profile = frequency_features_3d(enhanced_mri, tumor_segmented) if freq_features: print("n肿瘤频域特征:") for feat_name, feat_value in freq_features.items(): print(f" {feat_name}: {feat_value:.4f}") # 显示径向分布 plt.figure(figsize=(10, 6)) plt.plot(radial_profile) plt.title('肿瘤频谱径向分布') plt.xlabel('频率半径') plt.ylabel('平均幅度') plt.grid(True) plt.show() # 4. 结果可视化 # 创建3D渲染的简化版本（使用多个切片） plt.figure(figsize=(15, 10)) for i in range(9): slice_idx = int(i * brain_mri.shape[2] / 8) plt.subplot(3, 3, i+1) plt.imshow(enhanced_mri[:, :, slice_idx].T, cmap='gray', origin='lower') plt.contour(tumor_segmented[:, :, slice_idx].T, levels=[0.5], colors='r', linewidths=1, origin='lower') plt.title(f'切片 {slice_idx}') plt.tight_layout() plt.suptitle('肿瘤分割结果（红色轮廓）', y=1.02) plt.show()

Scipy与其他库的集成

Scipy虽然功能强大，但在医学图像处理中，通常需要与其他库结合使用，以构建完整的处理流程。

与NumPy的集成

NumPy是Python科学计算的基础库，Scipy构建在NumPy之上，两者无缝集成。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt # 创建一个示例图像 x = np.linspace(-5, 5, 500) y = np.linspace(-5, 5, 500) xx, yy = np.meshgrid(x, y) image = np.sin(xx**2 + yy**2) / (xx**2 + yy**2 + 1e-6) # 使用Scipy进行滤波 filtered_image = ndi.gaussian_filter(image, sigma=1.0) # 计算梯度 grad_y, grad_x = np.gradient(image) grad_mag = np.sqrt(grad_x**2 + grad_y**2) # 显示结果 plt.figure(figsize=(15, 5)) plt.subplot(131) plt.imshow(image, cmap='viridis') plt.title('原始图像') plt.subplot(132) plt.imshow(filtered_image, cmap='viridis') plt.title('高斯滤波后') plt.subplot(133) plt.imshow(grad_mag, cmap='viridis') plt.title('梯度幅值') plt.tight_layout() plt.show()

与Matplotlib的集成

Matplotlib是Python中最常用的绘图库，与Scipy结合使用可以可视化处理结果。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 创建一个示例3D图像 x, y, z = np.ogrid[-5:5:50j, -5:5:50j, -5:5:50j] sphere = x**2 + y**2 + z**2 <= 4**2 ellipsoid = x**2/9 + y**2/4 + z**2 <= 1 # 创建3D图像 image = sphere.astype(float) + 0.5 * ellipsoid.astype(float) image += 0.1 * np.random.randn(*image.shape) # 使用Scipy进行3D滤波 filtered_image = ndi.gaussian_filter(image, sigma=1.0) # 提取等值面 def plot_isosurface(ax, data, level, color='b', alpha=0.5): """绘制等值面""" verts, faces, _, _ = measure.marching_cubes(data, level=level) mesh = Poly3DCollection(verts[faces], alpha=alpha) mesh.set_facecolor(color) ax.add_collection3d(mesh) # 可视化3D图像 fig = plt.figure(figsize=(15, 5)) # 原始图像的切片 ax1 = fig.add_subplot(131) ax1.imshow(image[:, :, 25], cmap='gray') ax1.set_title('原始图像中间切片') # 滤波后图像的切片 ax2 = fig.add_subplot(132) ax2.imshow(filtered_image[:, :, 25], cmap='gray') ax2.set_title('滤波后图像中间切片') # 3D渲染 ax3 = fig.add_subplot(133, projection='3d') ax3.set_title('3D渲染') # 由于marching_cubes在较新版本的skimage中位置不同，这里使用简化的可视化 # 创建3D点云 x, y, z = np.where(image > 0.5) ax3.scatter(x, y, z, c=image[image > 0.5], cmap='viridis', alpha=0.1) ax3.set_xlim(0, image.shape[0]) ax3.set_ylim(0, image.shape[1]) ax3.set_zlim(0, image.shape[2]) plt.tight_layout() plt.show()

与Scikit-image的集成

Scikit-image是一个专门用于图像处理的库，与Scipy互补，提供了更多高级图像处理功能。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt from skimage import data, filters, feature, measure, segmentation, color # 加载图像 image = data.brain() # 如果图像是彩色的，转换为灰度 if len(image.shape) > 2: image_gray = color.rgb2gray(image) else: image_gray = image.copy() # 使用Scipy进行预处理 denoised = ndi.median_filter(image_gray, size=3) enhanced = denoised - ndi.gaussian_filter(denoised, sigma=3) # 使用Scikit-image进行边缘检测 edges = feature.canny(enhanced, sigma=1.0) # 使用Scikit-image进行分割 # 标记分水岭 markers = np.zeros_like(image_gray) markers[image_gray < 0.3] = 1 markers[image_gray > 0.7] = 2 # 应用分水岭算法 segmentation_result = segmentation.watershed(enhanced, markers) # 使用Scipy进行后处理 # 移除小对象 label_image = measure.label(segmentation_result) region_sizes = np.bincount(label_image.ravel()) small_regions = region_sizes < 100 small_regions[0] = False # 不移除背景 cleaned = small_regions[label_image] segmentation_result[cleaned] = 0 # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(image_gray, cmap='gray') plt.title('原始图像') plt.subplot(222) plt.imshow(enhanced, cmap='gray') plt.title('增强后图像') plt.subplot(223) plt.imshow(edges, cmap='gray') plt.title('边缘检测') plt.subplot(224) plt.imshow(segmentation_result, cmap='jet') plt.title('分割结果') plt.tight_layout() plt.show()

与SimpleITK的集成

SimpleITK是一个专门用于医学图像处理的库，与Scipy结合使用可以处理DICOM等医学图像格式。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt # 注意：以下代码需要安装SimpleITK库 # 由于可能没有安装，这里只提供示例代码结构 try: import SimpleITK as sitk # 读取DICOM图像 # 注意：需要替换为实际的DICOM文件路径 # dicom_image = sitk.ReadImage("path_to_dicom_file.dcm") # 转换为NumPy数组 # image_array = sitk.GetArrayFromImage(dicom_image) # 使用Scipy进行处理 # processed_array = ndi.gaussian_filter(image_array, sigma=1.0) # 转换回SimpleITK图像 # processed_image = sitk.GetImageFromArray(processed_array) # processed_image.CopyInformation(dicom_image) # 保存结果 # sitk.WriteImage(processed_image, "processed_image.dcm") print("SimpleITK集成示例代码（需要实际DICOM文件）") except ImportError: print("SimpleITK未安装，无法运行示例代码") # 创建模拟数据 image_array = np.random.rand(100, 100) # 使用Scipy进行处理 processed_array = ndi.gaussian_filter(image_array, sigma=1.0) # 显示结果 plt.figure(figsize=(10, 5)) plt.subplot(121) plt.imshow(image_array, cmap='gray') plt.title('模拟原始图像') plt.subplot(122) plt.imshow(processed_array, cmap='gray') plt.title('处理后图像') plt.tight_layout() plt.show()

创新应用与前沿发展

Scipy在医学图像处理领域的应用不断拓展，以下是一些创新应用和前沿发展方向。

深度学习与医学图像处理

深度学习技术在医学图像处理中取得了突破性进展，Scipy可以与深度学习框架（如TensorFlow、PyTorch）结合使用，构建端到端的医学图像分析系统。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt # 注意：以下代码需要安装TensorFlow或PyTorch # 由于可能没有安装，这里只提供示例代码结构 try: import tensorflow as tf from tensorflow.keras import layers, models # 创建一个简单的U-Net模型用于图像分割 def unet_model(input_size=(128, 128, 1)): inputs = layers.Input(input_size) # 编码器 conv1 = layers.Conv2D(32, 3, activation='relu', padding='same')(inputs) conv1 = layers.Conv2D(32, 3, activation='relu', padding='same')(conv1) pool1 = layers.MaxPooling2D(pool_size=(2, 2))(conv1) conv2 = layers.Conv2D(64, 3, activation='relu', padding='same')(pool1) conv2 = layers.Conv2D(64, 3, activation='relu', padding='same')(conv2) pool2 = layers.MaxPooling2D(pool_size=(2, 2))(conv2) # 中间层 conv3 = layers.Conv2D(128, 3, activation='relu', padding='same')(pool2) conv3 = layers.Conv2D(128, 3, activation='relu', padding='same')(conv3) # 解码器 up4 = layers.UpSampling2D(size=(2, 2))(conv3) up4 = layers.concatenate([up4, conv2], axis=-1) conv4 = layers.Conv2D(64, 3, activation='relu', padding='same')(up4) conv4 = layers.Conv2D(64, 3, activation='relu', padding='same')(conv4) up5 = layers.UpSampling2D(size=(2, 2))(conv4) up5 = layers.concatenate([up5, conv1], axis=-1) conv5 = layers.Conv2D(32, 3, activation='relu', padding='same')(up5) conv5 = layers.Conv2D(32, 3, activation='relu', padding='same')(conv5) outputs = layers.Conv2D(1, 1, activation='sigmoid')(conv5) model = models.Model(inputs=inputs, outputs=outputs) return model # 创建模型 model = unet_model() model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy']) model.summary() print("U-Net模型创建成功（需要TensorFlow）") except ImportError: print("TensorFlow未安装，无法创建深度学习模型") # 创建模拟数据展示传统方法 # 创建一个简单的图像分割示例 image = np.zeros((128, 128)) # 添加一些形状 x, y = np.ogrid[:128, :128] mask1 = (x - 40)**2 + (y - 40)**2 < 20**2 mask2 = (x - 80)**2 + (y - 80)**2 < 15**2 mask3 = ((x - 60) > 30) & ((x - 60) < 70) & ((y - 60) > 30) & ((y - 60) < 70) image[mask1] = 1 image[mask2] = 1 image[mask3] = 1 # 添加噪声 noisy_image = image + 0.1 * np.random.randn(*image.shape) noisy_image = np.clip(noisy_image, 0, 1) # 使用Scipy进行分割 # 阈值分割 threshold = 0.5 binary = noisy_image > threshold # 形态学操作 eroded = ndi.binary_erosion(binary) dilated = ndi.binary_dilation(binary) opened = ndi.binary_opening(binary) closed = ndi.binary_closing(binary) # 显示结果 plt.figure(figsize=(15, 10)) plt.subplot(231) plt.imshow(image, cmap='gray') plt.title('原始图像') plt.subplot(232) plt.imshow(noisy_image, cmap='gray') plt.title('添加噪声') plt.subplot(233) plt.imshow(binary, cmap='gray') plt.title('阈值分割') plt.subplot(234) plt.imshow(eroded, cmap='gray') plt.title('腐蚀') plt.subplot(235) plt.imshow(dilated, cmap='gray') plt.title('膨胀') plt.subplot(236) plt.imshow(opened, cmap='gray') plt.title('开运算') plt.tight_layout() plt.show()

多模态医学图像融合

多模态医学图像融合是将不同成像方式（如CT、MRI、PET等）的图像结合，以提供更全面的诊断信息。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt from skimage import exposure, transform # 创建模拟的多模态医学图像 def create_simulated_ct_mri_pet(size=128): """创建模拟的CT、MRI和PET图像""" x, y, z = np.ogrid[-size//2:size//2, -size//2:size//2, -size//2:size//2] # 创建基本形状 # 外部轮廓（头部） head = (x**2 / (size/2.5)**2 + y**2 / (size/2.5)**2 + z**2 / (size/3)**2) <= 1 # 脑部 brain = (x**2 / (size/3)**2 + y**2 / (size/3)**2 + z**2 / (size/3.5)**2) <= 1 # 脑室 ventricle = (x**2 / (size/8)**2 + y**2 / (size/8)**2 + z**2 / (size/6)**2) <= 1 # 肿瘤 tumor = ((x-size/4)**2 / (size/10)**2 + (y-size/4)**2 / (size/10)**2 + (z-size/4)**2 / (size/10)**2) <= 1 # 骨骼（头部边缘） skull = head & ~brain # 创建CT图像（骨骼高密度） ct = np.zeros((size, size, size)) ct[skull] = 200 # 骨骼 ct[brain] = 40 # 脑组织 ct[ventricle] = 10 # 脑室 ct[tumor] = 60 # 肿瘤 # 创建MRI图像（软组织对比度高） mri = np.zeros((size, size, size)) mri[skull] = 20 # 骨骼（在MRI中信号低） mri[brain] = 100 # 脑组织 mri[ventricle] = 30 # 脑室（脑脊液） mri[tumor] = 150 # 肿瘤 # 创建PET图像（代谢活动） pet = np.zeros((size, size, size)) pet[skull] = 5 # 骨骼（低代谢） pet[brain] = 20 # 脑组织（基础代谢） pet[ventricle] = 5 # 脑室（低代谢） pet[tumor] = 100 # 肿瘤（高代谢） # 添加噪声 ct += np.random.normal(0, 5, ct.shape) mri += np.random.normal(0, 5, mri.shape) pet += np.random.normal(0, 2, pet.shape) # 应用平滑 ct = ndi.gaussian_filter(ct, sigma=1.0) mri = ndi.gaussian_filter(mri, sigma=1.0) pet = ndi.gaussian_filter(pet, sigma=2.0) return ct, mri, pet # 创建模拟数据 ct_image, mri_image, pet_image = create_simulated_ct_mri_pet() # 显示中间切片 mid_slice = ct_image.shape[2] // 2 plt.figure(figsize=(15, 5)) plt.subplot(131) plt.imshow(ct_image[:, :, mid_slice].T, cmap='gray', origin='lower') plt.title('CT图像') plt.colorbar() plt.subplot(132) plt.imshow(mri_image[:, :, mid_slice].T, cmap='gray', origin='lower') plt.title('MRI图像') plt.colorbar() plt.subplot(133) plt.imshow(pet_image[:, :, mid_slice].T, cmap='hot', origin='lower') plt.title('PET图像') plt.colorbar() plt.tight_layout() plt.show() # 图像配准（这里简化处理，假设图像已经对齐） # 在实际应用中，需要使用更复杂的配准算法 # 多模态图像融合 def image_fusion(image1, image2, alpha=0.5): """简单的图像融合""" # 归一化图像 norm1 = (image1 - np.min(image1)) / (np.max(image1) - np.min(image1)) norm2 = (image2 - np.min(image2)) / (np.max(image2) - np.min(image2)) # 融合 fused = alpha * norm1 + (1 - alpha) * norm2 return fused # 融合CT和MRI ct_mri_fused = image_fusion(ct_image, mri_image, 0.5) # 融合MRI和PET mri_pet_fused = image_fusion(mri_image, pet_image, 0.7) # 融合CT、MRI和PET ct_mri_pet_fused = (ct_image/np.max(ct_image) + mri_image/np.max(mri_image) + pet_image/np.max(pet_image)) / 3 # 显示融合结果 plt.figure(figsize=(15, 10)) plt.subplot(221) plt.imshow(ct_mri_fused[:, :, mid_slice].T, cmap='gray', origin='lower') plt.title('CT-MRI融合') plt.colorbar() plt.subplot(222) plt.imshow(mri_pet_fused[:, :, mid_slice].T, cmap='hot', origin='lower') plt.title('MRI-PET融合') plt.colorbar() plt.subplot(223) plt.imshow(ct_mri_pet_fused[:, :, mid_slice].T, cmap='viridis', origin='lower') plt.title('CT-MRI-PET融合') plt.colorbar() plt.subplot(224) plt.imshow(mri_image[:, :, mid_slice].T, cmap='gray', origin='lower') plt.contour(pet_image[:, :, mid_slice].T, levels=[np.max(pet_image)*0.5], colors='r', linewidths=1, origin='lower') plt.contour(ct_image[:, :, mid_slice].T, levels=[np.max(ct_image)*0.7], colors='g', linewidths=1, origin='lower') plt.title('MRI with PET (red) and CT (green) contours') plt.tight_layout() plt.show() # 伪彩色融合 def false_color_fusion(ct, mri, pet): """创建伪彩色融合图像""" # 归一化 ct_norm = (ct - np.min(ct)) / (np.max(ct) - np.min(ct)) mri_norm = (mri - np.min(mri)) / (np.max(mri) - np.min(mri)) pet_norm = (pet - np.min(pet)) / (np.max(pet) - np.min(pet)) # 创建RGB图像 # 红色通道：CT # 绿色通道：MRI # 蓝色通道：PET rgb_image = np.zeros((ct.shape[0], ct.shape[1], ct.shape[2], 3)) rgb_image[:, :, :, 0] = ct_norm rgb_image[:, :, :, 1] = mri_norm rgb_image[:, :, :, 2] = pet_norm return rgb_image # 创建伪彩色融合图像 false_color = false_color_fusion(ct_image, mri_image, pet_image) # 显示伪彩色融合结果 plt.figure(figsize=(10, 8)) plt.imshow(false_color[:, :, mid_slice].transpose(1, 0, 2), origin='lower') plt.title('伪彩色融合 (R:CT, G:MRI, B:PET)') plt.show()

医学图像可视化与交互分析

医学图像的可视化和交互分析是医学图像处理的重要组成部分，Scipy可以与可视化库结合使用，创建交互式分析工具。

import numpy as np import scipy.ndimage as ndi import matplotlib.pyplot as plt from matplotlib.widgets import Slider, Button from mpl_toolkits.mplot3d import Axes3D # 创建模拟的3D医学图像数据 def create_simulated_medical_volume(size=64): """创建模拟的3D医学图像数据""" x, y, z = np.ogrid[-size//2:size//2, -size//2:size//2, -size//2:size//2] # 创建多个结构 # 外部轮廓 outer = (x**2 / (size/2.5)**2 + y**2 / (size/2.5)**2 + z**2 / (size/3)**2) <= 1 # 内部结构1 inner1 = (x**2 / (size/4)**2 + y**2 / (size/4)**2 + z**2 / (size/5)**2) <= 1 # 内部结构2 inner2 = ((x-size/4)**2 / (size/8)**2 + (y-size/4)**2 / (size/8)**2 + (z-size/4)**2 / (size/8)**2) <= 1 # 创建体积数据 volume = np.zeros((size, size, size)) volume[outer] = 30 volume[inner1] = 80 volume[inner2] = 120 # 添加噪声 volume += np.random.normal(0, 3, volume.shape) # 应用平滑 volume = ndi.gaussian_filter(volume, sigma=1.0) return volume # 创建模拟数据 volume = create_simulated_medical_volume() # 创建交互式切片查看器 class SliceViewer: def __init__(self, volume): self.volume = volume self.current_slice = volume.shape[2] // 2 self.fig, self.ax = plt.subplots(figsize=(8, 8)) plt.subplots_adjust(bottom=0.25) # 显示初始切片 self.im = self.ax.imshow(volume[:, :, self.current_slice].T, cmap='gray', origin='lower') self.ax.set_title(f'切片 {self.current_slice}') # 添加滑块 ax_slice = plt.axes([0.25, 0.1, 0.65, 0.03]) self.slice_slider = Slider( ax=ax_slice, label='切片', valmin=0, valmax=volume.shape[2]-1, valinit=self.current_slice, valstep=1 ) # 添加按钮 ax_prev = plt.axes([0.1, 0.1, 0.1, 0.04]) ax_next = plt.axes([0.8, 0.1, 0.1, 0.04]) self.btn_prev = Button(ax_prev, '上一个') self.btn_next = Button(ax_next, '下一个') # 连接事件 self.slice_slider.on_changed(self.update_slice) self.btn_prev.on_clicked(self.prev_slice) self.btn_next.on_clicked(self.next_slice) def update_slice(self, val): self.current_slice = int(self.slice_slider.val) self.im.set_data(self.volume[:, :, self.current_slice].T) self.ax.set_title(f'切片 {self.current_slice}') self.fig.canvas.draw_idle() def prev_slice(self, event): if self.current_slice > 0: self.current_slice -= 1 self.slice_slider.set_val(self.current_slice) def next_slice(self, event): if self.current_slice < self.volume.shape[2] - 1: self.current_slice += 1 self.slice_slider.set_val(self.current_slice) # 创建并显示切片查看器 viewer = SliceViewer(volume) plt.show() # 创建多平面重建（MPR）视图 fig = plt.figure(figsize=(15, 5)) # 轴位视图 ax1 = fig.add_subplot(131) ax1.imshow(volume[:, :, volume.shape[2]//2].T, cmap='gray', origin='lower') ax1.set_title('轴位视图') # 矢状位视图 ax2 = fig.add_subplot(132) ax2.imshow(volume[:, volume.shape[1]//2, :].T, cmap='gray', origin='lower') ax2.set_title('矢状位视图') # 冠状位视图 ax3 = fig.add_subplot(133) ax3.imshow(volume[volume.shape[0]//2, :, :].T, cmap='gray', origin='lower') ax3.set_title('冠状位视图') plt.tight_layout() plt.show() # 创建3D体积渲染 fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 创建等值面 def plot_isosurface(ax, data, level, color='b', alpha=0.5): """绘制等值面""" # 简化的等值面绘制（使用点云） x, y, z = np.where(data > level) ax.scatter(x, y, z, c=data[data > level], cmap='viridis', alpha=alpha) # 绘制不同强度的等值面 plot_isosurface(ax, volume, level=50, color='b', alpha=0.1) plot_isosurface(ax, volume, level=100, color='g', alpha=0.2) ax.set_xlim(0, volume.shape[0]) ax.set_ylim(0, volume.shape[1]) ax.set_zlim(0, volume.shape[2]) ax.set_title('3D体积渲染') plt.show() # 创建最大强度投影（MIP） fig = plt.figure(figsize=(15, 5)) # 轴位MIP ax1 = fig.add_subplot(131) mip_axial = np.max(volume, axis=2) ax1.imshow(mip_axial.T, cmap='gray', origin='lower') ax1.set_title('轴位MIP') # 矢状位MIP ax2 = fig.add_subplot(132) mip_sagittal = np.max(volume, axis=1) ax2.imshow(mip_sagittal.T, cmap='gray', origin='lower') ax2.set_title('矢状位MIP') # 冠状位MIP ax3 = fig.add_subplot(133) mip_coronal = np.max(volume, axis=0) ax3.imshow(mip_coronal.T, cmap='gray', origin='lower') ax3.set_title('冠状位MIP') plt.tight_layout() plt.show()

总结与展望

本文详细介绍了使用Scipy进行医学图像处理的完整工作流程，从图像预处理到特征提取，并探讨了Python科学计算库在医学影像领域的创新应用。通过丰富的代码示例和实际案例，我们展示了Scipy在医学图像处理中的强大功能和灵活性。

主要成果总结

图像预处理技术：我们介绍了多种图像预处理方法，包括图像读取与显示、对比度调整、图像滤波和图像配准。这些技术为后续的图像分析奠定了基础。
图像分割技术：我们详细讨论了多种分割方法，包括阈值分割、区域生长、边缘检测和基于机器学习的分割。这些方法可以有效地从医学图像中提取感兴趣的区域。
特征提取技术：我们介绍了多种特征提取方法，包括形态学特征、纹理特征、强度特征和频域特征。这些特征可以用于医学图像的分类、检测和诊断。
实际案例分析：通过一个完整的脑部MRI图像处理案例，我们展示了从预处理到特征提取的完整工作流程。
库集成与创新应用：我们讨论了Scipy与其他库（如NumPy、Matplotlib、Scikit-image等）的集成方法，并探讨了深度学习、多模态图像融合和交互式可视化等前沿应用。