1. 自编码器图像去噪  https://blog.csdn.net/baidu_28610773/article/details/84666918

2. TensorFlow实现卷积自编码器对图像进行去噪  https://blog.csdn.net/qq_34464926/article/details/80936150

3. 基于深度卷积神经网络的图像去噪方法 https://blog.csdn.net/qq_40716944/article/details/80098005

4. 基于深度学习的图像去噪（论文总结）https://blog.csdn.net/geekmanong/article/details/50572148

5. 深度学习与图像去噪总结 https://blog.csdn.net/weixin_41923961/article/details/80384499

6. 深度学习算法-图像去噪Noise2noise https://blog.csdn.net/weixin_40300585/article/details/90946856

#深度学习图像去噪发展概述
由于深度学习，特别是卷积神经网络(CNN)在图像识别等领域取得了较好的成果，近年来，基于深度学习的图像去噪方法也被提出并得到了发展。2008年，Viren Jain等提出用CNN处理自然图像的去噪问题[21]，得到了与常规方法(如小波变换和马尔可夫随机场)相近或更优的结果。并说明了特定形式的CNN可以被视为图像去噪的马尔可夫模型推断的结果的一种近似，但是神经网络模型可以避免了马尔可夫模型在概率学习和推断过程中的计算困难，从而降低了计算的复杂度。在训练神经网络的过程中为了更加快速和准确的收敛，采取了逐层训练的方法。该方法网络结构如下图所示：
图2.1 Viren Jain等的图像去噪卷积神经网络模型(CN1和CN2区别于训练集的不同)
该网络由四个隐含层组成，每个隐含层含有24个特征通道，每一个特征通道都与8个随机选择的上一层的特征通道相连接。该网络的卷积核大小为5×5，网络的损失函数为带噪图像经过神经网络去噪后的恢复图像和真实图像的误差平方和，即：。除了对噪声已知的情况进行训练和去噪，作者还测试了噪声未知情况下的去噪情况，将在σ=[0,100]的范围内随机变化的不同水平的噪声加入真实图像，并进行训练，此时将隐含层设置为5层，相对于传统方法也取得了较好的效果。
2012年，Xie等利用栈式去噪自编码器(stacked denoising auto-encoder)进行图像去噪以及图像修复(image inpainting)[22]，其网络结构为多层的全连接网络，如图2.2所示。
图2.2 栈式稀疏去噪自编码器(SSDA)
为了将稀疏编码的思想和深度神经网络结合，作者提出了栈式稀疏去噪自编码器(stacked sparse denoising auto-encoder, SSDA)，其损失函数如下：
其中KL散度中的ρ^ 代表隐含层激活函数输出的平均值，如果我们设定较小的ρ，则可以通过优化损失函数，使得每一层中被激活的神经元尽可能少，以此达到稀疏约束的目的。训练过程为逐层训练，每训练完一层后，将训练用的数据(带噪图像)和标签(真实图像)都通过该层，将输出送至下一层作为数据和标签进行训练。
Burger等人提出了用多层感知机(Multi Layer Perceptron，MLP)的方法实现神经网络去噪[23]，并说明如果训练集足够，MLP模型可以达到图像去噪的当前最优水平。MLP模型对噪声类型不敏感，对于不同类型的噪声，只要对带有该类噪声的图像进行处理，即可得到较好的结果。作者认为，MLP相对于CNN的优点在于，MLP可以近似所有函数，而CNN可以学习到的函数相对有限。为了使MLP网络达到最优水平，训练网络应足够大，分割图像得到的区块大小应足够大，以及训练集也应足够大。该方法的一个不足之处为它无法适应不同强度的噪声，若将不同水平的噪声图像作为输入进行训练则无法达到对特定噪声进行训练时的结果。
2016年，Mao等提出了深度的卷积编解码网络应用于图像去噪[24]。编码(Encoding)和解码(Decoding)结构对称，且每隔较少层就存在一条由编码卷积层到与其对称的解码层的跨层连接。该网络结构如图2.3所示。该网络共采用了对称的10层卷积层，其中前5层为编码层，后5层为解码层。在该CNN网络的编码层中，每隔2层就用一条跳线将编码层连接到对称的解码层，因此可以直接进行正向和反向传播。从该网络开始，图像去噪神经网络层数开始向更深层发展。以前的神经网络层数普遍较浅，原因在于：一方面，与图像分类等问题只需要提取关键特征的任务不同，图像的复原要求图像的细节在网络传播中不能丢失，而网络层数变深可能会导致图像细节的丢失；另一方面，层数较深容易在训练过程中导致梯度弥散，不容易收敛到最优解。由于该网络采用了卷积和解卷积层间的跳线连接，因此可以较好的解决上述的问题，因此网络相对较深。
图2.3 卷积编解码去噪网络结构(示跳线)
2017年，Zhang等提出用较深层的CNN网络，称为DnCNN[25]，实现去噪，为了解决网络层数加深导致的梯度弥散效应，DnCNN并不对图像进行学习，而是以输出与噪声的l2范数为损失函数来训练网络。网络结构如图2.4所示。DnCNN网络可以视为一个残差学习的过程，这样可以较好的训练。在该网络中利用了BN层(Batch Normalization)，实验表明BN层与残差学习共同使用可以提高模型的性能，DnCNN在不同噪声水平上训练，得到的结果要优于现在的最优结果，如BM3D等。
图2.4 DnCNN 网络结构

1. BM3D模型简介
BM3D模型是一个两阶段图像去噪方法，主要包含两个步骤: 
(1) 在噪声图像上，利用局部区域搜索相似块，并进行堆叠，在变换域(DCT域、FFT域）利用硬阈值去噪方法对堆叠的图像块进行去噪，获得堆叠相似块的估计值，最后，根据均值权重进行聚合； 
(2) 通过步骤(1) 获取初步估计的图像，在初步估计的图像上进行相似块的聚合； 然后，利用维纳协同滤波进行图像去噪，从而，获取最后的去噪结果

2. 模型实现(代码参考网络实现)： 
% BM3D_Color_Demo 
% BM3D 在彩色图像上去噪
% Author: HSW
% Date: 2018-05-06 
%

clc; 
close all; 
clear all; 

img_org = imread('timg.png'); 

figure(1); 
imshow(img_org); 
title('原图像'); 

% 加噪声
sigma = 25; 
img_noise = double(img_org)+sigma * randn(size(img_org));

figure; 
imshow(img_noise / 255, []); 
title('噪声图像'); 

img_denoise = BM3D_Color(img_noise, 0, sigma, 0, 1); 

figure; 
imshow(img_denoise / 255, []); 
title('去噪图像'); 

% BM3D_Gray_Demo
% BM3D 在灰度图像上去噪
% Author: HSW
% Date: 2018-05-06 
% 

clc; 
close all; 
clear all; 

img_org = imread('timg.png');

img_gray = rgb2gray(img_org); 

figure(1); 
imshow(img_gray); 
title('原图像'); 

% 加噪声
sigma = 25; 
img_noise = double(img_gray)+sigma * randn(size(img_gray));

figure; 
imshow(img_noise / 255, []); 
title('噪声图像'); 

img_denoise = BM3D_Gray(img_noise, 0, sigma, 1); 

figure; 
imshow(img_denoise / 255, []); 
title('去噪图像'); 

function img_denoise = BM3D_Color(img_noise, tran_mode, sigma, color_mode, isDisplay)
% BM3D实现去噪
% Inputs:
%       img_noise: 噪声图像
%       tran_mode: 变换方法： 默认值为0, tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1
%       sigma: 噪声水平，默认值为10
%       color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
%  Ouputs:
%       img_out: 去噪图像
% 参考文献：An Analysis and Implementation of the BM3D Image Denoising Method
% Inputs:
%        img_in: 噪声图像，必须为矩形方阵
%        tran_mode: = 0, FFT; = 1, DCT; = 2, DWT, = 3, db1
% Outputs:
%        img_denoise: 去噪图像
%
%
if ~exist('isDisplay', 'var')
    isDisplay = 0;
end

if ~exist('color_mode', 'var')
    color_mode = 0;
end

if ~exist('sigma', 'var')
    sigma = 10;
end

if ~exist('tran_mode', 'var')
    tran_mode = 0;
end

[row, col, dims] = size(img_noise);

img_trans = rgb2other(img_noise, color_mode);

% First Step 参数
kHard           = 8;          % 块大小
pHard           = 4;          % 块移动间隔
lambda_distHard = 0;          % 求相似的距离时，变换后，收缩的阈值
nHard           = 40;         % 搜索窗口大小
NHard           = 28;         % 最多相似块个数
tauHard         = 5000;       % 最大的相似距离for fft

% kaiser窗口的参数，实际上并没有特别大的影响
beta=2;
Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)';

% Second Step参数
kWien           = kHard;
pWien           = pHard;
lambda_distWien = lambda_distHard;
nWien           = nHard;
NWien           = NHard;
tauWien         = tauHard;
sigma2          = sigma*sigma;

if tran_mode == 0
    % FFT
    lambda2d=400;
    lambda1d=500;
    lambda2d_wie=50;
    lambda1d_wie=500;
elseif tran_mode == 1
    % DCT
    lambda2d=50;
    lambda1d=80;
    lambda2d_wie=20;
    lambda1d_wie=60;
elseif tran_mode == 2
    % DWT
    lambda2d=50;
    lambda1d=80;
    lambda2d_wie=20;
    lambda1d_wie=60;
end

fprintf('BM3D: First Stage Start...\n');
%block为原始图像块， tran_block为FFT变换且硬阈值截断后的频域系数(频域， 计算距离的时候采用的是变换块)
[block_ch1, tran_block_ch1, block2row_idx_ch1, block2col_idx_ch1] = im2block(img_trans(:,:,1), kHard, pHard, lambda_distHard, 0);
[block_ch2, tran_block_ch2, block2row_idx_ch2, block2col_idx_ch2] = im2block(img_trans(:,:,2), kHard, pHard, lambda_distHard, 0); 
[block_ch3, tran_block_ch3, block2row_idx_ch3, block2col_idx_ch3] = im2block(img_trans(:,:,3), kHard, pHard, lambda_distHard, 0); 

%bn_r和bn_c为行和列上的图像块个数
bn_r = floor((row - kHard) / pHard) + 1;
bn_c = floor((col - kHard) / pHard) + 1;
%基础估计的图像
img_basic_sum = zeros(row, col, 3);
img_basic_weight = zeros(row, col, 3);
%对每个块遍历
for i=1:bn_r
    for j=1:bn_c
        
        % 利用亮度通道进行相似块搜索
        [sim_blk_ch1, sim_num, sim_blk_idx] = search_similar_block(i, j, block_ch1, tran_block_ch1, floor(nHard/pHard), bn_r, bn_c, tauHard, NHard);
        
        % 进行亮度通道处理
        % 协同滤波： 公式(2)
        tran3d_blk_shrink_ch1 = transform_3d(sim_blk_ch1, tran_mode, lambda2d, lambda1d);
        tran3d_blk_shrink_ch2 = transform_3d(block_ch2(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); 
        tran3d_blk_shrink_ch3 = transform_3d(block_ch3(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); 
        % 聚合： 公式(3)中的说明
        NHard_P_ch1 = nnz(tran3d_blk_shrink_ch1);
        NHard_P_ch2 = nnz(tran3d_blk_shrink_ch2); 
        NHard_P_ch3 = nnz(tran3d_blk_shrink_ch3); 
        
        if NHard_P_ch1 > 1
            wHard_P_ch1 = 1 / NHard_P_ch1;
        else
            wHard_P_ch1 = 1;
        end
        
        if NHard_P_ch2 > 1
            wHard_P_ch2 = 1 / NHard_P_ch2; 
        else
            wHard_P_ch2 = 1; 
        end 
        
        if NHard_P_ch3 > 1
            wHard_P_ch3 = 1 / NHard_P_ch3; 
        else
            wHard_P_ch3 = 1; 
        end 
        
        blk_est_ch1 = inv_transform_3d(tran3d_blk_shrink_ch1,tran_mode);
        blk_est_ch1 = real(blk_est_ch1);
        
        blk_est_ch2 = inv_transform_3d(tran3d_blk_shrink_ch2, tran_mode); 
        blk_est_ch2 = real(blk_est_ch2); 
        
        blk_est_ch3 = inv_transform_3d(tran3d_blk_shrink_ch3, tran_mode); 
        blk_est_ch3 = real(blk_est_ch3); 
        
        % 公式(3): 对亮度通道，即第1个通道
        for k=1:sim_num
            idx = sim_blk_idx(k);
            ir = block2row_idx_ch1(idx);
            jr = block2col_idx_ch1(idx);
            img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1 * blk_est_ch1(:, :, k);
            img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1;
            
            img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2 * blk_est_ch2(:, :, k);
            img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2;
            
            img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3 * blk_est_ch3(:, :, k);
            img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3;
        end
        
    end
end

img_basic = img_basic_sum ./ img_basic_weight;

if isDisplay
    figure;
    img_rgb = other2rgb(img_basic, color_mode); 
    imshow(img_rgb / 255.0 ,[]);
    title('BM3D:Fist Stage Result');
end

fprintf('BM3D: First Stage End...\n');

fprintf('BM3D: Second Stage Start...\n');

[block_basic_ch1,tran_block_basic_ch1,block2row_idx_basic_ch1,block2col_idx_basic_ch1] = im2block(img_basic(:, :, 1), kWien, pWien, lambda_distWien, 0);
[block_basic_ch2,tran_block_basic_ch2,block2row_idx_basic_ch3,block2col_idx_basic_ch2] = im2block(img_basic(:, :, 2), kWien, pWien, lambda_distWien, 0);
[block_basic_ch3,tran_block_basic_ch3,block2row_idx_basic_ch3,block2col_idx_basic_ch3] = im2block(img_basic(:, :, 3), kWien, pWien, lambda_distWien, 0);

bn_r = floor((row - kWien) / pWien) + 1;
bn_c = floor((col - kWien) / pWien) + 1;

img_wien_sum = zeros(row, col, 3);
img_wien_weight = zeros(row, col, 3);

for i=1:1:bn_r
    for j=1:1:bn_c
        % 公式(5), 利用亮度进行相似性搜索
        [sim_blk_basic_ch1, sim_num, sim_blk_basic_idx] = search_similar_block(i, j, block_basic_ch1, tran_block_basic_ch1, floor(nWien/pWien), bn_r, bn_c, tauWien, NWien);
        
        % 公式(6)
        tran3d_blk_basic_ch1 = transform_3d(sim_blk_basic_ch1, tran_mode, lambda2d_wie, lambda1d_wie);
        tran3d_blk_basic_ch2 = transform_3d(block_basic_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
        tran3d_blk_basic_ch3 = transform_3d(block_basic_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
        
        omega_P_ch1 = (tran3d_blk_basic_ch1.^2) ./ ((tran3d_blk_basic_ch1.^2) + sigma2);
        omega_P_ch2 = (tran3d_blk_basic_ch2.^2) ./ ((tran3d_blk_basic_ch2.^2) + sigma2); 
        omega_P_ch3 = (tran3d_blk_basic_ch3.^2) ./ ((tran3d_blk_basic_ch3.^2) + sigma2); 
        
        % 公式(7)
        tran3d_blk_ch1 = transform_3d(block_ch1(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie);
        tran3d_blk_ch2 = transform_3d(block_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
        tran3d_blk_ch3 = transform_3d(block_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
        
        blk_est_ch1 = inv_transform_3d(omega_P_ch1 .* tran3d_blk_ch1, tran_mode);
        blk_est_ch2 = inv_transform_3d(omega_P_ch2 .* tran3d_blk_ch2, tran_mode); 
        blk_est_ch3 = inv_transform_3d(omega_P_ch3 .* tran3d_blk_ch3, tran_mode); 
        
        blk_est_ch1 = real(blk_est_ch1);
        blk_est_ch2 = real(blk_est_ch2); 
        blk_est_ch3 = real(blk_est_ch3); 
        
        NWien_P_ch1 = nnz(omega_P_ch1);
        NWien_P_ch2 = nnz(omega_P_ch2); 
        NWien_P_ch3 = nnz(omega_P_ch3); 
        
        if NWien_P_ch1 > 1
            wWien_P_ch1 = 1 / (NWien_P_ch1);
        else
            wWien_P_ch1 = 1;
        end
        
        if NWien_P_ch2 > 1
            wWien_P_ch2 = 1/(NWien_P_ch2);
        else
            wWien_P_ch2 = 1;
        end
        
        if NWien_P_ch3 > 1
            wWien_P_ch3 = 1 / (NWien_P_ch3);
        else
            wWien_P_ch3 = 1;
        end
        % 公式(8)
        for k=1:sim_num
            idx=sim_blk_basic_idx(k);
            ir=block2row_idx_basic_ch1(idx);
            jr=block2col_idx_basic_ch1(idx);
            img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1 * blk_est_ch1(:, :, k);
            img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1;
            
            img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2 * blk_est_ch2(:, :, k);
            img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2;
            
            img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3 * blk_est_ch3(:, :, k);
            img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3;
        end
    end
end

img_other = img_wien_sum ./ img_wien_weight; 
img_denoise = other2rgb(img_other, color_mode);

fprintf('BM3D: Second Stage End\n');

function img_denoise = BM3D_Gray(img_noise, tran_mode, sigma, isDisplay)
% 参考文献：An Analysis and Implementation of the BM3D Image Denoising Method
% Inputs:
%        img_noise: 灰度噪声图像，必须为矩形方阵
%        tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1
% Outputs:
%        img_denoise: 去噪图像
%
if ~exist('tran_mode', 'var')
    tran_mode = 0;
end

if ~exist('sigma', 'var')
    sigma = 10;
end

if ~exist('isDisplay', 'var')
    isDisplay = 0;
end

[row,col] = size(img_noise);

% First Step 参数
kHard           = 8;          % 块大小
pHard           = 4;          % 块移动间隔
lambda_distHard = 0;          % 求相似的距离时，变换后，收缩的阈值
nHard           = 40;         % 搜索窗口大小
NHard           = 28;         % 最多相似块个数
tauHard         = 5000;       % 最大的相似距离for fft

% kaiser窗口的参数，实际上并没有特别大的影响
beta=2;
Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)';

% Second Step参数
kWien           = kHard;
pWien           = pHard;
lambda_distWien = lambda_distHard;
nWien           = nHard;
NWien           = NHard;
tauWien         = tauHard;
sigma2          = sigma*sigma;

if(tran_mode==0)        %fft
    lambda2d=400;
    lambda1d=500;
    lambda2d_wie=50;
    lambda1d_wie=500;
elseif(tran_mode == 1)  %dct
    lambda2d=50;
    lambda1d=80;
    lambda2d_wie=20;
    lambda1d_wie=60;
elseif(tran_mode == 2)  %dwt
    lambda2d=50;
    lambda1d=80;
    lambda2d_wie=20;
    lambda1d_wie=60;
end

%block为原始图像块， tran_block为FFT变换且硬阈值截断后的频域系数(频域， 计算距离的时候采用的是变换块)
[block,tran_block,block2row_idx,block2col_idx]=im2block(img_noise,kHard,pHard,lambda_distHard,0);

%bn_r和bn_c为行和列上的图像块个数
bn_r=floor((row-kHard)/pHard)+1;
bn_c=floor((col-kHard)/pHard)+1;
%基础估计的图像
img_basic_sum=zeros(row,col);
img_basic_weight=zeros(row,col);
%basic处理
fprintf('BM3D: First Stage Start...\n');
%对每个块遍历
for i=1:bn_r
    for j=1:bn_c
        [sim_blk,sim_num,sim_blk_idx]=search_similar_block(i,j,block,tran_block,floor(nHard/pHard),bn_r,bn_c,tauHard,NHard);
        % 协同滤波： 公式(2)
        tran3d_blk_shrink=transform_3d(sim_blk,tran_mode,lambda2d,lambda1d);
        % 聚合： 公式(3)中的说明
        NHard_P=nnz(tran3d_blk_shrink);
        if(NHard_P >1)
            wHard_P=1/NHard_P;
        else
            wHard_P=1;
        end
        blk_est =inv_transform_3d(tran3d_blk_shrink,tran_mode);
        blk_est=real(blk_est);
        
        % 公式(3)
        for k=1:sim_num
            idx=sim_blk_idx(k);
            ir=block2row_idx(idx);
            jr=block2col_idx(idx);
            img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1) = img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1) + wHard_P*blk_est(:,:,k);
            img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1) = img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1) + wHard_P;
        end
    end
end
fprintf('BM3D: First Stage End...\n');
img_basic=img_basic_sum./img_basic_weight;

if isDisplay
    figure;
    imshow(img_basic,[]);
    title('BM3D:Fist Stage Result');
end


[block_basic,tran_block_basic,block2row_idx_basic,block2col_idx_basic] = im2block(img_basic,kWien,pWien,lambda_distWien,0);
bn_r=floor((row-kWien)/pWien)+1;
bn_c=floor((col-kWien)/pWien)+1;
img_wien_sum=zeros(row,col);
img_wien_weight=zeros(row,col);

fprintf('BM3D: Second Stage Start...\n');
for i=1:1:bn_r
    for j=1:1:bn_c
        % 公式(5)
        [sim_blk_basic,sim_num,sim_blk_basic_idx] = search_similar_block(i,j,block_basic,tran_block_basic,floor(nWien/pWien),bn_r,bn_c,tauWien,NWien);
        % 公式(6)
        tran3d_blk_basic = transform_3d(sim_blk_basic,tran_mode,lambda2d_wie,lambda1d_wie);
        omega_P=(tran3d_blk_basic.^2)./((tran3d_blk_basic.^2)+sigma2);
        % 公式(7)
        tran3d_blk = transform_3d(block(:,:,sim_blk_basic_idx),tran_mode,lambda2d_wie,lambda1d_wie);
        blk_est=inv_transform_3d(omega_P.*tran3d_blk,tran_mode);
        blk_est=real(blk_est);
        NWien_P=nnz(omega_P);
        if(NWien_P >1)
            wWien_P=1/(NWien_P);
        else
            wWien_P=1;
        end
        % 公式(8)
        for k=1:sim_num
            idx=sim_blk_basic_idx(k);
            ir=block2row_idx_basic(idx);
            jr=block2col_idx_basic(idx);
            img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1) = img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1) + wWien_P*blk_est(:,:,k);
            img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1) = img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1) + wWien_P;
        end
    end
end
fprintf('BM3D: Second Stage End\n');

img_denoise = img_wien_sum./img_wien_weight;


function [block,transform_block,block2row_idx,block2col_idx] =im2block(img,k,p,lambda2D,delta)
% 实现图像分块
% Inputs：
%        k: 块大小
%        p: 块移动步长
%        lambda_2D: 收缩阈值
%        delta: 收缩阈值
%  Outputs:
%        block: 返回的块
%        transform_block: 变换后的块
%        block2row_idx: 块索引与图像块的左上角行坐标对应关系
%        block2col_idx: 块索引与图像块的左上角列坐标对应关系
%
[row,col] = size(img);
% 频域去噪中的硬阈值，实际上原文中，对于噪声方差小于40时thres = 0, 具体见公式(1)的说明第2点(即距离计算)
thres = lambda2D*delta*sqrt(2*log(row*col));
% r_num 和 c_num分别表示行和列上可以采集的块的数目
r_num = floor((row-k)/p)+1;
c_num = floor((col-k)/p)+1;
block = zeros(k,k,r_num*c_num);
block2row_idx = [];
block2col_idx = [];
cnt = 1;
for i = 0:r_num-1
    rs = 1+i*p;
    for j = 0:c_num-1
        cs = 1+j*p;
        block(:,:,cnt) = img(rs:rs+k-1,cs:cs+k-1);
        block2row_idx(cnt) = rs;
        block2col_idx(cnt) = cs;
        tr_b = fft2(block(:,:,cnt));
        idx = find(abs(tr_b)<thres);
        tr_b(idx) = 0;
        transform_block(:,:,cnt) = tr_b;
        cnt = cnt+1;
    end
end
end

function [blk_est]=inv_transform_3d(blk_tran3d,tran_mode)
% 3D 逆变换
% Inputs:
%       blk_tran3d: 在频域中，硬阈值滤波的图像块
%       tran_mode: 变换方法
% Outputs:
%       blk_est:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=size(blk_tran3d);

blk_invtran1d=zeros(m,n,blk_num);
blk_est=zeros(m,n,blk_num);

if(tran_mode==0)    %fft
    for i=1:1:m
        for j=1:1:n
            blk_invtran1d(i,j,:)=ifft(blk_tran3d(i,j,:));
        end
    end
    for i=1:1:blk_num
        blk_est(:,:,i)=ifft2(blk_invtran1d(:,:,i));
    end
elseif(tran_mode==1)  %dct
    for i=1:1:m
        for j=1:1:n
            blk_invtran1d(i,j,:)=idct(blk_tran3d(i,j,:));
        end
    end
    for i=1:1:blk_num
        blk_est(:,:,i)=idct2(blk_invtran1d(:,:,i));
    end
elseif(tran_mode==2)    %dwt
    blk_num=length(blk_2d_s);
    blk_c=waverec2(blk_tran3d,blk_tran1d_s,'haar');
    blk_est=[];
    for i=1:1:blk_num
        blk_est(:,:,i)=waverec2(blk_c(:,i),blk_2d_s{i},'Bior1.5');
    end
    
else
    error('tran_mode error');
end

end

function img_trans = other2rgb(img_in, color_mode)
% 将RGB颜色空间转为其他颜色空间
% Inputs:
%        img_in: RGB颜色空间图像
%        color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
%        img_trans: 其他颜色空间
%
% Author: HSW
% Date: 2018-05-06

img_trans = zeros(size(img_in));
[row, col, dims] = size(img_in);

if color_mode == 0
    color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
    color_tran_inv = inv(color_tran);
    for i = 1:row
        for j = 1:col
            other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
            img_trans(i, j, :) = color_tran_inv * other;
        end
    end
elseif color_mode == 1
    color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
    color_tran_inv = inv(color_tran);
    for i = 1:row
        for j = 1:col
            other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
            img_trans(i, j, :) = color_tran_inv * other;
        end
    end
elseif color_mode == 2
    color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
    color_tran_inv = inv(color_tran);
    for i = 1:row
        for j = 1:col
            other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
            img_trans(i, j, :) = color_tran_inv * other;
        end
    end
end

end

function img_trans = rgb2other(img_in, color_mode)
% 将RGB颜色空间转为其他颜色空间
% Inputs:
%        img_in: RGB颜色空间图像
%        color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
%        img_trans: 其他颜色空间
%
% Author: HSW
% Date: 2018-05-06

img_trans = zeros(size(img_in));
[row, col, dims] = size(img_in);
if color_mode == 0
    color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
    for i = 1:row
        for j = 1:col
            rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)]; 
            img_trans(i, j, :) = (color_tran * rgb)';
        end
    end
elseif color_mode == 1
    color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
    for i = 1:row
        for j = 1:col
            rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)]; 
            img_trans(i, j, :) = (color_tran * rgb)';
        end
    end
elseif color_mode == 2
    color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
    for i = 1:row
        for j = 1:col
            rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)]; 
            img_trans(i, j, :) = (color_tran * rgb)';
        end
    end
end

end

function [sim_blk,sim_num,sim_blk_idx]=search_similar_block(ik,jk,block,tran_block,np,bn_r,bn_c,tau,max_sim_num)
% 搜索相似块
% Inputs:
%       ik, jk： 待搜索相似块的索引
%       block: 图像块集合
%       tran_block： 图像块FFT硬阈值过滤后的FFT系数
%       k: 图像块大小
%       np: floor(nHard / pHard)， 其中nHard表示图像的搜索区域大小， pHard表示块的移动步长
%       bn_r, bn_c: 图像总的行/列可以采集图像块的数目
%       tau: 图像块相似性判断阈值，见公式(1)
%       max_sim_num: 最多保留相似块的数目
% Ouputs:
%       sim_blk:
%       sim_num:
%       sim_blk_idx:
%
% 搜索窗口的左上角，右下角的块索引
in_s = max(ik-floor(np/2),1);
jn_s = max(jk-floor(np/2),1);
in_e = min(ik+floor(np/2),bn_r);
jn_e = min(jk+floor(np/2),bn_c);
% 当前参考块
ref_blk = tran_block(:,:,((ik-1)*bn_c+jk));
ii = in_s:1:in_e;
jj = jn_s:1:jn_e;
[II,JJ] = meshgrid(ii,jj);
IDX = (II-1)*bn_c+JJ;
blk_idx=IDX(:);
% 收缩范围内的全部图像块
cur_blk=tran_block(:,:,blk_idx);
cnt=size(cur_blk,3);
ref_blk_mat=repmat(ref_blk,[1,1,cnt]);
delta_blk=cur_blk-ref_blk_mat;
dist=sum(sum(delta_blk.*delta_blk,1),2);
[dist_sort,dist_idx]=sort(dist);
% 最大相似块是真实相似块和目标参数相似块的最小值
max_num=min(cnt,max_sim_num);
if(dist_sort(max_num)<tau)
    sim_num=max_num;
else
    sim_num=sum(dist_sort(1:max_num)<tau);
end
cnt_idx=dist_idx(1:sim_num);
sim_blk_idx=blk_idx(cnt_idx);
sim_blk=block(:,:,sim_blk_idx);
end

function [val]=thres_shrink(data,thres)
% 进行阈值截断： 即 data(i) < thres ? data(i) = 0 : data(i) = data(i)
% Inputs:
%       data: 阈值截断前的数据
%       thres: 阈值
% Ouputs:
%       val: 阈值截断后的数据
% 
val=data;
idx=find(abs(data)<thres);
val(idx)=0;
end

function blk_tran3d = transform_3d(blk_3d,tran_mode,lambda2d,lambda1d)
% 进行3D变换，即Collaborative Filtering: 在图像块内进行2D变换，在图像块间进行1D变换
% 公式(2)
% Inputs:
%        blk_3d:
%        tran_mode:
% Ouputs:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=size(blk_3d);

%变换不同时，可能需要修改？？
blk_2d_shrink=zeros(m,n,blk_num);
blk_1d_shrink=zeros(m,n,blk_num);

if(tran_mode==0)    %fft
    for i=1:1:blk_num
        blk_tran2d = fft2(blk_3d(:,:,i));
        blk_2d_shrink(:,:,i) = thres_shrink(blk_tran2d,lambda2d);
    end
    for i=1:1:m
        for j=1:1:n
            blk_tran1d = fft(blk_2d_shrink(i,j,:));
            blk_1d_shrink(i,j,:) = thres_shrink(blk_tran1d,lambda1d);
        end
    end
    blk_tran3d=blk_1d_shrink;
    
elseif(tran_mode==1)  %dct
    for i=1:1:blk_num
        blk_tran2d=dct2(blk_3d(:,:,i));
        blk_2d_shrink(:,:,i)=thres_shrink(blk_tran2d,lambda2d);
    end
    for i=1:1:m
        for j=1:1:n
            blk_tran1d=dct(blk_2d_shrink(i,j,:));
            blk_1d_shrink(i,j,:)=thres_shrink(blk_tran1d,lambda1d);
        end
    end
    blk_tran3d=blk_1d_shrink;
    
elseif(tran_mode==2)    %dwt
    blk_2d_s={};
    blk_2d_shrink=[];%zeros()
    for i=1:1:blk_num
        [blk_tran2d_c,blk_tran2d_s]=wavedec2(blk_3d(:,:,i),2,'Bior1.5');
        blk_2d_shrink(:,i)=thres_shrink(blk_tran2d_c,lambda2d);
        blk_2d_s{i}=blk_tran2d_s;
    end
    %这里应该用 wavedec.因为是对1维？？
    [blk_tran1d_c,blk_tran1d_s]=wavedec2(blk_2d_shrink,1,'haar');
    blk_tran3d=thres_shrink(blk_tran1d_c,lambda1d);
    %   elseif(strcmp(tran_mode,'db1')) %还未实现
    %       blk_2d_s={};
    %       blk_2d_shrink=[];%zeros()
    %       for i=1:1:blk_num
    %           [blk_tran2d_cA,blk_tran2d_cH,blk_tran2d_cV,blk_tran2d_cD]=...
    %               dwt2(blk_3d(:,:,i),'db1');
    %           blk_2d_shrink(:,i)=thres_shrink(blk_tran2d_c,lambda2d);
    %           blk_2d_s{i}=blk_tran2d_s;
    %       end
    %       [blk_tran1d_c,blk_tran1d_s]=wavedec2(blk_2d_shrink,1,'haar');
    %       blk_tran3d=thres_shrink(blk_tran1d_c,lambda1d);
else
    error('tran_mode error');
end
end

3. 模型效果：
3.1 灰度图像

3.2 彩色图像

Github个人博客：https://joeyos.github.io
图像去噪常用方法
图像去噪处理方法可分为空间域法和变换域法两大类。
基于离散余弦变换的图像去噪
一般而言，我们认为图像的噪声在离散余弦变换结果中处在其高频部分，而高频部分的幅值一般很小，利用这一性质，就可以实现去噪。然而，同时会失去图像的部分细节。
%读取图像
X=imread('wangshi.jpg'); 
X=rgb2gray(X);
%读取图像尺寸
[m,n]=size(X); 
%给图像加噪
Xnoised=imnoise(X,'speckle',0.01); 
%输出加噪图像
subplot(121); 
imshow(Xnoised);
%DCT变换
Y=dct2(Xnoised); 
I=zeros(m,n);
%高频屏蔽
I(1:m/3,1:n/3)=1; 
Ydct=Y.*I;
%逆DCT变换
Y=uint8(idct2(Ydct)); 
%结果输出
subplot(122);
imshow(Y);



基于小波变换的图像去噪
小波去噪是小波变换较为成功的一类应用，其去噪的基本思路为：含噪图像-小波分解-分尺度去噪-小波逆变换-恢复图像。含噪信号经过预处理，然后利用小波变换把信号分解到各尺度中，在每一尺度下把属于噪声的小波系数去掉，保留并增强属于信号的小波系数，最后再经过小波逆变换恢复检测信号。比基于傅里叶变换的去噪方法好。
clear;                 
X=imread('life.jpg');            
X=rgb2gray(X);
subplot(221);          
imshow(X);             
title('原始图像');                  
% 生成含噪图像并图示
init=2055615866;       
randn('seed',init);      
X=double(X);
% 添加随机噪声
XX=X+8*randn(size(X));  
subplot(222);             
imshow(uint8(XX));              
title(' 含噪图像 ');       
%用小波函数coif2对图像XX进行2层
% 分解
[c,l]=wavedec2(XX,2,'coif2'); 
% 设置尺度向量
n=[1,2];                  
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08]; 
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');   
subplot(223);              
imshow(uint8(X1));                
%colormap(map);            
title(' 第一次消噪后的图像 '); 
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
X2=waverec2(mc,l,'coif2');  
subplot(224);             
imshow(uint8(X2));               
title(' 第二次消噪后的图像 ');