基于DCT字典图像稀疏去噪算法学习

基于DCT字典图像稀疏去噪算法学习

理论基础：

评价一副图像质量的指标（MSE和PSNR）：

1.MSE（均方误差）：



其中，f'(i,j)和f(i,j)分别表示的是待评价图像和原始图像，M,N分别表示图像的长与宽。

2.PSNR（峰值信噪比）:

PSNR本质上与MSE相同，但它与图像的量化灰度级相联系，其表达式为：

主函数：main.m

%               demo2 - denoise an image
% this is a run_file the demonstrate how to denoise an image,
% using dictionaries. The methods implemented here are the same
% one as described in "Image Denoising Via Sparse and Redundant
% representations over Learned Dictionaries", (appeared in the
% IEEE Trans. on Image Processing, Vol. 15, no. 12, December 2006).
% ============================================================

clear
bb=8; % block size
RR=4; % redundancy factor  冗余因素
K=RR*bb^2; % number of atoms in the dictionary

sigma = 50;
%pathForImages ='';
%imageName = 'barbara.png';
% [IMin0,pp]=imread('cameraman.tif');
[IMin0,pp]=imread('w.jpg');
IMin0=im2double(IMin0);
if (length(size(IMin0))>2)
IMin0 = rgb2gray(IMin0);
end
if (max(IMin0(:))<2)
IMin0 = IMin0*255;
end

IMin=IMin0+sigma*randn(size(IMin0));%%%%%%此处有随机函数
PSNRIn = 20*log10(255/sqrt(mean((IMin(:)-IMin0(:)).^2)));

tic
[IoutDCT,output] = denoiseImageDCT(IMin, sigma, K);

PSNROut = 20*log10(255/sqrt(mean((IoutDCT(:)-IMin0(:)).^2)));
figure;
subplot(1,3,1); imshow(IMin0,[]); title('Original clean image');
subplot(1,3,2); imshow(IMin,[]); title(strcat(['Noisy image, ',num2str(PSNRIn),'dB']));
subplot(1,3,3);
imshow(IoutDCT,[]); title(strcat(['Clean Image by DCT dictionary, ',num2str(PSNROut),'dB']));
figure;
I = displayDictionaryElementsAsImage(output.D, floor(sqrt(K)), floor(size(output.D,2)/floor(sqrt(K))),bb,bb,0);
title('The DCT dictionary');
toc

核心子函数：denoiseImageDCT.m

function [IOut,output] = denoiseImageDCT(Image,sigma,K,varargin)
Reduce_DC = 1;
[NN1,NN2] = size(Image);
C = 1.15;%%%%%%%%%
waitBarOn = 1;
maxBlocksToConsider = 260000;
slidingDis = 1;
bb = 8;
errT = C*sigma;
% Create an initial dictionary from the DCT frame
Pn=ceil(sqrt(K));%%%Pn=16  bb=8  产生64*256的字典
DCT=zeros(bb,Pn);
for k=0:1:Pn-1,
V=cos([0:1:bb-1]'*k*pi/Pn);
if k>0,
V=V-mean(V);
end;
DCT(:,k+1)=V/norm(V);%norm(V)表示的是欧式距离
end;
%产生64*256的字典
DCT=kron(DCT,DCT);%http://blog.sina.com.cn/s/blog_7671b3eb0101132y.html
%%%%%%%%例如： a=[1,2;3,4]   kron(a,a)
while (prod(floor((size(Image)-bb)/slidingDis)+1)>maxBlocksToConsider)
slidingDis = slidingDis+1;
% Default value is 1. However, if the image is
% large, this number increases automatically (because of
% memory requirements)
end

[blocks,idx] = my_im2col(Image,[bb,bb],slidingDis);

if (waitBarOn)
counterForWaitBar = size(blocks,2);
h = waitbar(0,'Denoising In Process ...');
end
% go with jumps of 10000
for jj = 1:10000:size(blocks,2)
if (waitBarOn)
waitbar(jj/counterForWaitBar);
end
jumpSize = min(jj+10000-1,size(blocks,2));
if (Reduce_DC)%Reduce_DC为1
vecOfMeans = mean(blocks(:,jj:jumpSize));%取出10000列
blocks(:,jj:jumpSize) = blocks(:,jj:jumpSize) - repmat(vecOfMeans,size(blocks,1),1);%repmat(vecOfMeans,size(blocks,1),1):将vecOfMeans重复64行1列  公式3.12
end
%Coefs = mexOMPerrIterative(blocks(:,jj:jumpSize),DCT,errT);
Coefs = OMPerr(DCT,blocks(:,jj:jumpSize),errT);

if (Reduce_DC)%Reduce_DC为1
%DCT为64*256  Coefs256*10000         blocks为64*62001       vecOfMeans为1*10000
blocks(:,jj:jumpSize)= DCT*Coefs + ones(size(blocks,1),1) * vecOfMeans;%%%%blocks矩阵每一列的数值等于  该列乘以稀疏系数加上该列的平均值

else
blocks(:,jj:jumpSize)= DCT*Coefs ;
end
end

count = 1;
Weight= zeros(NN1,NN2);
IMout = zeros(NN1,NN2);
[rows,cols] = ind2sub(size(Image)-bb+1,idx);
for i  = 1:length(cols)
col = cols(i);
row = rows(i);
block =reshape(blocks(:,count),[bb,bb]);%block为值均相同的8*8的矩阵
%block
IMout(row:row+bb-1,col:col+bb-1)=IMout(row:row+bb-1,col:col+bb-1)+block;
Weight(row:row+bb-1,col:col+bb-1)=Weight(row:row+bb-1,col:col+bb-1)+ones(bb);%%%累计加了多少个block，到最后除以累加的次数即可
count = count+1;
end;
if (waitBarOn)
close(h);
end
C=100000;
IOut =(Image+C*sigma*IMout)./(1+C*sigma*Weight);%%%%%%%%%0.034,为什么要这样？？？  我改成C了，个人理解觉得这个值有，但是只要大于0就能得到去噪的效果
output.D = DCT;
% IOut =  IMout./Weight;  %%%其实只要这样就可以还原了，但是考虑到随机sigma的影响，后面加了一些变成(Image+C*sigma*IMout)./(1+C*sigma*Weight)
% IMout(100)
% Weight
%   blocks

子函数：my_im2col.m

function [blocks,idx] = my_im2col(I,blkSize,slidingDis);

if (slidingDis==1)
blocks = im2col(I,blkSize,'sliding');%行为blksize元素的总个数，列为(m-bb+1) x (n-bb+1)=62001
% http://fuda641.blog.163.com/blog/static/20751421620135483846711/ idx = [1:size(blocks,2)];
return
end

idxMat = zeros(size(I)-blkSize+1);
idxMat([[1:slidingDis:end-1],end],[[1:slidingDis:end-1],end]) = 1; % take blocks in distances of 'slidingDix', but always take the first and last one (in each row and column).
idx = find(idxMat);
[rows,cols] = ind2sub(size(idxMat),idx);
blocks = zeros(prod(blkSize),length(idx));
for i = 1:length(idx)
currBlock = I(rows(i):rows(i)+blkSize(1)-1,cols(i):cols(i)+blkSize(2)-1);
blocks(:,i) = currBlock(:);
End

子函数：OMPerr.m

function [A]=OMPerr(D,X,errorGoal);
%=============================================
% Sparse coding of a group of signals based on a given
% dictionary and specified number of atoms to use.
% input arguments: D - the dictionary
%                  X - the signals to represent
%                  errorGoal - the maximal allowed representation error for
%                  each siganl.
% output arguments: A - sparse coefficient matrix.
%=============================================
[n,P]=size(X);
[n,K]=size(D);
E2 = errorGoal^2*n;
maxNumCoef = n/2;%%%%%%32
A = sparse(size(D,2),size(X,2));%参考稀疏矩阵的帮助256*10000
for k=1:1:P,
a=[];
x=X(:,k);
residual=x;
indx = [];
a = [];
currResNorm2 = sum(residual.^2);
j = 0;

while currResNorm2>E2 & j < maxNumCoef,
j = j+1;
proj=D'*residual;%参考pinv函数的帮助
pos=find(abs(proj)==max(abs(proj)));%看看D（256列）中哪一列的值最大
pos=pos(1);
indx(j)=pos;%%%index的值为1到256
%c++的opm优化速度的算法     http://blog.csdn.net/pi9nc/article/details/26593003 a=pinv(D(:,indx(1:j)))*x;%j*64  *64*1=j*1
residual=x-D(:,indx(1:j))*a;
currResNorm2 = sum(residual.^2);
end;
if (length(indx)>0)
A(indx,k)=a;%%%a是j*1的矩阵,其中j=maxNumCoef
end
end;
return;

function [A]=OMPerr(D,X,errorGoal);

%=============================================

% Sparse coding of a group of signals based on a given

% dictionary and specified number of atoms to use.

% input arguments: D - the dictionary

% X - the signals to represent

% errorGoal - the maximal allowed representation error for

% each siganl.

% output arguments: A - sparse coefficient matrix.

%=============================================

[n,P]=size(X);

[n,K]=size(D);

E2 = errorGoal^2*n;

maxNumCoef = n/2;%%%%%%32

A = sparse(size(D,2),size(X,2));%参考稀疏矩阵的帮助256*10000

for k=1:1:P,

a=[];

x=X(:,k);

residual=x;

indx = [];

a = [];

currResNorm2 = sum(residual.^2);

j = 0;

while currResNorm2>E2 & j < maxNumCoef,

j = j+1;

proj=D'*residual;%参考pinv函数的帮助

pos=find(abs(proj)==max(abs(proj)));%看看D（256列）中哪一列的值最大

pos=pos(1);

indx(j)=pos;%%%index的值为1到256

%c++的opm优化速度的算法 http://blog.csdn.net/pi9nc/article/details/26593003
a=pinv(D(:,indx(1:j)))*x;%j*64 *64*1=j*1

residual=x-D(:,indx(1:j))*a;

currResNorm2 = sum(residual.^2);

end;

if (length(indx)>0)

A(indx,k)=a;%%%a是j*1的矩阵,其中j=maxNumCoef

end

end;

return;

运行结果如下：









整个代码中我个人觉得最难的部分就是那个OMP算法，我一直有一些疑问，为什么在这里用OMP算法，它到底的能解决什么样的问题？我尝试的根据csdn大牛的博客了解了一下算法的步骤，然后根据算法步骤，看懂了matlab代码。

参考资料：图像稀疏去噪算法的并行改进研究_孙黎明.