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

# ==========================================
# 1. 核心矩阵运算函数 (Matrix Operator)
# ==========================================
def matrix_rotate(img, theta, method='nearest'):
    """
    使用矩阵运算实现图像的反向映射旋转
    """
    H, W = img.shape
    cx, cy = (W - 1) / 2.0, (H - 1) / 2.0 

    # 生成目标图像的所有坐标点网格
    x = np.arange(W)
    y = np.arange(H)
    xx, yy = np.meshgrid(x, y)

    # 平移坐标系，使中心为原点
    xx_c = xx - cx
    yy_c = yy - cy

    # 展平坐标并堆叠成 2xN 矩阵
    coords = np.vstack((xx_c.flatten(), yy_c.flatten()))

    # 构建反向映射的旋转矩阵 R(-theta)
    cos_t = np.cos(theta)
    sin_t = np.sin(theta)
    R_inv = np.array([[cos_t, sin_t],
                      [-sin_t, cos_t]])

    # 矩阵乘法计算原图坐标
    src_coords = R_inv @ coords

    # 平移回原坐标系并重塑为图像形状
    src_x = (src_coords[0, :] + cx).reshape(H, W)
    src_y = (src_coords[1, :] + cy).reshape(H, W)

    out_img = np.zeros_like(img, dtype=np.float32)

    # 标记在原图范围内的有效坐标点 (防止越界)
    valid = (src_x >= 0) & (src_x <= W - 1) & (src_y >= 0) & (src_y <= H - 1)

    if method == 'nearest':
        # 最近邻插值
        src_x_round = np.round(src_x).astype(int)
        src_y_round = np.round(src_y).astype(int)
        out_img[valid] = img[src_y_round[valid], src_x_round[valid]]

    elif method == 'bilinear':
        # 双线性插值
        x0 = np.floor(src_x).astype(int)
        x1 = x0 + 1
        y0 = np.floor(src_y).astype(int)
        y1 = y0 + 1

        x0 = np.clip(x0, 0, W - 1)
        x1 = np.clip(x1, 0, W - 1)
        y0 = np.clip(y0, 0, H - 1)
        y1 = np.clip(y1, 0, H - 1)

        wa = (x1 - src_x) * (y1 - src_y)
        wb = (src_x - x0) * (y1 - src_y)
        wc = (x1 - src_x) * (src_y - y0)
        wd = (src_x - x0) * (src_y - y0)

        pixel_vals = (img[y0, x0] * wa + img[y0, x1] * wb + 
                      img[y1, x0] * wc + img[y1, x1] * wd)
        
        out_img[valid] = pixel_vals[valid]

    return out_img

# ==========================================
# 2. 辅助分析函数 (误差计算 & 动态ROI寻找)
# ==========================================
def calculate_mse(original, recovered):
    """精确计算 MSE (去除黑色背景干扰)"""
    mask = original > 1e-5
    return np.mean((original[mask] - recovered[mask])**2)

def find_max_diff_roi(img1, img2, roi_size):
    """通过滑动窗口（均值滤波）自动寻找两张图像差异最大的区域"""
    diff = np.abs(img1 - img2)
    local_mean_diff = cv2.blur(diff, (roi_size, roi_size))
    _, _, _, max_loc = cv2.minMaxLoc(local_mean_diff)
    
    x_center, y_center = max_loc
    x1 = max(0, x_center - roi_size // 2)
    y1 = max(0, y_center - roi_size // 2)
    
    x1 = min(x1, img1.shape[1] - roi_size)
    y1 = min(y1, img1.shape[0] - roi_size)
    return x1, y1

# ==========================================
# 3. 实验主体执行：加载数据与变换
# ==========================================
file_path = os.path.join('lab2_data', 'lab2_MRI.npy')
img = np.load(file_path)
theta = np.pi / 5

# 执行 1 次旋转循环
img_fwd_nn = matrix_rotate(img, theta, method='nearest')
img_rec_nn = matrix_rotate(img_fwd_nn, -theta, method='nearest')

img_fwd_bi = matrix_rotate(img, theta, method='bilinear')
img_rec_bi = matrix_rotate(img_fwd_bi, -theta, method='bilinear')

mse_nn = calculate_mse(img, img_rec_nn)
mse_bi = calculate_mse(img, img_rec_bi)

print("=" * 45)
print("量化指标对比 (1次旋转循环):")
print(f"最近邻插值 -> MSE: {mse_nn:.4f}")
print(f"双线性插值 -> MSE: {mse_bi:.4f}")
print("=" * 45)

# ==========================================
# 4. 自动锁定差异并可视化
# ==========================================
H, W = img.shape
w_roi = int(W * 0.15)
x1, y1 = find_max_diff_roi(img, img_rec_nn, w_roi)
x2, y2 = x1 + w_roi, y1 + w_roi

print(f"已自动锁定最大差异区域: 左上角坐标 ({x1}, {y1}), 窗口大小 {w_roi}x{w_roi}")

fig, axes = plt.subplots(2, 3, figsize=(16, 10))
fig.suptitle('Project 1: Spatial Transform Analysis with Dynamic ROI', fontsize=16)

axes[0, 0].imshow(img, cmap='gray')
axes[0, 0].set_title("Original Image")

axes[0, 1].imshow(img_rec_nn, cmap='gray')
axes[0, 1].set_title(f"Nearest (MSE: {mse_nn:.2f})")

axes[0, 2].imshow(img_rec_bi, cmap='gray')
axes[0, 2].set_title(f"Bilinear (MSE: {mse_bi:.2f})")

for ax in axes[0, :3]:
    rect = patches.Rectangle((x1, y1), w_roi, w_roi, linewidth=1.5, edgecolor='r', facecolor='none')
    ax.add_patch(rect)
    ax.axis('off')

axes[1, 0].imshow(img[y1:y2, x1:x2], cmap='gray')
axes[1, 0].set_title("Original (Zoomed ROI)")

axes[1, 1].imshow(img_rec_nn[y1:y2, x1:x2], cmap='gray')
axes[1, 1].set_title("Nearest (Zoomed ROI)")
axes[1, 1].annotate('Jagged / Artifacts', xy=(w_roi*0.4, w_roi*0.5), xytext=(w_roi*0.1, w_roi*0.8),
                    arrowprops=dict(facecolor='red', shrink=0.05, width=2, headwidth=8), color='red', weight='bold')

axes[1, 2].imshow(img_rec_bi[y1:y2, x1:x2], cmap='gray')
axes[1, 2].set_title("Bilinear (Zoomed ROI)")
axes[1, 2].annotate('Blurred Details', xy=(w_roi*0.4, w_roi*0.5), xytext=(w_roi*0.1, w_roi*0.8),
                    arrowprops=dict(facecolor='yellow', shrink=0.05, width=2, headwidth=8), color='yellow', weight='bold')

for ax in axes[1, :]:
    ax.axis('off')

plt.tight_layout()
plt.show()

=============================================
量化指标对比 (1次旋转循环):
最近邻插值 -> MSE: 166.7169
双线性插值 -> MSE: 60.3928
=============================================
已自动锁定最大差异区域: 左上角坐标 (64, 128), 窗口大小 76x76

# ==========================================
# Bonus: 5次连续旋转循环信息损失对比
# ==========================================
img_bonus_nn = img.copy()
img_bonus_bi = img.copy()

print("正在计算 5 次连续旋转循环，请稍候...")
for i in range(5):
    img_bonus_nn = matrix_rotate(img_bonus_nn, theta, method='nearest')
    img_bonus_nn = matrix_rotate(img_bonus_nn, -theta, method='nearest')
    
    img_bonus_bi = matrix_rotate(img_bonus_bi, theta, method='bilinear')
    img_bonus_bi = matrix_rotate(img_bonus_bi, -theta, method='bilinear')

mse_nn_5 = calculate_mse(img, img_bonus_nn)
mse_bi_5 = calculate_mse(img, img_bonus_bi)

print("-" * 40)
print(f"信息损失对比 (1次 vs 5次循环):")
print(f"最近邻 (Nearest): 1次 MSE = {mse_nn:.4f}  |  5次 MSE = {mse_nn_5:.4f}")
print(f"双线性 (Bilinear): 1次 MSE = {mse_bi:.4f}  |  5次 MSE = {mse_bi_5:.4f}")
print("-" * 40)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle('Bonus: 5 Continuous Rotation Cycles', fontsize=16)

axes[0].imshow(img, cmap='gray')
axes[0].set_title("Original Image")

axes[1].imshow(img_bonus_nn, cmap='gray')
axes[1].set_title(f"Nearest (5 Cycles)\nMSE: {mse_nn_5:.2f}")

axes[2].imshow(img_bonus_bi, cmap='gray')
axes[2].set_title(f"Bilinear (5 Cycles)\nMSE: {mse_bi_5:.2f}")

for ax in axes:
    ax.axis('off')

plt.tight_layout()
plt.show()

正在计算 5 次连续旋转循环，请稍候...
----------------------------------------
信息损失对比 (1次 vs 5次循环):
最近邻 (Nearest): 1次 MSE = 166.7169  |  5次 MSE = 166.7173
双线性 (Bilinear): 1次 MSE = 60.3928  |  5次 MSE = 446.6971
----------------------------------------

import numpy as np
import matplotlib.pyplot as plt

# 读取并归一化原始 CT 图像到 [0, 1] 区间
# 注意：文件路径已更新为 lab2_data 文件夹下
img_ct = np.load('lab2_data/lab2_CT.npy')
img_ct_norm = (img_ct - np.min(img_ct)) / (np.max(img_ct) - np.min(img_ct))

# 1. 探究对数变换参数 v 的影响
v_values = [1, 10, 100]  # v 越大，对数曲线在低灰度区的拉伸越剧烈
log_imgs = []
for v in v_values:
    # 计算归一化系数 c
    c_log = 1 / np.log(1 + v * np.max(img_ct_norm))
    # 矩阵运算实现参数化对数变换
    img_log = c_log * np.log(1 + v * img_ct_norm)
    log_imgs.append(img_log)

# 2. 探究幂律变换参数 gamma 的影响
gammas = [0.3, 0.6, 2.0, 4.0]
power_imgs = []
for gamma in gammas:
    c_power = 1 / np.power(np.max(img_ct_norm), gamma)
    img_power = c_power * np.power(img_ct_norm, gamma)
    power_imgs.append(img_power)

# 绘图展示所有参数对比
fig, axes = plt.subplots(2, 4, figsize=(20, 10))

# 第一行：原图与对数变换参数对比
axes[0, 0].imshow(img_ct_norm, cmap='gray')
axes[0, 0].set_title("Original CT Image")
axes[0, 1].imshow(log_imgs[0], cmap='gray')
axes[0, 1].set_title(f"Log Transform (v={v_values[0]})")
axes[0, 2].imshow(log_imgs[1], cmap='gray')
axes[0, 2].set_title(f"Log Transform (v={v_values[1]})")
axes[0, 3].imshow(log_imgs[2], cmap='gray')
axes[0, 3].set_title(f"Log Transform (v={v_values[2]})")

# 第二行：幂律变换参数对比
axes[1, 0].imshow(power_imgs[0], cmap='gray')
axes[1, 0].set_title(f"Power Transform (γ={gammas[0]})")
axes[1, 1].imshow(power_imgs[1], cmap='gray')
axes[1, 1].set_title(f"Power Transform (γ={gammas[1]})")
axes[1, 2].imshow(power_imgs[2], cmap='gray')
axes[1, 2].set_title(f"Power Transform (γ={gammas[2]})")
axes[1, 3].imshow(power_imgs[3], cmap='gray')
axes[1, 3].set_title(f"Power Transform (γ={gammas[3]})")

for ax in axes.flatten():
    ax.axis('off')
plt.tight_layout()
plt.show()

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

# ==========================================
# 1. 核心变换函数设计 (使用矩阵运算)
# ==========================================
def log_transform_v(img, v=1.0):
    """
    参数化对数变换: s = c * log(1 + v * r)
    参数 v 越大，对数曲线曲率越高，对暗部的拉伸越剧烈。
    """
    img_norm = (img - img.min()) / (img.max() - img.min())
    c = 1.0 / np.log1p(v) 
    img_log = c * np.log1p(v * img_norm)
    return np.clip(img_log, 0, 1)

def power_transform(img, gamma=1.0):
    """
    幂律(Gamma)变换: s = r^gamma
    """
    img_norm = (img - img.min()) / (img.max() - img.min())
    img_power = np.power(img_norm, gamma)
    return np.clip(img_power, 0, 1)

# ==========================================
# 2. 数据加载
# ==========================================
# 确保你的 lab2_CT.npy 文件存放在当前目录的 lab2_data 文件夹下
file_path = os.path.join('lab2_data', 'lab2_CT.npy')
img_ct = np.load(file_path)

# ==========================================
# 3. 特定组织的定向增强执行代码
# ==========================================

# 目标 A：肺部组织 (暗区拉伸) -> 使用高曲率对数变换
# v=50 能极大提升极暗区域的亮度，暴露出肺部微小纹理
img_lung_enhanced = log_transform_v(img_ct, v=50) 

# 目标 B：骨骼结构 (亮区拉伸) -> 使用高 Gamma 幂律变换
# gamma=3.0 的下凹曲线能压暗软组织背景，锐化高光骨骼区域的对比度
img_bone_enhanced = power_transform(img_ct, gamma=3.0) 

# ==========================================
# 4. 可视化结果展示
# ==========================================
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

axes[0].imshow(img_ct, cmap='gray')
axes[0].set_title('1. Original CT Image', fontsize=14)

axes[1].imshow(img_lung_enhanced, cmap='gray')
axes[1].set_title('2. Target: Lung Tissue\nLog Transform (v=50)', fontsize=14)

axes[2].imshow(img_bone_enhanced, cmap='gray')
axes[2].set_title('3. Target: Bone Structures\nPower Transform (γ=3.0)', fontsize=14)

for ax in axes:
    ax.axis('off')

plt.tight_layout()
plt.show()

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

# ==========================================
# 1. 加载数据与预处理
# ==========================================
file_path = os.path.join('lab2_data', 'lab2_filter.npy')
img = np.load(file_path)

if img.ndim == 3:
    img = img[:, :, 0]
img = (img - img.min()) / (img.max() - img.min())

# ==========================================
# 2. 算子定义
# ==========================================
laplacian_kernel = np.array([[-1, -1, -1],
                             [-1,  8, -1],
                             [-1, -1, -1]])

sobel_x_kernel = np.array([[-1,  0,  1],
                           [-2,  0,  2],
                           [-1,  0,  1]])
sobel_y_kernel = np.array([[-1, -2, -1],
                           [ 0,  0,  0],
                           [ 1,  2,  1]])

# ==========================================
# 3. 核心滤波与增强管线 (修复后的经典管线)
# ==========================================
# [Step A] Laplacian (提取细节，包含正负值和噪声)
laplacian_img = cv2.filter2D(img, -1, laplacian_kernel)

# [Step B] Simple Sharpening (初级锐化，仅做反面教材展示)
sharpened_img = np.clip(img + laplacian_img, 0, 1)

# [Step C] Sobel (提取真实边缘)
sobel_x = cv2.filter2D(img, -1, sobel_x_kernel)
sobel_y = cv2.filter2D(img, -1, sobel_y_kernel)
sobel_mag = np.clip(np.abs(sobel_x) + np.abs(sobel_y), 0, 1)

# [Step D] Smoothed Sobel (平滑掩膜)
box_kernel = np.ones((5, 5)) / 25
sobel_smoothed = cv2.filter2D(sobel_mag, -1, box_kernel)

# [Step E] 掩膜相乘 (应该用纯细节去乘掩膜，过滤掉背景细节噪声)
masked_laplacian = laplacian_img * sobel_smoothed

# [Step F] 最终融合 (将过滤后的干净细节加回原图)
enhanced_img = np.clip(img + masked_laplacian, 0, 1)

# [Step G] Gamma 对比度增强
gamma = 0.5
final_result = np.power(enhanced_img, gamma)

# ==========================================
# 4. 可视化分析展示
# ==========================================
def norm_for_display(image):
    # 将含有负数的矩阵居中归一化，0会变成0.5(灰色)
    return (image - image.min()) / (image.max() - image.min() + 1e-8)

fig, axes = plt.subplots(2, 4, figsize=(18, 10))
fig.suptitle('Project 3: Spatial Filtering Enhancement Pipeline', fontsize=18)

axes[0, 0].imshow(img, cmap='gray', vmin=0, vmax=1)
axes[0, 0].set_title("1. Original Image")

axes[0, 1].imshow(norm_for_display(laplacian_img), cmap='gray')
axes[0, 1].set_title("2. Laplacian (High Freq Details)")

axes[0, 2].imshow(sharpened_img, cmap='gray', vmin=0, vmax=1)
axes[0, 2].set_title("3. Simple Sharpening (Noisy)")

axes[0, 3].imshow(norm_for_display(sobel_mag), cmap='gray')
axes[0, 3].set_title("4. Sobel Magnitude (Edges)")

axes[1, 0].imshow(norm_for_display(sobel_smoothed), cmap='gray')
axes[1, 0].set_title("5. Smoothed Sobel (Mask)")

# 过滤后的纯净细节
axes[1, 1].imshow(norm_for_display(masked_laplacian), cmap='gray')
axes[1, 1].set_title("6. Masked Laplacian")

axes[1, 2].imshow(enhanced_img, cmap='gray', vmin=0, vmax=1)
axes[1, 2].set_title("7. Mask Added to Original")

axes[1, 3].imshow(final_result, cmap='gray', vmin=0, vmax=1)
axes[1, 3].set_title("8. Final: Gamma Transformed")

for ax in axes.flatten():
    ax.axis('off')

plt.tight_layout()
plt.show()

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

Project 1: Spatial Transform¶

实验要求¶

实验原理¶

实验结果与分析¶

Project1 Bonus¶

实验结果与分析¶

Project2 Intensity Transform¶

实验要求¶

实验原理¶

实验结果与分析¶

对数变换法和幂律变换方法的差异¶

对数变换参数 ($v$) 的影响¶

幂律变换参数 ($\gamma$) 的影响¶

最佳参数选择¶

特定临床目标的参数定向选择¶

肺部组织¶

骨骼结构¶

Project 3——Spatial Filter: High pass filter¶

实验内容¶

实验原理¶

8 邻域拉普拉斯算子 (8-neighbor Laplacian)¶

3x3 Sobel 算子¶

设计思路¶

算法设计流程图¶

实验结果与分析¶

图 1. Original Image（原始图像）¶

图 2. Laplacian (High Freq Details)¶

图 3. Simple Sharpening (Noisy)¶

图 4. Sobel Magnitude (Edges)¶

图 5. Smoothed Sobel (Mask)¶

图 6. Masked Laplacian¶

图 7. Mask Added to Original¶

图 8. Final: Gamma Transformed¶

实验心得¶

重采样的不可逆¶

理解临床导向的参数设计¶

领略增强图像中设计的智慧¶