www.pudn.com > denoiseImageKSVD.zip > denoiseImageKSVD.m, change:2015-12-23,size:10900b


function [IOut,output] = denoiseImageKSVD(Image,sigma,K,varargin)%[字典 稀疏系数] (加噪的图像 字典大小 原子个数 可变参数列表) 
%========================================================================== 
%   P E R F O R M   D E N O I S I N G   U S I N G   A  D I C T  I O N A R Y 
%                  T R A I N E D   O N   N O I S Y   I M A G E 
%   perform denoising using a dictionary trained on noisy image  
%========================================================================== 
% function IOut = denoiseImageKSVD(Image,sigma,K,varargin)  
% denoise an image by sparsely representing each block with the 稀疏表示去噪 每一块通过过完备的训练字典 
% already overcomplete trained Dictionary, and averaging the represented parts. 使被表示的部分 平均 
% Detailed description can be found 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). 
% This function may take some time to process. Possible factor that effect 
% the processing time are: 这个功能需要运行时间 可能影响运行时间的因素如下所示 
%  1. number of KSVD iterations - the default number of iterations is 10. 迭代次数10 
%  However, fewer iterations may, in most cases, result an acceleration in 
%  用尽可能少的迭代次数,使得这个处理过程花的时间很少,但是对结果的影响又不能太大 
%  the process, without effecting  the result too much. Therefore, when 
%  required, this parameter may be re-set.这个参数通常会被重置 
%  2. max Blocks To Consider - The maximal number of blocks to train 
%  on.训练样本的最大个数 
%  If this  number is larger the number of blocks in the image, random blocks 
%  from the image will be selected for training.如果大于待处理图像块,则从待处理图片中随机取块 进行训练  
% =================================================================== 
% INPUT ARGUMENTS(输入元素) : Image - the noisy image (gray-level scale) 噪声图像 灰度等级 
%                   sigma - the s.d. of the noise (assume to be white Gaussian).standard deviation 噪声的标准差 方差的平方根 
%                   K - the number of atoms in the trained dictionary. 训练字典中原子的数目 
%    Optional arguments:可选参数       
%                       'blockSize' - the size of the blocks the algorithm 算法 
%                        块的大小 
%                       works. All blocks are squares,所有的块都是方形的 therefore the given 
%                       parameter should be one number (width or 
%                       height).给定的参数应该是一个数字 长度和宽度应该是一样的 
%                         
%                       Default value: 8. 默认值  
%                       'errorFactor' - a factor that multiplies sigma in order 乘以sigma 
%                       to set the allowed representation error. In the 设定允许的表示误差 
%                       experiments presented in the paper, it was set to 1.15 
%                       (which is also the default  value here). 
%                  'maxBlocksToConsider' - maximal number of blocks that 
%                       can be processed. This number is dependent on the memory 
%                       capabilities of the machine, and performances? 机器性能  
%                       considerations. If the number of available blocks in the 
%                       image is larger than 'maxBlocksToConsider', the sliding 
%                       distance between the blocks increases. The default value变化的长度 
%                       is: 250000.   
%                  'slidingFactor' - the sliding distance between processed 
%                       blocks. Default value is 1. However, if the image is 
%                       large, this number increases automatically (because of 
%                       memory requirements). Larger values result faster 
%                       performances (because of fewer processed blocks). 
%                  'numKSVDIters' - the number of KSVD iterations processed对噪声图像块进行处理的迭代次数 
%                       blocks from the noisy image. If the number of 
%                       blocks in the image is larger than this number,如果图像块的数目大于迭代次数 
%                       random blocks from all available blocks will be 
%                       selected. The default value for this parameter is: 
%                       10 if sigma > 5, and 5 otherwise. 
%                  'maxNumBlocksToTrainOn' - the maximal number of blocks训练样本的最大数目 
%                       to train on. The default value for this parameter is 
%                       65000. However, it might not be enough for very large 
%                       images. 
%                  'displayFlag' - if this flag is switched on,这个标志被转换 
%                       announcement after finishing each iteration will appear,在完成每次迭代都会出现 
%                       as also a measure concerning the progress of the也作为算法的进展的标志 
%                       algorithm (the average number of required coefficients表示的要求系数的平均值 
%                       for representation). The default value is 1 (on). 
%                  'waitBarOn' - can be set to either 1 or 0. If 
%                       waitBarOn==1 a waitbar, presenting the progress of the 
%                       algorithm will be displayed.显示等待 将被延时 
% OUTPUT ARGUMENTS 输出内容: Iout - a 2-dimensional array in the same size of the二维数组 和输入图像的大小一致 
%                       input image, that contains the cleaned image.包括去噪后的图片 
%                       output.D - the trained dictionary.训练字典 
% ========================================================================= 
  
% first, train a dictionary on the noisy image 对噪声图片进行训练得到字典 
  
reduceDC = 1; 
[NN1,NN2] = size(Image); 
waitBarOn = 1;%进度为1 代表这个算法的进程会被延时 
if (sigma > 5)%sigma=50   numIterOfKsvd = 10; 
    numIterOfKsvd = 10;%迭代次数 
else 
    numIterOfKsvd = 5; 
end 
C = 1.15; 
maxBlocksToConsider = 260000;%可以处理的的最大的块的数目 
slidingDis = 1;%变化距离 
bb = 8; 
maxNumBlocksToTrainOn = 65000;%训练样本的最大数目 
displayFlag = 1; 
hh=length(varargin)%测试一下能不能进入下面的for循环中去。 可变参数 
% for argI = 1:2:length(varargin) 
%     if (strcmp(varargin{argI}, 'slidingFactor')) 
%         slidingDis = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'errorFactor')) 
%         C = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'maxBlocksToConsider')) 
%         maxBlocksToConsider = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'numKSVDIters')) 
%         numIterOfKsvd = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'blockSize')) 
%         bb = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'maxNumBlocksToTrainOn')) 
%         maxNumBlocksToTrainOn = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'displayFlag')) 
%         displayFlag = varargin{argI+1}; 
%     end 
%     if (strcmp(varargin{argI}, 'waitBarOn')) 
%         waitBarOn = varargin{argI+1}; 
%     end 
% end 
  
if (sigma <= 5) 
    numIterOfKsvd = 5;%迭代次数 
end 
  
% first, train a dictionary on blocks from the noisy image 
  
if(prod([NN1,NN2]-bb+1)> maxNumBlocksToTrainOn) 
    randPermutation =  randperm(prod([NN1,NN2]-bb+1));%随机排列 
    selectedBlocks = randPermutation(1:maxNumBlocksToTrainOn); 
  
    blkMatrix = zeros(bb^2,maxNumBlocksToTrainOn); 
    for i = 1:maxNumBlocksToTrainOn 
        [row,col] = ind2sub(size(Image)-bb+1,selectedBlocks(i)); 
        currBlock = Image(row:row+bb-1,col:col+bb-1); 
        blkMatrix(:,i) = currBlock(:); 
    end 
else 
    blkMatrix = im2col(Image,[bb,bb],'sliding');%8*8=64   所以blkMatrix矩阵大小为:64*[(NN1-bb+1)*(NN2-bb+1)] 
end 
  
param.K = K;%K=256 4*8*8=256 训练字典中 原子的数量 
param.numIteration = numIterOfKsvd ;%sigma=50   所以numIterOfKsvd = 10; 
  
param.errorFlag = 1; % decompose signals稀疏分解信号 until a certain error is reached.直到一个确定的错误出现  
param.errorGoal = sigma*C;% do not use fix number of coefficients.不要使用固定的系数 
param.preserveDCAtom = 0; 
  
Pn=ceil(sqrt(K));%%Pn=16 
DCT=zeros(bb,Pn);%%bb=8 
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); 
end; 
DCT=kron(DCT,DCT);%%%%%跟DCT中的代码一样的     64*256的矩阵 
  
param.initialDictionary = DCT(:,1:param.K );%%%% 取了256列。也就是全部都取了 
param.InitializationMethod =  'GivenMatrix'; 
  
if (reduceDC)%%reduceDC=1 
    vecOfMeans = mean(blkMatrix); 
    blkMatrix = blkMatrix-ones(size(blkMatrix,1),1)*vecOfMeans;%%%减去平均数  blkMatrix矩阵大小为:64*[(NN1-bb+1)*(NN2-bb+1)] 
end 
  
if (waitBarOn)%waitBarOn=1 
    counterForWaitBar = param.numIteration+1;%param.numIteration = numIterOfKsvd ;  =10 
    h = waitbar(0,'Denoising In Process ...'); 
    param.waitBarHandle = h; 
    param.counterForWaitBar = counterForWaitBar; 
end 
  
  
param.displayProgress = displayFlag;%displayFlag = 1; 
[Dictionary,output] = KSVD(blkMatrix,param);%%%%%%%最核心的函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
output.D = Dictionary; 
  
if (displayFlag)%displayFlag = 1; 
    disp('finished Trainning dictionary'); 
end 
  
  
  
% denoise the image using the resulted dictionary 
errT = sigma*C; 
IMout=zeros(NN1,NN2); 
Weight=zeros(NN1,NN2); 
%blocks = im2col(Image,[NN1,NN2],[bb,bb],'sliding'); 
while (prod(floor((size(Image)-bb)/slidingDis)+1)>maxBlocksToConsider) 
    slidingDis = slidingDis+1; 
end 
[blocks,idx] = my_im2col(Image,[bb,bb],slidingDis); 
  
if (waitBarOn) 
    newCounterForWaitBar = (param.numIteration+1)*size(blocks,2); 
end 
  
  
% go with jumps of 30000 
for jj = 1:30000:size(blocks,2) 
    if (waitBarOn) 
        waitbar(((param.numIteration*size(blocks,2))+jj)/newCounterForWaitBar); 
    end 
    jumpSize = min(jj+30000-1,size(blocks,2)); 
    if (reduceDC) 
        vecOfMeans = mean(blocks(:,jj:jumpSize)); 
        blocks(:,jj:jumpSize) = blocks(:,jj:jumpSize) - repmat(vecOfMeans,size(blocks,1),1); 
    end 
     
    %Coefs = mexOMPerrIterative(blocks(:,jj:jumpSize),Dictionary,errT); 
    Coefs = OMPerr(Dictionary,blocks(:,jj:jumpSize),errT); 
    if (reduceDC) 
        blocks(:,jj:jumpSize)= Dictionary*Coefs + ones(size(blocks,1),1) * vecOfMeans; 
    else 
        blocks(:,jj:jumpSize)= Dictionary*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]); 
    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); 
    count = count+1; 
end; 
  
if (waitBarOn) 
    close(h); 
end 
IOut = (Image+0.034*sigma*IMout)./(1+0.034*sigma*Weight);