4. 凸集投影超分辨率重建

任务说明

本案例实现凸集投影（Projection Onto Convex Sets, POCS）超分辨率重建算法。给定一张低分辨率输入图像，通过迭代交替投影到两个凸集约束上——数据一致性约束（保证重建结果下采样后与观测低分辨率图像一致）和平滑性约束（抑制噪声与振铃伪影），逐步恢复高分辨率细节。最终输出一张放大4倍的高分辨率图像，并给出可量化的收敛判定标准。

算法原理与实现

POCS的核心思想是将先验信息建模为希尔伯特空间中的闭凸集，通过交替投影逐步逼近所有约束集的交集。本案例定义两个凸集：C1（数据一致性）要求重建图像经下采样后与观测图像误差不超过阈值ε；C2（平滑性）限制局部像素变化幅度，抑制高频噪声。每次迭代中，先对当前估计执行C1投影（修正低频误差），再执行C2投影（平滑高频伪影），直至相邻两次迭代的均方误差变化小于阈值。

实现采用双线性插值作为降质模型（下采样），投影运算通过梯度下降完成。代码包含完整的初始化、迭代循环和收敛判断，不依赖任何外部图像库，仅使用Python标准库与NumPy。

代码实现

import numpy as np
from scipy.ndimage import zoom, gaussian_filter

def downsample(img, factor=4):
    """双线性下采样：将高分辨率图像降为低分辨率"""
    h, w = img.shape
    return zoom(img, (1/factor, 1/factor), order=1)

def upsample(img, factor=4):
    """双线性上采样：将低分辨率图像放大"""
    return zoom(img, (factor, factor), order=1)

def data_consistency_projection(x, y, factor=4, eps=1.0, lr=0.5):
    """C1投影：修正下采样后与观测图像的误差"""
    x_down = downsample(x, factor)
    diff = x_down - y
    # 仅对误差超过阈值的像素进行修正
    mask = np.abs(diff) > eps
    if np.any(mask):
        # 将误差上采样回高分辨率空间并更新
        correction = upsample(diff * mask, factor)
        x = x - lr * correction
    return x

def smoothness_projection(x, mu=0.1):
    """C2投影：高斯平滑抑制高频伪影"""
    # 计算局部方差，在平坦区域加强平滑
    local_mean = gaussian_filter(x, sigma=1.0)
    local_var = gaussian_filter((x - local_mean)**2, sigma=1.0)
    # 自适应平滑强度：方差越小平滑越强
    alpha = mu * np.exp(-local_var / (np.mean(local_var) + 1e-8))
    return x - alpha * (x - local_mean)

def pocs_super_resolution(y, factor=4, max_iter=50, tol=1e-4):
    """
    凸集投影超分辨率重建
    参数:
        y: 低分辨率输入图像 (M×N, float64, 范围[0,255])
        factor: 放大倍数
        max_iter: 最大迭代次数
        tol: 收敛阈值
    返回:
        x: 重建的高分辨率图像 (factor*M × factor*N)
        history: 每次迭代的MSE变化记录
    """
    # 初始化：双线性插值放大
    x = upsample(y, factor).astype(np.float64)
    history = []
    
    for k in range(max_iter):
        x_prev = x.copy()
        
        # C1投影：数据一致性
        x = data_consistency_projection(x, y, factor, eps=0.5, lr=0.3)
        
        # C2投影：平滑性约束
        x = smoothness_projection(x, mu=0.15)
        
        # 计算收敛指标
        mse_change = np.mean((x - x_prev)**2)
        history.append(mse_change)
        
        if mse_change < tol:
            break
    
    return x, history

# ========== 示例运行 ==========
if __name__ == "__main__":
    # 生成测试图像：8×8低分辨率（模拟真实场景）
    np.random.seed(42)
    lr_img = np.random.randint(50, 200, size=(8, 8)).astype(np.float64)
    
    # 执行POCS重建
    hr_img, hist = pocs_super_resolution(lr_img, factor=4, max_iter=30)
    
    print(f"输入尺寸: {lr_img.shape} → 输出尺寸: {hr_img.shape}")
    print(f"迭代次数: {len(hist)}")
    print(f"最终收敛值: {hist[-1]:.6f}")
    print(f"重建像素范围: [{hr_img.min():.1f}, {hr_img.max():.1f}]")
    
    # 验证数据一致性：重建结果下采样后应与输入接近
    verify = downsample(hr_img, factor=4)
    consistency_error = np.mean(np.abs(verify - lr_img))
    print(f"数据一致性平均误差: {consistency_error:.3f} (应小于1.0)")

运行与验证

依赖安装：需要NumPy和SciPy，执行pip install numpy scipy。

运行方式：将代码保存为pocs_sr.py，在终端执行python pocs_sr.py。

验证方法：程序输出应满足以下判定标准：

• 输出图像尺寸为输入尺寸的4倍（8×8 → 32×32）；
• 迭代次数在1～30之间，且最终收敛值小于1e-4；
• 数据一致性误差小于1.0，表明重建结果与观测数据吻合；
• 重建像素值保持在合理灰度范围[0,255]内。

若需处理真实图像，可将lr_img替换为plt.imread()读取的灰度图（需先转为float64并归一化到[0,255]）。