附录A
clear
clc
F=checkerboard(2); %原图像
D=imnoise(F,'gaussian',0,0.02); %退化图像=参考信号
nhood=[3 3];
% Estimate the local mean of f.
localMean = filter2(ones(nhood), D) / prod(nhood);
H=D-localMean;％仿照wiener2函数里的求取类似的“均值”
[h k]=size(D);
L=h*k;
D=reshape(D,L,1); %将图像矩阵变为矢量形式
f=zeros(size(D)); %系统输出＝误差信号
y=f; %噪声逼近信号初始化
E=f; %声明误差矢量
%设置输入噪声信号＝基本输入信号
%X=randn(h,k);
X=zeros(h,k);
X=imnoise(X,'gaussian',0,0.02);
X=reshape(X,L,1);
W=zeros(L,L); %设置权矢量初值
lmsMSE=f;
a=0;
%核心算法
u=0.000005;
for i=1:L
for n=1:L
for m=1:i
if n-m+1>0
y(n)=y(n)+W(m,n)*X(n-m+1);％滤波器输出
end
end
end
E=D-y;％计算误差
a=a+1;
lmsMSE(a)=mean(E.^2);％根据目标函数计算MSE
for n=1:L-1
W(i,n+1)=W(i,n)+u*E(n)*X(n);％更新权值
end
end
a=linspace(1,L,L);%设置画MSE变化曲线的横坐标
f=reshape(E,h,k)+H;％重构图像矩阵
figure,subplot(2,2,1),imshow(F);
title('原图像')
D=reshape(D,h,k);
subplot(2,2,2),imshow(D);
title('退化图像')
myMSE=mean2((F-f).^2)
subplot(2,2,3),imshow(f);
title('复原图像')
subplot(2,2,4)
plot(a,lmsMSE,'r-');
title('当u＝0.000001时的 MSE 变化曲线')
xlabel('迭代次数'),ylabel('MSE');
legend('lmsMSE(a)',0);
附录B
function [f,noise] = wiener2(varargin)
%WIENER2 Perform 2-D adaptive noise-removal filtering.
% WIENER2 lowpass filters an intensity image that has been
% degraded by constant power additive noise. WIENER2 uses a
% pixel-wise adaptive Wiener method based on statistics
% estimated from a local neighborhood of each pixel.
%
% J = WIENER2(I,[M N],NOISE) filters the image I using
% pixel-wise adaptive Wiener filtering, using neighborhoods of
% size M-by-N to estimate the local image mean and standard
% deviation. If you omit the [M N] argument, M and N default to
% 3. The additive noise (Gaussian white noise) power is assumed
% to be NOISE.
%
% [J,NOISE] = WIENER2(I,[M N]) also estimates the additive
% noise power before doing the filtering. WIENER2 returns this
% estimate as NOISE.
%
% Class Support
% -------------
% The input image I can be of class uint8, uint16, or double.
% The output image J is of the same class as I.
%
% Example
% -------
% I = imread('saturn.tif');
% J = imnoise(I,'gaussian',0,0.005);
% K = wiener2(J,[5 5]);
% imshow(J), figure, imshow(K)
%
% See also FILTER2, MEDFILT2.
% The following syntax is grandfathered:
%
% J = WIENER2(I,[M N],[MBLOCK NBLOCK],NOISE) or [J,NOISE] =
% WIENER2(I,[M N],[MBLOCK NBLOCK]) processes the intensity
% image I as above but in blocks of size MBLOCK-by-NBLOCK. Use
% J = WIENER2(I,[M N],SIZE(I),NOISE) to process the matrix all
% at once.
% Copyright 1993-2002 The MathWorks, Inc.
% $Revision: 5.18 $ $Date: 2002/03/15 15:29:28 $
% Uses algorithm developed by Lee (1980).
% Reference: "Two-Dimensional Signal and Image Processing" by
% Jae S. Lim, pp.536-540.
[g, nhood, block, noise, msg] = ParseInputs(varargin{:});%调用ParseInput()处理不同的输入参数的情况,且将输入参数以数组的形式表示
if (~isempty(msg))
error(msg);
end
classin = class(g);
classChanged = 0;
if ~isa(g, 'double')
classChanged = 1;
g = im2double(g);
end
% Estimate the local mean of f.
localMean = filter2(ones(nhood), g) / prod(nhood);
%用一个所有的值全相等的小矩阵中的二维FIR滤波器对二维数据g进行滤波
%其实质是利用小矩阵对g进行二维卷积运算，最终结果的维数大小
%等于g的维数，取的是卷积结果的中心部分。类似于求均值
% Estimate of the local variance of f.
localVar = filter2(ones(nhood), g.^2) / prod(nhood) - localMean.^2;
%求取类似于方差的localVar。这样的算法正好体现了图像矩阵中 %元素间的相关性
% Estimate the noise power if necessary.
if (isempty(noise))
noise = mean2(localVar);
%若没有输入noise，则取localVar的均值作为noise
end
% Compute result
% f = localMean + (max(0, localVar - noise) ./ ...
% max(localVar, noise)) .* (g - localMean);
%
% Computation is split up to minimize use of memory
% for temp arrays.
f = g - localMean;
%去除图像信号中的直流分量，使得图像信号成为零均值的信号
g = localVar - noise;
g = max(g, 0);
localVar = max(localVar, noise);
f = f ./ localVar;
f = f .* g;
f = f + localMean;
%对处理后的信号再重新加上直流分量
if classChanged,
f = changeclass(classin, f);
%使图像矩阵的数据类型与输入时相同
end
%%%
%%% Subfunction ParseInputs
%%%
function [g, nhood, block, noise, msg] = ParseInputs(varargin)
g = [];
nhood = [3 3];
block = [];
noise = [];
msg = '';
switch nargin%表示输入参数个数的系统参数
case 0
msg = 'Too few input arguments.';
return;
case 1
% wiener2(I)
g = varargin{1};%把输入的第一个参数赋值给g
case 2
g = varargin{1};
switch prod(size(varargin{2}))
case 1
% wiener2(I,noise)
noise = varargin{2};
case 2
% wiener2(I,[m n])
nhood = varargin{2};
otherwise
msg = 'Invalid input syntax';
return;
end
case 3
g = varargin{1};
if (prod(size(varargin{3})) == 2)
% wiener2(I,[m n],[mblock nblock]) OBSOLETE
warning(['WIENER2(I,[m n],[mblock nblock]) is an obsolete syntax.',...
'Omit the block size, the image matrix is processed all at once.']);
nhood = varargin{2};
block = varargin{3};
else
% wiener2(I,[m n],noise)
nhood = varargin{2};
noise = varargin{3};
end
case 4
% wiener2(I,[m n],[mblock nblock],noise) OBSOLETE
warning(['WIENER2(I,[m n],[mblock nblock],noise) is an obsolete syntax.',...
'Omit the block size, the image matrix is processed all at once.']);
g = varargin{1};
nhood = varargin{2};
block = varargin{3};
noise = varargin{4};
otherwise
msg = 'Too many input arguments.';
return;
end
% checking if input image is a truecolor image-not supported by WIENER2
if (ndims(g) == 3)%判断是否为三维图像
msg = 'WIENER2 does not support 3D truecolor images as an input.';
return;
end;
if (isempty(block))
block = bestblk(size(g));
附录C
% iterative wiener filter
clear all;
tic
% orignal image
I=checkerboard(2);
% num of iterations
num_iter=128;
for snr =20% 20:10:40 %db
res = size(I,1) % image size
res2 = res^2 % image size square
%%% degradation %%%
h=[6;1;1;1;1;1;zeros(res2-11,1);1;1;1;1;1]; % blur
h=h/sum(h); % normalize
Dh= fft(h); % diagonal
for c=1:1 % process RGB components individually
f=I; % orignal
f_vec=reshape(f, res2,1); % vector image
mn_f=mean(f_vec); % mean
amp_f = norm(f_vec)/res; % amplitude
noise_vec = randn(res2, 1)*amp_f/10^(snr/20); % Guassian noise given SNR
Dn = abs(fft(noise_vec )).^2/res2; % psd of noise
%%% observation %%%
g_vec = noise_vec;
for i=1:res2 % degraded image
g_vec(i) = g_vec(i) + [h(res2-i+2: res2)', h(1:res2-i+1)']*f_vec; % H %is circulant matrix
end
mn_g = mean(g_vec);
g_vec = g_vec -mn_g; % forced to be zero-mean
Dg = abs(fft(g_vec )).^2/res2; % psd of degraded image
%%%%%%%%%%%% Wiener %%%%%%%%%%%%%%%
Df = Dg; % basic iterative
Df_add = Dg; % additive
Dh_sq = abs(Dh).*abs(Dh);
for i=1:num_iter
i
B = real(ifft((Df.*conj(Dh))./(Dh_sq.*Df + Dn))) ;
B_add = real(ifft((Df_add.*conj(Dh))./(Dh_sq.*Df_add + Dn) ));
for k=1:res2
