图像复原
4.1图像退化/复原处理的模型
用退化函数对退化过程建模，它和附加噪声选项一起，作用于输入图像f(x,y),产生一幅退化的图像g(x,y):
复原的目标就是得到原始图像的估计。我们要使这个估计尽可能接近原始的输入图像。通常，我们对H和η(x,y)知道得越多，图像的估计就越接近f(x,y)。
4.2噪声模型
能够模拟噪声的行为和效果的能力是图像复原的核心。
4.2.1用imnoise函数为图像添加噪声
图像处理工具箱釆用 imnoise 函数，使噪声污染一幅图像。这个函数有如下基本语法:
g = imnoise(f, type, parameters)
函数 imnoise 在为图像添加噪声之前，将图像转换为范围在[0，1]的 double 类。
1.将均值为m，方差为var的高斯噪声加到图像 f 上。默认为均值是 0、方差是 0.01 的噪声。
g = imnoise(f, ’gaussian', m,var)
>> f = im2double(rgb2gray(imread('d.jpg')));
>> g1 = imnoise(f, 'gaussian', 0.1,0.01);
>> g2 = imnoise(f, 'gaussian', 0.1,0.05);
>> g3 = imnoise(f, 'gaussian', 0.2,0.01);
>> subplot(2, 2, 1), imshow('d.jpg'), title('原图像');
>> subplot(2, 2, 2), imshow(g1,[]), title('均值为0.1方差为0.01的高斯噪声加到图像');
>> subplot(2, 2, 3), imshow(g2,[]), title('均值为0.1方差为0.05的高斯噪声加到图像');
>> subplot(2, 2, 4), imshow(g3,[]), title('均值为0.2方差为0.01的高斯噪声加到图像');
2将均值为 0、局部方差为 V的高斯噪声添加到图像 f 上。其中，V 是与 f 大小相同的数组，其中包含了每个点的理想方差值。
g=imnoise (f,localvar,V)
>> f = im2double(rgb2gray(imread('car.jpg')));
>> [n,m] = size(f);
>> v = rand(n,m);
>> g = imnoise(f,'localvar',v);
>> subplot(1, 2, 1), imshow('car.jpg'), title('原图像');
>> subplot(1, 2, 2), imshow(g), title('均值为0、局部方差为V的高斯噪声加到图像');
3.用椒盐噪声污染图像 f。其中，d是噪声密度(也就是包含噪声值的图像区域的百分比)。因此，大约d*nurael (f) 个像素受到影响。 默认的噪声密度是 0.05。
g=imnoise (f, 'salt & pepper ',d)
>> f = im2double(rgb2gray(imread('car.jpg')));
>> g = imnoise (f, 'salt & pepper',0.12);
>> subplot(1, 2, 1), imshow('car.jpg'), title('原图像');
>> subplot(1, 2, 2), imshow(g), title('椒盐噪声加到图像');
- 用方程式 g=f+n.*f 将乘性噪声添加到图像 f 上。其中，n是均值为 0、方差为 v 的均匀分布的随机噪声。var 的默认值是0.04。
g = imnoise(f,'speckle',v)
f = im2double(rgb2gray(imread('car.jpg')));
>> g = imnoise(f,'speckle',0.2) ;
>> subplot(1, 2, 1), imshow('car.jpg'), title('原图像');
>> subplot(1, 2, 2), imshow(g), title('均值为0方差为V的均匀分布随机噪声加到图像');
5.从数据中生成泊松噪声来代替人造的噪声并添加到数据中。
g=imnoise (f, 'poisson' )
>> f = im2double(rgb2gray(imread('car.jpg')));
>> g = imnoise(f,'Poisson');
>> subplot(1, 2, 1), imshow('car.jpg'), title('原图像');
>> subplot(1, 2, 2), imshow(g), title('泊松噪声加到图像');
实验分析：均值增大，图像越明亮，大部分颜色细节丢失，只有原图像中颜色较深的部分显示出来，因此这种噪声污染比较严重。方差变大，噪声越明显，同时也会丢失细节，会有点状噪声产生，影响辨识。上述的泊松噪声似乎看不出什么变化，所以换了张图片还是没看出明显区别。（此为网上查询：由于光具有量子特效，到达光电检测器表面的量子数目存在统计涨落，因此，图像监测具有颗粒性，这种颗粒性造成了图像对比度的变小以及对图像细节信息的遮盖，我们对这种因为光量子而造成的测量不确定性成为图像的泊松噪声。）
4.2.2用给定分布产生空间随机噪声
在函数imnoise中，能够产生可用类型和参数的噪声是很有必要的。空间噪声值是随机数，虽然可以用 imnoise 产生这两种类型的噪声，但是从目前的上下文来看，更简单的是使用 MATLAB函数 rand 来产生均匀随机数，并使用函数 randn 产生正态(高斯)随机数。
假设拥有在区间(0， 1)内均匀分布的随机数产生器 W，并假设要用它来产生具有累积分布函数(CDF )的随机数z，形式如下：
其中，b>0。为得到z，解下面的方程:
由于平方根是非负的，因此我们可以确定：不会产生小于 a 值的 z。正如瑞利 CDF的定义中要求的那样。因此，来自产生器的均匀随机数 w可以被用在前面的方程式中，从而产生参数为 a 和 b的瑞利分布的随机变量 z。
随机函数rand:
M =10;N=10;
A = rand(M,N)
A =
1 至 7 列
0.8147 0.1576 0.6557 0.7060 0.4387 0.2760 0.7513
0.9058 0.9706 0.0357 0.0318 0.3816 0.6797 0.2551
0.1270 0.9572 0.8491 0.2769 0.7655 0.6551 0.5060
0.9134 0.4854 0.9340 0.0462 0.7952 0.1626 0.6991
0.6324 0.8003 0.6787 0.0971 0.1869 0.1190 0.8909
0.0975 0.1419 0.7577 0.8235 0.4898 0.4984 0.9593
0.2785 0.4218 0.7431 0.6948 0.4456 0.9597 0.5472
0.5469 0.9157 0.3922 0.3171 0.6463 0.3404 0.1386
0.9575 0.7922 0.6555 0.9502 0.7094 0.5853 0.1493
0.9649 0.9595 0.1712 0.0344 0.7547 0.2238 0.2575
8 至 10 列
0.8407 0.3517 0.0759
0.2543 0.8308 0.0540
0.8143 0.5853 0.5308
0.2435 0.5497 0.7792
0.9293 0.9172 0.9340
0.3500 0.2858 0.1299
0.1966 0.7572 0.5688
产生一个M*N大小的矩阵，其中数值在区间（0，1）之间。概率分布为均匀分布。应用：
>> f = rgb2gray(imread('b.jpg'));
>> I = find(f<128);
>> f(I) = 0;
>> subplot(121),imshow('b.jpg'),title('原图')
>> subplot(122),imshow(f),title('灰度值小于128的改为0')
与 imnoise 不同，下面的 M-函数产生 MxN 大小的噪声数组 R， 它不以任何方式缩放。另一个主要的不同是：imnoise 输出有噪声的图像， 而 inmoise2 产生噪声模式本身。用户可以直接指定希望的噪声参数值。
imnoise2方法：
function R=imnoise2(type,M,N,a,b)%type是函数类型，M*N是噪声数组的大小，a,b为两个参数
%设置默认值
if nargin==1%如果函数的输入参数为1，则默认a=0;b=1;M=1;N=1
a=0;b=1;
M=1;N=1;
elseif nargin==3%如果函数的输入参数为3，则默认a=0;b=1
a=0;b=1;
end
%开始运行程序
switch lower(type)
case 'uniform'
R = a + (b - a) * rand(M, N);
case 'gaussian'%如果是高斯类型，执行下面方程
R=a+b*randn(M,N);
case 'salt & pepper'%如果是焦盐类型，当输入参数小于等于3，a=0.05,b=0.05
if nargin<=3
a=0.05;b=0.05;
end
%检验Pa+Pb是否大于1
if(a+b)>1
error('The sum Pa+Pb must be not exceed >1')
end
R(1:M,1:N)=0.5;
X=rand(M,N);%(0,1)范围内产生一个M*N大小的均匀随机数组
c=find(X<=a);%寻找X中小于等于a的数，并赋值为0
R(c)=0;
u=a+b;
c=find(X>a & X<=u);%寻找X中大于a并小于等于u的数，并赋值为1
R(c)=1;
case 'lognormal'%对数正态类型，当输入参数小于等于3，a=1,b=0.25,执行下面方程
if nargin<=3
a=1;b=0.25;
end
R=a*exp(b*randn(M,N));
case 'rayleigh'%瑞利类型，执行下面方程
R=a+(-b*log(1-rand(M,N))).^0.5;
case 'exponential'%指数类型，执行下面方程
if nargin<=3%如果输入参数小于等于3，a=1
a=1;
end
if a<=0%如果a=0,错误类型
error('Parameter a must be positive for exponential type.')
end
k=-1/a;
R=k*log(1-rand(M,N));
case 'erlang'%厄兰类型，如果输入参数小于等于3，a=2，b=5
if nargin<=3
a=2;b=5;
end
if(b~=round(b)|b<=0)%如果b=0,错误类型
error('Param b must a positive integer for Erlang.')
end
k=-1/a;
R=zeros(M,N);
for j=1:b
R=R+k*log(1-rand(M,N));
end
otherwise%如果不是以上类型，输出未知分配类型
error('Unknown distribution type.')
end
观察imnoise2噪声直方图：
r = imnoise2('gaussian',1000,1,0,1);
p = imnoise2('uniform',10000,1,0,1);
q = imnoise2('lognormal',10000,1,1,0.25);
a = imnoise2('rayleigh',10000,1,0,1);
b = imnoise2('exponential',10000,1,1);
c = imnoise2('erlang',10000,1,2,5);
subplot(2, 3, 1), hist(r,50), title('高斯随机数直方图');
subplot(2, 3, 2), hist(p,50), title('均匀随机数直方图');
subplot(2, 3, 3), hist(q,50), title('对数正态随机数直方图');
subplot(2, 3, 4), hist(a,50), title('瑞利随机数直方图');
subplot(2, 3, 5), hist(b,50), title('指数随机数直方图');
subplot(2, 3, 6), hist(c,50), title('厄兰随机数直方图');
4.2.3周期噪声
一幅图像的周期噪声典型地产生于图像获取过程中的电气和/或电动机械的干扰。图像的周期噪声通常通过频域滤波 来处理。周期噪声的模型是 2D正弦波，它有如下所示的方程式：
其中，x=0,1,2，…，M-1且y=0,1,2,…,N-1。A 是振幅，
u
0
u_0
u0和
v
0
v_0
v0分别确定了关于x 轴和y 轴的正弦频率。
B
x
B_x
Bx和
B
y
B_y
By分别是关于原点的相移。
周期噪声imnoise3 方法：
function [r,R,S]=imnoise3(M,N,C,A,B)
%产生一个大小为M*N的正弦噪声模型r,R为傅里叶变换，S为正弦噪声模型的傅里叶的频谱
%C为冲击位置的坐标，A是1*k维向量包含振幅的冲击对，B是由k*2矩阵组成的冲击对
%处理输入参数
[K,n]=size(C);%矩阵C的行数返回给K，矩阵C的列数返回给n
if nargin==3
A(1:K)=1.0;
B(1:K,1:2)=0;
elseif nargin==4
B(1:K,1:2)=0;
end
R=zeros(M,N);%构造R
for j=1:K%j从1到K取遍，步长为1
u1=M/2+1+C(j,1);v1=N/2+1+C(j,2);
R(u1,v1)=i*(A(j)/2)*exp(i*2*pi*C(j,1)*B(j,1)/M);
%复共轭
u2=M/2+1-C(j,1);v2=N/2+1-C(j,2);
R(u2,v2)=-i*(A(j)/2)*exp(i*2*pi*C(j,2)*B(j,2)/N);
end
S=abs(R);%计算光谱
r=real(ifft2(ifftshift(R)));%产生空间正弦图形
空间正弦噪声模式：
>> C = [0 64; 0 128; 32 32; 64 0; 128 0;-32 32];
[r,R,S] = imnoise3(512,512,C);
subplot(1, 3, 1), imshow(S,[]), title('指定脉冲的谱');
subplot(1, 3, 2), imshow(r,[]), title('在空间域中相应的正弦噪声模式');
其他噪声模式：
C=[6 32;-2 2];
[r,R,S] = imnoise3(512,512,C);
imshow(r,[ ]), title('(e)其他噪声模式');
4.2.4估计噪声参数
估计周期噪声参数的典型方法是分析图像的傅里叶谱。周期噪声往往产生频率尖峰，频率尖峰可以通过目测来检测。当噪声尖峰足够明显时，或在干扰频率的一些知识可用的情况下，自动分析是有可能的。对于空间域噪声，PDF的参数可以通过传感器的技术来说明，但是通过样本图像来估计它们是很有必要的。
函数statmoments计算平均值和n阶中心距：
B = roipoly(f,c,r) %f是感兴趣图像，c和r是相应多边形的顶点列坐标和行坐标
函数计算ROI的直方图：
function [p, npix] = histroi(f, c, r)
% 用于计算一幅图像在多边形区域内的直方图，多边形的顶点由c和r指定，通过函数roipoly复制所定义的多边形区域
% Generate the binary mask image.
B = roipoly(f, c, r);
% Compute the histogram of the pixels in the ROI.
p = imhist(f(B));
% Obtain the number of pixels in the ROI if requested in the output.
if nargout > 1
npix = sum(B(:));
end
示例：
f = rgb2gray(imread('scene.jpg'));
figure,imshow(f);
c=[222 272 300 270 221 194];
r=[21 21 75 121 121 75];
[B,c,r] = roipoly(g); %产生交互式模板B
[h,npix] = histroi(g,c,r);
figure,bar(h,1);
4.3仅有噪声的复原——空间滤波
如果出现的退化仅仅是噪声，那么退化就遵循模型：
在这种情况下，选择的降低噪声的方法是空间滤波。
4.3.1空间噪声滤波器
spfilt自定义函数：
function f = spfilt(g, type, m, n, parameter)
%SPFILT Performs linear and nonlinear spatial filtering.
% F = SPFILT(G, TYPE, M, N, PARAMETER) performs spatial filtering
% of image G using a TYPE filter of size M-by-N. Valid calls to
% SPFILT are as follows:
%
% F = SPFILT(G, 'amean', M, N) Arithmetic mean filtering.
% F = SPFILT(G, 'gmean', M, N) Geometric mean filtering.
% F = SPFILT(G, 'hmean', M, N) Harmonic mean filtering.
% F = SPFILT(G, 'chmean', M, N, Q) Contraharmonic mean
% filtering of order Q. The
% default is Q = 1.5.
% F = SPFILT(G, 'median', M, N) Median filtering.
% F = SPFILT(G, 'max', M, N) Max filtering.
% F = SPFILT(G, 'min', M, N) Min filtering.
% F = SPFILT(G, 'midpoint', M, N) Midpoint filtering.
% F = SPFILT(G, 'atrimmed', M, N, D) Alpha-trimmed mean filtering.
% Parameter D must be a
% nonnegative even integer;
% its default value is D = 2.
%
% The default values when only G and TYPE are input are M = N = 3,
% Q = 1.5, and D = 2.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2003/10/27 20:07:00 $
% Process inputs.
if nargin == 2
m = 3; n = 3; Q = 1.5; d = 2;
elseif nargin == 5
Q = parameter; d = parameter;
elseif nargin == 4
Q = 1.5; d = 2;
else
error('Wrong number of inputs.');
end
% Do the filtering.
switch type
case 'amean'
w = fspecial('average', [m n]);
f = imfilter(g, w, 'replicate');
case 'gmean'
f = gmean(g, m, n);
case 'hmean'
f = harmean(g, m, n);
case 'chmean'
f = charmean(g, m, n, Q);
case 'median'
f = medfilt2(g, [m n], 'symmetric');
case 'max'
f = ordfilt2(g, m*n, ones(m, n), 'symmetric');
case 'min'
f = ordfilt2(g, 1, ones(m, n), 'symmetric');
case 'midpoint'
f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
f = imlincomb(0.5, f1, 0.5, f2);
case 'atrimmed'
if (d <= 0) | (d/2 ~= round(d/2))
error('d must be a positive, even integer.')
end
f = alphatrim(g, m, n, d);
otherwise
error('Unknown filter type.')
end
%-------------------------------------------------------------------%
function f = gmean(g, m, n)
% Implements a geometric mean filter.
inclass = class(g);
g = im2double(g);
% Disable log(0) warning.
warning off;
f = exp(imfilter(log(g), ones(m, n), 'replicate')).^(1 / m / n);
warning on;
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = harmean(g, m, n)
% Implements a harmonic mean filter.
inclass = class(g);
g = im2double(g);
f = m * n ./ imfilter(1./(g + eps),ones(m, n), 'replicate');
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = charmean(g, m, n, q)
% Implements a contraharmonic mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g.^(q+1), ones(m, n), 'replicate');
f = f ./ (imfilter(g.^q, ones(m, n), 'replicate') + eps);
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = alphatrim(g, m, n, d)
% Implements an alpha-trimmed mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
for k = (m*n - (d/2) + 1):m*n
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
f = f / (m*n - d);
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function image = changeclass(class, varargin)
switch class
case 'uint8'
image = im2uint8(varargin{:});
case 'uint16'
image = im2uint16(varargin{:});
case 'double'
image = im2double(varargin{:});
otherwise
error('Unsupported IPT data class.');
end
示例：
>>g = rgb2gray(imread('car.jpg'));
>> [M,N] = size(g);
>>
>> R1 = imnoise2('salt & pepper',M,N,0.1,0);
>> gp = g;
gp(R1 == 0) = 0;
R2 = imnoise2('salt & pepper',M,N,0,0.1);
gs = g;
gs(R2 == 1) = 255;
>> fp = spfilt(gp,'chmean',3,3,1.5); %用正值的反调和均值滤波
>> fs = spfilt(gs,'chmean',3,3,1.5); %用负值的反调和均值滤波
>> fpmax = spfilt(gp,'max',3,3);
fsmin = spfilt(gs,'min',3,3);
>> subplot(321);imshow(R1,[]);title('被胡椒噪声污染的图像');
subplot(322);imshow(R2,[ ]);title('被盐粒噪声污染的图像');
subplot(323);imshow(fp,[]);title('滤除胡椒噪声后的图像');
subplot(324);imshow(fs,[ ]);title('滤除被盐粒噪声后的图像');
subplot(325);imshow(fpmax,[]);title('用最大滤波后的结果');
subplot(326);imshow(fsmin,[ ]);title('用最小滤波后的结果');
分析：被胡椒噪声污染的图像会产生白色背景和黑色颗粒，被盐粒噪声污染的图像会产生黑色背景和白色颗粒。用正值的反调和均值法滤除胡椒噪声，基本可以滤除噪声并还原图像，虽然会产生模糊；但用负值的反调和均值法滤除盐粒噪声的效果较差，只能还原主体图像而背景图像仍然没有被还原。最大值滤波器和最小值滤波器能够比较好的滤除胡椒噪声和盐粒噪声，反调和均值则会模糊图像。
4.3.2自适应空间滤波
前面讨论的滤波器应用于不考虑图像特性在不同位置之间存在差异的那些图像。在有些应用中，可以通过使用能够被滤波的区域的图像特性自适应的滤波器来改进结果。主要考虑自适应中值滤波器。自适应中值滤波算法工作在两个层面，表示为 levelA和 levelB
Level A：如 果 Zmin<Zmed<Zmax ，转 向 levelB; 否则增大窗口尺寸 ，如果窗口尺寸小于等于Smax , 重复 levelA；否则输出Zmed
Level B：如 果 Zmin<Zxy<Zmax ，输出 Zxy; 否则输出Zmed。
Zmin=Sxy中的最小亮度值
Zmax=Sxy中的最大亮度值
Zmed=Sxy中的中值
Zxy=坐标(x,y)处的亮度值
Smax表示自适应滤波器窗口允许的最大尺寸。
名为 adpmedian的 M-函数可实现这个算法，语法为：
f = adpmedian(g, Smax） %g是被滤波的图像
function f = adpmedian(g,Smax)
if (Smax<=1)||(Smax/2==round(Smax/2))||(Smax~=round(Smax))
error('Smax must be an odd integer >1.')
end
f=g;
f(:)=0;
alreadyProcessed=false(size(g));
for k=3:2:Smax
zmin=ordfilt2(g,1,ones(k,k),'symmetric');
zmax=ordfilt2(g,k*k,ones(k,k),'symmetric');
zmed=medfilt2(g,[k k],'symmetric');
processUsingLevelB=(zmed>zmin)&(zmax>zmed)&...
~alreadyProcessed;
zB=(g>zmin)&(zmax>g);
outputZxy=processUsingLevelB & zB;
outputZmed=processUsingLevelB&~zB;
f(outputZxy)=g(outputZxy);
f(outputZmed)=zmed(outputZmed);
alreadyProcessed=alreadyProcessed | processUsingLevelB;
if all(alreadyProcessed(:))
break;
end
end
f(~alreadyProcessed)=zmed(~alreadyProcessed);
end
示例：
>> f = imread('scene.jpg');
>> g1 = rgb2gray(f);
g = imnoise(g1,'salt & pepper',0.25);
f1 = medfilt2(g,[7 7],'symmetric');
f2 = adpmedian(g,7);
subplot(131);imshow(g);title('被椒盐噪声污染的图像');
subplot(132);imshow(f1);title('中值滤波器得到的图像');
subplot(133);imshow(f2);title('自适应中值滤波器得到的图像');
结论： 经过对比，自适应中值滤波器比中值滤波器复原效果更加好。中值滤波器会模糊图像，而自适应中值滤波器会更好的还原图像，细节被复原得很清晰，基本与原图像无差别。总体来说使用自适应中值滤波得到的结果比使用大小为7x7的中值滤波得到的结果图还原效果要好，可以仔细观察下图
车的纹路使用7*7中值明显看出纹路模糊了，而自适应的清晰很多。
4.4 退化函数建模
在图像复原问题中，遇到的主要退化是图像模糊。场景和传感器产生的模糊可以用空间域 或频域的低通滤波器来建模。另一个重要的退化模型是由于在图像获取时传感器和场景之间的均匀线性运动而产生的图像模糊。图像模糊可以使用工具箱函数 fspecial 来建模:
PSF= fspecial('motion 1 , len, theta)
调用fspecial 将返回PSF,近似于由len像素的摄像机线性移动产生的效果。参数 theta 以度为单位，以顺时针方向对正水平轴进行量度。Len 和 theta 的默认值分别是 9 和 0, 这些设置对应在水平方向上移动9 个像素。使用函数 imfilter 来创建已知PSF 或用刚才描述的方法计算得到的PSF 的退化图像:
g = imfilter(f, PSF, 'circular');
'circular’用来减少边缘效应。然后通过添加适当的噪声来完成退化图像模拟:
g = g+noise;
noise 是与 g 大小相同的随机噪声图像。由函数 checkerboard 产生的测试模式对于实现这个目的非常有 用，因为大小可以缩放，但却不会影响主要特征。语法为：
C = checkerboard(NP, M, N)
NP 是每个正方形一边的像素数，M 是行数，N 是列数。如果省略 N，默认为 M。 如果 M 和 N 都省略，将产生一边为 8 个正方形的方形棋盘板。另外，如果省略NP, 将默认为 10 像素。在棋盘板上，左半部分的亮正方形是白色的，右半部分的亮正方形是灰色的。用下面 的命令可产生亮正方形全是白色的棋盘板：
K = checkerboard(NP, M, N) > 0.5;
函数 checkerboard 产生的图像属于 double 类，值在区间[0,1]内。 由于有些复原算法对于大图像来说很慢，因此好的办法是用小图像做实验，从而减少计算时间。在这种情况下，如果目的是显示，那么通过像素复制来放大图像是很有用的。下面的函数做了这项工作
B = pixeldup(A, m, n)
这个函数将 A的每个像素在垂直方向上总共复制了 m 次，在水平方向上总共复制了n 次。 如果省略n, 将默认为m。
示例：
f = checkerboard(8); %产生棋盘板图像
figure,subplot(221),imshow(f);title('原始图像');
PSF = fspecial('motion',7,45);
gb = imfilter(f,PSF,'circular'); %产生退化图像
subplot(222),imshow(gb);title('使用motion滤波器模糊后图像');
noise=imnoise(zeros(size(f)),'gaussian',0,0.002);%产生高斯噪声
subplot(223),imshow(noise,[]);title('高斯纯噪声');
g=gb+noise;
subplot(224),imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
4.5 直接逆滤波
逆滤波，就是简单将退化函数去除，直接的逆滤波没有什么意义，只处理了靠近直流分量的部分，其他不做处理。由退化函数H退化的图像复原的最简单的方法是直接做逆滤波，设图像退化前的傅里叶变换为F(u,v)，退化后的傅里叶变换为G(u,v)，系统函数即退化函数的傅里叶变换为H(u,v)。
所谓直接逆滤波，就是用退化函数除退化图像的傅里叶变换，得到退化前图像的傅里叶变换的估计。公式如下所示：
由该式可知，即使知道退化函数，也不能准确的复原图像，因为N(u,v)未知，甚至有更糟的情况是如果退化函数是零或是非常小的值，则N(u,v)/H(u,v)的值比较大，很容易支配F(u,v)的估计值。在无噪声的情况下，逆滤波是完美的。但实际中，所有都会噪声，因此这个类型逆滤波的方法很少使用到。
4.6 维纳滤波
维纳滤波也称为最小均方差误差滤波，它是一种最早、也是最为熟知的线性图像复原方法。这种滤波就是使图像尽可能的平滑，可以消除很严重的噪声，恢复图像。
维纳滤波器寻找统计误差函数最小的估计f ̂
其中，E是期望值算子，f是未退化图像。这个表达式在频域中的表达是:
维纳滤波是通过函数deconvwnr实现的，语法如下：（g是退化图像，frest是复原图像）
frest = deconvwnr(g,PSF); %信噪比为0
frest = deconvwnr(g,PSF,NAPR); %信噪比已知或是常量或是数组，NAPR是标量输入
frest = deconvwnr(g,PSF,NACORR,FACORR); %噪声和未退化的图像的自相关函数NACORR,FACORR是已知的
示例：
>> f = checkerboard(8); %产生棋盘板图像
gb = imfilter(f,PSF,'circular'); %产生退化图像
noise=imnoise(zeros(size(f)),'gaussian',0,0.002);%产生高斯噪声
g=gb+noise;
>> subplot(221),imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
frest1=deconvwnr(g,PSF); %维纳滤波
subplot(222),imshow(frest1,[]),title('维纳滤波后的结果');
Sn = abs(fft2(noise)).^2; %噪声功率谱
nA=sum(Sn(:))/prod(size(noise)); %平均噪声功率谱
Sf = abs(fft2(f)).^2; %信号功率谱
fA=sum(Sf(:))/prod(size(f)); %平均信号功率谱
R = nA/fA; %求信噪比
frest2=deconvwnr(g,PSF,R);
subplot(223),imshow(frest2,[]),title('使用常数比例的维纳滤波');
NACORR = fftshift(real(ifft(Sn)));
FACORR = fftshift(real(ifft(Sf)));
frest3=deconvwnr(g,PSF,NACORR,FACORR); %自相关后的维纳滤波
subplot(224),imshow(frest3,[]),title('使用自相关函数的维纳滤波');
结论：自相关函数的维纳滤波效果比较好，简单的维纳滤波效果最差。自相关函数的维纳滤波并没有很好地完全还原图像，而且还有一些噪声存在；相比之下，简单的维纳滤波存在的噪声更加多，效果比较差，亮度也没有得到还原。
当模糊但是不含噪声的图像使用维纳滤波还原：
>> clear
>> f = checkerboard(8);
PSF = fspecial('motion', 7, 45);
gb = imfilter(f, PSF,'circular');
g =gb;
fr1 = deconvwnr(g,PSF);
subplot(1, 2, 1), imshow(g,[]);title('(a)模糊的图像');
subplot(1, 2, 2), imshow(fr1,[]);title('(b)简单的维纳滤波(逆滤波)后的结果');
实验结论：不添加噪声后，简单的维纳滤波（逆滤波）后的结果，已经比添加噪声后的使用自相关函数的维纳滤波后的结果要好。
4.7 约束的最小二乘法(规则化)滤波
约束的最小二乘法滤波在工具箱中是通过函数 deconvreg 实现的，语法如下：
frest = deconvreg{gf PSF, NOISEPOWER, RANGE)
示例：
f = checkerboard(8); % 产生一个一面为8个正方形的测试板
PSF = fspecial('motion',7,45); % 运动模糊，PSF刚好为空间滤波器
gb = imfilter(f,PSF,'circular'); % 减少边界效应
noise = imnoise(zeros(size(f)),'gaussian',0,0.001); % 高斯噪声
g = gb + noise; % 添加高斯噪声构造退化的图像模型
fr1 = deconvreg(g,PSF,4); % 正则滤波器
fr2 = deconvreg(g,PSF,0.4,[1e-7 1e7]);
subplot(1,3,1);imshow(g);title('(a)运动模糊+高斯噪声污染后图像');
subplot(1,3,2);imshow(fr1);title('(b)用NOISEPOWER等于4的规则化滤波器复原');
subplot(1,3,3);imshow(fr2);title('(c)用NOISEPOWER等于0.4且RANGE为[1e-7 1e7]的规则化滤波器复原');
实验分析：NOISEPOWER的初始值为4是由图像大小是64*64，噪声有0.001的方差和0均值。
（
64
)
2
（64)^2
（64)2(0.001+0)≈4，可以发现使用这个估计值效果并不是很好。须把 NOISEPOWER 的值的量级下调一个等级之后，可以看出RANGE 比默认值更加紧缩。
4.8利用露西-理查德森算法的迭代非线性复原
在工具箱中，L-R算法用 deconvlucy 函数实现的，基本语法如下：
f = deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT)
f 代表复原的图像，g 代表退化图像，PSF是点扩散函数，NUMIT 为迭代的次数(默认为 10次)， DAMPAR 是标量，指定了结果图像与原始图像 g 之间的偏离阈值。当像素值偏离原值的范 围在 DAMPAR 内时，就不用再迭代。这既抑制了这些像素上的噪声，又保留了图像细节。默认值为 0(无衰减)。 WEIGHT是数组，大小与 g 相同，它为每个像素都施以权重以反映像素的质量。
示例：
>> f = checkerboard(8); % 产生64×64像素的方形图像
PSF = fspecial('gaussian',7,10); % 产生一个大小为7×7且标准偏差为10的高斯PSF
SD = 0.01; % 标准偏差
g = imnoise(imfilter(f,PSF),'gaussian',0,SD^2); % 添加均值为0、标准偏差为0.01的高斯噪声
DAMPAR = 10*SD; % 结果图像与原图像的偏离阈值
LIM = ceil(size(PSF,1)/2);
WEIGHT = zeros(size(g)); % WEIGHT数组大小64×64
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1; % WEIGHT数组有值为0的4像素宽的边界，其余像素为1
f1 = deconvlucy(g,PSF,5,DAMPAR,WEIGHT); % 迭代5次
f2 = deconvlucy(g,PSF,10,DAMPAR,WEIGHT); % 迭代10次
f3 = deconvlucy(g,PSF,20,DAMPAR,WEIGHT); % 迭代20次
f4 = deconvlucy(g,PSF,100,DAMPAR,WEIGHT); % 迭代100次
subplot(2,3,1);imshow(pixeldup(f,8));title('(a)原图像');
subplot(2,3,2);imshow(g);title('(b)由高斯噪声污染和模糊的图像');
subplot(2,3,3);imshow(f1);title('(c)使用L-R算法迭代5次复原图像');
subplot(2,3,4);imshow(f2);title('(d)使用L-R算法迭代10次复原图像'); % 图像虽然改进，但依然模糊
subplot(2,3,5);imshow(f3);title('(e)使用L-R算法迭代20次复原图像'); % 迭代20次为合理复原
subplot(2,3,6);imshow(f4);title('(f)使用L-R算法迭代100次复原图像'); % 除了稍微清晰明亮，并无显著改进
>>
实验分析：关于什么时候停止 L-R 算法通常很难回答。一种方法是观察输出，当在给定的应用中，可接受的结果已经得到时停止算法。从上面的实验看出，当迭代了20次的时候图像已经基本清晰.观察下图：
当迭代100次时（右图），图像并没有变得更清晰，而且还出现了些不希望的结果。
4.9 盲去卷积
在图像复原过程中，最困难的问题之一是前面讨论过的对这些复原算法 PSF 的恰当估计。正如先前表明过的那样，那些不以PSF 知识为基础的图像复原方法统称为盲去卷积算法。 盲去卷积的基本方法以最大似然估计为基础，也就是对随机噪声污染的量进行估计时采用的最佳策略。
工具箱用函数 deconvblind 来执行盲去卷积，基本语法如下：
[fr,PSF] = deconvblind{gr INITPSF)
g 是退化图像，INITPSF 是点扩散函数的初始估计。PSF 是这个函数最终计算得到 的估计值，fr 是利用估计的 PSF 复原图像。
示例：
>> f = checkerboard(8);
>> PSF = fspecial('gaussian',7,10);
SD = 0.01;
g = imnoise(imfilter(f,PSF),'gaussian',0,SD^2);
INITPSF = ones(size(PSF)); % 点扩散函数的初始估计
DAMPAR = 10*SD; % 结果图像与原图像的偏离阈值
LIM = ceil(size(PSF,1)/2);
WEIGHT = zeros(size(g)); % WEIGHT数组大小64×64
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1; % WEIGHT数组有值为0的4像素宽的边界，其余像素为1
[fr1,PSFe1] = deconvblind(g,INITPSF,5,DAMPAR,WEIGHT); % 盲去卷积迭代5次
[fr2,PSFe2] = deconvblind(g,INITPSF,10,DAMPAR,WEIGHT); % 盲去卷积迭代10次
[fr3,PSFe3] = deconvblind(g,INITPSF,20,DAMPAR,WEIGHT); % 盲去卷积迭代20次
subplot(2,2,1);imshow(pixeldup(PSF,73),[ ]);title('(a)原始PSF');
subplot(2,2,2);imshow(pixeldup(PSFe1,73),[ ]);title('盲去卷积迭代5次时PSF估计值');
subplot(2,2,3);imshow(pixeldup(PSFe2,73),[ ]);title('盲去卷积迭代10次时PSF估计值');
subplot(2,2,4);imshow(pixeldup(PSFe3,73),[ ]);title('盲去卷积迭代20次时PSF估计值'); % 最接近真正PSF
>>
4.10 来自投影的图像重建
4.10.1来自投影的图像重建
函数radon是用来为给定的二维矩阵列产生一组平行射线投影，基本语法是：
R =radon(I,theta); %I是二维阵列，theta是角度值的一维阵列
用来产生众所周知的图像的某个有用函数语法为：
P = phantom(def,n); %def是指定产生头部模型类型的字符串，n是行数和列数（默认值为256）
其中，字符串def的可用值为：‘Shepp-Logan’（对比度很低）和’ModifiedShepp-Logan’（改进了对比度）。
f = imread('b.jpg ');
g1=rgb2gray(f);
g1(10:50,50:150) = 1;
g2 = phantom('Modified Shepp-Logan',600);
subplot(221),imshow(g1),title('图一');
subplot(222),imshow(g2),title('图二');
theta = 0:0.5:179.5;
[R1,xp1] = radon(g1,theta); %执行雷登变换
[R2,xp2] = radon(g2,theta);
R1 = flipud(R1'); %投影图像
R2 = flipud(R2');
subplot(223),imshow(R1,[],'XData',xp1([1 end]),'YData',[179.5 0]) %显示投影图像
axis xy
axis on
xlabel('\rho'),ylabel('\theta')
subplot(224),imshow(R2,[],'XData',xp2([1 end]),'YData',[179.5 0])
axis xy
axis on
xlabel('\rho'),ylabel('\theta')
扫一扫
1230
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
