医学图像处理 期末速成

知识点清单

一、基础概念与系统（30分钟）

（一）数字图像基础

• 定义：离散化的图像，由像素矩阵构成，可被计算机存储与处理。
• 图像格式：
  • bmp：无损存储，文件容量大，保留原始图像完整信息。
  • jpg：有损压缩，日常通用，压缩后容量小但会损失部分细节。
  • png：支持透明通道，无损压缩，适合需保留透明背景的图像。
  • tif：专业图像格式，支持分层存储，常用于印刷、遥感等领域。
  • Dicom：医学专用，除图像外，还绑定患者信息（如姓名、检查号等 ），保障医疗数据合规性。
• 图像类型：
  • 二值图像：像素值仅为 0（黑 ）或 255（白 ），用于简单黑白对比场景（如文字识别 ）。
  • 灰度图像：单通道，像素值范围 0（黑 ）- 255（白 ），呈现从黑到白的渐变层次。
  • 彩色图像：多通道（如 RGB 三通道 ），通过不同通道色彩混合显示丰富色彩。

（二）医学影像与 PACS 系统

• 医学影像设备：
  • CT：通过 X 射线断层扫描，擅长显示骨骼、肺部等结构，图像有较高密度分辨率。
  • MRI：利用磁场和射频脉冲成像，对软组织（如脑部、肌肉 ）成像清晰，无辐射危害。
  • PETCT：结合 PET（正电子发射断层扫描 ）与 CT，既能显示代谢信息（找肿瘤等 ），又能呈现结构细节。
  • B 超：基于超声原理，实时成像，常用于腹部、妇科等部位检查。
• PACS 系统：
  • 定义：医学图像存储与传输系统，实现影像数字化管理。
  • 用途：统一存储、调阅医学影像，便于多科室协同诊断、远程会诊。
  • 构成：涵盖影像采集设备、存储服务器、诊断工作站等。
  • 图像格式：采用 Dicom 格式，确保图像与患者信息关联，符合医疗数据标准。

二、灰度变换与空间滤波（60分钟）

（一）灰度变换

• 线性变换：公式 g = a * f + b（g 为变换后灰度，f 为原灰度 ）。
  • a（对比度系数 ）：a > 1 增强对比度，图像明暗差异更显著；1 > a > 0 减弱对比度，图像趋于平淡；a < 0是图像求补，暗部变亮，亮部变暗。
  • b（亮度系数 ）：b > 0 整体提亮图像；b < 0 整体压暗图像。
• 分段线性变换：选取关键灰度区间（如医学图像中病灶对应的灰度范围 ），对区间内灰度拉伸、区间外灰度压缩，突出感兴趣区域。
• 灰级窗（窗口变换 ）：设定“感兴趣灰度范围 [T1, T2]”，将窗内灰度拉伸至 0 - 255 显示，窗外灰度压缩，便于观察特定组织（如 CT 肺窗聚焦肺部组织、骨窗聚焦骨骼 ）。

（二）直方图增强——直方图均衡化

• 核心目的：使图像灰度分布更均匀，提升对比度，解决医学图像中病灶与正常组织灰度黏连问题。
• 步骤：
  1. 统计原始图像各灰度值（0 - 255 ）出现的次数。
  2. 计算各灰度值的累计概率（累计出现次数 / 总像素数 ）。
  3. 将累计概率映射到 0 - 255 新灰度值（新灰度 = 累计概率 * 255 取整 ），实现灰度重新分布。

（三）图像平滑（去噪）

• 均值滤波：
  • 操作：以 3*3 窗口为例，计算窗口内 9 个像素的平均值（新像素 = （窗口内像素值之和）/ 9 ），替换中心像素。
  • 适用场景：有效去除高斯噪声（如医学图像中因设备产生的均匀分布噪声 ），但会使图像边缘、细节模糊。
• 中值滤波：
  • 操作：以 3*3 窗口为例，将窗口内 9 个像素值排序，取中间值替换中心像素。
  • 适用场景：对椒盐噪声（图像中随机出现的黑白斑点 ）去除效果极佳，且能较好保留图像边缘。

（四）图像锐化（突出边缘）

• 一阶微分算子（找边缘 ）：
  • Roberts 算子：
    • 模板（水平、垂直方向 ）：$\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$（水平 ）、$\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]$（垂直 ）。
    • 原理：计算邻域像素灰度差，快速定位边缘，对陡峭边缘敏感，但受噪声影响大。
  • Prewitt 算子：
    • 模板（水平、垂直 ）：水平 $\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]$ 、垂直 $\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]$ 。
    • 原理：通过邻域像素差分找边缘，相比 Roberts 算子，对噪声有一定平滑作用。
  • Sobel 算子：
    • 模板（水平、垂直 ）：水平 $\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$ 、垂直 $\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$ 。
    • 原理：加权计算邻域像素差分，突出边缘同时抑制噪声，是医学图像（如找器官轮廓 ）常用算子。
• 二阶微分算子（Laplacian ）：
  • 模板（4 邻域、8 邻域 ）：4 邻域 $\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]$ 、8 邻域 $\left[\begin{array}{ccc}1& 1& 1\\ 1& -8& 1\\ 1& 1& 1\end{array}\right]$ 。
  • 性质：各向同性（全方位找边缘 ），边缘定位精准，但会放大噪声；常与原始图像结合（g = f - ∇²f ，g 为处理后图像，f 为原图像 ），增强边缘与背景对比。

三、图像分割与特征提取（40分钟）

（一）图像分割算法

• 阈值法：
  • 原理：选定灰度阈值 T，将像素灰度值与 T 比较，大于 T 归为目标（或背景 ），小于 T 归为另一类。
  • 常用方法：双峰法，寻找图像灰度直方图两个峰值间的谷值作为 T，适合医学图像中组织与背景灰度差异明显场景（如 CT 骨头分割 ）。
• 区域生长法：
  • 步骤：
    1. 选择种子点（如医学图像中病灶内的像素 ）。
    2. 定义相似性准则（通常为灰度差，如相邻像素灰度差小于某阈值 ）。
    3. 从种子点开始，逐步合并满足准则的相邻像素，直到无符合条件像素，适合分割肿瘤等有连续边界的区域。
• 分水岭法：
  • 原理：将图像视为“地形”，灰度低的区域为山谷，灰度高的为山峰；模拟注水过程，不同“流域”（区域 ）被分割开。
  • 问题：易受噪声干扰导致过度分割（噪声被识别为山峰 ），一般先对图像滤波再使用。

（二）形态学操作

• 膨胀：用结构元（如 3*3 矩形 ）扩张亮区域，填补小空洞，使医学图像中病灶边界更完整。
• 腐蚀：用结构元收缩亮区域，去除小凸起，可分离粘连目标（如细胞粘连 ）。
• 开运算：先腐蚀后膨胀，能去除噪声、断开小连接，常用于医学图像背景噪声清理。
• 闭运算：先膨胀后腐蚀，可填补小空洞、连接相近目标，让目标区域更平滑。

（三）特征提取

• 连通性：
  • 4 连通：仅上下左右相邻的像素视为连通。
  • 8 连通：上下左右及对角线相邻的像素均视为连通，用于识别图像中连续的区域（如数医学图像中细胞团 ）。
• 目标特征计算：
  • 面积：统计连通分量包含的像素个数。
  • 周长：统计连通分量边缘像素数量。
  • 目标数量：通过识别不同连通分量，统计医学图像中目标（如 CT 里的结节 ）个数。

四、综合应用题答题套路（30分钟）

（一）算法选择

根据需求匹配算法：

• 去噪需求：椒盐噪声选“中值滤波”；高斯噪声选“均值滤波”。
• 增强需求：提升整体对比度选“直方图均衡化”；突出边缘选“图像锐化算子（如 Sobel ）”。
• 分割需求：组织与背景灰度差异大选“阈值法”；分割连续区域（如肿瘤 ）选“区域生长法”。

（二）步骤撰写

分点清晰罗列：

1. 输入图像：说明读取的图像格式、通道（如读取 Dicom 医学图像为灰度图 ）。
2. 算法 1 执行：写清算法名称、参数（如“中值滤波，3*3 窗口” ）、作用（如“去除椒盐噪声” ）。
3. 算法 2 执行（若有 ）：同上述格式，说明后续算法操作。
4. 输出结果：说明处理后图像的存储、应用（如“保存增强后图像用于病灶观察” ）。

（三）代码框架（Python + OpenCV ）

核心函数调用模板：

import cv2
# 1. 读取图像（以灰度图为例 ）
img = cv2.imread('medical_image.dcm', 0)  
# 2. 算法执行（如中值滤波去噪 ）
denoised_img = cv2.medianBlur(img, 3)  # 3*3 窗口中值滤波
# 3. 算法执行（如直方图均衡化增强 ）
enhanced_img = cv2.equalizeHist(denoised_img)  
# 4. 保存结果
cv2.imwrite('processed_image.png', enhanced_img)  

五、应试急救技巧（20分钟）

（一）选择、判断题

• 抓关键词判断：
  • 正确表述：“Dicom 是医学专用图像格式，包含患者信息”“中值滤波对椒盐噪声去除效果好” 。
  • 错误表述：“均值滤波最适合处理椒盐噪声”（应为中值滤波 ）“灰度图像是三通道图像”（灰度图像为单通道 ） 。

（二）简答题

• 背核心步骤与特点：
  • 直方图均衡化步骤：“统计灰度频次→计算累计概率→映射新灰度值” 。
  • Sobel 算子优点：“加权差分突出边缘，同时抑制噪声，适合医学图像轮廓提取” 。

（三）综合题

• 没时间写完整代码时：
  1. 写算法流程：如“先使用 cv2.medianBlur 函数进行 3*3 中值滤波去噪，再通过 cv2.equalizeHist 函数做直方图均衡化增强” 。
  2. 标注关键函数：即便不写完整代码，列出核心函数名，也可争取步骤分。

实验

一、基础

求平均降噪

import cv2
import numpy as np
from skimage.util import random_noise

# 读取原始图像
img = cv2.imread("headCT_ori.tif", -1)
# 获取图像的高度和宽度
height, width = img.shape[:2]

# 初始化一个数组用于存储8幅带噪声的图像，这里将每个通道（维度）存储一幅噪声图像
noise = np.zeros((height, width, 8), dtype=np.float64)

# 生成8幅带有高斯噪声的图像
for i in range(8):
    # 为每一幅图像添加高斯噪声，mode='gaussian'表示高斯噪声，mean为均值，var为方差
    noisy_img = random_noise(img, mode='gaussian', seed=None, clip=True, mean=0, var=0.1)
    # 将添加噪声后的图像转换为合适的数据类型并存储到noise数组中
    noise[:, :, i] = (noisy_img * 255).astype(np.float64)  # 还原到0-255范围

# 对8幅带噪声的图像对应像素求平均
img_aver = np.mean(noise, axis=2)  # axis=2表示沿着第三个维度（8幅图像这个维度）求平均
# 将结果转换为uint8类型以便显示
img_aver = img_aver.astype(np.uint8)

# 显示处理后的图像
cv2.imshow("result", img_aver)
cv2.waitKey(0)
cv2.destroyAllWindows()

数字剪影

import cv2
import numpy as np

# 1. 读取图像（蒙片与造影片）
img1 = cv2.imread("angiography_mask_image.tif")  
img2 = cv2.imread("angiography_live_image.tif")  

# 2. 图像减法（关键步骤：背景抵消）
# 转换为int32避免uint8溢出（减法可能出现负数）
img = np.int32(img2) - np.int32(img1)  

# 3. 归一化（将结果映射到0-255，保证视觉可读性）
# NORM_MINMAX：按最小-最大范围归一化，使结果覆盖0-255
cv2.normalize(img, img, 0, 255, norm_type=cv2.NORM_MINMAX)  

# 4. 类型转换（还原为uint8，OpenCV显示要求）
img = img.astype(np.uint8)  

# 5. 结果显示
cv2.imshow("result", img)  
cv2.waitKey(0)  
cv2.destroyAllWindows()  # 补全窗口销毁，避免程序残留

格式转换

import imageio
import os
import pydicom  # 需导入pydicom库

path1 = r"./BrainMRIDICOM_1/"
path2 = r"./BrainMRIDICOM_2/"
i = 1

for filename in os.listdir(path1):
    readpath = os.path.join(path1, filename)
    # 读取DICOM图像
    img = pydicom.dcmread(readpath)  
    writepath = os.path.join(path2, f"{i}.jpg")  
    # 将DICOM像素数组写入为JPG图像
    imageio.imwrite(writepath, img.pixel_array)  
    i += 1

二、图像增强 1

X射线图像做负片和gamma变换

import cv2
import numpy as np

# 读取乳腺X光图像（医学场景：乳腺癌筛查，钙化点检测）
img = cv2.imread('breast_xray.tif')

# 1. 负片转换（图像反转）
# - 数学原理：O = 255 - I
# - 医学价值：凸显被高密度区域掩盖的细节（如导管阴影）
img_neg = 255 - img
cv2.imshow("Negative", img_neg)

# 2. 伽马变换（γ=2）
# - 数学公式：O = 255 × (I/255)^γ
# - 医学价值：增强高灰度区域对比度（钙化点/肿块边缘）
img_gamma = np.power(img/255.0, 2) * 255
cv2.imshow("Gamma (γ=2)", img_gamma.astype(np.uint8))

# 窗口控制
cv2.waitKey(0)

直方图均衡化

import matplotlib.pyplot as plt
import cv2

# 读取胸部X光图像（医学场景：肺炎、肺结节等疾病筛查）
img = cv2.imread('chest.tif', -1)

# 直方图均衡化（核心算法）
# - 数学原理：通过累积分布函数(CDF)重新分配图像灰度
# - 医学价值：增强肺纹理(暗区)与骨骼/病灶(亮区)的对比度
imghis = cv2.equalizeHist(img)

# 可视化对比
cv2.imshow('Original', img)          # 原始图像
cv2.imshow('Equalized', imghis)      # 均衡化后图像

# 直方图分析（诊断价值）
plt.figure("Original Histogram")     # 原始图像灰度分布
plt.hist(img.flatten(), bins=256)
plt.figure("Equalized Histogram")    # 均衡化后灰度分布
plt.hist(imghis.flatten(), bins=256)
plt.show()

# 窗口控制
cv2.waitKey(0)
cv2.destroyAllWindows()

灰级窗处理

import pydicom
import cv2
import numpy as np

# 读取CT影像（DICOM格式）
ds = pydicom.dcmread("82821227.dcm")
im = ds.pixel_array

# 初始化结果数组
result = np.zeros_like(im, dtype=np.int16)

# 骨窗技术核心算法：
# - 骨组织CT值范围：1400-3000 HU
# - 低于下限/高于上限设为黑色
# - 范围内线性映射到0-255增强对比度
for i in range(im.shape[0]):
    for j in range(im.shape[1]):
        if im[i, j] <= 1400:
            result[i, j] = 0  # 软组织区域
        elif 1400 < im[i, j] <= 3000:
            result[i, j] = 255 * (im[i, j] - 1400) / 1600  # 骨骼区域
        else:
            result[i, j] = 0  # 金属植入物等超高密度区域

# 归一化与显示
cv2.normalize(result, result, 0, 255, cv2.NORM_MINMAX)
cv2.imshow("Bone Window", result.astype(np.uint8))
cv2.waitKey(0)

三、图像增强 2

平滑滤波

import cv2

# 读取头部CT图像（医学场景：脑出血、肿瘤等疾病诊断）
img = cv2.imread('headCT.bmp', -1)

# 滤波参数：3x3 窗口
n = 3  

# 1. 均值滤波
# 原理：邻域像素平均 → 模糊噪声
# 医学价值：平滑高斯噪声，保留整体轮廓（如脑组织边界）
img_avg = cv2.blur(img, (n, n))  

# 2. 中值滤波
# 原理：邻域像素排序取中值 → 抑制椒盐噪声
# 医学价值：保护边缘（如骨窗、病灶边界），适合含脉冲噪声的CT图像
img_med = cv2.medianBlur(img, n)  

# 可视化对比
cv2.imshow('Original', img)        # 原始图像
cv2.imshow('Average Filter', img_avg)  # 均值滤波（模糊噪声）
cv2.imshow('Median Filter', img_med)   # 中值滤波（保护边缘）

cv2.waitKey(0)

边缘锐化滤波

import cv2

# 读取脊柱MRI图像（医学场景：椎间盘突出、脊柱肿瘤等诊断）
img = cv2.imread('MRIspine.tif', -1)

# 1. Sobel 算子边缘检测
# 原理：计算 x/y 方向梯度 → 提取边缘信息
# 医学价值：突出脊柱轮廓、椎间盘边界、肿瘤边缘
x = cv2.Sobel(img, cv2.CV_16S, 1, 0)  # x 方向梯度（水平边缘）
y = cv2.Sobel(img, cv2.CV_16S, 0, 1)  # y 方向梯度（垂直边缘）

# 2. 梯度绝对值转换
# 作用：将有符号梯度转换为无符号（0-255），便于显示
ax = cv2.convertScaleAbs(x)
ay = cv2.convertScaleAbs(y)

# 3. 梯度融合（加权求和）
# 公式：dst = 0.5*ax + 0.5*ay → 平衡 x/y 方向边缘
dst = cv2.addWeighted(ax, 0.5, ay, 0.5, 0)

# 可视化对比
cv2.imshow('X Direction Edges', ax)   # x 方向边缘（水平结构：椎间盘间隙）
cv2.imshow('Y Direction Edges', ay)   # y 方向边缘（垂直结构：脊柱椎体）
cv2.imshow('Combined Edges', dst)     # 融合后边缘（完整脊柱轮廓）

cv2.waitKey(0)

拉普拉斯滤波

import cv2
import skimage as sk

# 读取肺部影像（医学场景：肺炎、肺结节、肺气肿诊断）
img = cv2.imread('lung.bmp', -1)
im = sk.img_as_float(img)  # 转换为浮点型，避免数值溢出

# 1. 拉普拉斯滤波（边缘增强）
# 原理：二阶导数计算边缘 → 突出高频细节
# 医学价值：增强肺纹理、结节边缘、支气管轮廓
laplacian = cv2.Laplacian(im, cv2.CV_64F, ksize=3)

# 2. 影像增强（原始影像 - 拉普拉斯边缘）
# 公式：result = 原始影像 - 拉普拉斯边缘 → 突出细节同时保留整体
result = cv2.subtract(im, laplacian)

# 可视化对比
cv2.imshow('Original Lung', im)        # 原始肺部影像
cv2.imshow('Enhanced Result', result)  # 增强后影像（肺纹理更清晰）

cv2.waitKey(0)

四、图像分割 1

阈值分割

import cv2
import numpy as np

# 读取聚合物囊泡图像（场景：微观结构分析，如材料科学/生物医学）
img = cv2.imread('polymersomes.tif', -1)

# 初始化迭代阈值（核心算法：迭代法阈值分割）
Tnext = np.mean(img)  # 第一步：用全局均值作为初始阈值
T = 0

# 迭代优化阈值的循环逻辑
# 原理：当阈值变化小于0.5时停止，保证收敛性
while (abs(T - Tnext) > 0.5):  
    T = Tnext  # 保存当前阈值
    
    # 分割前景（<=T）和背景（>T）并计算统计量
    part_1 = 0; n = 0  # 前景总和、前景像素数
    part_2 = 0; m = 0  # 背景总和、背景像素数
    
    # 遍历像素分类统计（双重循环实现）
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if img[i,j] <= T:
                part_1 += img[i,j]
                n += 1
            else:
                part_2 += img[i,j]
                m += 1
    
    # 计算新阈值（前景与背景均值的平均）
    ave_1 = part_1 / n
    ave_2 = part_2 / m
    Tnext = 0.5 * (ave_1 + ave_2)

# 两种阈值分割方法对比
# 1. 自定义迭代阈值分割
thres, result = cv2.threshold(img, T, 255, cv2.THRESH_BINARY)
# 2. Otsu自动阈值分割（经典对比算法）
thresotsu, resultotsu = cv2.threshold(img, 0, 255, cv2.THRESH_OTSU)

# 结果可视化
cv2.imshow('source', img)             # 原始图像
cv2.imshow('iteration result', result)  # 迭代法结果
cv2.imshow('otsu result', resultotsu)   # Otsu法结果

cv2.waitKey(0)

边缘检测

import cv2
import numpy as np

# 1. 定义 8 个方向的 Kirsch 卷积核（覆盖 360° 全方向边缘检测）
k1 = np.array([[5,5,5],[-3,0,-3],[-3,-3,-3]])  # 0° 水平边缘
k2 = np.array([[-3,5,5],[-3,0,5],[-3,-3,-3]])  # 45° 对角线边缘
k3 = np.array([[-3,-3,5],[-3,0,5],[-3,-3,5]])  # 90° 垂直边缘
k4 = np.array([[-3,-3,-3],[-3,0,5],[-3,5,5]])  # 135° 对角线边缘
k5 = np.array([[-3,-3,-3],[-3,0,-3],[5,5,5]])  # 180° 水平边缘
k6 = np.array([[-3,-3,-3],[5,0,-3],[5,5,-3]])  # 225° 对角线边缘
k7 = np.array([[5,-3,-3],[5,0,-3],[5,-3,-3]])  # 270° 垂直边缘
k8 = np.array([[5,5,-3],[5,0,-3],[-3,-3,-3]])  # 315° 对角线边缘

# 2. 读取医学图像（头部 CT 用于疾病诊断）
img = cv2.imread('headCT.tif', -1)  # -1 保持原始格式（灰度/彩色）

# 3. 图像补边（避免边缘像素因邻域不足导致信息丢失）
h, w = k1.shape  # 卷积核尺寸（3x3）
img_border = np.zeros((img.shape[0] + h - 1, img.shape[1] + w - 1))
img_border[(h-1)//2 : img.shape[0] + (h-1)//2, 
           (w-1)//2 : img.shape[1] + (w-1)//2] = img  # 原图居中放置

# 4. 初始化结果数组（存储边缘检测结果）
result = np.zeros(img.shape, np.int16)  # int16 避免卷积计算溢出

# 5. 逐像素执行 Kirsch 边缘检测（核心算法）
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        # 提取当前像素的 3x3 邻域（含补边）
        temp = img_border[i:i+h, j:j+w]
        
        # 计算 8 个方向的卷积响应（绝对值）
        k = [
            abs(np.sum(np.multiply(temp, k1))),  # 0° 方向边缘强度
            abs(np.sum(np.multiply(temp, k2))),  # 45° 方向边缘强度
            abs(np.sum(np.multiply(temp, k3))),  # 90° 方向边缘强度
            abs(np.sum(np.multiply(temp, k4))),  # 135° 方向边缘强度
            abs(np.sum(np.multiply(temp, k5))),  # 180° 方向边缘强度
            abs(np.sum(np.multiply(temp, k6))),  # 225° 方向边缘强度
            abs(np.sum(np.multiply(temp, k7))),  # 270° 方向边缘强度
            abs(np.sum(np.multiply(temp, k8))),  # 315° 方向边缘强度
        ]
        
        # 取最大响应作为最终边缘强度（保留最显著边缘方向）
        result[i, j] = max(k)

# 6. 结果归一化与类型转换（便于显示和后续处理）
cv2.normalize(result, result, 0, 255, cv2.NORM_MINMAX)  # 映射到 0-255 范围
result = result.astype(np.uint8)  # 转为 uint8 类型（OpenCV 显示要求）

# 7. 显示原始图像与边缘检测结果（医学诊断辅助）
cv2.imshow('Original Head CT', img)          # 原始头部 CT 图像
cv2.imshow('Kirsch Edge Detection', result)  # 边缘检测结果
cv2.waitKey(0)  # 等待按键退出

五、图像分割 2

区域生长

import cv2
import numpy as np
import pylab
import matplotlib.pyplot as plt

# 1. 读取医学图像（腹部 CT）
img = cv2.imread('AbdominalCT.bmp', -1)  # -1: 保持原始灰度格式加载

# 2. 初始化结果存储数组
result = np.zeros(img.shape, np.uint8)  # 存储最终分割结果
temp = np.zeros(img.shape, np.uint8)    # 存储当前生长区域

# 3. 显示初始空结果（调试用）
cv2.imshow('result', result)

# 4. 可视化原始图像（用于选择种子点）
pylab.imshow(img, cmap='Greys_r')  # 灰度反转显示（更符合视觉习惯）
pylab.title('Original Abdominal CT')
pylab.axis('off')
pylab.show(block=False)  # 非阻塞显示，方便后续交互

# 5. 手动选择种子点（核心交互步骤）
print("请在图像窗口点击选择种子点，按任意键继续...")
seed = pylab.ginput(1)  # 获取1个鼠标点击坐标 (x,y)
print(f"选择的种子点坐标: (x={seed[0][1]}, y={seed[0][0]})")  # 注意坐标转换

# 6. 种子点坐标转换（适配图像数组索引）
x = int(seed[0][1])  # pylab 坐标是 (列, 行)，转换为数组的 (行, 列)
y = int(seed[0][0])

# 7. 初始化种子点区域
result[x, y] = img[x, y]  # 标记种子点到结果
temp[x, y] = img[x, y]    # 标记种子点到临时区域

# 8. 区域生长参数初始化
flag = 1                  # 生长循环标志（1: 继续生长，0: 停止）
seed_ave = img[x, y]      # 初始种子点灰度均值
num = int(img[x, y])      # 灰度总和（初始为种子点灰度）
amount = 1                # 像素数量（初始为1个种子点）

# 9. 区域生长主循环
while flag:
    flag = 0  # 重置标志，若本次无生长则退出
    plt.pause(0.0000000001)  # 强制刷新图像显示
    
    # 临时数组用于存储本次生长的像素
    t = np.zeros(img.shape, np.uint8)
    
    # 遍历图像（跳过边缘，避免越界）
    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            # 如果当前像素是已生长区域
            if temp[i, j] != 0:
                # 遍历 3x3 邻域
                for k in range(i-1, i+2):
                    for l in range(j-1, j+2):
                        # 如果邻域像素未被标记
                        if result[k, l] == 0:
                            # 灰度差阈值判断（核心生长条件）
                            if abs(int(img[k, l]) - int(img[x, y])) < 33:
                                # 标记为生长区域
                                result[k, l] = img[k, l]
                                cv2.imshow('result', result)  # 实时显示生长过程
                                
                                # 更新统计量
                                num += int(img[k, l])
                                amount += 1
                                t[k, l] = img[k, l]  # 标记到临时数组
                                flag = 1  # 标记本次有生长，继续循环
    
    # 更新种子区域（本次生长的像素加入临时区域）
    temp = t
    # 重新计算种子区域灰度均值
    seed_ave = num / amount

# 10. 最终结果显示
cv2.imshow('result', result)  # 显示最终分割结果
cv2.waitKey(0)                # 等待按键退出
cv2.destroyAllWindows()       # 释放窗口资源

分水岭

import cv2
import numpy as np

# 1. 读取图像
img = cv2.imread('blob.tif')  # 读取彩色图像（BGR 格式）
cv2.imshow('img', img)        # 显示原始图像

# 2. 灰度转换
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转换为灰度图像

# 3. 二值化（反转）
# 阈值 110，最大值 255，反转二值化（背景为黑，前景为白）
ret, thresh = cv2.threshold(gray, 110, 255, cv2.THRESH_BINARY_INV)

# 4. 形态学操作（膨胀）
kernel = np.ones((3, 3), np.uint8)  # 3x3 结构元
# 膨胀操作：扩大背景区域，填充小空洞
background = cv2.dilate(thresh, kernel)

# 5. 距离变换
# 计算前景像素到背景的距离（DIST_L2: 欧氏距离）
dist_transform = cv2.distanceTransform(thresh, cv2.DIST_L2, 3)

# 6. 前景提取
# 阈值为距离变换最大值的 40%，提取前景区域
ret, foreground = cv2.threshold(dist_transform, 0.4 * dist_transform.max(), 255, 0)
foreground = np.uint8(foreground)  # 转换为 uint8 类型

# 7. 计算边界
# 背景减去前景，得到边界区域
boundary = cv2.subtract(background, foreground)

# 8. 标记连通组件
# 计算前景的连通组件（ret: 组件数量，marker: 标记结果）
ret, marker = cv2.connectedComponents(foreground)

# 9. 调整标记
marker = marker + 1          # 背景标记从 1 开始（0 保留给边界）
marker[boundary == 255] = 0  # 边界区域标记为 0

# 10. 应用分水岭算法
# 传入原始图像和标记，执行分水岭分割
marker = cv2.watershed(img, marker)

# 11. 标记边界
# 分水岭算法中，边界标记为 -1，这里用青色（0, 255, 255）标记
img[marker == -1] = [0, 255, 255]

# 12. 显示结果
cv2.imshow('result', img)  # 显示分割结果
cv2.waitKey(0)             # 等待按键退出
cv2.destroyAllWindows()    # 释放窗口资源

六、形态学

腐蚀、膨胀、开闭操作

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 1. 读取医学图像（细胞图像）
img = cv2.imread('cell.tif', -1)  # -1: 保持原始灰度格式加载

# 2. 定义形态学操作的核（3x3 结构元）
kernel = np.ones((3, 3), np.uint8)  # 常用 3x3 核，可调整为更大尺寸

# 3. 形态学操作
result_erode = cv2.erode(img, kernel)        # 腐蚀
result_dilate = cv2.dilate(img, kernel)      # 膨胀
result_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)    # 开运算
result_close = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel)  # 闭运算

# 4. 可视化结果（Matplotlib 子图布局）
plt.figure(figsize=(10, 6))  # 设置图像大小

# 4.1 原始图像
plt.subplot(2, 3, 1)  # 2行3列，第1个子图
plt.title('Original Cell Image')
plt.imshow(img, cmap='gray')  # 灰度显示
plt.axis('off')  # 隐藏坐标轴

# 4.2 腐蚀结果
plt.subplot(2, 3, 2)  # 第2个子图
plt.title('Erosion Result')
plt.imshow(result_erode, cmap='gray')
plt.axis('off')

# 4.3 膨胀结果
plt.subplot(2, 3, 3)  # 第3个子图
plt.title('Dilation Result')
plt.imshow(result_dilate, cmap='gray')
plt.axis('off')

# 4.4 开运算结果
plt.subplot(2, 3, 4)  # 第4个子图
plt.title('Opening Result')
plt.imshow(result_open, cmap='gray')
plt.axis('off')

# 4.5 闭运算结果
plt.subplot(2, 3, 5)  # 第5个子图
plt.title('Closing Result')
plt.imshow(result_close, cmap='gray')
plt.axis('off')

# 5. 显示图像
plt.tight_layout()  # 自动调整子图间距
plt.show()

提取骨架

import cv2
import numpy as np
from skimage import morphology

# 1. 读取二值图像（血管图像）
# -1: 保持原始灰度格式加载（图像应为二值化结果，如血管为255，背景为0）
img = cv2.imread('manual.tif', -1)  

# 2. 显示原始二值图像
cv2.imshow("binary vessel image", img)  

# 3. 二值化处理（确保图像为 0-1 矩阵，适配骨架提取算法）
# 将值为 255 的像素转为 1（skimage 骨架提取要求输入为 0-1 二值图像）
img[img == 255] = 1  

# 4. 提取骨架（核心操作）
# 使用 skimage 的 skeletonize 函数，提取图像的骨架
skeleton0 = morphology.skeletonize(img)  

# 5. 骨架结果转换（适配 OpenCV 显示）
# 将骨架结果转为 uint8 类型，并恢复为 0-255 范围
skeleton = skeleton0.astype(np.uint8) * 255  

# 6. 显示骨架图像
cv2.imshow("vessel skeleton image", skeleton)  

# 7. 等待按键退出
cv2.waitKey()  
cv2.destroyAllWindows()