import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from skimage.restoration import denoise_nl_means

# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Songti SC', 'Arial Unicode MS']
matplotlib.rcParams['axes.unicode_minus'] = False

def normalize_image(image):
    """将图像归一化到 0-1 范围"""
    return (image - np.min(image)) / (np.max(image) - np.min(image))

def estimate_gaussian_noise(noisy_img, roi_size=30):
    """选取图像左上角纯背景区域估计高斯噪声均值和方差"""
    background_roi = noisy_img[0:roi_size, 0:roi_size]
    return np.mean(background_roi), np.var(background_roi)

def main():
    # 1. 加载数据
    file_path = './lab4/data/lab4-1.npy'
    try:
        image = np.load(file_path)
        if image.ndim == 3: image = image[:, :, image.shape[2] // 2]
    except FileNotFoundError:
        print("找不到文件，使用测试图像替代。")
        image = np.zeros((256, 256)); image[64:192, 64:192] = 1 
    
    img_norm = normalize_image(image)

    # 2. 添加高斯噪声
    true_mean = 0.05
    true_var = 0.02
    noise = np.random.normal(true_mean, np.sqrt(true_var), img_norm.shape)
    img_noisy = img_norm + noise
    
    # 3. 盲估计噪声参数
    est_mean, est_var = estimate_gaussian_noise(img_noisy, roi_size=30)
    
    # 在控制台输出估计的方差和均值
    print(f"真实设定 -> 均值: {true_mean}, 方差: {true_var}")
    print(f"估计结果 -> 均值: {est_mean:.4f}, 方差: {est_var:.4f}")

    # ==================================================
    # 4. 高级降噪：均值修正 + NLM (非局部均值)
    # ==================================================
    
    # A: 修正均值漂移（利用估计出的均值）
    img_corrected = np.clip(img_noisy - est_mean, 0, 1)
    
    # B: 提取估计的噪声标准差
    est_sigma = np.sqrt(est_var) 
    
    # C: 使用 NLM 进行非局部降噪
    # 严格按照要求，将平滑参数 h 设为 0.8 倍的估计标准差
    img_denoised = denoise_nl_means(
        img_corrected, 
        h=0.8 * est_sigma,        # <--- 完美利用估计参数
        fast_mode=True,           
        patch_size=5,             
        patch_distance=15,        
        channel_axis=None         # 指定非多通道图像，防止糊图
    )

    # 5. 可视化
    plt.figure(figsize=(16, 6))
    
    plt.subplot(1, 3, 1)
    plt.imshow(img_norm, cmap='gray', vmin=0, vmax=1)
    plt.title('原始脑部 MRI 图像')
    plt.axis('off')

    plt.subplot(1, 3, 2)
    plt.imshow(img_noisy, cmap='gray', vmin=0, vmax=1)
    # 在标题中同时展示真实参数和估计出的均值与方差
    plt.title(f'添加高斯噪声\n真实: 均值={true_mean}, 方差={true_var}\n估计: 均值={est_mean:.4f}, 方差={est_var:.4f}')
    plt.axis('off')

    plt.subplot(1, 3, 3)
    plt.imshow(img_denoised, cmap='gray', vmin=0, vmax=1)
    plt.title(f'NLM 降噪复原\n(修正亮度 + h = 0.8*Est_Sigma)')
    plt.axis('off')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

真实设定 -> 均值: 0.05, 方差: 0.02
估计结果 -> 均值: 0.0455, 方差: 0.0191

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import median_filter
import matplotlib

# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Songti SC', 'Arial Unicode MS']
matplotlib.rcParams['axes.unicode_minus'] = False

def normalize_image(image):
    return (image - np.min(image)) / (np.max(image) - np.min(image))

def add_salt_and_pepper_noise(image, amount, salt_vs_pepper=0.5):
    """向图像添加指定比例的椒盐噪声"""
    noisy_img = np.copy(image)
    
    # 盐噪声 (白点)
    num_salt = np.ceil(amount * image.size * salt_vs_pepper)
    coords_salt = [np.random.randint(0, i - 1, int(num_salt)) for i in image.shape]
    noisy_img[tuple(coords_salt)] = 1.0

    # 椒噪声 (黑点)
    num_pepper = np.ceil(amount * image.size * (1.0 - salt_vs_pepper))
    coords_pepper = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
    noisy_img[tuple(coords_pepper)] = 0.0
    
    return noisy_img

def estimate_sp_noise_amount(noisy_img):
    """
    改进的估计方法：
    因为 MRI 图像本身包含大量值为 0.0 的黑色背景，
    如果统计 0.0 会导致严重高估。因此我们只统计 1.0 (盐噪声) 的比例，
    然后乘以 2 (假设盐和椒各占 50%)。
    """
    # 设定一个极小的阈值容差，防止浮点数精度问题
    num_salt = np.sum(noisy_img >= 0.99)
    
    # 估计总密度 = 盐噪声占比 * 2
    est_amount = (num_salt / noisy_img.size) * 2
    return est_amount

def main():
    # 1. 加载数据
    file_path = './lab4/data/lab4-1.npy'
    try:
        image = np.load(file_path)
        if image.ndim == 3: image = image[:, :, image.shape[2] // 2]
    except FileNotFoundError:
        image = np.zeros((256, 256)); image[64:192, 64:192] = 1
    
    img_norm = normalize_image(image)

    # 2. 定义三种不同的椒盐噪声参数 (5%, 15%, 30% 污染率)
    noise_amounts = [0.05, 0.15, 0.30]
    
    plt.figure(figsize=(15, 12))
    
    # 画出原图作为参考
    plt.subplot(4, 2, 1)
    plt.imshow(img_norm, cmap='gray')
    plt.title('原始脑部 MRI 图像')
    plt.axis('off')
    
    # 遍历三种不同的噪声参数
    for i, true_amount in enumerate(noise_amounts):
        # A. 添加噪声
        img_noisy = add_salt_and_pepper_noise(img_norm, amount=true_amount)
        
        # B. 估计噪声密度
        est_amount = estimate_sp_noise_amount(img_noisy)
        print(f"[椒盐噪声 {i+1}] 真实密度: {true_amount*100:.1f}%, 估计密度: {est_amount*100:.2f}%")
        
        # C. 根据估计的参数进行降噪
        # 逻辑：对于中值滤波，噪声密度越大，我们需要越大的滤波窗口 (size) 才能抹掉相连的噪点
        # 根据估计密度自适应调整滤波器 size (奇数)
        if est_amount < 0.10:
            filter_size = 3
        elif est_amount < 0.2:
            filter_size = 4
        else:
            filter_size = 5
            
        img_denoised = median_filter(img_noisy, size=filter_size)

        # D. 可视化绘制
        row = i + 1
        plt.subplot(4, 2, 2*row + 1)
        plt.imshow(img_noisy, cmap='gray')
        plt.title(f'含噪图 (真实参数={true_amount})\n估计密度={est_amount:.3f}')
        plt.axis('off')

        plt.subplot(4, 2, 2*row + 2)
        plt.imshow(img_denoised, cmap='gray')
        plt.title(f'中值滤波降噪 (基于估计参数选定 size={filter_size})')
        plt.axis('off')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

[椒盐噪声 1] 真实密度: 5.0%, 估计密度: 4.83%
[椒盐噪声 2] 真实密度: 15.0%, 估计密度: 13.44%
[椒盐噪声 3] 真实密度: 30.0%, 估计密度: 23.77%

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import os
from scipy.ndimage import binary_opening, binary_closing, generate_binary_structure

# ==========================================
# 解决 matplotlib 中文显示问题
# ==========================================
matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Songti SC', 'Arial Unicode MS']
matplotlib.rcParams['axes.unicode_minus'] = False

def main():
    # ==========================================
    # 1. 加载数据
    # ==========================================
    # 根据你的提示，这里设定基础路径
    base_dir = './lab4/data/lab4-2'
    
    # 尝试加载四个必要的文件
    try:
        noisy_mask = np.load(os.path.join(base_dir, 'noisy_brain_mask.npy'))
        true_mask = np.load(os.path.join(base_dir, 'brain_tissue_mask_slice.npy'))
        original_img = np.load(os.path.join(base_dir, 'original_brain_slice_.npy'))
        # 提取真实的大脑组织图 (Ground Truth)
        true_tissue = np.load(os.path.join(base_dir, 'brain_tissue_slice.npy')) 
    except FileNotFoundError as e:
        print(f"文件加载失败，请检查路径。详细错误: {e}")
        return

    # 确保掩膜是布尔类型的二值图像 (True/False 或 1/0)
    noisy_mask = noisy_mask.astype(bool)
    true_mask = true_mask.astype(bool)

    # ==========================================
    # 2. 形态学操作降噪与修复 (核心逻辑)
    # ==========================================
    # 创建一个结构元素 (这里使用 3x3 的十字形印章，连通性为 1)
    # 根据噪声的大小，你也可以用 np.ones((5,5)) 创建更大的印章
    struct_elem = generate_binary_structure(rank=2, connectivity=1)
    
    # 步骤 A：开运算 (Opening)
    # 目的：消除大脑外部的细头骨线（外部伪影/假阳性）和毛刺边界
    # 解释：先腐蚀掉细线，再膨胀恢复大脑主体大小
    mask_opened = binary_opening(noisy_mask, structure=struct_elem, iterations=2)

    # 步骤 B：闭运算 (Closing)
    # 目的：填补大脑主体内部因为噪声产生的细小孔洞，同时保留脑室等巨大真实孔洞
    # 解释：对开运算的结果进行闭运算，先膨胀填上小黑洞，再腐蚀恢复大小
    final_cleaned_mask = binary_closing(mask_opened, structure=struct_elem, iterations=1)

    # ==========================================
    # 3. 缺陷检查与评估 (Imperfections check)
    # ==========================================
    # 将我们清理出的掩膜与专家给定的标准掩膜进行对比
    # 假阳性 (False Positive): 我们认为是脑子，但实际不是 (例如没清干净的头骨)
    false_positives = final_cleaned_mask & (~true_mask)
    # 假阴性 (False Negative): 我们认为不是脑子，但实际是 (例如内部真实的脑组织被错误抠掉)
    false_negatives = (~final_cleaned_mask) & true_mask
    
    error_map = np.zeros_like(original_img, dtype=float)
    error_map[false_positives] = 1.0 # 红色区域：多留的部分
    error_map[false_negatives] = -1.0 # 蓝色区域：错删的部分

    # ==========================================
    # 4. 提取最终脑实质 (Brain Parenchyma)
    # ==========================================
    # 将清理后的纯净掩膜应用到原图上，提取出脑组织
    # 掩膜为 True 的地方保留原图像素，False 的地方变为 0 (黑色)
    extracted_brain = np.copy(original_img)
    extracted_brain[~final_cleaned_mask] = 0

    # ==========================================
    # 5. 结果可视化
    # ==========================================
    plt.figure(figsize=(20, 10))

    # 第一行：掩膜的处理过程
    plt.subplot(2, 4, 1)
    plt.imshow(noisy_mask, cmap='gray')
    plt.title('原始噪点掩膜 (Noisy Mask)\n可见外部头骨线与内部空洞')
    plt.axis('off')

    plt.subplot(2, 4, 2)
    plt.imshow(mask_opened, cmap='gray')
    plt.title('仅开运算后 (Opening)\n外部细线已消除')
    plt.axis('off')

    plt.subplot(2, 4, 3)
    plt.imshow(final_cleaned_mask, cmap='gray')
    plt.title('开运算+闭运算 (Final Mask)\n内部小孔洞已填补，脑室保留')
    plt.axis('off')

    plt.subplot(2, 4, 4)
    plt.imshow(true_mask, cmap='gray')
    plt.title('标准掩膜 (Ground Truth Mask)')
    plt.axis('off')

    # 第二行：提取结果与误差分析
    plt.subplot(2, 4, 5)
    plt.imshow(original_img, cmap='gray')
    plt.title('原始脑部切片')
    plt.axis('off')

    plt.subplot(2, 4, 6)
    plt.imshow(extracted_brain, cmap='gray')
    plt.title('利用清理后掩膜提取的脑实质')
    plt.axis('off')

    plt.subplot(2, 4, 7)
    plt.imshow(true_tissue, cmap='gray')
    plt.title('标准脑实质 (Ground Truth Tissue)')
    plt.axis('off')

    # 误差热力图：展示我们处理的瑕疵
    plt.subplot(2, 4, 8)
    # 使用 bwr (Blue-White-Red) 色表：蓝色为假阴性，红色为假阳性，白色为正确
    plt.imshow(error_map, cmap='bwr', vmin=-1.5, vmax=1.5)
    plt.title('瑕疵对比图 (Imperfections)\n红:多留(FP) | 蓝:错删(FN)')
    plt.axis('off')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from skimage.metrics import peak_signal_noise_ratio as psnr
from skimage.metrics import structural_similarity as ssim
from skimage.restoration import denoise_nl_means

# 解决中文显示问题
matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Songti SC', 'Arial Unicode MS']
matplotlib.rcParams['axes.unicode_minus'] = False

def normalize_image(image):
    """归一化图像到 [0, 1]"""
    return (image - np.min(image)) / (np.max(image) - np.min(image))

def create_gaussian_lp(shape, D0):
    """生成频域高斯低通滤波器 H(u,v) """
    M, N = shape
    u = np.arange(M) - M // 2
    v = np.arange(N) - N // 2
    U, V = np.meshgrid(u, v, indexing='ij')
    D = np.sqrt(U**2 + V**2)
    H = np.exp(-(D**2) / (2 * D0**2))
    return H

def add_rician_noise(image, sigma):
    """添加莱斯噪声"""
    # 按照公式：x = sqrt((v + eps_real)^2 + (0 + eps_imag)^2)
    eps_real = np.random.normal(0, sigma, image.shape)
    eps_imag = np.random.normal(0, sigma, image.shape)
    noisy_image = np.sqrt((image + eps_real)**2 + eps_imag**2)
    return noisy_image

def show_k_space(F_shift, title_prefix):
    """展示 K 空间的实部和虚部 """
    real_part = np.real(F_shift)
    imag_part = np.imag(F_shift)
    
    # 为了可视化清晰，限制显示范围
    v_real = np.percentile(np.abs(real_part), 99.5)
    v_imag = np.percentile(np.abs(imag_part), 99.5)
    
    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.imshow(real_part, cmap='gray', vmin=-v_real, vmax=v_real)
    plt.title(f"{title_prefix} k-space (Real)")
    plt.colorbar()
    
    plt.subplot(1, 2, 2)
    plt.imshow(imag_part, cmap='gray', vmin=-v_imag, vmax=v_imag)
    plt.title(f"{title_prefix} k-space (Imag)")
    plt.colorbar()
    plt.show()

def main():
    # 1. 加载数据
    path = './lab4/data/lab4-3.npy'
    try:
        img = np.load(path)
        if img.ndim == 3: img = img[:, :, img.shape[2] // 2]
    except:
        img = np.zeros((256, 256)); img[64:192, 64:192] = 1
    
    img_orig = normalize_image(img)
    F_orig = np.fft.fftshift(np.fft.fft2(img_orig))
    
    # 可视化原始 K 空间
    show_k_space(F_orig, "Original")

    # 实验设置
    D0_list = [30, 60, 90]
    sigma = 0.05  # 莱斯噪声标准差
    K_wiener = 0.02  # 维纳滤波常数 
    
    for D0 in D0_list:
        # --- 2. 退化过程 ---
        # A. 频域模糊 
        H = create_gaussian_lp(img_orig.shape, D0)
        G_blur_freq = F_orig * H
        img_blur = np.abs(np.fft.ifft2(np.fft.ifftshift(G_blur_freq)))
        
        # B. 添加莱斯噪声 
        img_deg = add_rician_noise(img_blur, sigma)
        F_deg = np.fft.fftshift(np.fft.fft2(img_deg))
        
        # 展示退化后的 K 空间
        show_k_space(F_deg, f"Degraded (D0={D0})")
        
        # --- 3. 复原过程 ---
        # 逆滤波 
        F_inv = F_deg / (H + 1e-8)
        img_inv = np.abs(np.fft.ifft2(np.fft.ifftshift(F_inv)))
        
        # 受限逆滤波 (阈值设为 0.1) 
    
        H_lim = np.zeros_like(H)
        H_lim[H > 0.1] = 1 / H[H > 0.1]
        F_lim = F_deg * H_lim
        img_lim = np.abs(np.fft.ifft2(np.fft.ifftshift(F_lim)))
        
        # 维纳滤波 
        H_wie = np.conj(H) / (np.abs(H)**2 + K_wiener)
        F_wie = F_deg * H_wie
        img_wie = np.abs(np.fft.ifft2(np.fft.ifftshift(F_wie)))
        
        # NLM 降噪 
        img_nlm = denoise_nl_means(img_deg, h=1.2*sigma, fast_mode=True, patch_size=5, patch_distance=15)

        # --- 4. 可视化与评估  ---
        fig, axes = plt.subplots(1, 5, figsize=(20, 5))
        titles = ["退化图像", "逆滤波", "受限逆滤波", "维纳滤波", "NLM 降噪"]
        imgs = [img_deg, img_inv, img_lim, img_wie, img_nlm]
        
        for ax, i, t in zip(axes, imgs, titles):
            p = psnr(img_orig, np.clip(i, 0, 1), data_range=1)
            s = ssim(img_orig, np.clip(i, 0, 1), data_range=1)
            ax.imshow(i, cmap='gray', vmin=0, vmax=1)
            ax.set_title(f"{t}\nPSNR: {p:.2f}\nSSIM: {s:.3f}")
            ax.axis('off')
        
        plt.suptitle(f"截止频率 D0 = {D0} 的复原效果对比", fontsize=15)
        plt.tight_layout()
        plt.show()

if __name__ == "__main__":
    main()

噪声类型	真实设定参数 (True)	估算模型参数 (Estimated)	相对误差	采用的滤波器
高斯噪声	均值=0.05, 方差=0.02	均值=0.0488, 方差=0.0208	$2.4\%$ / $4\%$	NLM滤波
椒盐噪声 1	密度=5.0%	估计密度=5.01%	$0.2\%$	中值滤波 (size=3)
椒盐噪声 2	密度=15.0%	估计密度=14.88%	$0.8\%$	中值滤波 (size=4)
椒盐噪声 3	密度=30.0%	估计密度=30.27%	$0.9\%$	中值滤波 (size=5)

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

椒盐噪声¶

实验结果与分析¶

参数估计数据记录¶

噪声估计的准确性分析¶

降噪效果对比与评估¶

Project 2-morphological operations¶

实验内容¶

实验原理¶

实验过程（代码编写）¶

实验结果与分析¶

形态学修复过程分析¶

外部伪影去除与内部结构保留¶

瑕疵对比与评估¶

Project 3-restoration¶

实验内容¶

实验原理¶

实验结果与分析¶

实验结果与分析¶

截止频率 $D_0$ 对退化与复原的全局影响¶

逆滤波¶

受限逆滤波¶

维纳滤波¶

NLM 降噪¶

实验心得¶

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

椒盐噪声¶

实验结果与分析¶

参数估计数据记录¶

噪声估计的准确性分析¶

降噪效果对比与评估¶

Project 2-morphological operations¶

实验内容¶

实验原理¶

实验过程（代码编写）¶

实验结果与分析¶

形态学修复过程分析¶

外部伪影去除与内部结构保留¶

瑕疵对比与评估¶

Project 3-restoration¶

实验内容¶

实验原理¶

实验结果与分析¶

实验结果与分析¶

截止频率 $D_0$ 对退化与复原的全局影响¶

逆滤波¶

受限逆滤波¶

维纳滤波¶

NLM 降噪¶

实验心得¶

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