《生物医学图像处理》实验7 实验报告¶
声明：本实验所有代码在gemini 3.1 pro下辅助完成，本人已充分理解实验内容
Project 1: Shape Feature & Fourier Descriptors¶
实验内容¶
实现肿瘤轮廓的傅里叶描述子（Fourier Descriptors）特征提取算法 。
逐步增加保留的描述子数量 $P$，重建肿瘤的几何轮廓并观察其视觉和定量变化 。
探讨和分析描述子数量 $P$ 与形状保真度（Shape Fidelity）之间的内在关系 。
实验原理¶
傅里叶描述子（Fourier Descriptors, FDs）是一种强有力的二维闭合轮廓特征表达方法。其基本原理是将空间域的二维几何轮廓映射到频域，利用离散傅里叶变换（DFT）的系数来描述物体的形状特征：
- 轮廓的复数化表达：
设在二维空间中提取出的肿瘤闭合轮廓由一系列按顺时针或逆时针排列的边界坐标点序列组成：$(x(n), y(n))$，其中 $n = 0, 1, 2, \dots, N-1$，$N$ 为轮廓点总数。 为了消除二维坐标的维度解耦限制，将这些坐标点投射到复平面上，表示为一个一维复数序列 $s(n)$：
$$s(n) = x(n) + j \cdot y(n)$$
- 离散傅里叶变换 (DFT)：
对复数序列 $s(n)$ 进行一维离散傅里叶变换，将其从空间域转换到频率域。变换后得到的复数系数 $a(k)$ 即为该形状的傅里叶描述子：
$$a(k) = \frac{1}{N} \sum_{n=0}^{N-1} s(n) e^{-j \frac{2\pi k n}{N}}, \quad k = 0, 1, 2, \dots, N-1$$
其中，$a(k)$ 的低频分量（靠近 $k=0$ 和 $k=N-1$ 的部分）代表了轮廓的宏观几何形状和整体轮廓（如位置、大小、大体长宽比等）；高频分量则代表了轮廓的微观精细结构、边缘凹凸细节以及像素级噪声。
- 逆傅里叶变换 (IDFT) 与轮廓重建：
在实际应用或进行形状重建时，我们可以选择只保留前 $P$ 个低频描述子（通常通过对变换后的频域系数进行中心平移 fftshift，截取中心区域的 $P$ 个主要能量系数，其余高频系数置零），然后再通过逆离散傅里叶变换（IDFT）还原出重建后的复数序列 $\hat{s}(n)$，从而解耦得到重建后的轮廓坐标 $(\hat{x}(n), \hat{y}(n))$：
$$\hat{s}(n) = \sum_{k=0}^{P-1} a(k) e^{j \frac{2\pi k n}{N}}$$
实验过程（代码编写）¶
import h5py
import numpy as np
import cv2
import matplotlib.pyplot as plt
import os
# ================= 1. 数据读取 (使用 h5py) =================
file_path = 'data/209.mat'
if not os.path.exists(file_path):
print(f"找不到文件：{file_path}，请检查 data 文件夹和文件是否存在！")
exit()
with h5py.File(file_path, 'r') as mat:
# 提取肿瘤掩膜 (Mask)
mask_data = np.array(mat['cjdata']['tumorMask'])
tumor_mask = mask_data.T.astype(np.uint8)
# 提取原始 MRI 脑部图像
image_data = np.array(mat['cjdata']['image'])
original_image = image_data.T
# ================= 2. 提取真实轮廓 =================
contours, _ = cv2.findContours(tumor_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
original_contour = max(contours, key=cv2.contourArea)
num_points = len(original_contour)
# ================= 3. 核心算法：傅里叶描述子 =================
complex_contour = original_contour[:, 0, 0] + 1j * original_contour[:, 0, 1]
fd = np.fft.fft(complex_contour)
# ================= 4. 实验观察：8个不同的 P 值 =================
# 【修改代码】：设置 8 个不同的 P 值进行对比
# 建议使用 2 的幂次方递增，这样能最好地展示低频到高频的过渡
P_values = [2, 4, 8, 16, 32, 64, 128, num_points // 2]
# 【修改代码】：调整画布大小，准备 2 行 4 列的布局
plt.figure(figsize=(20, 10))
plt.suptitle(f'Project 1: Fourier Descriptors Reconstruction (1.mat)\nTotal Points: {num_points}', fontsize=20, y=0.98)
for i, P in enumerate(P_values):
# --- 截断高频，保留 P 个低频 ---
fd_shifted = np.fft.fftshift(fd)
truncated_fd = np.zeros_like(fd_shifted)
center = len(fd_shifted) // 2
half_P = P // 2
start_idx = max(0, center - half_P)
end_idx = min(len(fd_shifted), center + half_P)
truncated_fd[start_idx:end_idx] = fd_shifted[start_idx:end_idx]
# --- 逆傅里叶变换，重建 ---
reconstructed_complex = np.fft.ifft(np.fft.ifftshift(truncated_fd))
recon_x = np.real(reconstructed_complex)
recon_y = np.imag(reconstructed_complex)
# --- 画图展示 (2行4列布局) ---
plt.subplot(2, 4, i + 1)
plt.imshow(original_image, cmap='gray')
plt.plot(recon_x, recon_y, 'r-', linewidth=2)
# 动态显示标题，标明当前的 P 值
plt.title(f'P = {P}', fontsize=16)
plt.axis('off')
# 调整子图之间的间距，防止重叠
plt.tight_layout()
plt.subplots_adjust(top=0.9) # 给总标题留出空间
plt.show()
实验结果与分析¶
在具体实验过程中，我们通过改变描述子数量 $P$（分别取 $P = 2, 4, 8, 16, 32, 64, 128$ 等），对真实脑肿瘤 MRI 轮廓进行重建处理，观察到的规律以及描述子数量 $P$ 与形状保真度（Shape Fidelity）的关系如下 ：
- 极低描述子区间 ($P = 2, 4, 8$)：宏观轮廓初步拟合
- 当 $P=2$ 时，重建的轮廓在几何形式上仅表现为一个标准的圆形，它仅携带了肿瘤病灶的大小和位置质心信息，无法呈现任何边缘边界特征。
- 随着 $P$ 增加到 4 或 8，轮廓开始拉伸，逐渐向肿瘤的主要伸展方向靠拢，形成了一个粗糙的块状多边形。这一阶段保真度较低，属于全局低频粗糙拟合。
- 中等描述子区间 ($P = 16, 32, 64$)：核心特征
- 伴随 $P$ 值的进一步增大，傅里叶描述子的能量开始向中高频渗透。重建轮廓在此时发生显著的质变，开始快速捕捉肿瘤轮廓上明显的凹陷、突起和病灶特有的不规则拐角。
- 在该区间内，形状保真度曲线通常呈现出一种极快的上升趋势。此时重建出的边界已经能较好地用于医学临床上对肿瘤整体体积和主侵袭方向的评估。
- 高描述子区间 ($P \ge 128$)：细节微调与过拟合/噪声收敛
- 当 $P$ 达到 128 或更高时，重建出的边界线条在肉眼上与经过算法提取的绝对原始边界几乎无法区分，形状保真度无限趋近于 100%。
- 继续增加 $P$ 值后，保真度曲线进入极度平缓的饱和平台期。这意味着剩余的高频分量在数学上捕捉的仅是像素级别的锯齿轮廓或分割时产生的微弱离散噪声。在实际医学特征工程中，适当丢弃这些超高频项不仅可以显著降低特征矩阵的维度，还能起到边缘平滑与抗噪的滤波作用。
结论： 描述子数量 $P$ 与形状保真度之间存在非线性的收益递减关系。较小的 $P$ 只能代表粗糙的宏观外形；中等大小的 $P$ 是医学图像特征提取在“数据压缩”与“几何精度”之间的平衡点；而过高的 $P$ 则会导致特征冗余并引入像素级别的细微噪声。
Project 2: First-Order Stats¶
实验内容¶
本实验使用 Lab 7 的真实脑肿瘤 MRI 数据 。首先需要定位肿瘤的质心（Center of mass），以此为中心提取感兴趣区域（ROI） 。
接着，计算该区域的一阶直方图统计特征，包括均值（Mean）、方差（Variance）、一致性（Uniformity）和熵（Entropy） 。
在此基础上，实现一个迭代的 $3 \times 3$ 均值滤波器，用以模拟渐进的高斯模糊过程，并记录在迭代过程中这些统计指标的动态变化。
最后，将分析统计指标随迭代次数的演变趋势，并将此模糊与特征提取过程在整张原始图像上重复进行，以对比局部 ROI 与全局图像在特征演变上的差异。
实验原理¶
1. 质心计算 (Center of Mass)
图像中肿瘤区域的质心可以通过图像矩（Image Moments）来确定。对于二值化后的肿瘤掩膜 $M(x,y)$，其零阶矩 $M_{00}$ 代表肿瘤的总面积（或像素总数），一阶矩 $M_{10}$ 和 $M_{01}$ 则分别代表 $x$ 和 $y$ 方向上的灰度加权坐标总和。质心坐标 $(\bar{x}, \bar{y})$ 的计算公式为：
$$\bar{x} = \frac{M_{10}}{M_{00}}, \quad \bar{y} = \frac{M_{01}}{M_{00}}$$
以此质心为基准点向外辐射，即可精准截取包含肿瘤主体的 ROI 矩形区域。
2. 一阶统计特征 (First-Order Stats)
一阶统计特征基于图像区域内像素灰度值的概率分布（即直方图）进行计算，而不考虑像素之间的空间相对位置关系。设图像总像素数为 $N$，$p(z_i)$ 为灰度级 $z_i$ 在区域内出现的概率，$L$ 为总灰度级数（本实验归一化为 256 级）：
- 均值 (Mean)：反映区域的整体平均亮度。
$$\mu = \sum_{i=0}^{L-1} z_i p(z_i)$$
- 方差 (Variance)：反映图像的灰度对比度与异质性。方差越大，说明灰度分布越分散，细节差异越明显。
$$\sigma^2 = \sum_{i=0}^{L-1} (z_i - \mu)^2 p(z_i)$$
- 一致性 (Uniformity)：衡量直方图分布的平坦度或集中度。图像灰度越集中于少数几个级别，该值越接近 1。
$$U = \sum_{i=0}^{L-1} p(z_i)^2$$
- 熵 (Entropy)：衡量图像纹理的复杂度和灰度分布的随机性。细节越丰富、分布越混乱的区域，熵值越高。
$$H = -\sum_{i=0}^{L-1} p(z_i) \log_2 p(z_i)$$
3. 迭代滤波模拟高斯模糊
基于中心极限定理，对图像进行多次连续的 $3 \times 3$ 均值滤波操作，其等效的空间卷积核会逐渐平滑，最终在数学上逼近高斯滤波核。这种迭代方式能够动态展示图像细节被逐步抹平的过程。
实验过程（代码编写）¶
import h5py
import numpy as np
import cv2
import matplotlib.pyplot as plt
import os
import warnings
# ================= 0. 环境设置与预处理 =================
# 1. 彻底屏蔽 UserWarning 及其他非致命警告
warnings.filterwarnings("ignore")
# 2. 设置 Matplotlib 的中文字体，确保图表中文和负号显示正常
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'PingFang SC']
plt.rcParams['axes.unicode_minus'] = False
# ================= 1. 数据读取 =================
file_path = 'data/1.mat'
if not os.path.exists(file_path):
print(f"找不到文件：{file_path}，请检查 data 文件夹和文件是否存在！")
exit()
with h5py.File(file_path, 'r') as mat:
mask_data = np.array(mat['cjdata']['tumorMask'])
tumor_mask = mask_data.T.astype(np.uint8)
image_data = np.array(mat['cjdata']['image'])
img_min, img_max = image_data.min(), image_data.max()
# 归一化到 0-255 范围，方便计算直方图
original_image = (((image_data - img_min) / (img_max - img_min)) * 255).T.astype(np.uint8)
# ================= 2. 计算肿瘤质心并确定 ROI 边界 =================
M = cv2.moments(tumor_mask)
cx = int(M['m10'] / M['m00'])
cy = int(M['m01'] / M['m00'])
print(f"--- 算法初始化 ---")
print(f"计算获得肿瘤质心坐标: (X={cx}, Y={cy})\n")
# 以质心为中心，切取 120x120 的矩形区域
half_size = 60
ymin, ymax = max(0, cy - half_size), min(original_image.shape[0], cy + half_size)
xmin, xmax = max(0, cx - half_size), min(original_image.shape[1], cx + half_size)
roi_image = original_image[ymin:ymax, xmin:xmax].copy()
# ================= 3. 统计指标计算函数 =================
def compute_metrics(img_block):
pixels = img_block.flatten()
total_pixels = len(pixels)
mean_val = np.mean(pixels)
var_val = np.var(pixels)
# 计算 256 级的灰度直方图概率分布
hist, _ = np.histogram(pixels, bins=256, range=(0, 256))
prob = hist / total_pixels
# 一致性 (Uniformity)
uniformity_val = np.sum(prob ** 2)
# 熵 (Entropy)
active_prob = prob[prob > 0]
entropy_val = -np.sum(active_prob * np.log2(active_prob))
return mean_val, var_val, uniformity_val, entropy_val
# ================= 4. 核心迭代、均值滤波与终端日志打印 =================
max_iters = 30
metrics_roi = {"mean": [], "var": [], "uniformity": [], "entropy": []}
metrics_full = {"mean": [], "var": [], "uniformity": [], "entropy": []}
current_roi = roi_image.copy()
current_full = original_image.copy()
# 打印工整的表头，包含全部 8 个核心指标
print(f"{'Iter':<4} | {'ROI Mean':<8} {'ROI Var':<10} {'ROI Unif':<8} {'ROI Ent':<8} | {'Full Mean':<9} {'Full Var':<11} {'Full Unif':<9} {'Full Ent':<9}")
print("-" * 105)
for i in range(max_iters + 1):
# 计算当前迭代下的 ROI 和 全图特征
m_r, v_r, u_r, e_r = compute_metrics(current_roi)
m_f, v_f, u_f, e_f = compute_metrics(current_full)
# 记录数据用于后续画图
metrics_roi["mean"].append(m_r)
metrics_roi["var"].append(v_r)
metrics_roi["uniformity"].append(u_r)
metrics_roi["entropy"].append(e_r)
metrics_full["mean"].append(m_f)
metrics_full["var"].append(v_f)
metrics_full["uniformity"].append(u_f)
metrics_full["entropy"].append(e_f)
# 实时向终端输出每一次迭代的全部具体数值（第0次即为原始图像特征）
print(f"{i:<4} | {m_r:<8.2f} {v_r:<10.2f} {u_r:<8.4f} {e_r:<8.4f} | {m_f:<9.2f} {v_f:<11.2f} {u_f:<9.4f} {e_f:<9.4f}")
# 执行 3x3 均值滤波（最后一次迭代计算完指标后无需再滤波）
if i < max_iters:
current_roi = cv2.blur(current_roi, (3, 3))
current_full = cv2.blur(current_full, (3, 3))
print("-" * 105)
print("计算完成，正在弹出可视化成果窗口...\n")
# ================= 5. 图像可视化：质心与模糊结果 (Figure 1) =================
plt.figure(figsize=(15, 5))
plt.suptitle('Project 2: 肿瘤质心定位与高斯模糊模拟结果对照', fontsize=16)
# 子图 1：原图 + 质心标记 + ROI 绿框
plt.subplot(1, 3, 1)
plt.imshow(original_image, cmap='gray')
# 【已修复】：将 speakersize 改为正确的 markersize
plt.plot(cx, cy, 'r+', markersize=15, markeredgewidth=2, label='肿瘤质心')
# 绘制绿色矩形框标出 ROI 区域
rect = plt.Rectangle((xmin, ymin), xmax - xmin, ymax - ymin, edgecolor='lime', facecolor='none', linewidth=2, label='ROI 边界')
plt.gca().add_patch(rect)
plt.title('1. 原始 MRI (含质心与ROI框)')
plt.legend(loc='upper right')
plt.axis('off')
# 子图 2：原始肿瘤 ROI（未模糊）
plt.subplot(1, 3, 2)
plt.imshow(roi_image, cmap='gray')
plt.title('2. 原始肿瘤 ROI (清晰)')
plt.axis('off')
# 子图 3：最终迭代 30 次后的模糊 ROI（模拟高斯模糊结果，对应 PPT Slide 7）
plt.subplot(1, 3, 3)
plt.imshow(current_roi, cmap='gray')
plt.title(f'3. 均值滤波迭代 {max_iters} 次后 (渐进高斯模糊)')
plt.axis('off')
plt.tight_layout()
# ================= 6. 指标动态变化折线图 (Figure 2) =================
iters = np.arange(max_iters + 1)
plt.figure(figsize=(14, 9))
plt.suptitle('Project 2: 一阶统计指标动态演变趋势 (ROI vs 全图)', fontsize=16)
keys = ["mean", "var", "uniformity", "entropy"]
titles = ["均值 (Mean) - 整体亮度", "方差 (Variance) - 图像对比度/异质性", "一致性 (Uniformity) - 直方图集中度", "熵 (Entropy) - 纹理复杂度/随机性"]
for idx, key in enumerate(keys):
plt.subplot(2, 2, idx + 1)
plt.plot(iters, metrics_roi[key], 'g-o', linewidth=2, label='肿瘤 ROI')
plt.plot(iters, metrics_full[key], 'r--s', linewidth=2, label='整张图像 (选做部分)')
plt.title(titles[idx], fontsize=12)
plt.xlabel('滤波迭代次数 (Iterations)', fontsize=10)
plt.grid(True, linestyle=':', alpha=0.6)
plt.legend(fontsize=10)
plt.tight_layout()
plt.subplots_adjust(top=0.9)
# 同时显示两个窗口
plt.show()
--- 算法初始化 --- 计算获得肿瘤质心坐标: (X=327, Y=216) Iter | ROI Mean ROI Var ROI Unif ROI Ent | Full Mean Full Var Full Unif Full Ent --------------------------------------------------------------------------------------------------------- 0 | 76.24 1407.54 0.0081 7.1227 | 33.09 2176.54 0.1085 5.1656 1 | 76.25 1350.03 0.0083 7.0888 | 33.09 2128.65 0.1097 5.1371 2 | 76.25 1311.75 0.0084 7.0676 | 33.08 2096.80 0.1099 5.1320 3 | 76.25 1282.01 0.0085 7.0518 | 33.08 2071.93 0.1100 5.1293 4 | 76.25 1257.73 0.0086 7.0393 | 33.08 2051.09 0.1099 5.1288 5 | 76.25 1236.88 0.0087 7.0288 | 33.08 2032.88 0.1098 5.1292 6 | 76.26 1218.51 0.0087 7.0184 | 33.07 2016.54 0.1097 5.1295 7 | 76.26 1202.16 0.0087 7.0111 | 33.07 2001.80 0.1096 5.1302 8 | 76.26 1187.43 0.0088 7.0046 | 33.07 1988.11 0.1095 5.1311 9 | 76.26 1173.90 0.0088 6.9975 | 33.07 1975.33 0.1094 5.1320 10 | 76.26 1161.43 0.0088 6.9926 | 33.07 1963.33 0.1093 5.1328 11 | 76.26 1149.86 0.0089 6.9856 | 33.07 1952.09 0.1092 5.1336 12 | 76.26 1139.16 0.0089 6.9806 | 33.07 1941.52 0.1091 5.1344 13 | 76.26 1129.32 0.0089 6.9759 | 33.07 1931.45 0.1090 5.1351 14 | 76.26 1119.86 0.0089 6.9711 | 33.07 1921.82 0.1088 5.1358 15 | 76.25 1111.16 0.0090 6.9658 | 33.07 1912.66 0.1087 5.1365 16 | 76.25 1102.99 0.0090 6.9611 | 33.07 1903.84 0.1086 5.1370 17 | 76.25 1095.05 0.0090 6.9548 | 33.07 1895.34 0.1085 5.1377 18 | 76.25 1087.62 0.0091 6.9493 | 33.07 1887.23 0.1084 5.1382 19 | 76.25 1080.75 0.0091 6.9453 | 33.07 1879.41 0.1083 5.1387 20 | 76.24 1074.18 0.0091 6.9400 | 33.07 1871.85 0.1083 5.1392 21 | 76.24 1067.81 0.0092 6.9336 | 33.07 1864.55 0.1082 5.1397 22 | 76.24 1061.63 0.0092 6.9285 | 33.07 1857.49 0.1081 5.1400 23 | 76.23 1055.69 0.0092 6.9221 | 33.07 1850.66 0.1080 5.1405 24 | 76.23 1050.29 0.0093 6.9178 | 33.07 1844.00 0.1079 5.1407 25 | 76.23 1045.15 0.0093 6.9117 | 33.07 1837.52 0.1079 5.1410 26 | 76.22 1040.36 0.0094 6.9075 | 33.07 1831.26 0.1078 5.1414 27 | 76.22 1035.80 0.0094 6.9035 | 33.07 1825.17 0.1077 5.1415 28 | 76.23 1031.45 0.0094 6.8996 | 33.07 1819.22 0.1076 5.1418 29 | 76.23 1026.96 0.0094 6.8960 | 33.06 1813.39 0.1076 5.1419 30 | 76.22 1022.82 0.0094 6.8927 | 33.06 1807.72 0.1075 5.1420 --------------------------------------------------------------------------------------------------------- 计算完成，正在弹出可视化成果窗口...
实验结果与分析¶
根据算法初始化结果，提取到的肿瘤质心坐标为 (X=327, Y=216)。随着 $3 \times 3$ 均值滤波的 30 次迭代，各项指标的动态演变趋势以及局部与全局的对比分析如下：
1. 肿瘤 ROI 区域的指标动态演变趋势
- 均值 (Mean)：在 30 次迭代中，ROI 的均值始终稳定在 76.24 左右。均值滤波本质上是像素值的局部平均与重新分配，它不会改变图像的总能量，因此整体亮度保持恒定。
- 方差 (Variance)：呈现出显著且持续的下降趋势（从 1407.54 降至 1022.82）。这是因为滤波操作不断抹平肿瘤内部的纹理细节和噪声，相邻像素之间的灰度差异被强制拉近，导致对比度和异质性大幅削弱。
- 一致性 (Uniformity)：呈现稳步上升趋势（从 0.0081 升至 0.0094）。随着高斯模糊的加深，像素的灰度值逐渐向均值靠拢，直方图从宽幅分布聚拢为较窄的峰，使得各灰度级概率的平方和显著增加。
- 熵 (Entropy)：呈现持续下降趋势（从 7.1227 降至 6.8927）。图像细节的丢失意味着纹理复杂度的降低和信息量的减少，灰度分布不再具有高随机性，从而导致熵值变小。
2. 全局图像与 ROI 的对比分析 (Entire Image vs. ROI)
基础特征量级差异：
- 全图的均值（33.09）远低于 ROI（76.24），且全图的方差（2176.54）远高于 ROI（1407.54）。这是因为整张 MRI 图像包含了大面积灰度值为 0 的黑色背景，以及极亮的高颅骨边缘。这种极端的明暗反差拉低了整体均值，同时极大地拉高了全局方差。
- 全图的一致性（0.1085）比 ROI 高出一个数量级，而熵（5.1656）低于 ROI。庞大的纯黑背景构成了直方图上的绝对主导峰，这种高度的同质化导致全图的一致性极高、随机性（熵）偏低；而 ROI 内部全是充满细节的肿瘤组织，因而表现出更低的集中度和更高的熵。
演变趋势的分化：
- 在迭代滤波过程中，全图的一致性并没有像 ROI 那样上升，反而出现了微弱的下降（0.1085 降至 0.1075）；同时，全图的熵在最初几次迭代中下降后，又出现了缓慢的回升波动。
- 原因分析：在局部 ROI 中，模糊主要发生在质地相近的组织内部，促使灰度集中；而在整张图像中，存在着大量“背景-头皮-脑组织”之间的极锐利边界。在迭代模糊时，高亮的头皮和脑组织会发生严重的灰度扩散，渗透进原本统一的纯黑背景中。这种跨越宏观边界的渗透，反而创造出了大量原本不存在的中间过渡灰度级，打破了原本绝对集中的背景峰，从而导致全局的一致性下降和熵值的微弱波动。
Project 3 : Texture Feature & GLCM Analysis¶
实验内容¶
本实验聚焦于医学图像的微观空间纹理分析。首先，提取 Lab 7 提供的真实脑肿瘤感兴趣区域（ROI），并对其进行灰度量化（Quantization）以降低计算复杂度并凸显主要纹理特征。
随后，基于特定的空间关系（不同的采样步长/距离），计算该区域的灰度共生矩阵（GLCM）。
接着，调用相关函数从 GLCM 中提取二阶统计特征，重点观测对比度（Contrast）、同质性（Homogeneity）以及能量（Energy）。
最后，结合定量数据与 GLCM 可视化图像，深入探讨采样步长与图像中可能存在的周期性结构纹理之间的交互关系 。
实验原理¶
1. 灰度共生矩阵 (GLCM)
一阶统计特征仅关注单个像素灰度的出现概率，而缺乏空间位置信息。GLCM 是一种二阶统计纹理分析方法，它通过统计图像中成对像素的灰度分布规律来描述纹理。 具体而言，对于设定好的空间方向 $\theta$ 和步长距离 $d$，GLCM 矩阵中的每一个元素 $P(i,j)$ 表示：在整幅图像中，灰度值为 $i$ 的像素与相距为 $d$、方向为 $\theta$ 的灰度值为 $j$ 的像素结对出现的相对频率。 由于原始 $256 \times 256$ 的矩阵计算量极大且分布稀疏，在计算特征前通常会将图像量化（Quantization）至更低的灰度级（如 16 级），从而使统计分布更加集中，纹理特征更具代表性。
2. 二阶统计特征 (Second-Order Stats)
基于计算出的 GLCM，可以提取多种二阶特征来定量描述纹理性质：
- 对比度 (Contrast)：衡量矩阵的值如何分布在主对角线周围，反映图像局部灰度变化的剧烈程度（即纹理沟槽的深浅）。公式中赋予了远离对角线的元素 $(i-j)^2$ 的高权重：$\sum_{i,j} (i-j)^2 P(i,j)$。
- 同质性 (Homogeneity / 逆差分矩)：衡量图像局部的均匀度。当 GLCM 中的非零元素高度集中在主对角线附近时（即相邻像素灰度差异小），同质性极高。
- 能量 (Energy / 角二阶矩)：衡量纹理的稳定性和规则度。矩阵元素分布越集中（存在占主导地位的像素对），能量值越大：$\sum_{i,j} P(i,j)^2$。
实验过程（代码编写）¶
import h5py
import numpy as np
import cv2
import matplotlib.pyplot as plt
from skimage.feature import graycomatrix, graycoprops
import os
import warnings
# ================= 0. 环境与画图配置 =================
warnings.filterwarnings("ignore")
# 配置中文字体，Jupyter 中同样适用
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'PingFang SC']
plt.rcParams['axes.unicode_minus'] = False
# ================= 1. 数据读取与 ROI 提取 =================
file_path = 'data/1.mat'
if not os.path.exists(file_path):
print(f"找不到文件：{file_path}，请检查当前工作目录下的 data 文件夹！")
else:
with h5py.File(file_path, 'r') as mat:
mask_data = np.array(mat['cjdata']['tumorMask'])
tumor_mask = mask_data.T.astype(np.uint8)
image_data = np.array(mat['cjdata']['image'])
img_min, img_max = image_data.min(), image_data.max()
original_image = (((image_data - img_min) / (img_max - img_min)) * 255).T.astype(np.uint8)
# 利用质心截取 120x120 的肿瘤局部 ROI
M = cv2.moments(tumor_mask)
cx, cy = int(M['m10'] / M['m00']), int(M['m01'] / M['m00'])
half_size = 60
ymin, ymax = max(0, cy - half_size), min(original_image.shape[0], cy + half_size)
xmin, xmax = max(0, cx - half_size), min(original_image.shape[1], cx + half_size)
roi_image = original_image[ymin:ymax, xmin:xmax].copy()
# ==========================================================
# 第一部分：二阶纹理特征分析 (16级量化)
# ==========================================================
# 1. 图像灰度量化
levels = 16
roi_quantized = (roi_image / (256 / levels)).astype(np.uint8)
# 2. GLCM 计算与提取 (16级)
distances_feat = [1, 2, 3, 5, 8, 12, 15, 20, 25, 30]
glcm_feat = graycomatrix(roi_quantized, distances=distances_feat, angles=[0], levels=levels, symmetric=True, normed=True)
contrast = graycoprops(glcm_feat, 'contrast').flatten()
homogeneity = graycoprops(glcm_feat, 'homogeneity').flatten()
energy = graycoprops(glcm_feat, 'energy').flatten()
# 3. 打印数据表
print(f"{'采样步长 (Distance)':<18} | {'对比度 (Contrast)':<15} | {'同质性 (Homogeneity)':<18} | {'能量 (Energy)'}")
print("-" * 80)
for i, d in enumerate(distances_feat):
print(f"{d:<22} | {contrast[i]:<19.2f} | {homogeneity[i]:<22.4f} | {energy[i]:.4f}")
# 4. 绘制特征分析图 (Figure 1)
plt.figure(figsize=(16, 5))
plt.suptitle('Project 3 Part A: GLCM 纹理特征与采样步长的关系', fontsize=16)
plt.subplot(1, 3, 1)
plt.imshow(roi_image, cmap='gray')
plt.title('1. 原始肿瘤 ROI (256 级)')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.imshow(roi_quantized, cmap='gray')
plt.title(f'2. 量化后的 ROI ({levels} 级)')
plt.axis('off')
ax1 = plt.subplot(1, 3, 3)
ax1.plot(distances_feat, contrast, 'r-o', linewidth=2, label='对比度 (Contrast)')
ax1.set_xlabel('采样步长 Step Size (像素)', fontsize=12)
ax1.set_ylabel('对比度', color='r', fontsize=12)
ax1.tick_params(axis='y', labelcolor='r')
ax1.grid(True, linestyle=':', alpha=0.6)
ax2 = ax1.twinx()
ax2.plot(distances_feat, homogeneity, 'g--s', linewidth=2, label='同质性 (Homogeneity)')
ax2.set_ylabel('同质性', color='g', fontsize=12)
ax2.tick_params(axis='y', labelcolor='g')
plt.title('3. 采样步长对二阶统计指标的影响')
lines_1, labels_1 = ax1.get_legend_handles_labels()
lines_2, labels_2 = ax2.get_legend_handles_labels()
ax1.legend(lines_1 + lines_2, labels_1 + labels_2, loc='center right')
plt.tight_layout()
plt.show() # 在 Jupyter 中渲染出第一张图
# ==========================================================
# 第二部分：GLCM 矩阵可视化 (256级原图)
# ==========================================================
# 1. 计算 256 级原图的高分辨率 GLCM
d1, d2 = 1, 15
# 注意这里用的是没有量化的 roi_image，levels=256，不归一化以便肉眼看清频数
glcm_v1 = graycomatrix(roi_image, distances=[d1], angles=[0], levels=256, symmetric=True, normed=False)
glcm_v2 = graycomatrix(roi_image, distances=[d2], angles=[0], levels=256, symmetric=True, normed=False)
# 2. 对数变换以突出暗部细节 (核心步骤)
visual_1 = np.log1p(glcm_v1[:, :, 0, 0])
visual_2 = np.log1p(glcm_v2[:, :, 0, 0])
# 3. 绘制矩阵可视化图 (Figure 2)
plt.figure(figsize=(14, 6))
plt.suptitle('Project 3 Part B: GLCM 灰度共生矩阵可视化对比', fontsize=16)
plt.subplot(1, 2, 1)
plt.imshow(visual_1, cmap='gray')
plt.title(f'GLCM (采样步长 d = {d1})', fontsize=14)
plt.xlabel('像素 j 的灰度级 (0-255)')
plt.ylabel('像素 i 的灰度级 (0-255)')
plt.colorbar(label='Log(频数)')
plt.subplot(1, 2, 2)
plt.imshow(visual_2, cmap='gray')
plt.title(f'GLCM (采样步长 d = {d2})', fontsize=14)
plt.xlabel('像素 j 的灰度级 (0-255)')
plt.ylabel('像素 i 的灰度级 (0-255)')
plt.colorbar(label='Log(频数)')
plt.tight_layout()
plt.show()
采样步长 (Distance) | 对比度 (Contrast) | 同质性 (Homogeneity) | 能量 (Energy) -------------------------------------------------------------------------------- 1 | 0.30 | 0.8695 | 0.2700 2 | 0.80 | 0.7692 | 0.2280 3 | 1.40 | 0.7018 | 0.2056 5 | 2.49 | 0.6212 | 0.1822 8 | 3.59 | 0.5579 | 0.1661 12 | 4.30 | 0.5204 | 0.1569 15 | 4.65 | 0.5019 | 0.1519 20 | 5.16 | 0.4760 | 0.1462 25 | 5.68 | 0.4477 | 0.1405 30 | 6.50 | 0.4179 | 0.1359
实验结果与分析¶
根据实验记录，随着采样步长（Distance）从 1 逐渐增加到 30，肿瘤 ROI 的二阶纹理特征表现出显著的单调变化：对比度从 0.30 升至 6.50；同质性从 0.8695 降至 0.4179；能量从 0.2700 降至 0.1359。
1. 采样步长与周期性结构纹理的交互关系探讨
在图像处理中，如果目标区域包含明显的周期性结构纹理（如超声波纹、肌肉纤维栅格，设其空间周期为 $T$），采样步长 $d$ 会与其产生强烈的空间共振交互：
- 同相采样：当 $d$ 刚好等于周期 $T$（或 $T$ 的整数倍）时，采样的两个像素会同时落在纹理的波峰或波谷，灰度高度一致。此时，对比度会骤降至局部极小值，同质性会飙升至局部极大值。
- 反相采样：当 $d$ 等于半个周期（如 $0.5T, 1.5T$）时，一个像素在波峰，另一个必然在波谷，灰度反差最大。此时，对比度会达到局部最高峰。
- 本实验肿瘤数据分析：在本实验的数据曲线中，对比度随着步长的增加呈现平滑的上升趋势，并未出现明显的波谷或突变振荡。这说明该脑肿瘤内部的组织纹理是高度无序和随机的，并不存在宏观的、规律性的周期结构纹理（或者其潜在的周期 $T$ 远大于我们测试的最大步长 30 像素）。
2. GLCM 矩阵可视化的动态演变分析
- 短步长 ($d=1$)：如可视化左图所示，此时矩阵呈现出一条极其明亮、纤细的主对角线，周围几乎是纯黑。这在物理上意味着：在极短的距离内，相邻细胞组织的灰度几乎完全相同（即 $i \approx j$）。因此，能量高度集中在对角线上，对应了数据表中极高的同质性（0.8695）和极低的对比度（0.30）。
- 长步长 ($d=15$)：如可视化右图所示，随着跨越的物理距离增加，采样点跨越了肿瘤内部不同的微观组织结构。原本集中的明亮对角线如同“散焦”一般，发散、弥漫到了矩阵的非对角线区域（即 $i$ 与 $j$ 的差异变大）。这种“能量向两侧发散”的视觉现象，在数学上直接导致了 $(i-j)^2$ 权重项的增加，完美解释了为何数据表中对比度骤增至 4.65，且同质性大幅稀释。
实验心得¶
在本次生物医学图像处理实验中，我完整经历了一套从底层数据读取到高级特征提取的医学图像分析流水线。这不仅加深了我对图像处理理论的理解，也极大提升了我的工程代码实现能力。通过这三个循序渐进的 Project，我获得了以下几点心得。
1. 频域视角的降维与重构
通过实现傅里叶描述子（Fourier Descriptors），我直观地认识到二维几何形状是可以转化为一维信号并投射到频率域进行分析的。实验中不断调整描述子数量 $P$ 的过程，让我深刻体会到了数据压缩与信息保真度之间的博弈关系。低频分量决定了肿瘤的宏观走向，高频分量勾勒了边缘的不规则突起。舍弃极高频的分量不仅能实现特征的降维，还能有效过滤掉分割过程中产生的边缘噪声，这对于后续构建机器学习分类器具有极大的现实意义。
2. 统计指标的物理意义与边界效应
在计算一阶统计特征并模拟高斯模糊的过程中，我看到了抽象数学公式在图像上的具象表达。方差（Variance）的下降和一致性（Uniformity）的上升，完美量化了图像细节被抹平、灰度趋于集中的过程。在选做部分（ROI vs 全图）的对比实验中：全图在迭代模糊时，发生了明显的灰度扩散，导致其统计演变曲线与单纯的肿瘤 ROI 产生了截然不同的趋势。这提醒我在未来的医学影像分析中，精准的 ROI 提取（如利用质心定位）是防止背景干扰、提取有效特征的绝对前提。
3. 空间纹理的二阶量化与共振
灰度共生矩阵（GLCM）的学习让我跳出了单一像素的局限，开始关注像素之间的空间协同关系。在实验中我认识到，灰度量化（Quantization）是进行 GLCM 分析前必不可少的步骤，否则稀疏的矩阵将掩盖真实的纹理特征。随着采样步长的增加，对比度逐渐上升而同质性不断下降，这从侧面证明了肿瘤内部纹理的无序性。通过对 GLCM 矩阵本身进行对数变换可视化，我更是直观地看到了“能量是如何随着距离增加而从主对角线向四周发散的”，这为理解空间结构与采样步长之间的“共振交互”提供了最直接的证据。
4. 算法落地与工程实践的挑战
除了理论上的收获，本次实验也是一次极佳的 Python 图像处理工程演练。从传统 MATLAB .mat (v7.3) 格式的解析，我学习到了使用 h5py 库并处理 Fortran-style 到 C-style 的矩阵转置问题。