import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import pearsonr
from sklearn.metrics import mutual_info_score

# ==========================================
# 1. 手动实现的相似度计算函数 (基于 Mask)
# ==========================================

def calc_ssd_manual(img1, img2, mask):
    """手动计算 SSD (Sum of Squared Differences)"""
    # 提取 mask 内的体素
    v1 = img1[mask]
    v2 = img2[mask]
    return np.sum((v1 - v2) ** 2)

def calc_riu_manual(img1, img2, mask, epsilon=1e-8):
    """手动计算 RIU (Ratio Image Uniformity)"""
    v1 = img1[mask]
    v2 = img2[mask]
    # 避免除以0
    ratio = v1 / (v2 + epsilon)
    return np.std(ratio) / (np.mean(ratio) + epsilon)

def calc_cc_manual(img1, img2, mask):
    """手动计算 CC (Cross-Correlation)"""
    v1 = img1[mask]
    v2 = img2[mask]
    
    mean1 = np.mean(v1)
    mean2 = np.mean(v2)
    
    numerator = np.sum((v1 - mean1) * (v2 - mean2))
    denominator = np.sqrt(np.sum((v1 - mean1)**2) * np.sum((v2 - mean2)**2))
    
    if denominator == 0:
        return 0
    return numerator / denominator

def calc_mi_manual(img1, img2, mask, bins=256):
    """手动计算 MI (Mutual Information)"""
    v1 = img1[mask]
    v2 = img2[mask]
    
    # 计算联合直方图
    hist_2d, _, _ = np.histogram2d(v1, v2, bins=bins)
    
    # 将直方图转换为概率 p(a,b)
    pxy = hist_2d / float(np.sum(hist_2d))
    
    # 计算边缘概率 p(a) 和 p(b)
    px = np.sum(pxy, axis=1) # 对 y 积分
    py = np.sum(pxy, axis=0) # 对 x 积分
    
    # 广播计算 p(a)p(b)
    px_py = px[:, None] * py[None, :]
    
    # 仅计算非零项的 MI 避免 log(0)
    non_zero = pxy > 0
    mi = np.sum(pxy[non_zero] * np.log(pxy[non_zero] / px_py[non_zero]))
    return mi

# ==========================================
# 2. 标准库实现的对比函数
# ==========================================

def calc_metrics_library(img1, img2, mask, bins=256):
    """使用标准库计算 SSD, CC, MI (RIU除外)"""
    v1 = img1[mask]
    v2 = img2[mask]
    
    # SSD: numpy
    ssd_lib = np.sum((v1 - v2) ** 2)
    
    # CC: scipy
    cc_lib, _ = pearsonr(v1, v2)
    
    # MI: scikit-learn
    # sklearn 需要一维直方图的输入，或者自己分箱
    # 离散化数据以匹配 sklearn 的行为
    c1 = np.digitize(v1, np.histogram_bin_edges(v1, bins=bins))
    c2 = np.digitize(v2, np.histogram_bin_edges(v2, bins=bins))
    mi_lib = mutual_info_score(c1, c2)
    
    return ssd_lib, cc_lib, mi_lib

# ==========================================
# 3. 主程序执行模块
# ==========================================

def main():
    # 1. 使用 nibabel 读取 .nii 文件 
    t1_img = nib.load('T1w.nii').get_fdata()
    t2_img = nib.load('T2w.nii').get_fdata()
    
    # 2. 提取并绘制一个轴状切片 (Axial Slice) 
    # 假设轴状切片在第三个维度 (Z轴)，我们取中间的一层
    slice_idx = t1_img.shape[2] // 2
    
    
    t1_slice = t1_img[:, slice_idx, :] 
    t2_slice = t2_img[:, slice_idx, :]
    
    # 3. 创建 Mask 
    # 简单的阈值法，过滤掉背景噪声 (可根据实际图像强度调整)
    threshold = 10 
    mask = (t1_slice > threshold) & (t2_slice > threshold)
    
    # 可视化切片和Mask
    plt.figure(figsize=(15, 5))
    plt.subplot(1, 3, 1)
    plt.title('T1w Slice')
    plt.imshow(t1_slice, cmap='gray')
    
    plt.subplot(1, 3, 2)
    plt.title('T2w Slice')
    plt.imshow(t2_slice, cmap='gray')
    
    plt.subplot(1, 3, 3)
    plt.title('Mask')
    plt.imshow(mask, cmap='gray')
    plt.show()
    
    # 4. 定义需要比较的组合
    # 组合 1: T1w 和它自己 (Identity) 
    # 组合 2: T1w 和 T2w 
    comparisons = [
        ("T1w vs T1w (Identity)", t1_slice, t1_slice),
        ("T1w vs T2w", t1_slice, t2_slice)
    ]
    
    for name, imgA, imgB in comparisons:
        print(f"\n--- {name} ---")
        
        # 手动计算 
        ssd_m = calc_ssd_manual(imgA, imgB, mask)
        riu_m = calc_riu_manual(imgA, imgB, mask)
        cc_m  = calc_cc_manual(imgA, imgB, mask)
        mi_m  = calc_mi_manual(imgA, imgB, mask)
        
        # 库计算 
        ssd_l, cc_l, mi_l = calc_metrics_library(imgA, imgB, mask)
        
        print(f"SSD: Manual = {ssd_m:.2f} | Library = {ssd_l:.2f}")
        print(f"RIU: Manual = {riu_m:.4f} | Library = (N/A)")
        print(f"CC : Manual = {cc_m:.4f}  | Library = {cc_l:.4f}")
        print(f"MI : Manual = {mi_m:.4f}  | Library = {mi_l:.4f}")

if __name__ == "__main__":
    main()

--- T1w vs T1w (Identity) ---
SSD: Manual = 0.00 | Library = 0.00
RIU: Manual = 0.0000 | Library = (N/A)
CC : Manual = 1.0000  | Library = 1.0000
MI : Manual = 4.8985  | Library = 4.8985

--- T1w vs T2w ---
SSD: Manual = 252153445.00 | Library = 252153445.00
RIU: Manual = 0.4371 | Library = (N/A)
CC : Manual = -0.1079  | Library = -0.1079
MI : Manual = 0.8779  | Library = 0.8780

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import warp, AffineTransform

# ==========================================
# 1. 基础变换与辅助函数
# ==========================================
def apply_rigid_transform(img, theta_deg, tx, ty):
    """施加刚性变换"""
    center_y, center_x = np.array(img.shape) / 2.0
    t_center = AffineTransform(translation=(-center_x, -center_y))
    rot = AffineTransform(rotation=np.deg2rad(theta_deg))
    t_uncenter_translate = AffineTransform(translation=(center_x + tx, center_y + ty))
    
    matrix = t_center.params
    matrix = rot.params @ matrix
    matrix = t_uncenter_translate.params @ matrix
    return warp(img, AffineTransform(matrix=matrix).inverse, preserve_range=True)

def calc_overlap_dice(mask1, mask2):
    intersection = np.logical_and(mask1, mask2)
    if mask1.sum() + mask2.sum() == 0: return 0
    return 2. * intersection.sum() / (mask1.sum() + mask2.sum())

def create_overlap_image(target, moving, threshold=10):
    h, w = target.shape
    rgb_img = np.zeros((h, w, 3))
    rgb_img[..., 0] = (target > threshold) * 1.0  
    rgb_img[..., 1] = (moving > threshold) * 1.0  
    return rgb_img

# ==========================================
# 2. 相似度指标定义 (全部转换为 Loss 函数，越小越好)
# ==========================================

def loss_ssd(target, moving, threshold=10):
    """SSD 误差平方和 (直接作为 Loss)"""
    mask = target > threshold
    if mask.sum() == 0: return float('inf')
    return np.mean((target[mask] - moving[mask]) ** 2)

def loss_riu(target, moving, threshold=10):
    """RIU 比值均匀度 (标准差/均值，直接作为 Loss)"""
    mask = target > threshold
    t_vals = target[mask]
    m_vals = moving[mask]
    
    # 避免除以 0 的情况
    valid = m_vals > 1e-5
    if valid.sum() < 500: return float('inf')
    
    ratio = t_vals[valid] / m_vals[valid]
    mu = np.mean(ratio)
    if mu == 0: return float('inf')
    
    return np.std(ratio) / mu

def loss_cc(target, moving, threshold=10):
    """CC 互相关 (取 -CC^2 作为 Loss，使得正/负强相关都趋向极小值)"""
    mask = target > threshold
    if mask.sum() < 500: return 0.0 
    
    v_t = target[mask]
    v_m = moving[mask]
    m_t, m_m = np.mean(v_t), np.mean(v_m)
    
    num = np.sum((v_t - m_t) * (v_m - m_m))
    den = np.sqrt(np.sum((v_t - m_t)**2) * np.sum((v_m - m_m)**2))
    
    if den == 0: return 0.0
    cc = num / den
    return -(cc ** 2)

# ==========================================
# 3. 归一化梯度的最速下降法
# ==========================================
def basic_steepest_descent(target, moving, metric_fn, iterations=150, initial_lr=2.0):
    """基础最速下降法"""
    current_params = np.array([0.0, 0.0, 0.0]) # [theta, tx, ty]
    delta = 0.5 # 数值求导的微元
    
    for i in range(iterations):
        theta, tx, ty = current_params
        lr = initial_lr * (1.0 - i / iterations) # 线性衰减学习率
        
        # 计算数值梯度
        l_theta_plus = metric_fn(target, apply_rigid_transform(moving, theta + delta, tx, ty))
        l_theta_minus = metric_fn(target, apply_rigid_transform(moving, theta - delta, tx, ty))
        grad_theta = (l_theta_plus - l_theta_minus) / (2 * delta)
        
        l_tx_plus = metric_fn(target, apply_rigid_transform(moving, theta, tx + delta, ty))
        l_tx_minus = metric_fn(target, apply_rigid_transform(moving, theta, tx - delta, ty))
        grad_tx = (l_tx_plus - l_tx_minus) / (2 * delta)
        
        l_ty_plus = metric_fn(target, apply_rigid_transform(moving, theta, tx, ty + delta))
        l_ty_minus = metric_fn(target, apply_rigid_transform(moving, theta, tx, ty - delta))
        grad_ty = (l_ty_plus - l_ty_minus) / (2 * delta)
        
        grad = np.array([grad_theta, grad_tx, grad_ty])
        
        # 【极其重要】：因为 SSD 的梯度可能上万，CC 的梯度只有零点几
        # 这里必须将梯度方向归一化，让步伐大小完全由 lr 决定，否则优化器必崩
        norm = np.linalg.norm(grad)
        if norm > 1e-8:
            grad = grad / norm
            
        current_params = current_params - lr * grad
        
    return current_params

# ==========================================
# 4. 主程序：批量遍历所有测试组合
# ==========================================
def main():
    print("加载并提取数据...")
    t1_img = nib.load('T1w.nii').get_fdata()
    t2_img = nib.load('T2w.nii').get_fdata()
    slice_idx = t1_img.shape[2] // 2
    
    t1_slice = t1_img[:, slice_idx, :] 
    t2_slice = t2_img[:, slice_idx, :]
    
    true_theta, true_tx, true_ty = 15.0, 10.0, -5.0
    print(f"\n[施加的真实变换] 角度: {true_theta}°, Tx: {true_tx}, Ty: {true_ty}")
    
    target_t1 = t1_slice
    moving_t1_warped = apply_rigid_transform(t1_slice, true_theta, true_tx, true_ty)
    moving_t2_warped = apply_rigid_transform(t2_slice, true_theta, true_tx, true_ty)

    # 实验要求：组合所有除了 MI 的指标
    metrics_to_test = {
        'SSD': loss_ssd,
        'RIU': loss_riu,
        'CC': loss_cc
    }
    
    tasks = [
        ("Intra-modality (T1w to T1w)", target_t1, moving_t1_warped),
        ("Inter-modality (T2w to T1w)", target_t1, moving_t2_warped)
    ]
    
    for task_name, target, moving in tasks:
        print(f"\n=======================================================")
        print(f"执行任务: {task_name}")
        print(f"=======================================================")
        
        for metric_name, metric_fn in metrics_to_test.items():
            print(f"\n---> 使用指标: {metric_name}")
            
            mask_target = target > 10
            overlap_before = calc_overlap_dice(mask_target, moving > 10)
            
            # 执行最速下降法
            final_params = basic_steepest_descent(target, moving, metric_fn, iterations=150, initial_lr=2.0)
            calc_theta, calc_tx, calc_ty = final_params
            
            registered_img = apply_rigid_transform(moving, calc_theta, calc_tx, calc_ty)
            overlap_after = calc_overlap_dice(mask_target, registered_img > 10)
            
            print(f"  配准前 Dice: {overlap_before:.4f} | 配准后 Dice: {overlap_after:.4f}")
            print(f"  恢复参数: 角度 {calc_theta:.2f}°, Tx {calc_tx:.2f}, Ty {calc_ty:.2f}")
            print(f"  残余误差: 角度 {abs(true_theta + calc_theta):.2f}°, Tx {abs(true_tx + calc_tx):.2f}, Ty {abs(true_ty + calc_ty):.2f}")
            
            # 画图
            plt.figure(figsize=(15, 4))
            plt.suptitle(f"{task_name} | Metric: {metric_name}", fontsize=14)
            plt.subplot(1, 4, 1); plt.title('Target'); plt.imshow(target, cmap='gray'); plt.axis('off')
            plt.subplot(1, 4, 2); plt.title('Moving (Before)'); plt.imshow(moving, cmap='gray'); plt.axis('off')
            plt.subplot(1, 4, 3); plt.title('Registered (After)'); plt.imshow(registered_img, cmap='gray'); plt.axis('off')
            plt.subplot(1, 4, 4); plt.title('Overlap'); plt.imshow(create_overlap_image(target, registered_img)); plt.axis('off')
            plt.tight_layout()
            plt.show()

if __name__ == "__main__":
    main()

加载并提取数据...

[施加的真实变换] 角度: 15.0°, Tx: 10.0, Ty: -5.0

=======================================================
执行任务: Intra-modality (T1w to T1w)
=======================================================

---> 使用指标: SSD
  配准前 Dice: 0.9262 | 配准后 Dice: 0.9893
  恢复参数: 角度 -15.01°, Tx -8.38, Ty 7.44
  残余误差: 角度 0.01°, Tx 1.62, Ty 2.44

---> 使用指标: RIU
  配准前 Dice: 0.9262 | 配准后 Dice: 0.8706
  恢复参数: 角度 -22.66°, Tx 1.33, Ty -4.43
  残余误差: 角度 7.66°, Tx 11.33, Ty 9.43

---> 使用指标: CC
  配准前 Dice: 0.9262 | 配准后 Dice: 0.9892
  恢复参数: 角度 -14.99°, Tx -8.37, Ty 7.44
  残余误差: 角度 0.01°, Tx 1.63, Ty 2.44

=======================================================
执行任务: Inter-modality (T2w to T1w)
=======================================================

---> 使用指标: SSD
  配准前 Dice: 0.9199 | 配准后 Dice: 0.9877
  恢复参数: 角度 -11.83°, Tx -8.01, Ty 8.33
  残余误差: 角度 3.17°, Tx 1.99, Ty 3.33

---> 使用指标: RIU
  配准前 Dice: 0.9199 | 配准后 Dice: 0.8879
  恢复参数: 角度 -3.83°, Tx 3.14, Ty -3.03
  残余误差: 角度 11.17°, Tx 13.14, Ty 8.03

---> 使用指标: CC
  配准前 Dice: 0.9199 | 配准后 Dice: 0.9735
  恢复参数: 角度 -7.89°, Tx -6.45, Ty 6.04
  残余误差: 角度 7.11°, Tx 3.55, Ty 1.04

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import warp, AffineTransform
from scipy.optimize import minimize

# ==========================================
# 1. 刚性变换工具与辅助函数
# ==========================================

def apply_rigid_transform(img, theta_deg, tx, ty):
    """施加刚性变换 (绕图像中心旋转，并平移)"""
    center_y, center_x = np.array(img.shape) / 2.0
    t_center = AffineTransform(translation=(-center_x, -center_y))
    rot = AffineTransform(rotation=np.deg2rad(theta_deg))
    t_uncenter_translate = AffineTransform(translation=(center_x + tx, center_y + ty))
    
    matrix = t_center.params
    matrix = rot.params @ matrix
    matrix = t_uncenter_translate.params @ matrix
    
    t_final = AffineTransform(matrix=matrix)
    return warp(img, t_final.inverse, preserve_range=True, order=1)

def calc_overlap_dice(mask1, mask2):
    intersection = np.logical_and(mask1, mask2)
    if mask1.sum() + mask2.sum() == 0: return 0
    return 2. * intersection.sum() / (mask1.sum() + mask2.sum())

def create_overlap_image(target, moving, threshold=10):
    h, w = target.shape
    rgb_img = np.zeros((h, w, 3))
    rgb_img[..., 0] = (target > threshold) * 1.0  
    rgb_img[..., 1] = (moving > threshold) * 1.0  
    return rgb_img

# ==========================================
# 2. 相似度度量（因地制宜）
# ==========================================

def loss_ssd(target, moving, threshold=10):
    """【专供模态内】固定以 target 为掩膜，防止旋转带来的边缘空洞干扰"""
    mask = target > threshold 
    if mask.sum() == 0: return float('inf')
    return np.mean((target[mask] - moving[mask]) ** 2)

def loss_mi(target, moving, threshold=10, bins=32):
    """【专供跨模态】互信息"""
    mask = (target > threshold) & (moving > threshold)
    if mask.sum() < 100: return float('inf')
    v_t = target[mask]
    v_m = moving[mask]
    hist_2d, _, _ = np.histogram2d(v_t, v_m, bins=bins)
    pxy = hist_2d / float(np.sum(hist_2d))
    px = np.sum(pxy, axis=1)
    py = np.sum(pxy, axis=0)
    px_py = px[:, None] * py[None, :]
    non_zero = pxy > 0
    mi = np.sum(pxy[non_zero] * np.log(pxy[non_zero] / (px_py[non_zero] + 1e-8)))
    return -mi  

# ==========================================
# 3. 不同的优化器实现（因地制宜）
# ==========================================

def run_powell_ssd(target, moving):
    """【方法 A 升级版】同模态也使用无梯度的 Powell 算法，追求极致精度"""
    print("    -> 正在执行 Powell 行搜索 (针对 SSD)...")
    def obj(p): 
        # Powell 寻找极小值，SSD 本身就是越小越好
        return loss_ssd(target, apply_rigid_transform(moving, p[0], p[1], p[2]))
    
    # 初始搜索步长
    init_directions = np.eye(3) * [5.0, 5.0, 5.0]
    result = minimize(obj, [0.0, 0.0, 0.0], method='Powell', 
                      options={'xtol': 1e-3, 'ftol': 1e-3, 'direc': init_directions, 'disp': False})
    return result.x

def run_powell_mi(target, moving):
    """【方法 B】专门为跨模态准备：无梯度 Powell 算法"""
    print("    -> 正在执行 Powell 行搜索...")
    def obj(p): return loss_mi(target, apply_rigid_transform(moving, p[0], p[1], p[2]))
    
    # 允许较大的初始搜索步长
    init_directions = np.eye(3) * [5.0, 5.0, 5.0]
    result = minimize(obj, [0.0, 0.0, 0.0], method='Powell', 
                      options={'xtol': 1e-2, 'ftol': 1e-2, 'direc': init_directions, 'disp': False})
    return result.x

# ==========================================
# 4. 主程序
# ==========================================

def main():
    print("读取数据...")
    t1_img = nib.load('T1w.nii').get_fdata()
    t2_img = nib.load('T2w.nii').get_fdata()
    
    # 严格采用指定的切片方式
    slice_idx = t1_img.shape[2] // 2
    t1_slice = t1_img[:, slice_idx, :] 
    t2_slice = t2_img[:, slice_idx, :]
    
    # 真实畸变参数
    true_theta, true_tx, true_ty = 15.0, 10.0, -5.0
    print(f"\n[施加的真实畸变] 角度: {true_theta}°, Tx: {true_tx}, Ty: {true_ty}")
    
    target_t1 = t1_slice
    moving_t1_warped = apply_rigid_transform(t1_slice, true_theta, true_tx, true_ty) 
    moving_t2_warped = apply_rigid_transform(t2_slice, true_theta, true_tx, true_ty) 

    # 核心：将不同的任务解耦，分派给最适合它们的算法！
    results = []
    
    print("\n========== 任务 1: 模态内 (T1w to T1w) ==========")
    print("策略: 误差平方和 (SSD) + 无梯度搜索 (Powell)")
    # 把原来的 run_gradient_descent_ssd 换成 run_powell_ssd
    params_intra = run_powell_ssd(target_t1, moving_t1_warped)
    results.append(("Intra-modality (T1w to T1w) [SSD + Powell]", target_t1, moving_t1_warped, params_intra))
    
    print("\n========== 任务 2: 跨模态 (T2w to T1w) ==========")
    print("策略: 互信息 (MI) + 无梯度搜索 (Powell)")
    params_inter = run_powell_mi(target_t1, moving_t2_warped)
    results.append(("Inter-modality (T2w to T1w) [MI + Powell]", target_t1, moving_t2_warped, params_inter))

    # 评估与可视化
    for name, target, moving, final_params in results:
        print(f"\n--- {name} 结果评估 ---")
        calc_theta, calc_tx, calc_ty = final_params
        registered_img = apply_rigid_transform(moving, calc_theta, calc_tx, calc_ty)
        
        print(f"配准前重叠度 (Dice): {calc_overlap_dice(target>10, moving>10):.4f}")
        print(f"配准后重叠度 (Dice): {calc_overlap_dice(target>10, registered_img>10):.4f}")
        print(f"[恢复变换] 角度: {calc_theta:.2f}°, Tx: {calc_tx:.2f}, Ty: {calc_ty:.2f}")
        print(f"残余误差: 角度 {abs(true_theta+calc_theta):.2f}°, Tx {abs(true_tx+calc_tx):.2f}, Ty {abs(true_ty+calc_ty):.2f}")
        
        # 画图
        plt.figure(figsize=(15, 4))
        plt.suptitle(name, fontsize=14)
        
        plt.subplot(1, 4, 1); plt.title('Target'); plt.imshow(target, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 2); plt.title('Moving'); plt.imshow(moving, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 3); plt.title('Registered'); plt.imshow(registered_img, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 4); plt.title('Overlap'); plt.imshow(create_overlap_image(target, registered_img)); plt.axis('off')
        
        plt.tight_layout()
        plt.show()

if __name__ == "__main__":
    main()

读取数据...

[施加的真实畸变] 角度: 15.0°, Tx: 10.0, Ty: -5.0

========== 任务 1: 模态内 (T1w to T1w) ==========
策略: 误差平方和 (SSD) + 无梯度搜索 (Powell)
    -> 正在执行 Powell 行搜索 (针对 SSD)...

========== 任务 2: 跨模态 (T2w to T1w) ==========
策略: 互信息 (MI) + 无梯度搜索 (Powell)
    -> 正在执行 Powell 行搜索...

--- Intra-modality (T1w to T1w) [SSD + Powell] 结果评估 ---
配准前重叠度 (Dice): 0.9262
配准后重叠度 (Dice): 0.9892
[恢复变换] 角度: -15.00°, Tx: -8.38, Ty: 7.43
残余误差: 角度 0.00°, Tx 1.62, Ty 2.43

--- Inter-modality (T2w to T1w) [MI + Powell] 结果评估 ---
配准前重叠度 (Dice): 0.9199
配准后重叠度 (Dice): 0.9950
[恢复变换] 角度: -14.72°, Tx: -8.44, Ty: 7.35
残余误差: 角度 0.28°, Tx 1.56, Ty 2.35

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import warp, AffineTransform
import SimpleITK as sitk

# ==========================================
# 1. 人为施加刚性畸变
# ==========================================
def apply_rigid_transform_skimage(img, theta_deg, tx, ty):
    center_y, center_x = np.array(img.shape) / 2.0
    t_center = AffineTransform(translation=(-center_x, -center_y))
    rot = AffineTransform(rotation=np.deg2rad(theta_deg))
    t_uncenter_translate = AffineTransform(translation=(center_x + tx, center_y + ty))
    
    matrix = t_center.params
    matrix = rot.params @ matrix
    matrix = t_uncenter_translate.params @ matrix
    t_final = AffineTransform(matrix=matrix)
    return warp(img, t_final.inverse, preserve_range=True, order=1)

# ==========================================
# 2. 评估辅助函数 (计算四个指标)
# ==========================================
def calc_overlap_dice(mask1, mask2):
    intersection = np.logical_and(mask1, mask2)
    if mask1.sum() + mask2.sum() == 0: return 0
    return 2. * intersection.sum() / (mask1.sum() + mask2.sum())

def create_overlap_image(target, moving, threshold=10):
    h, w = target.shape
    rgb_img = np.zeros((h, w, 3))
    rgb_img[..., 0] = (target > threshold) * 1.0  
    rgb_img[..., 1] = (moving > threshold) * 1.0  
    return rgb_img

def calc_ssd(target, moving, threshold=10):
    """计算 SSD (Sum of Squared Differences)"""
    mask = target > threshold 
    if mask.sum() == 0: return float('inf')
    return np.mean((target[mask] - moving[mask]) ** 2)

def calc_cc(target, moving, threshold=10):
    """计算 CC (Cross-Correlation)"""
    mask = target > threshold
    if mask.sum() == 0: return 0.0
    t_val = target[mask]
    m_val = moving[mask]
    
    t_mean = np.mean(t_val)
    m_mean = np.mean(m_val)
    
    numerator = np.sum((t_val - t_mean) * (m_val - m_mean))
    denominator = np.sqrt(np.sum((t_val - t_mean)**2) * np.sum((m_val - m_mean)**2))
    
    if denominator == 0: return 0.0
    return numerator / denominator

def calc_riu(target, moving, threshold=10):
    """计算 RIU (Ratio Image Uniformity)"""
    mask = target > threshold
    if mask.sum() == 0: return float('inf')
    t_val = target[mask]
    m_val = moving[mask]
    
    ratio = t_val / (m_val + 1e-8)
    std_r = np.std(ratio)
    mean_r = np.mean(ratio)
    
    if mean_r == 0: return float('inf')
    return std_r / mean_r

# ==========================================
# 3. 核心工具包封装：SimpleITK 自动配准框架
# ==========================================
def register_simpleitk(fixed_arr, moving_arr, metric_type='ssd'):
    fixed_image = sitk.Cast(sitk.GetImageFromArray(fixed_arr), sitk.sitkFloat32)
    moving_image = sitk.Cast(sitk.GetImageFromArray(moving_arr), sitk.sitkFloat32)

    R = sitk.ImageRegistrationMethod()

    if metric_type == 'ssd':
        R.SetMetricAsMeanSquares()
    elif metric_type == 'cc':
        R.SetMetricAsCorrelation()
    else:
        raise ValueError("不支持该指标作为 SITK 优化器目标")

    R.SetOptimizerAsRegularStepGradientDescent(
        learningRate=1.0, 
        minStep=1e-5,
        numberOfIterations=250,
        gradientMagnitudeTolerance=1e-8
    )
    
    R.SetOptimizerScalesFromPhysicalShift()

    R.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
    R.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
    R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    initial_transform = sitk.Euler2DTransform()
    center = [c / 2.0 for c in fixed_image.GetSize()]
    initial_transform.SetCenter(center)
    R.SetInitialTransform(initial_transform, inPlace=True)

    mask_arr = (fixed_arr > 10).astype(np.uint8)
    fixed_mask = sitk.GetImageFromArray(mask_arr)
    R.SetMetricFixedMask(fixed_mask)

    R.SetInterpolator(sitk.sitkLinear)

    final_transform = R.Execute(fixed_image, moving_image)

    resampled = sitk.Resample(
        moving_image, fixed_image, final_transform,
        sitk.sitkLinear, 0.0, moving_image.GetPixelID()
    )
    registered_arr = sitk.GetArrayFromImage(resampled)

    theta_deg = np.rad2deg(initial_transform.GetAngle())
    tx, ty = initial_transform.GetTranslation()

    return registered_arr, theta_deg, tx, ty

# ==========================================
# 4. 主程序
# ==========================================
def main():
    try:
        t1_img = nib.load('T1w.nii').get_fdata()
        t2_img = nib.load('T2w.nii').get_fdata()
    except FileNotFoundError:
        print("未找到 NIfTI 文件，请确保 T1w.nii 和 T2w.nii 存在于同一目录下。")
        return
    
    slice_idx = t1_img.shape[2] // 2
    t1_slice = t1_img[:, slice_idx, :] 
    t2_slice = t2_img[:, slice_idx, :]
    
    true_theta, true_tx, true_ty = 15.0, 10.0, -5.0
    print(f"\n[人为施加的畸变] 角度: {true_theta}°, Tx: {true_tx}, Ty: {true_ty}")
    
    target_t1 = t1_slice
    moving_t1_warped = apply_rigid_transform_skimage(t1_slice, true_theta, true_tx, true_ty) 
    moving_t2_warped = apply_rigid_transform_skimage(t2_slice, true_theta, true_tx, true_ty) 

    # 任务列表：移除了所有主观描述后缀
    tasks = [
        ("Intra-modality (T1w to T1w) [Optimizer: SSD]", target_t1, moving_t1_warped, 'ssd'),
        ("Intra-modality (T1w to T1w) [Optimizer: CC]", target_t1, moving_t1_warped, 'cc'),
        ("Inter-modality (T2w to T1w) [Optimizer: SSD]", target_t1, moving_t2_warped, 'ssd'), 
        ("Inter-modality (T2w to T1w) [Optimizer: CC]", target_t1, moving_t2_warped, 'cc') 
    ]
    
    for name, target, moving, optim_metric in tasks:
        print(f"\n========== {name} ==========")
        
        dice_b = calc_overlap_dice(target>10, moving>10)
        print("【配准前指标】")
        print(f" ┣ Dice 重叠度 : {dice_b:.4f}")
        print(f" ┣ SSD : {calc_ssd(target, moving):.2f}")
        print(f" ┣ CC  : {calc_cc(target, moving):.4f}")
        print(f" ┗ RIU : {calc_riu(target, moving):.4f}")
        
        registered_img, calc_theta, calc_tx, calc_ty = register_simpleitk(target, moving, metric_type=optim_metric)
        
        dice_a = calc_overlap_dice(target>10, registered_img>10)
        print("\n【配准后指标】")
        print(f" ┣ Dice 重叠度 : {dice_a:.4f}")
        print(f" ┣ SSD : {calc_ssd(target, registered_img):.2f}")
        print(f" ┣ CC  : {calc_cc(target, registered_img):.4f}")
        print(f" ┗ RIU : {calc_riu(target, registered_img):.4f}")
        
        print(f"\n[SITK求得参数] 角度: {calc_theta:.2f}°, Tx: {calc_tx:.2f}, Ty: {calc_ty:.2f}")
        print(f"[参数残差误差] 角度差: {abs(true_theta - calc_theta):.2f}°, Tx差: {abs(true_tx - calc_tx):.2f}, Ty差: {abs(true_ty - calc_ty):.2f}")
        
        plt.figure(figsize=(15, 4))
        plt.suptitle(name, fontsize=14)
        
        plt.subplot(1, 4, 1); plt.title('Target'); plt.imshow(target, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 2); plt.title('Moving'); plt.imshow(moving, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 3); plt.title('Registered'); plt.imshow(registered_img, cmap='gray'); plt.axis('off')
        plt.subplot(1, 4, 4); plt.title('Overlap'); plt.imshow(create_overlap_image(target, registered_img)); plt.axis('off')
        
        plt.tight_layout()
        plt.show()

if __name__ == "__main__":
    main()

[人为施加的畸变] 角度: 15.0°, Tx: 10.0, Ty: -5.0

========== Intra-modality (T1w to T1w) [Optimizer: SSD] ==========
【配准前指标】
 ┣ Dice 重叠度 : 0.9262
 ┣ SSD : 5229.49
 ┣ CC  : 0.1342
 ┗ RIU : 3.8503

【配准后指标】
 ┣ Dice 重叠度 : 0.9893
 ┣ SSD : 64.53
 ┣ CC  : 0.9790
 ┗ RIU : 0.0521

[SITK求得参数] 角度: 15.00°, Tx: 10.02, Ty: -5.01
[参数残差误差] 角度差: 0.00°, Tx差: 0.02, Ty差: 0.01

========== Intra-modality (T1w to T1w) [Optimizer: CC] ==========
【配准前指标】
 ┣ Dice 重叠度 : 0.9262
 ┣ SSD : 5229.49
 ┣ CC  : 0.1342
 ┗ RIU : 3.8503

【配准后指标】
 ┣ Dice 重叠度 : 0.9893
 ┣ SSD : 64.53
 ┣ CC  : 0.9790
 ┗ RIU : 0.0520

[SITK求得参数] 角度: 15.00°, Tx: 10.02, Ty: -5.01
[参数残差误差] 角度差: 0.00°, Tx差: 0.02, Ty差: 0.01

========== Inter-modality (T2w to T1w) [Optimizer: SSD] ==========
【配准前指标】
 ┣ Dice 重叠度 : 0.9199
 ┣ SSD : 14473.04
 ┣ CC  : 0.1533
 ┗ RIU : 3.3976

【配准后指标】
 ┣ Dice 重叠度 : 0.9861
 ┣ SSD : 11806.66
 ┣ CC  : 0.1737
 ┗ RIU : 9.9333

[SITK求得参数] 角度: 18.28°, Tx: 11.66, Ty: -5.65
[参数残差误差] 角度差: 3.28°, Tx差: 1.66, Ty差: 0.65

========== Inter-modality (T2w to T1w) [Optimizer: CC] ==========
【配准前指标】
 ┣ Dice 重叠度 : 0.9199
 ┣ SSD : 14473.04
 ┣ CC  : 0.1533
 ┗ RIU : 3.3976

【配准后指标】
 ┣ Dice 重叠度 : 0.9721
 ┣ SSD : 12272.27
 ┣ CC  : 0.2171
 ┗ RIU : 6.9594

[SITK求得参数] 角度: 7.76°, Tx: 6.95, Ty: -5.04
[参数残差误差] 角度差: 7.24°, Tx差: 3.05, Ty差: 0.04

对比维度	最速下降法 (SGD)	Powell法 (本次Bonus采用)
梯度依赖	强。需通过有限差分计算数值梯度，在离散的图像空间中，梯度极易受噪声和插值伪影干扰而失真。	无。直接比较目标函数值（SSD/MI），鲁棒性强。
超参数敏感度	极高。步长（学习率）过大易震荡，过小易陷入局部最优或迭代缓慢。	低。自适应的线搜索机制，只需给定初始方向和容差阈值即可。
对MI的适应性	不可用。互信息的计算涉及直方图分箱（Binning），函数极不平滑，梯度下降在此处失效。	极佳。不关心函数的平滑度，非常适合求解互信息这类复杂、非凸的目标函数。

任务类型	相似度指标	优化器框架	角度残差 ($^{\circ}$)	Tx 残差 (像素)	Ty 残差 (像素)
同模态 (T1w-T1w)	SSD	手动 SGD	0.01	1.62	2.44
		SITK (Project 3)	0.00	0.02	0.01
	CC	手动 SGD	0.01	1.63	2.44
		SITK (Project 3)	0.00	0.02	0.01
跨模态 (T2w-T1w)	SSD	手动 SGD	3.17	1.99	3.33
		SITK (Project 3)	3.28	1.66	0.65
	CC	手动 SGD	7.11	3.55	1.04
		SITK (Project 3)	7.24	3.05	0.04

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

Project 1: Similarity Measure¶

实验内容¶

实验原理¶

1. SSD (Sum of Squared Differences, 误差平方和)¶

2. RIU (Ratio Image Uniformity, 比值图像均匀度)¶

3. CC (Cross-Correlation, 互相关系数)¶

4. MI (Mutual Information, 互信息)¶

实验过程（代码编写）¶

实验结果与分析¶

模态内比较：T1w vs T1w (恒等图像)¶

跨模态比较：T1w vs T2w¶

Project 2:Image Registration¶

实验内容¶

实验原理¶

1. 最速下降法 (Steepest Descent Method)¶

2. Dice 系数 (Dice Coefficient)¶

3. 讨论题：不能选择MI作为相似度指标的原因¶

实验过程（代码编写）¶

实验结果与分析¶

模态内配准分析 (Intra-modality: T1w to T1w)¶

跨模态配准分析 (Inter-modality: T2w to T1w)¶

Bonus¶

Bonus实验结果与分析¶

优化策略¶

实验结果评价与配准效果优异的原因分析¶

与SGD的对比¶

Project 3: Tool-based Co-registration¶

实验内容¶

实验原理¶

面向对象的配准框架 (Image Registration Method)¶

Mattes 互信息 (Mattes Mutual Information)¶

正则化步长梯度下降 (Regular Step Gradient Descent)¶

多分辨率金字塔策略 (Multi-Resolution Pyramid)¶

物理空间参数缩放 (Physical Shift Scaling)¶

实验过程（代码编写）¶

实验结果与分析¶

SITK与手动SGD的对比¶

实验结果分析¶

同模态配准¶

跨模态配准¶

实验心得¶

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