@dragonfive 2016-03-10T07:44:45.000000Z 字数 18876 阅读 1237

数字图像处理课程笔记

计算机视觉

一些图像的像素范围映射变化

灰度平方运算

% 下面是在做平方运算;
figure;
img=imread('E:\资料\onedrive\code\test\image\lenargb.jpg');
img=rgb2gray(img);             %matlab把真彩图转化为灰度图;
img=mat2gray(img);  %读入的图片刚开始是uint8的，然后转化为double类型之后才能用mat2gray转化为0-1的double
 i=1; 
for c=[0.04,0.1,0.2,0.4,0.67,1,1.5,2.5,5.0,10.0,25.0]
     y=img.^c;
     subplot(3,4,i);           %把多个图片放在一起的方法;
     imshow(y);
     i=i+1;
     %hold on
end

灰度区间映射

matlab中的imadjust函数可以将灰度在一定范围内的图像到另一个区间里面

g=imadjust(f,[low_in,high_in],[low_out,high_out],gamma)

在[0-1]范围内做平方，能调整灰度范围

%下面是测试matlab画曲线
x=0:0.01:1;
for c=[0.04,0.1,0.2,0.4,0.67,1,1.5,2.5,5.0,10.0,25.0]
    y=x.^c;
    plot(x,y);
    hold on;
end

% 下面测试用imadjust方法做灰度调整
img=imread('E:\资料\onedrive\code\test\image\result.png');
class(img)
figure;
subplot(1,2,1);imshow(img);
img=imadjust(mat2gray(img),[0,1],[0.4,0.6]); %这个函数对于输入图像不要求其类型为[0,1]的double型
subplot(1,2,2);imshow(img);

灰度直方图均衡化

% 下面做灰度直方图均衡和分配
f=imread('E:\资料\onedrive\code\test\image\Fig0206(a)(rose-original).tif');
imshow(f);figure;imhist(f);
%ylim('auto');
g=histeq(f,256);  %进行直方图均衡;
figure;imshow(g);
figure;imhist(g);
%ylim('atuo'); %sets the axis limit mode to auto. 
hnorm=imhist(f)./numel(f); %归一化直方图
cdf=cumsum(hnorm);         %累加
x=linspace(0,1,256);   %在0和1之间产生256个等距离的点;
figure;
plot(x,cdf);
axis([0 1 0 1]);

log变换

% 下面进行log变换
%Fig0206(a)(rose-original).tif
%Fig0305(a)(spectrum).tif
f=imread('E:\资料\onedrive\code\test\image\Fig0305(a)(spectrum).tif');
subplot(2,2,1);imshow(f);
g=im2uint8(mat2gray(log(1+double(f))));%进行log变换
%g=mat2gray(log(1+double(f)));
%h=histeq(f,256);  %进行直方图均衡;
%j=histeq(g,256);
subplot(2,2,3);imhist(f);%画出直方图
subplot(2,2,4);imhist(g);

相对拉伸变换

subplot(2,2,2);imshow(g);

% 下面进行相对伸展变换
%Fig0206(a)(rose-original).tif
%Fig0305(a)(spectrum).tif
%Fig0308(a)(pollen).tif
%下面先画出拉伸变换的图像;
m=0.5;
A=[0:1:10];
B=[10:10:100];
r=[0:0.01:1];
for E=[A B]
    s=1./(1+(m./r).^E);
    subplot(2,2,1);plot(r,s);
    axis([0 1 0 1])
    set(gca,'xtick',0:0.1:1);
    set(gca,'ytick',0:0.1:1);
    xlabel('x');
    ylabel('y');
    hold on
end
%s=1./(1+(m./r).^E)
%plot(r,s);
%下面对图像进行操作
f=imread('E:\资料\onedrive\code\test\image\Fig0308(a)(pollen).tif');
E=2;
s=1./(1+(m./r).^E);
subplot(2,2,2);plot(r,s);
 title('函数图像')     %添加图像标题
 text(0.5,0.5,'E=2')    %将cosx这个注解加到坐标中的某个位置
 %gtext('s=1./(1+(m./r).^E)')      % 用鼠标的光标定位,将sinx这个注解放在你鼠标点击的地方
s=1./(1+(m./double(f)+eps).^E);
subplot(2,2,3);imshow(f);title('原图') 
subplot(2,2,4);imshow(mat2gray(s));title('扩展后的图')

画出自适应的坐标轴

% 画出自适应的坐标轴
function myplot(x,y)
% xyplot Plot 2D axes through the origin
%   Example 1:
%   t = linspace(0,2*pi,500); 
%   y1 = 80*sin(t); 
%   y2 = 100*cos(t);
%   xyplot(t,[y1;y2])
%
%   Example 2:
%   x = -2*pi:pi/10:2*pi;
%   y = sin(x);
%   plot(x,y)
%   xyplot
% PLOT
if nargin>0
    if nargin == 2
        plot(x,y);
    else
        display('   Not 2D Data set !')
    end
end
hold on;
% GET TICKS
X=get(gca,'Xtick');
Y=get(gca,'Ytick');
% GET LABELS
XL=get(gca,'XtickLabel');
YL=get(gca,'YtickLabel');
% GET OFFSETS
Xoff=diff(get(gca,'XLim'))./40;
Yoff=diff(get(gca,'YLim'))./40;
% DRAW AXIS LINEs
plot(get(gca,'XLim'),[0 0],'k');
plot([0 0],get(gca,'YLim'),'k');
% Plot new ticks
for i=1:length(X)
    plot([X(i) X(i)],[0 Yoff],'-k');
end;
for i=1:length(Y)
    plot([Xoff, 0],[Y(i) Y(i)],'-k');
end;
% ADD LABELS
text(X,zeros(size(X))-2.*Yoff,XL);
text(zeros(size(Y))-3.*Xoff,Y,YL);
box off;
% axis square;
axis off;
set(gcf,'color','w');
end

空域滤波　　

基础知识

输出图像是对输入图像的像素做简单的运算得到的，邻域像素也会参与运算
可以是线性运算，也可以是非线性运算

线性滤波函数

g = imfilter(f, w, ltering mode,boundary options, size options)

计算过程中使用的是double, 但是imfilter的结果要保证输出图像和输入图像的数据类型一致，这样可能会在最后转换的时候出现截断的情况，所以正确的做法应该是在调用函数之前将图像转化为double的格式。

 fspecial('type', parameters)

这个函数可以产生线性滤波器，
type取值：

拉普拉斯滤波器

用拉普拉斯滤波器检测边缘来做图像增强,图像锐化。

%下面做拉普拉斯空域滤波
f=imread('E:\资料\onedrive\code\test\image\Fig0316(a)(moon).tif');
w4=fspecial('laplacian',0);
w8=[1 1 1;1 -8 1;1 1 1] ;%自己生成一个矩阵做算子;
g1=imfilter(f,w4,'replicate');%拉普拉斯边缘检测结果;
g2=imfilter(f,w8,'replicate');
douf=im2double(f);%转换为double类型
g3=imfilter(douf,w4,'replicate');
g4=imfilter(douf,w8,'replicate');
g5=f-g1;g6=f-g2;g7=douf-g3;g8=douf-g4;
subplot(2,5,1);imshow(f);
subplot(2,5,2);imshow(douf);
for i=3:10
    subplot(2,5,i),imshow(eval(['g',num2str(i-2)]));%把字符串当作变量名
end

拉普拉斯变换算子的矩阵的所有元素的和为1，而梯度算子的所有元素的和为0

非线性空域滤波

colfilt方法比nlfilter方法更快，在使用cofit方法之前，我们需要对图像进行装填。
在图像处理中非空域里面使用最多的滤波器是统计顺序滤波器，使用最多的函数是ordfilt2

g=ordfilt2(f,order,domain)
通过使用在邻域中的从小到大排序在第order个的元素来替换原图中的每一个元素，邻域的大小由domain确定

g=ordfilt2(f,median(1:m*n),ones(m,n))

这里的median函数用来产生中位数；比如：

g=ordfilt2(f,1,ones(m,n));              %用邻域中最小的元素替代
g=ordfilt2(f,m*n,ones(m,n));            %用邻域中最大的元素替代
g=ordfilt2(f,median(1:m*n),ones(m,n));  %用邻域中中位数元素替代

在统计顺序滤波中最常使用的是中位数滤波，可以用来进行降噪处理;

g=medfilt2(f,[m,n],padopt);

在m*n的邻域中寻找中位数来替代中间元素，padopt是图像边缘的填充方式；

频域滤波

DFT:discrete Fouriter Transform离散傅立叶变换
傅立叶表换具有唯一性。傅立叶变换揭示了信号的时域特性和频域特性之间的确定的内在联系。
在图像里面时域特性表现为空域特性。

一个例子

下面是用傅立叶与逆傅立叶变换操作图像

% 下面求图像的傅立叶变换频谱
f=imread('E:\资料\onedrive\code\test\image\Fig0403(a)(image).tif');
F=fft2(f);
S=abs(F);     % 去F的频谱
w=2;h=3;i=0;
i=i+1;subplot(w,h,i);imshow(f,[]);
i=i+1;subplot(w,h,i);imshow(S,[]);
% imshow(K)与imshow(K,[])的区别;
% imshow(K)直接显示K;
% imshow(K,[])显示K，并将K的最大值和最小值分别作为纯白(255)和纯黑(0)，中间的K值映射为0到255之间的标准灰度值。
%fftshift移动零频点到频谱中间，重新排列fft,fft2和fftn的输出结果。
Fc = fftshift(F);
i=i+1;subplot(w,h,i);imshow(abs(Fc),[]);
%log变换将像素范围压缩;
S2=log(1+abs(Fc));
i=i+1;subplot(w,h,i);imshow(S2,[]);
f=real(ifft2(F));% 逆傅立叶变换;逆变换后取得实部就好
%这次变换中傅立叶变换后没有做什么操作，所以逆傅立叶变换后还是原来的结果;
i=i+1;subplot(w,h,i);imshow(f,[]);

显示的优化

从上面的代码我们看到，在一个窗口输出多个图片，需要自己设置位置什么的，但是刚开始我们不知道最终要显示多少行多少列的图片，而且随着图片的增加需要经常修改显示的参数。同时如果调换显示顺序和插入显示都需要调整很多代码，基于这些不便，我自己写了个imshow函数。程序会一直保存你要输出的图像，知道输入0，把图像全部输出出来，同时能够实现自适应。

function myImshow(img,argTitle)
persistent myNumOfImg %记录最终需要的图片个数;
persistent M %记录所有的图片
persistent Title ; % 保存标题;
if isempty(myNumOfImg)
    myNumOfImg=0;
    M=cell(0);
    Title=cell(0);
end
if isequal(img,0) %表示结束
    h=ceil(sqrt(myNumOfImg)); %宽
    w=ceil(myNumOfImg*1.0/h);  %高
    for i=1:myNumOfImg
        subplot(w,h,i);
        imshow(M{i},[]);
        title(Title{i});
    end
    clear myNumOfImg;clear Title;clear M;
else
    myNumOfImg=myNumOfImg+1;
    M{myNumOfImg}=img;
    if nargin==2
        Title{myNumOfImg}=argTitle;
    else
        Title{myNumOfImg}='';
    end
end

运行示例

clear all;
close all;
f=imread('E:\资料\onedrive\code\test\image\Fig0403(a)(image).tif');
F=fft2(f);
S=abs(F);     % 去F的频谱
myImshow(f);
myImshow(S);
% w=2;h=3;i=0;
% i=i+1;subplot(w,h,i);imshow(f,[]);
% i=i+1;subplot(w,h,i);imshow(S,[]);
% imshow(K)与imshow(K,[])的区别;
% imshow(K)直接显示K;
% imshow(K,[])显示K，并将K的最大值和最小值分别作为纯白(255)和纯黑(0)，中间的K值映射为0到255之间的标准灰度值。
%fftshift移动零频点到频谱中间，重新排列fft,fft2和fftn的输出结果。
Fc = fftshift(F);
% i=i+1;subplot(w,h,i);imshow(abs(Fc),[]);
myImshow(abs(Fc));
%log变换将像素范围压缩;
S2=log(1+abs(Fc));
% i=i+1;subplot(w,h,i);imshow(S2,[]);
myImshow(S2);
f=real(ifft2(F));% 逆傅立叶变换;
%这次变换中傅立叶变换后没有做什么操作，所以逆傅立叶变换后还是原来的结果;
% i=i+1;subplot(w,h,i);imshow(f,[]);
myImshow(f);
myImshow(0);

运行结果

傅立叶变换滤波

滤波步骤

获得填充大小：PQ=paddedSize(size(f))
获得填充后的傅立叶变换 F=fft2(f,PQ(1),PQ(2))
生成一个滤波器，大小为PQ(1)*PQ(2)，需要使用fftshift把重心的波峰移到边角处
频域点乘 G=H.*F
获得逆傅立叶变换出的图像的实部g=real(ifft2(G))
恢复图像原始大小 g=g(1:size(f,1),1:size(f,2))

从空域滤波生成频域滤波

下面测试几个基础函数的使用效果

% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 测试几个频域常用的函数的用法
% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
clear all;
close all;
% 测试gscale 测试freqz2 测试fftshift与ifftshift
f=imread('E:\资料\onedrive\code\test\image\Fig0409(a)(bld).tif');
F=fft2(f);
S=fftshift(log(1+abs(F)));
%subplot(1,3,1);imshow(S);
%subplot(1,3,2);imshow(S,[]);
S=gscale(S);   % 测试得这个函数的作用是让频谱自适应地显示出来;
%subplot(1,3,3);imshow(S);
w=1;h=3;i=0;
% 由空域滤波器生成频域滤波器 对比fftshift与ifftshift
H=fspecial('sobel'); %梯度算子
H3=freqz2(H);         %由空域滤波器生成频域滤波器;
i=i+1;subplot(w,h,i);mesh(abs(H3));%生成函数图像？
% figure;
H1=ifftshift(H3);
i=i+1;subplot(w,h,i);mesh(abs(H1));
H2=fftshift(H3);
i=i+1;subplot(w,h,i);mesh(abs(H2));
%view(45,30); %设置看图像的视角;第一个是左右倾斜度，第二个是上下倾斜度
%我就感觉这个ifftshift和fftshift的效果是一样的

测试后发现矩阵经过gscale后的效果与imshow(f,[])的效果一致。但是明显gscale适用范围更广，可以先把这个矩阵标准化，进行相应操作而不是简简单单的输出。
freqz2:frequcy表示生成频域滤波，用空域滤波生成频域滤波
fftshift:表示把位于中间的波峰移动边上，大概是这个意思？
mesh：二维显示频域函数
view：调整观测角度

paddedsize函数实现

计算傅立叶滤波需要的填充大小

function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering. 
%   PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
%   computes the two-element size vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
%   vectors, computes the two-element size vector PQ.  The elements
%   of PQ are the smallest even integers greater than or equal to 
%   AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]). 
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $
if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 & ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.
   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
elseif nargin == 3
   m = max([AB CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
else 
   error('Wrong number of inputs.')
end

gscale函数的实现

代码来自冈萨雷斯，输出和输入的图像格式一样。效果是规格化图像元素的像素值

function g = gscale(f, varargin)
%GSCALE Scales the intensity of the input image.
%   G = GSCALE(F, 'full8') scales the intensities of F to the full
%   8-bit intensity range [0, 255].  This is the default if there is
%   only one input argument.
%
%   G = GSCALE(F, 'full16') scales the intensities of F to the full
%   16-bit intensity range [0, 65535]. 
%
%   G = GSCALE(F, 'minmax', LOW, HIGH) scales the intensities of F to
%   the range [LOW, HIGH]. These values must be provided, and they
%   must be in the range [0, 1], independently of the class of the
%   input. GSCALE performs any necessary scaling. If the input is of
%   class double, and its values are not in the range [0, 1], then
%   GSCALE scales it to this range before processing.
%
%   The class of the output is the same as the class of the input.
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:36:09 $
if length(varargin) == 0 % If only one argument it must be f.
   method = 'full8';
else 
   method = varargin{1};
end
if strcmp(class(f), 'double') & (max(f(:)) > 1 | min(f(:)) < 0)
   f = mat2gray(f);
end
% Perform the specified scaling.
switch method
case 'full8'
   g = im2uint8(mat2gray(double(f)));
case 'full16'
   g = im2uint16(mat2gray(double(f)));
case 'minmax'
   low = varargin{2}; high = varargin{3};
   if low > 1 | low < 0 | high > 1 | high < 0
      error('Parameters low and high must be in the range [0, 1].')
   end
   if strcmp(class(f), 'double')
      low_in = min(f(:));
      high_in = max(f(:));
   elseif strcmp(class(f), 'uint8')
      low_in = double(min(f(:)))./255;
      high_in = double(max(f(:)))./255;
   elseif strcmp(class(f), 'uint16')
      low_in = double(min(f(:)))./65535;
      high_in = double(max(f(:)))./65535;    
   end
   % imadjust automatically matches the class of the input.
   g = imadjust(f, [low_in high_in], [low high]);   
otherwise
   error('Unknown method.')
end

dftfilt函数

函数的效果是合并了傅立叶变换滤波中的
2. 填充原始图像并进行傅立叶变换
4. 进行滤波
5. 傅立叶反变换并获得实部
6. 剪切成原始尺寸

所以剩下需要做的就是生成一个合适尺寸的频域滤波器并移动规范化

function g = dftfilt(f, H)
%DFTFILT Performs frequency domain filtering.
%   G = DFTFILT(F, H) filters F in the frequency domain using the
%   filter transfer function H. The output, G, is the filtered
%   image, which has the same size as F.  DFTFILT automatically pads
%   F to be the same size as H.  Function PADDEDSIZE can be used to
%   determine an appropriate size for H.
%
%   DFTFILT assumes that F is real and that H is a real, uncentered
%   circularly-symmetric filter function. 
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $
% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));
% Perform filtering. 
g = real(ifft2(H.*F));
% Crop to original size.
g = g(1:size(f, 1), 1:size(f, 2));

空域生成频域的例子

clear all;
close all;
f=imread('E:\资料\onedrive\code\test\image\Fig0409(a)(bld).tif');
w=1;h=3;i=0;
h=fspecial('sobel'); %梯度算子
PQ=paddedsize(size(f));
H3=freqz2(h,PQ(1),PQ(2));         %由空域滤波器生成频域滤波器;
H1=fftshift(H3);
g1=dftfilt(f,H1);
g2=imfilter(double(f),h);
myImshow(g1);myImshow(g2);myImshow(0);

运行结果

从结果可以看到空域滤波和频域滤波效果差不多；
但是时间上频域滤波用时：0.044662 seconds.
空域滤波用时.046712 seconds.
这便是频域滤波的优势所在。

直接生成频域滤波

理想低通滤波器

巴特沃兹低通滤波器

高斯低通滤波器

dftuv的实现

dftuv函数的作用是生成频域的距离向量组。这个距离是距离频域中心的距离

function [U, V] = dftuv(M, N)
%DFTUV Computes meshgrid frequency matrices.
%   [U, V] = DFTUV(M, N) computes meshgrid frequency matrices U and
%   V.  U and V are useful for computing frequency-domain filter
%   functions that can be used with DFTFILT.  U and V are both
%   M-by-N.
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.3 $  $Date: 2003/04/16 22:30:34 $
% Set up range of variables.
u = 0:(M - 1);
v = 0:(N - 1);
% Compute the indices for use in meshgrid.
idx = find(u > M/2);
u(idx) = u(idx) - M;
idy = find(v > N/2);
v(idy) = v(idy) - N;
% Compute the meshgrid arrays.
[V, U] = meshgrid(v, u);

直接用频域高斯低通滤波的例子

clear all;
close all;
f=imread('E:\资料\onedrive\code\test\image\Fig0413(a)(original_test_pattern).tif');
PQ=paddedsize(size(f));
% 生成低通滤波器(注意大小) 移动
[U,V]=dftuv(PQ(1),PQ(2));
D0=0.05*PQ(2);
H=exp(-1*(U.^2+V.^2)/(2*D0.^2));
H1=fftshift(H);
% 滤波输出结果
g=dftfilt(f,H);
f1=fft2(f);
myImshow(f);myImshow(H1); myImshow(g);myImshow(0);

运行结果

直接生成频域高通滤波

高通滤波器可以直接从低通滤波器获得

理想高通滤波器

波特沃兹高通滤波器

高斯高通滤波器
　　

用频域波特沃兹高通滤波的例子

% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 一个频域高斯低通滤波的例子
% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 读取图片

clear all;
close all;
f=imread('E:\资料\onedrive\code\test\image\Fig0413(a)(original_test_pattern).tif');
% 生成高通滤波器
PQ=paddedsize(size(f));
[U,V]=dftuv(PQ(1),PQ(2));
D0=[15,30,80];

for i=1:numel(D0)
H=1-exp(-1*(U.^2+V.^2)/(2*D0(i).^2));
%myImshow(fftshift(H));
g=dftfilt(f,H);
myImshow(g);
end
myImshow(0);

运行结果

用高通滤波进行图像锐化

用图像的原始频率同时在高频部分进行增强。(拉普拉斯检测边缘也是检测的高频部分其实)
% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 用高通滤波做频域增强
% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 读取图片

clear all;
close all;
f=imread('E:\资料\onedrive\code\test\image\Fig0419(a)(chestXray_original).tif');
% 生成波特沃兹滤波器
PQ=paddedsize(size(f));
[U,V]=dftuv(PQ(1),PQ(2));
D0=0.05*PQ(2);
n=[1,5,15,50];
% D0=[15,30,80];
myImshow(gscale(f));
for i=1:numel(n);
b=((U.^2+V.^2)/(D0.^2)).^n(i);
H=b/(b+1);
%myImshow(H);
% 按照公式进行运算
Hh=0.5+2*H;
g=dftfilt(f,Hh);
g=histeq(gscale(g),256);
myImshow(g);
end
myImshow(0);

运行结果

图像恢复

图像降指与恢复过程

噪声模型

imnoise函数生成噪声

用之前需要先将f规范化

g = imnoise(f, type, parameters)

type:

'gussian' 'salt & pepper' 'motion'

parameters

生成特定分布的空域随机噪声

使用累积分布的逆函数，可以生成任意你想要的原始分布；这里假设w是均匀分布，可以得到z的分布；

CDF是从PDF积分得到的

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 模拟噪声
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 
close all;
clear all;
r={};
r{1}=imnoise2('gaussian',100000,2,0,1);
r{2}=imnoise2('uniform',100000,2,0,1);
r{3}=imnoise2('rayleigh',100000,2,0,1);
% for i=1:3
%     subplot(2,2,i);imshow(r{i},50);
% end
% figure;
for i=1:numel(r)
   subplot(2,2,i);mesh(r{i}); 
end
figure;
for i=1:numel(r)
   subplot(2,2,i);hist(r{i});
   size(r{i})
end

只存在空域噪声的恢复

与图像增强的技术一致

均值滤波器

统计排序与调和均值滤波器

空域滤波基础函数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);

逆谐波函数用于消除椒盐噪声

% 逆调和均值滤波的例子
gp  = imread('E:\资料\onedrive\code\test\image\Fig0505(a)(ckt_pepper_only).tif');
gme = spfilt(gp,'gmean',3,3);
gch = spfilt(gp,'chmean',3,3,1.5); % 椒噪声用正的Q;
gar = spfilt(gp,'amean',3,3); 
gme = spfilt(gp,'gmean',3,3);
myImshow(gp,'椒噪');myImshow(gch,'逆调和除噪');myImshow(gar,'均值滤波');%gme
myImshow(gme,'几何均值滤波');
myImshow(0);
figure;
gs  = imread('E:\资料\onedrive\code\test\image\Fig0505(b)(ckt_salt_only).tif');
gch = spfilt(gs,'chmean',3,3,-1.5); % 椒噪声用正的Q;
gar = spfilt(gs,'amean',3,3); 
gme = spfilt(gs,'gmean',3,3);
myImshow(gs,'盐噪');myImshow(gch,'逆调和除噪');myImshow(gar,'均值滤波');
myImshow(gme,'几何均值滤波');
myImshow(0);

自适应噪声滤波器

% 中值滤波与自适应中值滤波
gsp = imread('E:\资料\onedrive\code\test\image\Fig0506(a)(ckt_salt_pep_pt25).tif');%Fig0506(a)(ckt_salt_pep_pt25).tif
gm2 = medfilt2(gsp);% 中值滤波
gch1 = spfilt(gsp,'chmean',3,3,1.5); % 椒噪声用正的Q;
gch2 = spfilt(gsp,'chmean',3,3,-1.5); % 盐噪声用负的Q;
gch3 = spfilt(gch1,'chmean',3,3,-1.5); % 除椒后除盐
gadp = adpmedian(gsp,7);
myImshow(gsp,'椒盐噪声');myImshow(gm2,'中值滤波');myImshow(gadp,'自适应中值');
myImshow(gch1,'除椒');myImshow(gch2,'除盐');myImshow(gch3,'除椒后除盐');
myImshow(0);

运动模型的退化函数

运动模型

在快门曝光的瞬间内，取景区与镜头发生了位移，使得原始图片f(x,y)变成了g(x,y)

基础函数

function B = pixeldup(A, m, n)
%PIXELDUP Duplicates pixels of an image in both directions.
%   B = PIXELDUP(A, M, N) duplicates each pixel of A M times in the
%   vertical direction and N times in the horizontal direction.
%   Parameters M and N must be integers.  If N is not included, it
%   defaults to M.
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/06/14 16:29:54 $
% Check inputs.
if nargin < 2 
   error('At least two inputs are required.'); 
end
if nargin == 2 
   n = m; 
end
% Generate a vector with elements 1:size(A, 1).
u = 1:size(A, 1);
% Duplicate each element of the vector m times.
m = round(m); % Protect against nonintergers.
u = u(ones(1, m), :);
u = u(:);
% Now repeat for the other direction.
v = 1:size(A, 2);
n = round(n);
v = v(ones(1, n), :);
v = v(:);
B = A(u, v);

u = u(ones(1, m), :);这代码可以推广到一般的形式A(a,b);其中A为矩阵，a、b为向量，则a表示从A中取元素时的所取的行下标，而b则表示列下标，比如
A=[1,2,3;3,4,5] 则A([1,2],[2,3])=[2,3;4,5]
故最终结果有d(a)*d(b)个;d(a)表示a的维数。

模拟运动模型

clear all;
f=checkerboard(8);%生成一个8格的板子
g=pixeldup(f,8);
PSF = fspecial('motion',7,45);
gb=imfilter(f,PSF,'circular');
noise=imnoise(zeros(size(f)),'gaussian',0,0.001);
gbn=gb+noise;
myImshow(f,'拓展图');myImshow(gb,'运动后图像');myImshow(noise,'噪声图像');
myImshow(gbn,'运动后噪声图像');myImshow(0);

维纳滤波（最小均方误差误差）

用来做图像复原。其中H(u,v)是降指模型的点扩散函数

　
　　
第一部分完

数字图像处理课程笔记

一些图像的像素范围映射变化

灰度平方运算

灰度区间映射

灰度直方图均衡化

log变换

相对拉伸变换

画出自适应的坐标轴

空域滤波

基础知识

线性滤波函数

非线性空域滤波

频域滤波

一个例子

显示的优化

运行示例

运行结果

相关参考

傅立叶变换滤波

滤波步骤

从空域滤波生成频域滤波

下面测试几个基础函数的使用效果

paddedsize函数实现

gscale函数的实现

dftfilt函数

空域生成频域的例子

运行结果

直接生成频域滤波

dftuv的实现

直接用频域高斯低通滤波的例子

运行结果

直接生成频域高通滤波

用频域波特沃兹高通滤波的例子

运行结果

用高通滤波进行图像锐化

运行结果

图像恢复

图像降指与恢复过程

噪声模型

imnoise函数生成噪声

type:

parameters

生成特定分布的空域随机噪声

只存在空域噪声的恢复

均值滤波器

统计排序与调和均值滤波器

空域滤波基础函数Spfilt

逆谐波函数用于消除椒盐噪声

自适应噪声滤波器

运动模型的退化函数

运动模型

基础函数

模拟运动模型

维纳滤波（最小均方误差误差）

内容目录

选择主题

空域滤波　　