matlab程序是压缩感知图像去噪的程序

追梦1991322 2013-05-28 10:34:22

该程序是压缩感知图像去噪代码，可是在程序中有一个函数没有完成Compute_Error_Stats(currResIm , trueIm);我也是初学者，不知道该函数是干什么的，自己也不知道怎么写，函数里参数trueIm也不知道是什么，请各位帮帮忙啊，谢谢你们啊！！

function ResIm = Image_Denoising_Global_Denoising(Im , trueIm)
% Denoise image by denoising its patches, using a pre-determined
% dictionary.去噪图像通过去噪补丁，使用预先确定的字典。
% Patches are denoised WITHOUT overlap.补丁都去噪不重叠
%
% Inputs :
% Im : image to denoise, double, [0 255]IM：图像降噪，双，[0 255]
% trueIm : This is used (unfairly) to optimize the parameters for the algorithm
%trueIm：这是用来（不公平）的算法中的参数进行优化
% Outputs :
% ResIm : Result Image for the best parameter ResIm：结果图像的最佳参数
%%
%
% $$e^{\pi i} + 1 = 0$$
%

%% Parameters 参数n = numel(A) 返回数组A中元素个数。
imLen = numel(Im);%numel所属函数，在MATLAB中，该函数用于计算数组中满足指定条件的元素个数。

%% Create the dictionary - Haar创建的词典 - 哈尔
[HaarDict , atomNorms] = Generate_Haar_Matrix(size(Im) , 2);
atomNorms = atomNorms(:);
invAtomNorms = 1 ./ atomNorms;
% Show_Haar_Dict(HaarDict , size(Im));
nAtoms = size(HaarDict , 2); % = 1 + 3 * nLevels

%% Prepare the projections of the image onto the dictionary准备图像预测到词典
projs = HaarDict' * Im(:);
absProjs = abs(projs);

%% Test for different values of T测试不同的T值
Tvalues = [0 : 200];
resPSNRs = zeros(1 , length(Tvalues));
fprintf('%d T values: ' , length(Tvalues));
for Tind = 1 : length(Tvalues)
%length是求某一个矩阵或者向量的长度。ones(length(t))指的是生成一个
% length(t)*length(t)全是1的矩阵。z=z=0*ones(length(t))
%指的是在执行完上面的语句后，将矩阵中的数都乘以0，因此，z是全为0的矩阵。

if mod(Tind , 10) == 0, fprintf('%d ' , Tind); end;

% Current threshold 电流阈值
T = Tvalues(Tind);

% Compute the hard-thresholding result计算硬阈值的结果
resProjs = projs .* (absProjs > (T .* invAtomNorms));

% Compute the result signal by multiplying by the
% dictionary计算结果的信号乘以由字典
currResIm = HaarDict * resProjs;
currResIm = reshape(currResIm , size(Im));

% Compute error result计算错误结果
resPSNRs(Tind) = Compute_Error_Stats(currResIm , trueIm);

end
fprintf('\n');

%% Find the best value and plot the results 寻找最好的价值，并绘制出结果
[~, maxPSNRind] = max(resPSNRs);
bestT = Tvalues(maxPSNRind);
if 1
figure;
plot(Tvalues , resPSNRs , '+r');
title('PSNR as a function of T');
end

%% Reconstruct the image for the best value of T重构图像的最佳T值
resProjs = projs .* (absProjs > (bestT .* invAtomNorms));
ResIm = HaarDict * resProjs;
ResIm = reshape(ResIm , size(Im));

%%
return;

% [n,m]=size(Haar);
% W=[0.25*ones(n*4,1); 0.5*ones(n*3,1)];

function [HaarDict , atomNorms] = Generate_Haar_Matrix(imSize , nLevels)

atomNorms = [];

% Prepare the first level bands 准备第一级带
basicH = [1 -1];
basicL = [1 1];
HH1 = kron(basicH' , basicH);
%得出结果HH1 =

% 1 -1
% -1 1
%矩阵的Kronecker乘法(kronecker product)对n×m阶矩阵A和p×q阶矩阵B，
%A和B的Kronecher乘法运算可定义为：
%由上面的式子可以看出，Kronecker乘积A B表示矩阵A的所有元素与B之间的乘积组合而成的较大的矩阵，
%B A则完全类似．A B和B A均为np×mq矩阵，但一般情况下A B B A．和普通矩阵的乘法不同，Kronecker
%乘法并不要求两个被乘矩阵满足任何维数匹配方面的要求
%．Kronecker乘法的Matlab命令为 C=kron(A,B)，例如给定两个矩阵A和B：
%则由以下命令可以求出A和B的Kronecker乘积C：
%A=[1 2; 3 4]; B=[1 3 2; 2 4 6]; C=kron(A,B)
%C =
%1 3 2 2 6 4
%2 4 6 4 8 12
%3 9 6 4 12 18
%6 12 18 8 16 24
%作为比较，可以计算B和A的Kronecker乘积D，可以看出C、D是不同的：
%A=[1 2; 3 4]; B=[1 3 2; 2 4 6]; D=kron(B,A)
%D =
%1 2 3 6 2 4
%3 4 9 12 6 8
%2 4 4 8 6 12
%6 8 12 16 18 24
%后者矩阵分别于前者矩阵每个元素相乘

HaarDict = Construct_One_Haar_Band(imSize , HH1);
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(HH1(:).^2))];

LH1 = kron(basicL' , basicH);
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , LH1)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(LH1(:).^2))];

HL1 = kron(basicH' , basicL);
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , HL1)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(HL1(:).^2))];

LL1 = kron(basicL' , basicL);
if nLevels == 1,
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , LL1)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(LL1(:).^2))];
end;

% Create the rest of the bands from the LL of the previous level and the
% basic ones创建其余的频段从以前的水平的LL和基本的
currLL = LL1;
for cLevel = 2 : nLevels
currHH = kron(HH1 , currLL);
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currHH)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currHH(:).^2))];

currLH = kron(LH1 , currLL);
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currLH)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currLH(:).^2))];

currHL = kron(HL1 , currLL);
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currHL)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currHL(:).^2))];%prod计算数组元素的连乘积

currLL = kron(LL1 , currLL);
if cLevel == nLevels,
HaarDict = [HaarDict Construct_One_Haar_Band(imSize , currLL)];
atomNorms = [atomNorms ones(1 , prod(imSize)) * sqrt(sum(currLL(:).^2))];
end;

end

function haarBand = Construct_One_Haar_Band(imSize , haarMatrix)

% Number of pixels 的像素数
nPix = prod(imSize); %计算数组元素的连乘积。prod（[1:5]）返回120
[c,r] = meshgrid(1:imSize(2) , 1:imSize(1));
%meshgrid是MATLAB中用于生成网格采样点的函数。在使用MATLAB进行3-D图形绘制方面有着广泛的应用。

c = c(:)'; r = r(:)';
nOnes = numel(haarMatrix);

% Prepare all the indices we will need to insert into the sparse
% matrix准备好所有的指标，我们将需要插入的稀疏矩阵
sparseMatRowInd = zeros(numel(haarMatrix) , nPix);
spraseMatColInd = repmat([1:nPix] , numel(haarMatrix) , 1 );
%repmat 即 Replicate Matrix ，复制和平铺矩阵，是 MATLAB 里面的一个函数。B = repmat(A,m,n)
%将矩阵 A 复制 m×n 块，即把 A 作为 B 的元素，B 由 m×n 个 A 平铺而成。B 的维数是 [size(A,1)*m, (size(A,2)*n] 。

sparseMatVal = zeros(numel(haarMatrix) , nPix);

% Run on the haarMatrix 运行在haarMatrix上
insertInd = 0;
for dc = 0 : size(haarMatrix , 2)-1
for dr = 0 : size(haarMatrix , 1)-1
insertInd = insertInd + 1;
sparseMatRowInd(insertInd , :) = sub2ind(imSize , mod((r + dr) - 1 , imSize(1)) + 1 , mod((c + dc) - 1 , imSize(2)) + 1);
sparseMatVal(insertInd , :) = haarMatrix(dr + 1 , dc + 1);
end
end
sparseMatVal = sparseMatVal / sum(abs(haarMatrix(:)),1);

% % Create the sparse matrix itself 创建稀疏矩阵本身
haarBand = sparse(sparseMatRowInd(:) , spraseMatColInd(:) , sparseMatVal(:) , nPix , nPix);

function Show_Haar_Dict(D , imSize)
figure;
multSize = max(min(floor(600 ./ imSize)) , 1);%[Y,I]=max(M,[],2), 在第2维方向上取最大值，也就是每行最大值，结果存在Y里，I里存的是每行最大值的列位置。
for c1 = 1 : size(D , 2)
currAtom = reshape(full(D(:,c1)) , imSize);%函数功能：重新调整矩阵的行数、列数、维数。在matlab命令窗口中键入doc reshape或help reshape即可获得该函数的帮助信息。
%B = reshape(A,m,n)返回一个m*n的矩阵B， B中元素是按列从A中得到的。如果A中元素个数不等于m*n，则会引发错误。

currAtom2 = imresize(currAtom , multSize , 'nearest');%该函数用于对图像做缩放处理。在matlab的命令窗口中输入doc imresize或者help imresize即可获得该函数的帮助信息。
%B = imresize(A, m)返回的图像B的长宽是图像A的长宽的m倍，即缩放图像。 m大于1，则放大图像； m小于1，缩小图像。B = imresize(A, [numrows numcols])
%numrows和numcols分别指定目标图像的高度和宽度。显而易见，由于这种格式允许图像缩放后长宽比例和源图像长宽比例不相同，因此所产生的图像有可能发生畸变。
%[Y newmap] = imresize(X, map, scale)[...] = imresize(..., method)method参数用于指定在改变图像尺寸时所使用的算法，可以为一下几种：'nearest'：这个参数也是默认的，即改变图像尺寸时采用最近邻插值算法；
%'bilinear'：采用双线性插值算法；'bicubic'：采用双三次插值算法；
cla;
imshow(currAtom2 , [-1 1]);
title(sprintf('Atom %d' , c1));
drawnow;
%drawnow用于刷新屏幕的。
pause(0.05);

end

...全文