【图像去噪】基于稀疏表示KSVD实现图像去噪matlab代码

1 简介

图像信号的稀疏表示是要将图像中复杂的信息用尽可能少的基本信号的线性组合来表示原始信号，这些基本信号被叫做原子，而这些原子组成了字典。通过这种变换使得模型复杂度降低，让原始信号更容易处理。由于图像信号一般是复杂的，因此，通过稀疏表示简化处理过程，使得稀疏表示技术在图像处理领域应用广泛。最初用于信号表示的字典出现在 20 世纪 80 年代左右，研究表明，通过傅立叶变换或小波变换获得的字典具有稀疏性，但是获得的字典简单，对细节丰富的图像效果较差，容易出现对信号低频特征保持良好，而高频特征保持效果较差的问题。近年来，基于冗余字典的稀疏表示方法受到众多学者关注，人们发现利用过完备字典能够获得更稀疏的表示系数，恢复的图像与原始图像有更好的结构特征匹配度，且这种字典具有自适应的特性。信号的稀疏表示过程主要是通过稀疏编码和字典学习完成，先把目标信号投影到一组基空间中，获得信号投影到这组非正交基上的系数的过程就是稀疏编码，这组系数就是稀疏系数。通过稀疏系数和过完备字典可以重构信号，达到图像去噪的目的

2 部分代码

function [Dictionary,output] = KSVD_NN(...
   Data,... % an nXN matrix that contins N signals (Y), each of dimension n.
   param)
% =========================================================================
%                       Non Negative K-SVD algorithm
% =========================================================================
% The NN-K-SVD algorithm finds a non-negative dictionary for linear representation of
% signals. Given a set of signals, it searches for the best dictionary that
% can sparsely represent each signal. Detailed discussion on the algorithm
% and possible applications can be found in "K-SVD and its non-negative 
% variant for dictionary design", written by M. Aharon, M. Elad, and A.M. Bruckstein 
% and appeared in the Proceedings of the SPIE conference wavelets, Vol.
% 5914, July 2005. 
% =========================================================================
% INPUT ARGUMENTS:
% Data                         an nXN matrix that contins N signals (Y), each of dimension n. 
% param                       structure that includes all required
%                                 parameters for the K-SVD execution.
%                                 Required fields are:
%   K, ...                   the number of dictionary elements to train
%   numIteration,...         number of iterations to perform.
%   L,...                     maximum coefficients to use in OMP coefficient calculations.
%   InitializationMethod,... mehtod to initialize the dictionary, can
%                                 be one of the following arguments: 
%                                 * 'DataElements' (initialization by the signals themselves), or: 
%                                 * 'GivenMatrix' (initialization by a given matrix param.initialDictionary).
%   (optional, see InitializationMethod) initialDictionary,...     % if the initialization method 
%                                 is 'GivenMatrix', this is the matrix that will be used.
%   (optional) TrueDictionary, ...       % if specified, in each
%                                 iteration the difference between this dictionary and the trained one
%                                 is measured and displayed.
%   displayProgress, ...     if =1 progress information is displyed. If param.errorFlag==0, 
%                                 the average repersentation error (RMSE) is displayed, while if 
%                                 param.errorFlag==1, the average number of required coefficients for 
%                                 representation of each signal is displayed.
% =========================================================================
% OUTPUT ARGUMENTS:
% Dictionary                 The extracted dictionary of size nX(param.K).
% output                     Struct that contains information about the current run. It may include the following fields:
%   CoefMatrix                 The final coefficients matrix (it should hold that Data equals approximately Dictionary*output.CoefMatrix.
%   ratio                       If the true dictionary was defined (in
%                               synthetic experiments), this parameter holds a vector of length
%                               param.numIteration that includes the detection ratios in each
%                               iteration).
%   totalerr                   The total representation error after each
%                               iteration (defined only if
%                               param.displayProgress=1)
% =========================================================================

if (~isfield(param,'displayProgress'))
   param.displayProgress = 0;
end
totalerr(1) = 99999;
if (isfield(param,'errorFlag')==0)
   param.errorFlag = 0;
end

%Data(Data<0) = 0;

if (isfield(param,'TrueDictionary'))
   displayErrorWithTrueDictionary = 1;
   ErrorBetweenDictionaries = zeros(param.numIteration+1,1);
   ratio = zeros(param.numIteration+1,1);
else
   displayErrorWithTrueDictionary = 0;
ratio = 0;
end
if (param.preserveDCAtom>0)
   FixedDictionaryElement(:,1) = 1/sqrt(size(Data,1));
else
   FixedDictionaryElement = [];
end
% coefficient calculation method is OMP with fixed number of coefficients

if (size(Data,2) < param.K)
   disp('Size of data is smaller than the dictionary size. Trivial solution...');
   Dictionary = Data(:,1:size(Data,2));
   return;
elseif (strcmp(param.InitializationMethod,'DataElements'))
   Dictionary(:,1:param.K) = Data(:,1:param.K-param.preserveDCAtom);
elseif (strcmp(param.InitializationMethod,'GivenMatrix'))
   Dictionary = param.initialDictionary(:,1:param.K-param.preserveDCAtom);
end
% reduce the components in Dictionary that are spanned by the fixed
% elements
if (param.preserveDCAtom)
   tmpMat = FixedDictionaryElement \ Dictionary;
   Dictionary = Dictionary - FixedDictionaryElement*tmpMat;
end
%normalize the dictionary.
Dictionary = Dictionary*diag(1./sqrt(sum(Dictionary.*Dictionary)));
totalErr = zeros(1,param.numIteration);

CoefMatrix = sparse(param.K,size(Data,2));

% the K-SVD algorithm starts here.
for iterNum = 1:param.numIteration
   % find the coefficients
CoefMatrix = NN_BP(Data, [FixedDictionaryElement,Dictionary],param.L,CoefMatrix);
   
   replacedVectorCounter = 0;
rPerm = randperm(size(Dictionary,2));
   for j = rPerm
      [betterDictionaryElement,CoefMatrix,addedNewVector] = I_findBetterDictionaryElement(Data,...
          [FixedDictionaryElement,Dictionary],j+size(FixedDictionaryElement,2),...
           CoefMatrix ,param.L);
       Dictionary(:,j) = betterDictionaryElement;
       if (param.preserveDCAtom)
           tmpCoef = FixedDictionaryElement\betterDictionaryElement;
           Dictionary(:,j) = betterDictionaryElement - FixedDictionaryElement*tmpCoef;
           Dictionary(:,j) = Dictionary(:,j)./sqrt(Dictionary(:,j)'*Dictionary(:,j));
       end
       replacedVectorCounter = replacedVectorCounter+addedNewVector;
   end

   if (iterNum>1 & param.displayProgress)
       if (param.errorFlag==0)
           output.totalerr(iterNum-1) = sqrt(sum(sum((Data-[FixedDictionaryElement,Dictionary]*CoefMatrix).^2))/prod(size(Data)));
           disp(['Iteration   ',num2str(iterNum),'   Total error is: ',num2str(output.totalerr(iterNum-1))]);
       else
           output.numCoef(iterNum-1) = length(find(CoefMatrix))/size(Data,2);
           disp(['Iteration   ',num2str(iterNum),'   Average number of coefficients: ',num2str(output.numCoef(iterNum-1))]);
       end
   end
   if (displayErrorWithTrueDictionary ) 
      [ratio(iterNum+1),ErrorBetweenDictionaries(iterNum+1)] = I_findDistanseBetweenDictionaries(param.TrueDictionary,Dictionary);
       disp(strcat(['Iteration ', num2str(iterNum),' ratio of restored elements: ',num2str(ratio(iterNum+1))]));
       output.ratio = ratio;
   end
   Dictionary = I_clearDictionary(Dictionary,CoefMatrix,Data);
   
   if (isfield(param,'waitBarHandle'))
       waitbar(iterNum/param.counterForWaitBar);
   end
end

output.CoefMatrix = CoefMatrix;
Dictionary = [FixedDictionaryElement,Dictionary];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% findBetterDictionaryElement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [betterDictionaryElement,CoefMatrix,NewVectorAdded] = I_findBetterDictionaryElement(Data,Dictionary,j,CoefMatrix,numCoefUsed)
if (length(who('numCoefUsed'))==0)
   numCoefUsed = 1;
end
relevantDataIndices = find(CoefMatrix(j,:)); % the data indices that uses the j'th dictionary element.
if (length(relevantDataIndices)<1) %(length(relevantDataIndices)==0)
   ErrorMat = Data-Dictionary*CoefMatrix;
   ErrorNormVec = sum(ErrorMat.^2);
  [d,i] = max(ErrorNormVec);
   betterDictionaryElement = Data(:,i);%ErrorMat(:,i); %
   betterDictionaryElement = betterDictionaryElement./sqrt(betterDictionaryElement'*betterDictionaryElement);
   betterDictionaryElement = betterDictionaryElement.*sign(betterDictionaryElement(1));
   CoefMatrix(j,:) = 0;
   NewVectorAdded = 1;
   return;
end

NewVectorAdded = 0;
reduced_coeff = CoefMatrix(:, relevantDataIndices);
reduced_Data = Data (:, relevantDataIndices);
saveDebugDict = Dictionary(:,j);
saveDebugCoef = reduced_coeff(j,:);
%debug1 = sum(sum((reduced_Data - Dictionary*reduced_coeff).^2));
reduced_coeff (j, :) = 0;    % all but the j-th element

err_mat = reduced_Data - Dictionary * reduced_coeff;
[U S V flag] = svds((err_mat), 1);


% check for sign, flip U and V's sign if negative.
Idx_U = find(U<0);
Idx_V = find(V<0);

u1 = U; u1(Idx_U) = 0; v1 = V; v1(Idx_V) = 0;approx1 = norm(err_mat- u1*v1'*S);
u1 = zeros(size(U)); u1(Idx_U) = -U(Idx_U); v1 = zeros(size(V)); v1(Idx_V) = -V(Idx_V);approx2 = norm(err_mat- u1*v1'*S);
if (approx1<= approx2)
betterDictionaryElement = U;
betterDictionaryElement(Idx_U) = 0;
coefs = V;
coefs(Idx_V) = 0;
else
betterDictionaryElement = zeros(size(U));
betterDictionaryElement(Idx_U) = -U(Idx_U);
coefs = zeros(size(V));
coefs(Idx_V) = -V(Idx_V);
end

newAtomNorm = sqrt(betterDictionaryElement'*betterDictionaryElement);
betterDictionaryElement = betterDictionaryElement/newAtomNorm;
coefs = coefs * newAtomNorm;
% coefs(coefs<0) = 0;

newE = sum(sum(((reduced_Data - Dictionary(:,[1:j-1,j+1:end])*reduced_coeff([1:j-1,j+1:end],:))-betterDictionaryElement*coefs').^2));
oldE = sum(sum(((reduced_Data - Dictionary(:,[1:j-1,j+1:end])*reduced_coeff([1:j-1,j+1:end],:))-saveDebugDict*saveDebugCoef).^2));
if (newE>oldE)
for iter = 1:30 % the number of iterations
betterDictionaryElement = err_mat*coefs/(coefs'*coefs);
betterDictionaryElement(betterDictionaryElement<0) = 0;
coefs = err_mat'*betterDictionaryElement/(betterDictionaryElement'*betterDictionaryElement);
coefs(coefs<0) = 0;
end
newAtomNorm = sqrt(betterDictionaryElement'*betterDictionaryElement);
betterDictionaryElement = betterDictionaryElement/newAtomNorm;
coefs = coefs * newAtomNorm;
reduced_coeff(j,:) =  coefs;
end
reduced_coeff(j,:) =  coefs;
CoefMatrix (:, relevantDataIndices) =reduced_coeff;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% findDistanseBetweenDictionaries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ratio,totalDistances] = I_findDistanseBetweenDictionaries(original,new)
% first, all the column in oiginal starts with positive values.
catchCounter = 0;
totalDistances = 0;
for i = 1:size(new,2)
   new(:,i) = new(:,i);
end
for i = 1:size(original,2)
   d = original(:,i);
   distances =sum ( (new-repmat(d,1,size(new,2))).^2);
  [minValue,index] = min(distances);
   errorOfElement = 1-abs(new(:,index)'*d);
   totalDistances = totalDistances+errorOfElement;
   catchCounter = catchCounter+(errorOfElement<0.01);
end
ratio = 100*catchCounter/size(original,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I_clearDictionary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Dictionary = I_clearDictionary(Dictionary,CoefMatrix,Data)
T2 = 0.999;
T1 = 3;
K=size(Dictionary,2);
Er=sum((Data-Dictionary*CoefMatrix).^2,1); % remove identical atoms
G=Dictionary'*Dictionary; G = G-diag(diag(G));
for jj=1:1:K,
   if max(G(jj,:))>T2 | length(find(abs(CoefMatrix(jj,:))>1e-7))<=T1 ,
      [val,pos]=max(Er);
       Er(pos(1))=0;
       Dictionary(:,jj)=Data(:,pos(1))/norm(Data(:,pos(1)));
       G=Dictionary'*Dictionary; G = G-diag(diag(G));
   end;
end;

3 仿真结果

4 参考文献

[1]沈琦, and 关云霞. "基于稀疏表示的图像去噪方法.", CN111242873A. 2020.