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

最新推荐文章于 2025-04-24 22:48:46 发布

happinesslz

最新推荐文章于 2025-04-24 22:48:46 发布

阅读量1.6w

点赞数 11

分类专栏：模式识别机器学习图像处理

本文链接：https://blog.youkuaiyun.com/ZHE123ZHE123ZHE123/article/details/49668557

版权

机器学习同时被 3 个专栏收录

18 篇文章

订阅专栏

图像处理

17 篇文章

订阅专栏

模式识别

9 篇文章

订阅专栏

基于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.youkuaiyun.com/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);