全球

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

# -------------------------- 核心工具函数 --------------------------
def gradient(x):
    """计算图像的水平和垂直梯度（前向差分），对应论文中的∇x算子"""
    dx = np.zeros_like(x)
    dy = np.zeros_like(x)
    dx[:, :-1] = x[:, 1:] - x[:, :-1]  # 水平梯度
    dy[:-1, :] = x[1:, :] - x[:-1, :]  # 垂直梯度
    return dx, dy

def divergence(dx, dy):
    """计算梯度的散度（梯度的转置运算），用于求解x子问题"""
    div = np.zeros_like(dx)
    div[:, 1:] += dx[:, 1:] - dx[:, :-1]  # 水平散度
    div[1:, :] += dy[1:, :] - dy[:-1, :]  # 垂直散度
    return div

def soft_threshold(z, tau):
    """软阈值算子，对应论文式(4-5)，用于求解d子问题"""
    return np.sign(z) * np.maximum(np.abs(z) - tau, 0)

# -------------------------- ADMM-TV凸优化去噪核心算法 --------------------------
def admm_tv_denoising(y, lam=0.15, rho=1.0, max_iter=200, tol=1e-4):
    """
    论文提出的"最小二乘保真项+TV正则项+ADMM求解"船舶图像去噪算法
    参数:
        y: 输入含噪图像（归一化到[0,1]）
        lam: 正则化系数λ，论文默认0.15（适配中等盐雾噪声）
        rho: ADMM惩罚参数ρ，论文默认1.0
        max_iter: 最大迭代次数，论文默认200
        tol: 收敛阈值，论文默认1e-4
    返回:
        x: 去噪后的图像
    """
    # 初始化参数（对应论文4.4.1节）
    x = y.copy()
    dx, dy = gradient(x)
    ux = np.zeros_like(dx)  # 拉格朗日乘子μ的水平分量
    uy = np.zeros_like(dy)  # 拉格朗日乘子μ的垂直分量
    
    # 预计算傅里叶变换核（FFT加速求解x子问题，大幅提升运算速度）
    H, W = y.shape
    I = np.fft.fft2(np.eye(H, W))
    kernel = 1 + rho * (np.abs(np.fft.fft2(np.array([[1, -1]]), (H, W)))**2 + 
                        np.abs(np.fft.fft2(np.array([[1], [-1]]), (H, W)))**2)
    
    for k in range(max_iter):
        # 步骤1：更新图像变量x（对应论文4.3.1子问题1，FFT快速求解闭式解）
        rhs = y + rho * divergence(dx - ux/rho, dy - uy/rho)
        x = np.fft.ifft2(np.fft.fft2(rhs) / kernel).real
        
        # 步骤2：更新辅助变量d（对应论文4.3.2子问题2，软阈值求解）
        grad_x, grad_y = gradient(x)
        dx = soft_threshold(grad_x + ux/rho, lam/rho)
        dy = soft_threshold(grad_y + uy/rho, lam/rho)
        
        # 步骤3：更新拉格朗日乘子μ（对应论文4.3.3子问题3）
        ux += rho * (grad_x - dx)
        uy += rho * (grad_y - dy)
        
        # 计算迭代残差，判断收敛（对应论文4.4.3终止条件）
        res = np.linalg.norm(np.concatenate([grad_x - dx, grad_y - dy]))
        if res < tol:
            print(f"算法在第{k+1}次迭代收敛")
            break
    
    return np.clip(x, 0, 1)  # 限制像素值在[0,1]区间

# -------------------------- 实验辅助函数 --------------------------
def add_salt_fog_noise(img, salt_density=0.02, gauss_var=0.01):
    """模拟海上盐雾+高斯混合噪声（对应论文5.2节退化样本制备）"""
    noisy = img.copy()
    # 盐雾噪声（模拟离散椒盐噪声）
    salt = np.random.choice([0, 1, 2], size=img.shape, p=[salt_density/2, salt_density/2, 1-salt_density])
    noisy[salt == 0] = 0
    noisy[salt == 1] = 1
    # 高斯噪声（模拟海面散射噪声）
    gauss = np.random.normal(0, np.sqrt(gauss_var), img.shape)
    noisy += gauss
    return np.clip(noisy, 0, 1)

def evaluate_image(clean, denoised):
    """计算PSNR和SSIM定量指标（对应论文5.3节）"""
    psnr_val = psnr(clean, denoised, data_range=1.0)
    ssim_val = ssim(clean, denoised, data_range=1.0)
    return psnr_val, ssim_val

# -------------------------- 主实验流程 --------------------------
if __name__ == "__main__":
    # 1. 读取并预处理图像（对应论文3.5节步骤一）
    img_path = "ship_test.jpg"  # 替换为你的船舶图像路径
    clean_img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    clean_img = cv2.resize(clean_img, (512, 512))  # 统一为512x512，与论文一致
    clean_img = clean_img.astype(np.float32) / 255.0  # 归一化到[0,1]
    
    # 2. 添加模拟盐雾噪声
    noisy_img = add_salt_fog_noise(clean_img, salt_density=0.02, gauss_var=0.01)
    
    # 3. 运行三种算法（对应论文5.4节对比方案）
    print("正在运行均值滤波...")
    mean_filtered = cv2.blur(noisy_img, (3, 3))
    
    print("正在运行中值滤波...")
    median_filtered = cv2.medianBlur((noisy_img * 255).astype(np.uint8), 3) / 255.0
    
    print("正在运行ADMM-TV凸优化去噪...")
    tv_denoised = admm_tv_denoising(noisy_img, lam=0.15, rho=1.0, max_iter=200, tol=1e-4)
    
    # 4. 计算定量指标
    metrics = {
        "均值滤波": evaluate_image(clean_img, mean_filtered),
        "中值滤波": evaluate_image(clean_img, median_filtered),
        "本文凸优化算法": evaluate_image(clean_img, tv_denoised)
    }
    
    # 5. 打印指标结果（对应论文表5-1）
    print("\n" + "="*50)
    print("算法性能对比（平均指标）")
    print("="*50)
    print(f"{'算法':<15} {'PSNR(dB)':<10} {'SSIM':<10}")
    for alg, (p, s) in metrics.items():
        print(f"{alg:<15} {p:<10.2f} {s:<10.3f}")
    print("="*50)
    
    # 6. 结果可视化
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    axes = axes.flatten()
    
    axes[0].imshow(clean_img, cmap='gray')
    axes[0].set_title("原始清晰图像")
    axes[0].axis('off')
    
    axes[1].imshow(noisy_img, cmap='gray')
    axes[1].set_title(f"含盐雾噪声图像\nPSNR={psnr(clean_img, noisy_img):.2f}dB, SSIM={ssim(clean_img, noisy_img):.3f}")
    axes[1].axis('off')
    
    axes[2].imshow(median_filtered, cmap='gray')
    axes[2].set_title(f"中值滤波结果\nPSNR={metrics['中值滤波'][0]:.2f}dB, SSIM={metrics['中值滤波'][1]:.3f}")
    axes[2].axis('off')
    
    axes[3].imshow(tv_denoised, cmap='gray')
    axes[3].set_title(f"本文凸优化算法结果\nPSNR={metrics['本文凸优化算法'][0]:.2f}dB, SSIM={metrics['本文凸优化算法'][1]:.3f}")
    axes[3].axis('off')
    
    plt.tight_layout()
    plt.savefig("ship_denoising_result.png", dpi=300, bbox_inches='tight')
    plt.show()