import numpy as np
import cv2
import matplotlib.pyplot as plt
from scipy.signal import convolve2d
import os

def apply_sobel_kernel(image):
    """
    实现 Sobel 梯度核计算
    """
    Kx = np.array([[-1, 0, 1], 
                   [-2, 0, 2], 
                   [-1, 0, 1]], dtype=np.float32)
    
    Ky = np.array([[-1, -2, -1], 
                   [ 0,  0,  0], 
                   [ 1,  2,  1]], dtype=np.float32)
    
    Gx = convolve2d(image, Kx, mode='same', boundary='symm')
    Gy = convolve2d(image, Ky, mode='same', boundary='symm')
    
    magnitude = np.sqrt(Gx**2 + Gy**2)
    
    if magnitude.max() > 0:
        magnitude = (magnitude / magnitude.max()) * 255
        
    return magnitude.astype(np.uint8)

def segment_lesion_edge_only(img, threshold_T):
    # 1. 5x5 平滑滤波器
    img_smooth = cv2.GaussianBlur(img, (5, 5), 0)
    
    # 2. Sobel 滤波器
    grad_mag = apply_sobel_kernel(img_smooth)
    
    # 3. 阈值 T
    _, edge_mask = cv2.threshold(grad_mag, threshold_T, 255, cv2.THRESH_BINARY)
    
    # 4. 去除大脑外围边缘
    # 获取大脑主体区域，用较大的核进行腐蚀，彻底排除头骨和皮层边缘
    _, brain_mask = cv2.threshold(img_smooth, 20, 255, cv2.THRESH_BINARY)
    kernel_erode = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (25, 25))
    brain_core = cv2.erode(brain_mask, kernel_erode, iterations=1)
    
    internal_edges = cv2.bitwise_and(edge_mask, brain_core)
    
    # 5. 闭运算
    # 注意：这里缩小了闭运算的核 (5x5 或 7x7)，为了连接断点，同时防止边缘线糊成一团厚厚的面
    kernel_close = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
    closed_edges = cv2.morphologyEx(internal_edges, cv2.MORPH_CLOSE, kernel_close)
    
    # 6. 连通域处理 (提取病灶边缘环)
    # 对线段进行连通域分析
    num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(closed_edges, connectivity=8)
    
    lesion_edge_mask = np.zeros_like(img)
    valid_labels = []
    
    for i in range(1, num_labels): # 跳过背景(0)
        area = stats[i, cv2.CC_STAT_AREA] # 这里的 area 实际上是这条“线”包含的像素个数
        # 过滤掉细碎的短线(噪点)，保留较长的边缘线
        if 150 < area < 3000:
            valid_labels.append(i)
            
    if valid_labels:
        # 在符合长度的边缘线中，取像素最多(最长/最完整)的那条作为病灶轮廓
        best_label = max(valid_labels, key=lambda x: stats[x, cv2.CC_STAT_AREA])
        lesion_edge_mask[labels == best_label] = 255
        
    return grad_mag, closed_edges, lesion_edge_mask

def main():
    if not (os.path.exists('T1.npy') and os.path.exists('T2.npy')):
        print("请确保 T1.npy 和 T2.npy 在当前目录下。")
        return

    img_t1 = np.load('T1.npy')
    img_t2 = np.load('T2.npy')
    
    img_t1 = cv2.normalize(img_t1, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
    img_t2 = cv2.normalize(img_t2, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

    # 阈值可根据实际效果微调
    T1_THRESHOLD = 45
    T2_THRESHOLD = 55 

    t1_grad, t1_edges, t1_lesion = segment_lesion_edge_only(img_t1, T1_THRESHOLD)
    t2_grad, t2_edges, t2_lesion = segment_lesion_edge_only(img_t2, T2_THRESHOLD)

    # ================= 结果可视化 =================
    fig, axes = plt.subplots(2, 4, figsize=(16, 8))
    plt.subplots_adjust(wspace=0.3, hspace=0.3)
    
    titles = ['Original', 'Sobel Gradient', 'Closed Edges', 'Final Lesion Edge']
    
    images_t1 = [img_t1, t1_grad, t1_edges, t1_lesion]
    for i in range(4):
        axes[0, i].imshow(images_t1[i], cmap='gray')
        axes[0, i].set_title(f'T1: {titles[i]}')
        axes[0, i].axis('off')
        
    images_t2 = [img_t2, t2_grad, t2_edges, t2_lesion]
    for i in range(4):
        axes[1, i].imshow(images_t2[i], cmap='gray')
        axes[1, i].set_title(f'T2: {titles[i]}')
        axes[1, i].axis('off')

    plt.show()

if __name__ == "__main__":
    main()

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

def iterative_threshold_manual(image, initial_T, delta=0.5):
    """
    实现迭代全局阈值分割
    """
    # 排除背景像素(<=10)，仅对脑组织进行计算
    pixels = image[image > 10].astype(np.float32)
    if len(pixels) == 0: return 0
    
    T = float(initial_T)
    while True:
        # 分为两类
        g1 = pixels[pixels > T]
        g2 = pixels[pixels <= T]
        
        # 计算均值
        mu1 = np.mean(g1) if len(g1) > 0 else 0
        mu2 = np.mean(g2) if len(g2) > 0 else 0
        
        T_new = (mu1 + mu2) / 2.0
        
        # 判断收敛
        if abs(T - T_new) < delta:
            break
        T = T_new
    return T

def otsu_threshold_manual(image):
    """
    手动实现 Otsu 最优阈值算法，并返回方差曲线
    """
    # 排除背景像素，只统计脑组织部分
    brain_pixels = image[image > 10]
    total_pixels = brain_pixels.size
    
    # 计算直方图
    hist = np.histogram(brain_pixels, bins=256, range=(0, 256))[0]
    
    variances = np.zeros(256)
    max_variance = -1
    best_T = 0
    
    # 遍历所有可能的阈值 t
    for t in range(1, 255):
        # 计算 w0, w1 (类概率权重)
        w0 = np.sum(hist[:t]) / total_pixels
        w1 = np.sum(hist[t:]) / total_pixels
        
        if w0 == 0 or w1 == 0:
            continue
            
        # 计算 mu0, mu1 (类平均灰度)
        mu0 = np.sum(np.arange(t) * hist[:t]) / (w0 * total_pixels)
        mu1 = np.sum(np.arange(t, 256) * hist[t:]) / (w1 * total_pixels)
        
        # 计算类间方差 sigma_b^2
        variance = w0 * w1 * ((mu0 - mu1) ** 2)
        variances[t] = variance
        
        if variance > max_variance:
            max_variance = variance
            best_T = t
            
    return best_T, variances

def main():
    # 检查数据
    if not (os.path.exists('T1.npy') and os.path.exists('T2.npy')):
        print("请确保数据文件 T1.npy 和 T2.npy 在当前目录下。")
        return

    # 加载并归一化到 0-255
    img_t1 = cv2.normalize(np.load('T1.npy'), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
    img_t2 = cv2.normalize(np.load('T2.npy'), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

    # 设定三个不同的初始 T 值
    initial_Ts = [60, 125, 190]
    
    # 绘图设置：2行(T1, T2) x 6列(原图, 3个迭代结果, Otsu结果, 方差曲线)
    fig, axes = plt.subplots(2, 6, figsize=(22, 8))
    plt.subplots_adjust(wspace=0.3, hspace=0.3)

    for row, img, label in zip([0, 1], [img_t1, img_t2], ['T1', 'T2']):
        # 1. 原始图像
        axes[row, 0].imshow(img, cmap='gray')
        axes[row, 0].set_title(f"{label} Original")
        
        # 2. 迭代法：三个不同的初始 T 值
        for i, t_init in enumerate(initial_Ts):
            t_final = iterative_threshold_manual(img, t_init)
            _, bin_img = cv2.threshold(img, t_final, 255, cv2.THRESH_BINARY)
            axes[row, i+1].imshow(bin_img, cmap='gray')
            axes[row, i+1].set_title(f"Iterative T_init={t_init}\nFinal T={t_final:.1f}")
        
        # 3. Otsu 法结果
        otsu_t, variances = otsu_threshold_manual(img)
        _, otsu_bin = cv2.threshold(img, otsu_t, 255, cv2.THRESH_BINARY)
        axes[row, 4].imshow(otsu_bin, cmap='gray')
        axes[row, 4].set_title(f"Otsu's Result\nOptimal T={otsu_t}")
        
        # 4. Otsu 类间方差曲线
        axes[row, 5].plot(variances, color='blue', linewidth=1.5)
        axes[row, 5].axvline(otsu_t, color='red', linestyle='--', label=f'Best T={otsu_t}')
        axes[row, 5].set_title(f"{label} Otsu Variance Curve")
        axes[row, 5].set_xlabel("Threshold t")
        axes[row, 5].set_ylabel("Variance")
        axes[row, 5].legend()
        axes[row, 5].grid(True, alpha=0.3)

    # 隐藏图像部分的坐标轴
    for ax in axes[:, :5].flatten():
        ax.axis('off')

    plt.suptitle("Project 2: Global and Optimal Thresholding Segmentation", fontsize=16)
    plt.show()

if __name__ == "__main__":
    main()

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

# 全局变量
seed_point = None
img_t2 = None

def get_kmeans_mask(img, k, w_spatial, seed_x, seed_y):
    """手动实现改进 K-means：灰度 + 空间坐标"""
    h, w = img.shape
    y_idx, x_idx = np.indices((h, w))
    # 构建特征向量：[灰度, 归一化x*权重, 归一化y*权重]
    data = np.stack([
        img.flatten() / 255.0, 
        (x_idx.flatten() / w) * w_spatial, 
        (y_idx.flatten() / h) * w_spatial
    ], axis=1).astype(np.float32)
    
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    _, labels, _ = cv2.kmeans(data, k, None, criteria, 10, cv2.KMEANS_PP_CENTERS)
    labels = labels.reshape((h, w))
    
    # 根据种子点坐标提取对应的聚类簇
    target_label = labels[seed_y, seed_x]
    mask = np.zeros_like(img)
    mask[labels == target_label] = 255
    return mask, labels

def main():
    global seed_point, img_t2
    if not os.path.exists('T2.npy'): 
        print("错误：未找到 T2.npy 文件")
        return

    # 1. 预处理
    raw = np.load('T2.npy')
    img_t2 = cv2.normalize(raw, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
    h, w = img_t2.shape

    # 2. 交互选取种子点
    print(">>> 请在弹出窗口点击病灶中心...")
    cv2.imshow("1. Click Lesion Center", img_t2)
    def on_click(event, x, y, flags, param):
        global seed_point
        if event == cv2.EVENT_LBUTTONDOWN:
            seed_point = (x, y)
            cv2.destroyWindow("1. Click Lesion Center")
    cv2.setMouseCallback("1. Click Lesion Center", on_click)
    cv2.waitKey(0)
    
    if seed_point is None: 
        print("未选取种子点，程序退出。")
        return
    sx, sy = seed_point

    # 3. 区域生长：带滑动条实时调优
    print(">>> 正在运行区域生长，请调整滑动条后按【回车键】确认...")
    cv2.namedWindow("2. Adjust Region Growing")
    cv2.createTrackbar("Tol", "2. Adjust Region Growing", 15, 50, lambda x: None)
    
    mask_rg = np.zeros_like(img_t2)
    final_rg_vol = 0
    while True:
        tol = cv2.getTrackbarPos("Tol", "2. Adjust Region Growing")
        ff_mask = np.zeros((h + 2, w + 2), np.uint8)
        flags = 4 | (255 << 8) | cv2.FLOODFILL_MASK_ONLY
        cv2.floodFill(img_t2, ff_mask, (sx, sy), 0, (tol,), (tol,), flags)
        mask_rg = ff_mask[1:-1, 1:-1]
        final_rg_vol = np.sum(mask_rg == 255)
        
        preview = cv2.cvtColor(img_t2, cv2.COLOR_GRAY2BGR)
        preview[mask_rg == 255] = [0, 0, 255] 
        cv2.putText(preview, f"Vol: {final_rg_vol}mm2", (10, 30), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,0,255), 2)
        cv2.imshow("2. Adjust Region Growing", preview)
        
        key = cv2.waitKey(1) & 0xFF
        if key == 13: # 回车确认
            cv2.destroyAllWindows()
            break

    # 4. 改进 K-Means 多参数对比测试
    print(">>> 正在进行 K-Means 多参数对比测试...")
    k_values = [3, 5, 7]
    weights = [0.0, 0.1, 0.2] 
    kmeans_results = [] # 用于最后输出体积

    km_fig, km_axes = plt.subplots(len(k_values), len(weights), figsize=(15, 12))
    plt.subplots_adjust(hspace=0.4)
    
    for i, k in enumerate(k_values):
        for j, w_s in enumerate(weights):
            mask, _ = get_kmeans_mask(img_t2, k, w_s, sx, sy)
            vol = np.sum(mask == 255)
            
            # 记录结果
            method_name = f"Standard (K={k})" if w_s == 0 else f"Improved (K={k}, W={w_s})"
            kmeans_results.append((method_name, vol))
            
            km_axes[i, j].imshow(mask, cmap='gray')
            km_axes[i, j].set_title(f"{method_name}\nVol: {vol}mm²", fontsize=10)
            km_axes[i, j].axis('off')

    # 5. 分水岭 (Watershed)
    _, brain_mask = cv2.threshold(img_t2, 20, 255, cv2.THRESH_BINARY)
    markers = np.zeros(img_t2.shape, dtype=np.int32)
    markers[brain_mask == 0] = 1 
    cv2.circle(markers, (sx, sy), 5, 2, -1) 
    
    gradient = cv2.morphologyEx(img_t2, cv2.MORPH_GRADIENT, np.ones((3,3)))
    cv2.watershed(cv2.cvtColor(gradient, cv2.COLOR_GRAY2BGR), markers)
    mask_water = np.zeros_like(img_t2)
    mask_water[markers == 2] = 255
    water_vol = np.sum(mask_water == 255)

    # 6. 最终展示图
    res_fig, res_axes = plt.subplots(1, 3, figsize=(18, 6))
    res_axes[0].imshow(img_t2, cmap='gray')
    res_axes[0].scatter(sx, sy, c='red', s=50, label='Seed')
    res_axes[0].set_title("Original T2 + Seed")
    
    res_axes[1].imshow(mask_rg, cmap='gray')
    res_axes[1].set_title(f"Region Growing Mask\nVol: {final_rg_vol}mm²")
    
    res_axes[2].imshow(mask_water, cmap='gray')
    res_axes[2].set_title(f"Watershed Mask\nVol: {water_vol}mm²")
    
    for ax in res_axes: ax.axis('off')
    
    # 7. 控制台输出体积统计结果
    print("\n" + "="*40)
    print("      病灶体积统计报告 (mm²)")
    print("="*40)
    print(f"{'分割方法':<30} | {'体积 (mm²)':<10}")
    print("-" * 45)
    print(f"{'区域生长 (Region Growing)':<30} | {final_rg_vol:<10}")
    print(f"{'分水岭 (Watershed)':<30} | {water_vol:<10}")
    print("-" * 45)
    for name, v in kmeans_results:
        print(f"{name:<30} | {v:<10}")
    print("="*40)

    print("\n[完成] 图像已生成。请查看对比图表并关闭窗口结束程序。")
    plt.show()

if __name__ == "__main__":
    main()

>>> 请在弹出窗口点击病灶中心...
>>> 正在运行区域生长，请调整滑动条后按【回车键】确认...
>>> 正在进行 K-Means 多参数对比测试...

========================================
      病灶体积统计报告 (mm²)
========================================
分割方法                           | 体积 (mm²)  
---------------------------------------------
区域生长 (Region Growing)          | 2469      
分水岭 (Watershed)                | 3095      
---------------------------------------------
Standard (K=3)                 | 4594      
Improved (K=3, W=0.1)          | 5060      
Improved (K=3, W=0.2)          | 5174      
Standard (K=5)                 | 3452      
Improved (K=5, W=0.1)          | 4164      
Improved (K=5, W=0.2)          | 5647      
Standard (K=7)                 | 3250      
Improved (K=7, W=0.1)          | 4195      
Improved (K=7, W=0.2)          | 4749      
========================================

[完成] 图像已生成。请查看对比图表并关闭窗口结束程序。

算法	优点	缺点
迭代阈值法	1. 原理简单，计算开销极小。 2. 逻辑直观，通过物理意义上的“中点”寻找平衡。 3. 对初始值不敏感，收敛速度快。	1. 属于启发式方法，缺乏严格的统计学证明其为最优解。 2. 在多噪声或灰度分布极度不均时，可能陷入局部震荡。
大津法 (Otsu)	1. 全自动，无需设置任何初值。 2. 数学严谨，基于类间方差最大化，具有理论上的最优性。 3. 适应性广，对直方图双峰明显的图像效果极佳。	1. 计算量相对略大（需遍历 0-255 灰度级）。 2. 当目标与背景的面积比例极其悬殊时，阈值会向面积大的一方偏移，导致分割失败。

算法	优点	缺点	应用场景
区域生长	保证结果的强连通性，适合提取单一组织。	对种子点极度敏感，易发生泄漏。	边界清晰、背景简单的病灶提取。
标准 K-means	无需手动设点，计算速度快。	结果破碎，受噪声影响大，无法区分不同位置的亮区。	图像初步分类、多组织普查。
改进 K-means	结果圆润连续，能剔除远处噪点干扰。	$W$ 过大会引入非病灶区，产生几何畸变。	存在信号不均、需平滑边缘的病灶分割。
分水岭算法	边缘定位最精准，能处理粘连。	极易过度分割，参数调整极其复杂。	对体积精确度要求极高的解剖结构分割。

《生物医学图像处理》实验五实验报告¶

Project 1: Edge Detection¶

实验内容¶

实验原理¶

实验结果与分析¶

T1 与 T2 提取效果对比¶

边缘提取效果不佳的原因分析¶

算法改进思路¶

Project 2: Thresholding¶

实验内容¶

实验原理¶

迭代全局阈值法¶

大津法 (Otsu's Method)¶

实验过程（代码编写）¶

实验结果与分析¶

初始阈值 $T_{init}$ 对迭代法的影响¶

T1 与 T2 模态的分割效果差异¶

Otsu 法与迭代法的效果比较¶

算法优缺点说明¶

Project 3: Segmentation¶

实验内容¶

实验原理¶

区域生长 (Region Growing)¶

改进 K-means 聚类¶

分水岭算法 (Watershed)¶

实验结果与分析¶

区域生长算法分析¶

分水岭算法分析¶

改进 K-means 算法参数测试与分析¶

各算法比较与应用场景¶

实验心得¶

《生物医学图像处理》实验五 实验报告¶

Project 1: Edge Detection¶

实验内容¶

实验原理¶

实验结果与分析¶

T1 与 T2 提取效果对比¶

边缘提取效果不佳的原因分析¶

算法改进思路¶

Project 2: Thresholding¶

实验内容¶

实验原理¶

迭代全局阈值法¶

大津法 (Otsu's Method)¶

实验过程（代码编写）¶

实验结果与分析¶

初始阈值 $T_{init}$ 对迭代法的影响¶

T1 与 T2 模态的分割效果差异¶

Otsu 法与迭代法的效果比较¶

算法优缺点说明¶

Project 3: Segmentation¶

实验内容¶

实验原理¶

区域生长 (Region Growing)¶

改进 K-means 聚类¶

分水岭算法 (Watershed)¶

实验结果与分析¶

区域生长算法分析¶

分水岭算法分析¶

改进 K-means 算法参数测试与分析¶

各算法比较与应用场景¶

实验心得¶

《生物医学图像处理》实验五实验报告¶