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 网络结构

齐艳丽
【摘 要】图像噪声会严重影响其视觉效果，采用适当的方法进行图像去噪是一项必不可少的预处理操作。本文对几种常用的图像去噪方法进行了分析，通过仿真高斯噪声和椒盐噪声，着重对基于均值滤波、中值滤波以及卷积神经网络的实际去噪效果进行了比较。实验结果表明三种方法的去噪效果各有优劣，均能实现一定程度上的图像去噪。
【关键词】图像噪声，均值滤波，中值滤波，卷积神经网络
中图分类号： TP391.41文献标识码： A文章编号： 2095-2457(2019)26-0024-002
DOI：10.19694/j.cnki.issn2095-2457.2019.26.010
【Abstract】Image noise affects its visual effects seriously.Appropriate method for image denoising is an essential pre-processing operation.Several common methods for  image denoising are analyzed in this paper.Comparing the actual denoising effects based on mean filtering，median filtering and convolutional neural network with simulating Gaussian noise and salt and pepper noise.The experimental results show that the denoising effects of the three methods have their own advantages and disadvantages， and can achieve a certain degree of image denoising.
【Key words】Image noise;On mean filtering;Median filtering;Convolutional neural network
0 引言
由于实际应用中图像不可避免地被噪声污染，这类噪声对于从图像中提取有效信息产生严重干扰[1]。因此，去除噪声是一项必不可少的图像预处理操作，而有效的去噪算法一直是图像处理领域的研究热点。
本文选取常见的三种图像去噪方法，即均值滤波、中值滤波以及卷积神經网络，进行研究。通过仿真高斯噪声和椒盐噪声，对这三种去噪方法的实际效果进行比较和分析。
1 去噪方法原理
中值滤波和均值滤波是经典的图像去噪方法，而基于卷积神经网络的图像去噪是一种随着人工智能以及机器学习发展而兴起的去噪方法[2]。
中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术[3]。其基本原理是以当前像素点(x，y，z)的邻域中各点灰度值的中值代替当前像素点的灰度值，让周围的像素值接近的真实值，从而消除孤立的噪声点，即：
不同于上述两种算法，基于卷积神经网络的图像去噪是随着深度学习特别是卷积神经网络(CNN)在图像识别等领域的研究和应用而被提出的[5]。通过对大量图像集的学习和训练，导出训练好并经验证的卷积神经网络来实现图像去噪。本文直接基于MATLAB所提供DnCNN进行去噪效果的分析和比较。
2 仿真结果分析
为了分析和比较三种去噪算法的效果，对原始图像分别加入噪声浓度为0.05的高斯噪声和椒盐噪声，如图1-3所示。可以看出，加入噪声后的图像对提取图像中的有效信息造成了严重的干扰，因此去噪显得尤为重要。
图4和图5为分别经过中值滤波和均值滤波处理后的含有噪声的图像。从图4中可以看出，使用中值滤波对高斯噪声图像进行去噪后图像颜色严重失真。另一方面，中值滤波对椒盐图像的去噪效果较好，基本还原了原图像。从图5可以看出，均值滤波可以实现对高斯噪声图像和椒盐噪声图像的较好的去噪效果。总体上来看，均值滤波好中值滤波对高斯噪声和椒盐噪声的去噪效果相当。
图6的结果中可以看出，本文所使用的DnCNN对于高斯白噪声图像的去噪效果较好，但是对于添加了椒盐噪声后的图像去噪效果一般。其主要原因可能是所使用的卷积神经网络主要针对高斯噪声图像进行的训练，所以并不完全适用于椒盐噪声图像的处理。因此，使用卷积神经网络进行图像去噪时需要进行有针对性的训练。
3 总结
本文对三种常用的图像去噪算法进行了比较和分析，结果表明，在三种方法中，基于卷积神经网络的图像去噪方法能够实现最好的高斯噪声去噪处理，而中值滤波和均值滤波可以对高斯噪声图像和椒盐噪声图像都能实现比较稳定均衡的去噪效果。值得注意的是，由于所使用的卷积神经网络主要针对高斯噪声训练，因此其对椒盐噪声图像的去噪效果在三种方法中表现最差，需要在下一步的工作中对椒盐噪声图像的处理进行有针对性的卷积神经网络训练。
【参考文献】
[1]刘祝华.图像去噪方法的研究[D].江西师范大学，2005.
[2]朱习军，隋思涟，张宾，刘尊年.MATLAB在信号与图像处理中的应用[M].电子工业出版社，2009.
[3]赵高长，张磊，武风波.改进的中值滤波算法在图像去噪中的应用[J].应用光学，2011，32(4)：678-682.
[4]张丽果.快速非局部均值滤波图像去噪[J].信号处理， 2013，29(8)：1043-1049.
[5]吕永标，赵建伟，曹飞龙.基于复合卷积神经网络的图像去噪算法[J].模式识别与人工智能，2017，30(2)：97-105.

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 彩色图像