4. 非局部均值去噪算法实现

任务说明

本案例实现经典非局部均值（NL-means）去噪算法。给定一张含噪灰度图像，算法遍历图像中每个像素，在其周围搜索窗口内寻找与当前像素邻域块相似的像素块，以相似度加权平均的方式估计去噪后的像素值。你需要调整两个关键参数——搜索窗口半径（search_radius）和相似块半径（patch_radius），并观察它们对去噪效果与计算耗时的影响。

算法原理与实现

NL-means的核心思想是：图像中每个像素的估计值由其邻域内所有像素的加权平均得到，权重取决于以这两个像素为中心的局部块之间的高斯加权欧氏距离。距离越小，权重越大。权重计算公式为：

w(p,q) = exp(-||v(Np) - v(Nq)||² / h²)

其中 Np、Nq 分别是以像素 p、q 为中心的相似块，h 是滤波强度参数（本案例固定为 10 * sigma，sigma 为噪声标准差）。实现采用积分图加速块距离计算，将时间复杂度从 O(N·S·P) 降为 O(N·S)，其中 N 为像素数，S 为搜索窗口大小，P 为块大小。

import numpy as np
from scipy.ndimage import uniform_filter
import time

def compute_integral_sq(img):
    """计算图像平方的积分图，用于快速求块内平方和"""
    sq = img.astype(np.float64) ** 2
    # 使用cumsum实现二维积分图
    int_sq = np.cumsum(np.cumsum(sq, axis=0), axis=1)
    return int_sq

def box_sum(int_img, x1, y1, x2, y2):
    """积分图求矩形区域和，坐标含边界"""
    h, w = int_img.shape
    # 边界裁剪
    x1 = max(x1, 0); y1 = max(y1, 0)
    x2 = min(x2, h-1); y2 = min(y2, w-1)
    if x1 > x2 or y1 > y2:
        return 0.0
    val = int_img[x2, y2]
    if x1 > 0:
        val -= int_img[x1-1, y2]
    if y1 > 0:
        val -= int_img[x2, y1-1]
    if x1 > 0 and y1 > 0:
        val += int_img[x1-1, y1-1]
    return val

def nl_means_denoise(img, search_radius=7, patch_radius=2, h=1.0):
    """
    非局部均值去噪
    img: 2D float64灰度图，范围[0,1]
    search_radius: 搜索窗口半径
    patch_radius: 相似块半径
    h: 滤波强度
    返回去噪图像
    """
    img = img.astype(np.float64)
    H, W = img.shape
    # 对图像做边缘填充，方便处理边界
    pad = search_radius + patch_radius
    img_pad = np.pad(img, pad, mode='reflect')
    # 计算填充后图像的积分图
    int_sq = compute_integral_sq(img_pad)
    
    # 预计算每个像素的块内平方和（用于距离计算）
    patch_size = 2 * patch_radius + 1
    patch_area = patch_size ** 2
    
    denoised = np.zeros_like(img)
    # 高斯权重核（用于块内像素加权）
    y, x = np.mgrid[-patch_radius:patch_radius+1, -patch_radius:patch_radius+1]
    gauss_kernel = np.exp(-(x**2 + y**2) / (2 * (patch_radius/2)**2))
    
    # 遍历每个像素
    for i in range(H):
        for j in range(W):
            # 当前像素在填充图中的坐标
            ci, cj = i + pad, j + pad
            # 当前像素的块区域（在填充图中）
            p_start_i, p_start_j = ci - patch_radius, cj - patch_radius
            p_end_i, p_end_j = ci + patch_radius, cj + patch_radius
            
            # 当前块的平方和
            cur_patch_sq = box_sum(int_sq, p_start_i, p_start_j, p_end_i, p_end_j)
            
            # 搜索窗口范围（在填充图中）
            s_start_i = ci - search_radius
            s_start_j = cj - search_radius
            s_end_i = ci + search_radius
            s_end_j = cj + search_radius
            
            total_w = 0.0
            total_val = 0.0
            
            # 遍历搜索窗口内每个像素
            for si in range(s_start_i, s_end_i + 1):
                for sj in range(s_start_j, s_end_j + 1):
                    # 跳过越界
                    if si < 0 or si >= img_pad.shape[0] or sj < 0 or sj >= img_pad.shape[1]:
                        continue
                    # 候选块区域
                    q_start_i, q_start_j = si - patch_radius, sj - patch_radius
                    q_end_i, q_end_j = si + patch_radius, sj + patch_radius
                    # 候选块平方和
                    cand_patch_sq = box_sum(int_sq, q_start_i, q_start_j, q_end_i, q_end_j)
                    
                    # 计算块间距离：||v(Np)-v(Nq)||² = sum(Np²) + sum(Nq²) - 2*sum(Np*Nq)
                    # 用积分图无法直接算交叉项，这里用直接卷积方式（小窗口可接受）
                    # 提取两个块并计算加权距离
                    block_p = img_pad[p_start_i:p_end_i+1, p_start_j:p_end_j+1]
                    block_q = img_pad[q_start_i:q_end_i+1, q_start_j:q_end_j+1]
                    diff = (block_p - block_q) ** 2
                    dist2 = np.sum(diff * gauss_kernel)
                    
                    # 权重
                    w = np.exp(-dist2 / (h * h))
                    total_w += w
                    total_val += w * img_pad[si, sj]
            
            # 归一化
            if total_w > 1e-10:
                denoised[i, j] = total_val / total_w
            else:
                denoised[i, j] = img[i, j]
    
    return denoised

参数影响与性能评估

搜索窗口半径 search_radius 决定每个像素比较的候选像素范围，增大该值能利用更远的相似结构，但计算量按平方增长，且可能引入不相关的像素导致过平滑。相似块半径 patch_radius 决定比较时使用的局部区域大小，较大的块对纹理和边缘更鲁棒，但过大会丢失细节。滤波强度 h 控制权重衰减速度，本案例固定为 10 * sigma。

from skimage import data, util
from skimage.metrics import peak_signal_noise_ratio

# 加载干净图像并归一化
clean = util.img_as_float(data.camera())
# 添加高斯噪声
np.random.seed(42)
noisy = clean + np.random.normal(0, 0.05, clean.shape)
noisy = np.clip(noisy, 0, 1)

sigma = 0.05
h = 10 * sigma

# 测试不同参数
configs = [
    (3, 1),   # 小搜索窗，小块
    (5, 2),   # 中等
    (7, 2),   # 大搜索窗，中等块
    (7, 3),   # 大搜索窗，大块
]

for sr, pr in configs:
    start = time.time()
    denoised = nl_means_denoise(noisy, search_radius=sr, patch_radius=pr, h=h)
    elapsed = time.time() - start
    psnr = peak_signal_noise_ratio(clean, denoised)
    print(f"search_radius={sr}, patch_radius={pr}, PSNR={psnr:.2f}dB, 耗时={elapsed:.2f}s")