import SimpleITK as sitk
from radiomics import featureextractor, setVerbosity, getFeatureClasses
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import os
import pandas as pd
import logging
import matplotlib.colors as mcolors
from matplotlib.colors import LinearSegmentedColormap
from skimage import measure
from mpl_toolkits.mplot3d.art3d import Poly3DCollection # 添加缺失的导入
def set_chinese_font():
from matplotlib import font_manager
font_path = "C:/Windows/WinSxS/amd64_microsoft-windows-font-truetype-simhei_31bf3856ad364e35_10.0.26100.1_none_f11b5ebedca17dd9/simhei.ttf" # 或 msyh.ttc、STSong.ttf 等
font_prop = font_manager.FontProperties(fname=font_path)
plt.rcParams['font.family'] = font_prop.get_name()
plt.rcParams['axes.unicode_minus'] = False
# 设置日志详细程度
setVerbosity(logging.INFO)
# 设置路径
base_path = r"C:\Users\53145\OneDrive\Desktop\1\radiomics\radiomic\p"
image_dir = os.path.join(base_path, "image")
mask_dir = os.path.join(base_path, "mask")
results_dir = os.path.join(base_path, "results")
os.makedirs(results_dir, exist_ok=True)
# 设定设置
settings = {
"binWidth": 25,
"resampledPixelSpacing": [1, 1, 1],
"interpolator": sitk.sitkBSpline,
"normalize": True,
"voxelBased": True # 启用基于体素的特征提取
}
# 初始化特征提取器
extractor = featureextractor.RadiomicsFeatureExtractor(**settings)
extractor.enableImageTypes(Original={}, LoG={}, Wavelet={})
# 打印启用的特征
print("启用的图像类型:", extractor.enabledImagetypes)
print("启用的特征类:", getFeatureClasses().keys())
# ************ 数据准备 ************
img_file = "3.nii.gz" # 这里替换成你要处理的图像文件名
base_name = img_file.replace(".nii.gz", "")
img_path = os.path.join(image_dir, img_file)
mask_file = f"{base_name}mask.nii.gz"
mask_path = os.path.join(mask_dir, mask_file)
try:
if not os.path.exists(mask_path):
raise FileNotFoundError(f"掩膜文件 {mask_file} 不存在")
# 读取图像和掩膜
image = sitk.ReadImage(img_path)
mask = sitk.ReadImage(mask_path)
# 打印图像信息
print(f"图像尺寸: {image.GetSize()}")
print(f"图像间距: {image.GetSpacing()}")
print(f"掩膜尺寸: {mask.GetSize()}")
print(f"掩膜间距: {mask.GetSpacing()}")
# 强制空间对齐
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(image)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
aligned_mask = resampler.Execute(mask)
# 检查掩膜有效性
mask_array = sitk.GetArrayFromImage(aligned_mask)
if np.sum(mask_array) < 10:
raise ValueError("掩膜区域太小,无法提取特征")
# ************ 特征提取 ************
print("正在提取特征图...")
# 提取特征
feature_vector = extractor.execute(image, aligned_mask, voxelBased=True)
# 检查特征提取结果
if not feature_vector:
raise ValueError("未提取到任何特征")
# 查找我们需要的特征图
target_feature_name = "wavelet-LHH_glszm_ZoneEntropy"
feature_map = None
# 尝试找到匹配的特征图
for feature_name, feature_value in feature_vector.items():
if target_feature_name in feature_name:
if isinstance(feature_value, sitk.Image):
feature_map = feature_value
target_feature_name = feature_name # 使用完整的特征名
break
if feature_map is None:
# 如果找不到特定特征,使用第一个特征图
for feature_name, feature_value in feature_vector.items():
if isinstance(feature_value, sitk.Image):
feature_map = feature_value
target_feature_name = feature_name
print(f"未找到 '{target_feature_name}',使用第一个特征图: {feature_name}")
break
if feature_map is None:
raise ValueError("未提取到任何特征图")
print(f"使用特征图: {target_feature_name}")
print(f"特征图尺寸: {feature_map.GetSize()}")
# 转换为NumPy数组
feature_array = sitk.GetArrayFromImage(feature_map)
# ************ 特征图可视化 ************
print("可视化特征图...")
# 创建自定义颜色映射
colors = ["darkblue", "blue", "cyan", "green", "yellow", "orange", "red"]
cmap = LinearSegmentedColormap.from_list("custom_cmap", colors)
# 确定要可视化的切片范围
z_start = max(0, feature_array.shape[0] // 3)
z_end = min(feature_array.shape[0], 2 * feature_array.shape[0] // 3)
# 创建多切片可视化
fig, axes = plt.subplots(3, 3, figsize=(15, 12))
fig.suptitle(f"特征图可视化: {target_feature_name}\n{base_name}", fontsize=16)
# 选择9个切片均匀分布在ROI范围内
slice_indices = np.linspace(z_start, z_end - 1, 9, dtype=int)
# 计算全局最小值和最大值用于颜色标准化
# 注意:背景为0,我们只考虑非零区域
non_zero_mask = (feature_array != 0)
if non_zero_mask.any():
vmin = np.percentile(feature_array[non_zero_mask], 5)
vmax = np.percentile(feature_array[non_zero_mask], 95)
else:
vmin = 0
vmax = 1
for i, slice_idx in enumerate(slice_indices):
ax = axes[i // 3, i % 3]
slice_data = feature_array[slice_idx]
# 显示热图
sns.heatmap(slice_data,
cmap=cmap,
cbar=i == 8, # 只在最后一个子图显示颜色条
ax=ax,
vmin=vmin,
vmax=vmax,
xticklabels=False,
yticklabels=False)
ax.set_title(f"切片 {slice_idx}")
set_chinese_font()
plt.tight_layout(rect=[0, 0, 1, 0.96])
output_path = os.path.join(results_dir, f"{base_name}_{target_feature_name.replace(':', '_')}_feature_map.png")
plt.savefig(output_path, dpi=300)
plt.close() # 关闭图形以释放内存
print(f"特征图已保存为: {output_path}")
except Exception as e:
print(f"处理失败:{img_file} | 错误:{str(e)}")
import traceback
traceback.print_exc()