影像组学实战5大纹理特征Python代码全解与乳腺癌影像分析如果你是一名医学影像处理的新手或者是一位希望将数据分析能力扩展到医学领域的Python开发者那么“纹理特征”这个词对你来说可能既熟悉又陌生。熟悉的是它在图像处理领域无处不在陌生的是当面对一张CT或MRI影像想要从中提取出有临床价值的量化信息时却不知从何下手。理论文献往往堆砌着复杂的公式和术语而实际操作的“最后一公里”——如何用代码实现这些特征的计算——却常常语焉不详。这篇文章就是为你准备的。我们不打算重复教科书上的定义而是直接切入核心如何用Python一行行代码从零开始计算并理解影像组学中最关键的5种纹理特征。我们将以乳腺癌影像分析为实战案例带你走过从数据加载、预处理、特征提取到结果解读的完整流程。你会发现那些听起来高深的GLCM、GLRLM其背后的计算逻辑清晰而优雅而Python正是实现它们的最佳工具。准备好你的开发环境我们开始吧。1. 环境搭建与数据准备迈出第一步在开始编写任何特征提取代码之前一个稳定、高效的开发环境是基石。对于医学影像处理我们不仅需要通用的科学计算库还需要专门处理医学图像格式和进行高级图像操作的库。首先我们通过pip安装核心依赖包。建议使用虚拟环境来管理依赖避免包冲突。pip install numpy pandas scikit-learn scikit-image matplotlib seaborn pip install SimpleITK # 用于读取DICOM、NIfTI等医学影像格式 pip install pyradiomics # 影像组学特征提取的权威库我们将参考其实现逻辑接下来我们准备数据。假设我们有一组乳腺癌患者的动态增强磁共振成像DCE-MRI数据存储为NIfTI格式.nii或.nii.gz同时每个患者都有一个对应的肿瘤区域分割掩膜Segmentation Mask通常由放射科医生勾画也存储为NIfTI格式。我们将使用SimpleITK来读取它们。import SimpleITK as sitk import numpy as np import matplotlib.pyplot as plt # 加载原始影像和分割掩膜 image_path “path/to/patient_image.nii.gz” mask_path “path/to/tumor_mask.nii.gz” image_sitk sitk.ReadImage(image_path) mask_sitk sitk.ReadImage(mask_path) # 转换为NumPy数组以便后续处理 image_array sitk.GetArrayFromImage(image_sitk) # 形状通常为 (深度, 高度, 宽度) mask_array sitk.GetArrayFromImage(mask_sitk) # 检查数据和掩膜是否对齐 print(f“影像形状: {image_array.shape}”) print(f“掩膜形状: {mask_array.shape}”) print(f“掩膜中肿瘤体素数量: {np.sum(mask_array 0)}”)注意医学影像数据通常涉及患者隐私公开数据集如TCIAThe Cancer Imaging Archive提供了匿名化的数据供研究使用。在实战中请确保你使用的数据已获得合规授权。数据加载后通常需要进行预处理例如灰度归一化。影像组学特征对灰度强度标度敏感不同扫描设备或协议产生的强度差异可能主导特征值掩盖真实的生物学纹理信息。一种常见的方法是将肿瘤区域ROI的灰度值缩放到一个标准范围比如0到2558位或0到409512位具体取决于原始数据的位深。def normalize_roi_intensity(roi_voxels): “”“将ROI内的体素强度进行归一化到[0, 255]区间。”“” roi_min roi_voxels.min() roi_max roi_voxels.max() # 避免除零错误 if roi_max roi_min: return np.zeros_like(roi_voxels) normalized (roi_voxels - roi_min) / (roi_max - roi_min) * 255 return normalized.astype(np.uint8) # 转换为8位整数 # 提取肿瘤区域体素 tumor_voxels image_array[mask_array 0] normalized_tumor_voxels normalize_roi_intensity(tumor_voxels) # 将归一化后的值放回原图像数组中仅用于特征提取 image_normalized image_array.copy().astype(np.float32) image_normalized[mask_array 0] normalized_tumor_voxels至此我们完成了环境和数据的准备工作。接下来我们将深入每一种纹理特征并附上可运行的Python实现。2. 灰度共生矩阵GLCM捕捉像素对的空间关系灰度共生矩阵是影像组学中最经典、使用最广泛的纹理特征之一。它的核心思想不是看单个像素的灰度而是研究具有特定空间关系的像素对例如右边1个像素下方1个像素的灰度组合出现的频率。这种“二阶统计”能有效描述纹理的粗糙度、对比度和规律性。计算GLCM本身是一个构建联合概率矩阵的过程。假设我们有一幅归一化后灰度级为8级0-7的图像我们关注像素与其右侧相邻像素距离d1方向θ0°的关系。我们遍历图像中所有这样的像素对统计灰度值i和j同时出现的次数填入矩阵P(i, j)中。最后通常将矩阵归一化使得所有元素之和为1表示概率。下面是一个简化但完整的GLCM计算函数并从中提取4个关键特征对比度、相关性、能量角二阶矩和同质性。import numpy as np from scipy import ndimage def compute_glcm_features(roi_2d_slice, distance1, angle0, gray_levels8): “”” 计算2D图像片段的GLCM及其特征。 参数: roi_2d_slice: 2D numpy数组包含肿瘤区域的切片。 distance: 像素对之间的距离。 angle: 方向角度弧度0表示水平向右。 gray_levels: 量化后的灰度级数。 返回: 包含对比度、相关性、能量、同质性的字典。 “”” # 1. 图像灰度级量化 min_val, max_val roi_2d_slice.min(), roi_2d_slice.max() if max_val min_val: return None quantized ((roi_2d_slice - min_val) / (max_val - min_val) * (gray_levels - 1)).astype(np.uint8) # 2. 初始化GLCM矩阵 glcm np.zeros((gray_levels, gray_levels), dtypenp.float64) rows, cols quantized.shape # 3. 根据角度计算偏移量 angle_rad np.deg2rad(angle) if not isinstance(angle, float) else angle offset_row int(round(distance * np.sin(angle_rad))) offset_col int(round(distance * np.cos(angle_rad))) # 4. 遍历图像填充GLCM for r in range(max(0, -offset_row), min(rows, rows - offset_row)): for c in range(max(0, -offset_col), min(cols, cols - offset_col)): i quantized[r, c] j quantized[r offset_row, c offset_col] glcm[i, j] 1 # 5. 归一化GLCM为概率矩阵 if glcm.sum() 0: return None glcm / glcm.sum() # 6. 计算特征 i_idx, j_idx np.indices(glcm.shape) # 对比度: 衡量局部灰度变化值大表示纹理沟壑深 contrast np.sum(glcm * ((i_idx - j_idx) ** 2)) # 相关性: 衡量行元素和列元素的线性依赖程度 mu_i np.sum(i_idx * glcm) mu_j np.sum(j_idx * glcm) sigma_i np.sqrt(np.sum(glcm * (i_idx - mu_i) ** 2)) sigma_j np.sqrt(np.sum(glcm * (j_idx - mu_j) ** 2)) if sigma_i * sigma_j 0: correlation np.sum(glcm * (i_idx - mu_i) * (j_idx - mu_j)) / (sigma_i * sigma_j) else: correlation 0 # 能量角二阶矩: 衡量图像纹理的均匀性值大表示纹理稳定 energy np.sum(glcm ** 2) # 同质性逆差矩: 衡量GLCM中元素分布对主对角线的靠近程度值大表示纹理局部均匀 homogeneity np.sum(glcm / (1 (i_idx - j_idx) ** 2)) return { ‘contrast’: contrast, ‘correlation’: correlation, ‘energy’: energy, ‘homogeneity’: homogeneity } # 实战应用提取某一切片的GLCM特征 # 假设我们有一个2D的肿瘤切片 tumor_slice_2d features compute_glcm_features(tumor_slice_2d, distance1, angle0, gray_levels16) print(f“GLCM特征: {features}”)在实际分析中我们通常会在多个方向如0°, 45°, 90°, 135°上计算GLCM然后对每个特征取各个方向上的平均值或最大值作为最终特征值以使其具有旋转不变性。3. 灰度游程矩阵GLRLM分析相同灰度的连续游程如果说GLCM关注的是“点对”关系那么灰度游程矩阵关注的则是“线”的关系。它统计的是图像中在特定方向上具有相同灰度值的连续像素即一个“游程”的长度分布。这对于描述纹理的粗糙度、线性结构或方向性特别有用。例如一个具有长游程的图像其纹理可能更粗糙或包含更多的线性图案。GLRLM矩阵P(i, j)的定义是在给定方向上灰度级为i、游程长度为j的游程出现的次数。其中i是灰度级索引j是游程长度。计算GLRLM需要沿着指定方向扫描图像识别并测量每一个连续的、灰度相同的像素序列。以下是其Python实现及几个核心特征的计算def compute_glrlm_features(roi_2d_slice, angle0, gray_levels8): “”” 计算2D图像片段的GLRLM及其特征。 参数: roi_2d_slice: 2D numpy数组。 angle: 扫描方向0度表示水平方向0°。 gray_levels: 量化灰度级数。 返回: 包含短游程强调、长游程强调、灰度不均匀性、游程长度不均匀性、游程百分比的字典。 “”” # 图像量化 min_val, max_val roi_2d_slice.min(), roi_2d_slice.max() if max_val min_val: return None quantized ((roi_2d_slice - min_val) / (max_val - min_val) * (gray_levels - 1)).astype(np.uint8) rows, cols quantized.shape # 根据角度确定扫描顺序 if angle 0: # 水平扫描 scan_array quantized primary_dim, secondary_dim rows, cols step_primary, step_secondary 0, 1 # 沿行固定沿列移动 elif angle 90: # 垂直扫描 scan_array quantized.T primary_dim, secondary_dim cols, rows step_primary, step_secondary 0, 1 else: # 对于对角线方向需要更复杂的遍历此处为简化以水平为例 scan_array quantized primary_dim, secondary_dim rows, cols step_primary, step_secondary 0, 1 # 初始化GLRLM矩阵 max_run_length secondary_dim # 最大可能游程长度不会超过该维度 glrlm np.zeros((gray_levels, max_run_length), dtypenp.float64) # 扫描每一行或列 for p in range(primary_dim): current_gray scan_array[p, 0] if secondary_dim 0 else None run_length 0 for s in range(secondary_dim): if scan_array[p, s] current_gray: run_length 1 else: if current_gray is not None and run_length 0: glrlm[current_gray, run_length-1] 1 # 索引从0开始 current_gray scan_array[p, s] run_length 1 # 处理行末的游程 if current_gray is not None and run_length 0: glrlm[current_gray, run_length-1] 1 # 归一化可选转为概率 total_runs glrlm.sum() if total_runs 0: return None glrlm_prob glrlm / total_runs # 计算特征 i_idx, j_idx np.indices(glrlm_prob.shape) j_idx j_idx 1 # 游程长度从1开始计数 # 短游程强调 (SRE): 对短游程给予更高权重 sre np.sum(glrlm_prob / (j_idx ** 2)) / total_runs if total_runs 0 else 0 # 长游程强调 (LRE): 对长游程给予更高权重 lre np.sum(glrlm_prob * (j_idx ** 2)) / total_runs if total_runs 0 else 0 # 灰度不均匀性 (GLN): 衡量不同灰度级游程分布的均匀性 sum_over_runs glrlm_prob.sum(axis1) gln np.sum(sum_over_runs ** 2) # 游程长度不均匀性 (RLN): 衡量不同长度游程分布的均匀性 sum_over_gray glrlm_prob.sum(axis0) rln np.sum(sum_over_gray ** 2) # 游程百分比 (RP): 游程总像素数占ROI总像素数的比例 rp np.sum(glrlm * j_idx) / (primary_dim * secondary_dim) if (primary_dim * secondary_dim) 0 else 0 return { ‘short_run_emphasis’: sre, ‘long_run_emphasis’: lre, ‘gray_level_nonuniformity’: gln, ‘run_length_nonuniformity’: rln, ‘run_percentage’: rp } # 使用示例 glrlm_feats compute_glrlm_features(tumor_slice_2d, angle0, gray_levels16) print(f“GLRLM特征 (水平方向): {glrlm_feats}”)提示在实际的3D影像分析中GLRLM可以在13个空间方向在3D中进行计算然后聚合特征。pyradiomics库提供了非常成熟且高效的3D GLRLM实现生产环境中推荐直接使用。4. 灰度区域大小矩阵GLSZM与灰度依赖矩阵GLDM除了GLCM和GLRLM还有几种强大的纹理描述符。灰度区域大小矩阵关注的是“面”或“体”的关系。它不统计线性的游程而是统计在2D或3D空间中相互连接且灰度值相同的像素/体素所构成的“区域”Zone的大小。一个由许多小区域组成的纹理其GLSZM矩阵会集中在短区域长度一侧而包含大块均匀区域的纹理则会有更长的区域长度分布。灰度依赖矩阵则从一个中心像素出发考察其邻域内灰度值与其相似的像素数量即“依赖”于该中心像素的像素数。它量化了局部区域的灰度一致性。一个高“小依赖强调”值的区域意味着图像中存在许多小的、均匀的斑块而高“大依赖强调”值则表明存在大的均匀区域。由于GLSZM和GLDM的算法实现相对复杂涉及连通域分析或邻域依赖判断代码量较大。为了清晰和实用性这里我们展示如何利用业界标准的pyradiomics库来一站式提取所有这些特征这在实际项目中是最可靠和高效的做法。import SimpleITK as sitk from radiomics import featureextractor # 1. 初始化特征提取器 extractor featureextractor.RadiomicsFeatureExtractor() # 2. 查看和修改提取参数可选 print(“默认提取的特征类型:”, extractor.enabledFeatures) # 我们可以只启用纹理特征 extractor.disableAllFeatures() extractor.enableFeatureClassByName(‘glcm’) extractor.enableFeatureClassByName(‘glrlm’) extractor.enableFeatureClassByName(‘glszm’) extractor.enableFeatureClassByName(‘gldm’) extractor.enableFeatureClassByName(‘ngtdm’) # 邻域灰度差矩阵 # 3. 设置特定参数例如GLCM的距离和方向 extractor.settings[‘distanceThreshold’] 5 # 用于GLCM等特征的最大距离 extractor.settings[‘force2D’] True # 如果是2D分析 extractor.settings[‘force2Ddimension’] 0 # 沿哪个维度提取2D切片 # 4. 执行特征提取 # 使用之前加载的SimpleITK图像和掩膜 result extractor.execute(image_sitk, mask_sitk) # 5. 结果是一个字典键为特征名值为特征值 print(f“共提取 {len(result)} 个特征”) # 过滤出我们关心的纹理特征并打印 texture_features {k: v for k, v in result.items() if any( ft in k.lower() for ft in [‘glcm’, ‘glrlm’, ‘glszm’, ‘gldm’, ‘ngtdm’] )} for feat_name, feat_value in texture_features.items(): print(f“{feat_name}: {feat_value:.6f}”)使用pyradiomics的优势在于其经过了严格的验证支持多种图像类型和滤波并且计算效率高。对于初学者理解其底层原理后直接使用这样的成熟工具是快速上手的捷径。5. 实战案例乳腺癌影像纹理特征分析与解读现在让我们将这些知识应用于一个模拟的实战场景。假设我们有一个小型数据集包含20例乳腺癌患者的MRI影像和对应的肿瘤分割其中10例为恶性10例为良性。我们的目标是提取每个肿瘤的多种纹理特征并初步探索哪些特征可能对区分良恶性有帮助。首先我们批量处理所有患者数据构建特征表格。import os import pandas as pd from tqdm import tqdm # 用于显示进度条 def extract_features_for_patient(image_path, mask_path, extractor): “”“为单个患者提取特征。”“” try: image sitk.ReadImage(image_path) mask sitk.ReadImage(mask_path) # 确保掩膜是二值化的 mask sitk.BinaryThreshold(mask, lowerThreshold1, upperThreshold255, insideValue1, outsideValue0) features extractor.execute(image, mask) # 移除一些非纹理的通用特征 features.pop(‘diagnostics_Configuration_EnabledImageTypes’, None) features.pop(‘diagnostics_Image-original_Hash’, None) # ... 可以移除其他diagnostics开头的特征 return features except Exception as e: print(f“处理 {image_path} 时出错: {e}”) return None # 假设数据组织在如下目录结构 data_dir “./breast_mri_data/” patients [‘P001’, ‘P002’, …, ‘P020’] # 患者ID列表 labels [‘malignant’, ‘benign’] # 对应的标签需要从临床数据中获取 extractor featureextractor.RadiomicsFeatureExtractor() # 配置提取器启用所有纹理特征类 extractor.enableAllFeatures() all_features [] for pid, label in zip(patients, labels): img_path os.path.join(data_dir, pid, ‘image.nii.gz’) msk_path os.path.join(data_dir, pid, ‘mask.nii.gz’) if os.path.exists(img_path) and os.path.exists(msk_path): feats extract_features_for_patient(img_path, msk_path, extractor) if feats: feats[‘patient_id’] pid feats[‘label’] label all_features.append(feats) # 转换为DataFrame df_features pd.DataFrame(all_features) print(f“特征表形状: {df_features.shape}”) # (患者数, 特征数2) # 保存到CSV df_features.to_csv(‘breast_cancer_texture_features.csv’, indexFalse)得到特征表格后我们可以进行简单的探索性数据分析。例如计算每个特征在良恶性两组间的差异并用箱线图可视化最有区分度的几个特征。import seaborn as sns import matplotlib.pyplot as plt from scipy.stats import mannwhitneyu # 非参数检验 # 分离特征和标签 X df_features.drop([‘patient_id’, ‘label’], axis1) y df_features[‘label’] # 初始化一个列表来存储有显著差异的特征 significant_features [] # 对每个特征进行Mann-Whitney U检验假设数据不服从正态分布 for col in X.columns: malignant_vals X[y ‘malignant’][col].dropna() benign_vals X[y ‘benign’][col].dropna() if len(malignant_vals) 1 and len(benign_vals) 1: stat, p_value mannwhitneyu(malignant_vals, benign_vals, alternative‘two-sided’) if p_value 0.05: # 使用0.05作为显著性阈值 significant_features.append((col, p_value)) # 按p值排序 significant_features.sort(keylambda x: x[1]) print(f“发现 {len(significant_features)} 个显著特征 (p0.05):”) for feat, p in significant_features[:10]: # 显示前10个最显著的 print(f“ {feat}: p {p:.4e}”) # 可视化前3个最显著的特征 top_n 3 fig, axes plt.subplots(1, top_n, figsize(15, 5)) for idx, (feat, p) in enumerate(significant_features[:top_n]): ax axes[idx] sns.boxplot(x‘label’, yfeat, datadf_features, axax, palette‘Set2’) ax.set_title(f‘{feat}\n(p{p:.2e})’) ax.set_xlabel(‘’) ax.set_ylabel(‘Feature Value’) plt.tight_layout() plt.show()通过这样的分析我们可能会发现例如GLCM的对比度original_glcm_Contrast在恶性肿瘤中显著更高这可能反映了恶性肿瘤内部细胞排列更紊乱、异质性更强导致影像纹理的局部对比更明显。而GLSZM的小区域强调original_glszm_SmallAreaEmphasis在良性肿瘤中更高或许暗示良性病变内部结构更均质由许多小均匀区域构成。当然这只是一个初步的探索。在实际的影像组学研究中后续步骤还包括特征标准化、使用LASSO等算法进行特征选择以消除冗余、构建机器学习模型如逻辑回归、支持向量机或随机森林进行预测并采用交叉验证严格评估模型性能。但无论如何扎实地理解并能够计算这些纹理特征是整个流程中不可或缺的第一步。当你亲手运行这些代码看到一个个数字从影像中浮现出来时你对影像组学的理解就不再停留在纸面上了。