前面曾经就介绍过图像处理的一些知识了，这里就MATLAB来介绍一些图像处理的案例。今天算是填一个坑。

图像都在https://gitee.com/lcl1026504480/MATLAB-imgprocess/tree/master/im/%E5%AE%9E%E9%AA%8C%E5%9B%BE%E5%83%8F/%E5%AE%9E%E9%AA%8C%E5%9B%BE%E5%83%8F

里面。

图像处理基本操作与灰度变换

1. 基础知识

1) Matlab支持的图像格式

TIFF, JEPG, GIF, BMP, PNG, XWD (X Window Dump)，其中GIF不支持写。

2) 与图像处理相关的最基本函数

读：imread; 写：imwrite; 显示：imshow; 信息查看：imfinfo;

3) Matlab支持的数据类型

double, unit8, int8, uint16, int16, uint32, int32, single, char (2 bytes per element), logical.

4) Matlab支持的图像类型

Intensity images, binary images, indexed images, RGB image

2. 图像信息显示与图像读写

>> imfinfo('parrots.bmp')%%图像文件信息显示

Filename: 'parrots.bmp'

FileModDate: '27-Sep-2010 19:42:59'

FileSize: 1179702

Format: 'bmp'

FormatVersion: 'Version 3 (Microsoft Windows 3.x)'

Width: 768

Height: 512

BitDepth: 24

ColorType: 'truecolor'

……

这是个彩色图像。

>> imfinfo('pollen.tif')%%图像文件信息显示

Filename: 'pollen.tif'

FileModDate: '27-Sep-2010 19:42:59'

FileSize: 250398

Format: 'tif'

FormatVersion: []

Width: 500

Height: 500

BitDepth: 8

ColorType: 'grayscale'

……

这个是灰度图像。

>> img=imread('parrots.bmp');%%读入图像文件到内存变量

>> imshow(img);%%图像显示

>> whos img%%显示变量信息

Name      Size                           Bytes  Class

img     512x768x3                      1179648  uint8 array

Grand total is 1179648 elements using 1179648 bytes

>> imwrite(img,'NewParrots.png');%%以其它格式另外保存图像

>> img2=imread('pollen.tif');

>> whos img2

Name       Size                    Bytes  Class

img2     500x500                  250000  uint8 array

Grand total is 250000 elements using 250000 bytes

>> figure, imshow(img2);

3. 直方图显示与均衡化

>> im=imread('pollen.tif');%%读入图像文件

>> size(im)%%观察2维还是3维

ans =%%如果是3维，使用im=im(:,:,1);转成2维

500   500

>> subplot(221), imshow(im)%%显示原图像

>> subplot(222), imhist(im)%%显示直方图

>> im2=histeq(im);%%直方图均衡化

>> subplot(223), imshow(im2)%%显示均衡化后的图像

>> subplot(224), imhist(im2)%%显示均衡化后的图像的直方图

4. 灰度变换

函数形式：J = imadjust(I,[low_in; high_in],[low_out; high_out], Gamma)

参数 Gamma用于指定变换函数形状，小于1时，输出更亮，大于1时，输出更暗。

>> srcimage=imread('chest_xray.tif');%%读入图像文件

>> g1=imadjust(srcimage, [0 1], [0 1], 0.5);%%灰度变换

>> g2=imadjust(srcimage, [0 1], [0 1], 2.5);

>> subplot(1, 3, 1); imshow(srcimage); title ('Source Image');

>> subplot(1, 3, 2); imshow(g1); title ('Gamma=0.5');

>> subplot(1, 3, 3); imshow(g2); title ('Gamma=2.5');

实验

1. 显示不同图像及其直方图，选择chest_xray.tif, pollen.tif和Einstein.png等图像进行直方图均衡化，并显示结果图像和直方图。记录其中一幅图像的原图像、直方图和均衡化后的结果图像与直方图。

2. 选择不同的参数应用灰度变换方法对图像Einstein. png进行增强，观察效果,记录下满意的参数和结果图像。

3. 对图像进行Gamma 变换增强，观察不同Gamma值对增强结果的影响，并分析原因。

原图像及直方图：

直方图均衡化之后的图像及直方图：

显然提高了对比度。

增强的Einstein.png图像参数和结果

命令语句及参数：

>> e=imread('Einstein.PNG');

>> e1=imadjust(e,[0 1],[0 1],0.5);

>> subplot(121),imshow(e);subplot(122),imshow(e1)

结果图像：

不同Gamma值对应的图像增强结果

总结灰度变换中Gamma值对于结果的影响和原因：

γ>1的时候，增强后的图像会比原图像更暗，可以使高灰度值区域分级更明显，低灰度之间的分层基本消失，反之，输出图像会更亮，可以更好的区分低灰度的层次，而使较高灰度的层次几乎消失。从下图可以很好的理解，γ>1的时候输出的灰度值在输入灰度的下方，所以更暗，还可以看到在输入灰度比较高的地方，输出灰度以很陡的斜率增长，而在输入灰度级比较低的地方，斜率很小，因此，低灰度会出现一定程度的输出合并，而高灰度会有输出分层更细的现象。

图像处理应用

读入图像 Problem1.bmp，输出由姓名的拼音和学号组成的结果图像，并把该结果图像缩小后做为水印加到一幅灰度图像的左上角

这个代码是https://gitee.com/lcl1026504480/MATLAB-imgprocess/blob/master/problem1.m

代码执行过程：

首先输入学号和姓名的小写字母，比如150210116和lcl。

然后选择要添加水印的图片：

然后选择添加水印的左上角，程序自动添加水印：

然后选择添加水印之后文件的保存位置，就可以保存了。

频域图像处理

1. 基础知识

1) 二维离散Fourier变换对

设图像f(x,y)的大小为(M,N),其二维离散Fourier变换对

2) 二维离散Fourier谱的计算与可视化

为了便于观察，一般取Fourier谱的一个整周期，并且平移到坐标轴的正向，如下图所示。

二维离散Fourier变换的FFT计算函数：

F = fft2 (f, P, Q)

F: 输出大小为P x Q, 通过对输入图像f添加需要数目的0实现

二维离散Fourier变换的Fourier谱计算

S = abs ( F )

二维离散Fourier变换的原点移动

Fc = fftshift (F)

Log变换

S2=log(1+abs(Fc))

2. Fourier谱的计算与显示

>>close all

>>clear all

>> im=imread('Xraychest.png');%% 读入图像文件

>> imdouble=double(im);%% fft2函数不处理uint8类型数据，转换成double

>> F=fft2(imdouble);%% 2D FFT计算

>> Fc=fftshift(F);%% 平移

>> S2=log(1+abs(Fc));%% Log变换

>> subplot(121);imshow(im);%% 显示原图像

>> subplot(122);imshow(S2,[]);%% 显示Fourier谱

3.Fourier逆变换计算

>> F = ifftshift(Fc);%% 对于Fc要使用平移复原， 对于F则不需要平移

>> f=real(ifft2(F));%% 2D FFT逆变换

>> imifft=uint8(f);%% 转换成uint8类型

>> figure, imshow(imifft);%% 显示图像

4.频域滤波

频域滤波基本步骤：

输入图像->图像预处理->Fourier变换->滤波-> Fourier反变换->图像后处理->输出图像

F = fft2(f, size(H, 1), size(H, 2));

g = real(ifft2(H.*F));

g = g(1:size(f, 1), 1:size(f, 2));

5. 空域和频域滤波转换

%显示原图像及其Fourier谱

close all

clear all

f = imread('bld.tif');

F = fft2(double(f));

S = fftshift(log(1+abs(F)));

figure (1)

subplot(221);imshow(f);title('Original image');

subplot(222);imshow(S, []);title(' Fourier spectrum');

%生成空域和频域滤波器

h = fspecial('sobel');%%使用Sobel滤波器

PQ = paddedsize(size(f));%%自编函数，计算填0后的图像大小

H = freqz2(h, PQ(1), PQ(2));%%生成Sobel算子的频率滤波器

H1 = ifftshift(H);%%频率滤波器原点移位到左上角

figure(2);

subplot(121);

mesh(abs(H(1:20:end, 1:20:end)));         %%绘制频率滤波器

axis off;

title('perspective plot of’ frequency filter');

subplot(122);

mesh(abs(H1(1:20:end, 1:20:end)));        %%绘制移位后的频率滤波器

axis off;

title('frequency filter with shift');

% 分别在空域和频域滤波

gs = imfilter(double(f), h);%%空域滤波

F = fft2(double(f), size(H1, 1), size(H1, 2));%%频域滤波

gf = real(ifft2(H1.*F));

gf = gf(1:size(f, 1), 1:size(f, 2));

figure(1);

subplot(223);imshow(gs, []);

title('Filtered image in spatial domain');

subplot(224);imshow(gf, []);

title('Filtered image in frequency domain');

6. Gaussian低通和高通频域滤波

Gaussian低通滤波器的传递函数为：

Gaussian低通频域滤波

close all

clear all

f = imread('FigMark.tif');

PQ = paddedsize(size(f));%这个函数就是把图像的长和宽都变成2倍，为了防止缠绕

[U, V] = dftuv(PQ(1), PQ(2));%%生成距离计算需要的网格

D0 = 0.05*PQ(2);                  %% D0取宽度的5%,根据需要调整

F = fft2(double(f), PQ(1), PQ(2));

H = exp(-(U.^2+V.^2) / (2*(D0^2)));  %%滤波器计算

g = real(ifft2(H.*F));               %%频域滤波

g = g(1:size(f, 1), 1:size(f, 2));

figure                            %%结果显示

subplot(221); imshow(f, []);

title('Original image');

subplot(222); imshow(log(1 + abs(fftshift(F))), []);

title(' Fourier spectrum of original image');

subplot(223); imshow(g, []);

title('Filtered image by Gaussian lowpass filter');

F2=fft2(double(g), PQ(1), PQ(2));

subplot(224); imshow(log(1 + abs(fftshift(F2))), []);

title(' Fourier spectrum of lowpass filtered image');

Gaussian高通频域滤波

H2 =1-H;%%高通滤波器计算

g2 = real(ifft2(H2.*F));               %%频域滤波

g2 = g2(1:size(f, 1), 1:size(f, 2));

figure%%结果显示

subplot(221); imshow(f, []);

title('Original image');

subplot(222); imshow(log(1 + abs(fftshift(F))), []);

title(' Fourier spectrum of original image');

subplot(223); imshow(g2, []);

title('Filtered image by Gaussian highpass filter');

FH2=fft2(double(g2), PQ(1), PQ(2));

subplot(224); imshow(log(1 + abs(fftshift(FH2))), []);

title(' Fourier spectrum ofhighpass filtered image');

原图像和Fourier谱显示结果：

傅里叶谱显示的结果是经过fftshift和对数变换之后的，为了方便看。

傅里叶再反变换之后：

虽然人看不出来什么不同，但是由于计算精度的原因，灰度值可能会有一些细微的改变。

2. 原图像、空域和频域的滤波结果图像

空域滤波：

频域滤波：

Sobel频域滤波器如下图：

频域滤波结果：

计算空域和频域滤波结果图像之差，并求最大的差值，分析产生差值的原因：

>> max(max(gf-gs))

ans =

4.5475e-13

这是因为计算误差，傅里叶反变换回来的值会有很小的虚部，我们舍弃了这些虚部。

3. 原图像和高通、低通滤波结果图像

参数

25.6

标准差越大，对低通滤波效果越差，这是因为标准差越大，高斯曲线越“胖”，高频衰减越少，低通效果越差，但是对高通滤波正好相反。

图像处理应用

条纹噪声在频域滤波的效果较好。读入含条纹噪声图像problem2.png,恢复得到不含条纹噪声的图像。

代码执行过程：

原图像：

二维傅里叶变换之后：

然后用鼠标点击四个非中心的亮点，代码会自动把这4个点的位置参数传进滤波函数。滤波之后结果：

看到边界还会有一些条纹，这个应该和边界扩展的方式有一定关系。而且图像变得比之前模糊，这很正常，因为用了滤波器。

图像的恢复与分割

1. 基础知识

1) 图像恢复

设真实图像为f(x,y)，其退化和恢复模型如下图所示。

这个过程也叫做维纳滤波，特例是逆滤波。

图像恢复就是利用某种先验知识对退化过程建模，并在某个质量准则下，应用其逆过程来得到原图像的最优估计。实际的恢复过程相当于设计一个滤波器，使其能从降质图像计算得到真实图像的估值，根据预先规定的误差准则，最大程度接近真实图像。

2) 噪声模型

对噪声的行为和效果进行仿真的能力是图像恢复处理中的重点。除了周期噪声，通常假定噪声与像素坐标无关。

Matlab中噪声仿真函数为：g = imnoise(f, type, parameters)，常见的高斯、椒盐噪声的仿真实例如下：

g = imnoise(f, ‘gaussian’, m, var);

g = imnoise(f,’salt & pepper’, d);  (其中 d 代表噪声密度)

Matlab中两个常用的随机分布产生函数是：rand(均匀分布噪声), randn(标准正态分布噪声)。

常用的噪声模型包括：

Gaussian噪声分布：

Rayleigh噪声分布：

脉冲噪声分布：

选择图像区域可以通过Matlab函数roipoly进行：

BW = roipoly(I,c,r): 返回图像I中由c, r确定的区域，结果是一个二值图像，其大小和图像I相同，区域内为1，区域外为0。参数c和r保存组成该区域多边形的各个顶点的坐标对。

[BW, x, y]=roipoly(I); 允许用户交互式地通过鼠标选择图像I的多边形区域，结果返回多边形坐标对x,y以及二值区域BW。

均值和方差的估计可以通过自编函数statmoments完成。

对于逆谐波均值滤波器, type设置为'chmean', 同时还应该加上参数Q, 用法如下

f=spfilt (img, type, m, n, Q)

2. 噪声的产生和分布的显示

噪声可以直接用自编函数imnoise2来产生，分布可以使用直方图进行观察。

g=imnoise2 (type, M, N, a, b )

Type — 'uniform', ‘gaussian‘, 'salt & pepper', ‘rayleigh', ‘exponential', ‘erlang‘

M, N — 数组大小

a, b — 对应type的参数

g — 返回的噪声数组

例如：

R=imnoise2('gaussian',500,500,0,1); %%产生Gaussian噪声 500x500像素, 均值为0, 方差为1

figure, hist(R(:), 50);  %% 显示直方图 50为统计的区间个数

title('Gaussian Noise');

R2=imnoise2('rayleigh', 100000, 1); %% 产生Rayleigh噪声

figure, hist(R2(:), 50);

title('Rayleigh Noise');

3. 仅有噪声影响图像的恢复

噪声可以直接用自编函数imnoise2来产生，空域降噪可以采用函数spfilt完成。

close all

clear all

im=imread('bld.tif');

subplot(231);imshow(im); title ('Source image');

%% 计算噪声图像

imnoise=imnoise2('gaussian',size(im,1), size(im,2),0,0.1); %% 噪声参数设置

im2=double(im)+255.0*imnoise;      %% 加性噪声影响

subplot(232);imshow(im2,[]);  title ('Noise image');

%% 进行空域滤波

f=spfilt (im2,'amean' , 3, 3);  %% 算术均值滤波

subplot(233);imshow(f,[]);  title ('Arithmetic mean');

f=spfilt (im2,'gmean' , 3, 3);  %% 几何均值滤波

subplot(234);imshow(f,[]);  title ('Geometric mean');

f=spfilt (uint8(im2),'hmean' , 3, 3);  %% :谐波均值滤波

subplot(235);imshow(f,[]);  title ('Harmonic mean');

f=spfilt (im2,'median' , 3, 3);  %% 中值滤波

subplot(236);imshow(f,[]);  title ('Median');

4. 含椒盐噪声图像的去噪

椒盐噪声中椒(黑色)和盐(白色)可以分别设置强度(即在图像中的占的比例，注意椒和盐的强度之和应该小于1)，也可以仅出现椒噪声或者盐噪声。中值滤波和逆谐波滤波都经常用于椒盐噪声的消除。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  产生椒盐噪声图像的方法：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all

clear all

f=imread('cktboard.tif');

%%产生椒盐噪声

[M, N]=size(f);

R=imnoise2('salt & pepper',M, N, 0.1, 0.1);  %% 噪声强度 20%--Peppers  20%--Salt

c = find(R==0);

g2 = f;

g2(c)=0;

c = find(R==1);

g2(c)=255;

%%产生椒噪声

R=imnoise2('salt & pepper',M, N, 0.2, 0);    %%噪声强度30%--Peppers

c = find(R==0);

gp = f;

gp(c)=0;

%%产生盐噪声

R=imnoise2('salt & pepper',M, N, 0,0.2);      %%噪声强度30%--Salt

c = find(R==1);

gs = f;

gs(c)=255;

figure(1)

subplot(241); imshow(f); title('source');

subplot(242); imshow(g2); title('salt & pepper');

subplot(243); imshow(gp); title('pepper');

subplot(244); imshow(gs); title('salt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  中值滤波降噪

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f1=spfilt (g2,'median' , 3, 3);  %% 中值滤波,窗口大小3x3

f2=spfilt (gp,'median' , 3, 3);

f3=spfilt (gs,'median' , 3, 3);

subplot(246); imshow(f1);  title ('Median Denoise Salt&Pepper');

subplot(247); imshow(f2);  title ('Median Denoise Pepper');

subplot(248); imshow(f3);  title ('Median Denoise Salt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  不同Q值的逆谐波均值滤波降噪

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)

subplot(241); imshow(f); title('source');

subplot(242); imshow(g2); title('salt & pepper');

subplot(243); imshow(gp); title('pepper');

subplot(244); imshow(gs); title('salt');

f4=spfilt (g2,'chmean' , 3, 3, 0.2);  %% 逆谐波滤波,窗口大小3x3,Q=0.2

f5=spfilt (gp,'chmean' , 3, 3, 0.2);

f6=spfilt (gs,'chmean' , 3, 3, 0.2);

subplot(246); imshow(f4);  title ('Chmean Q=0.2 Denoise Salt&Pepper');

subplot(247); imshow(f5);  title ('Chmean Q=0.2 Denoise Pepper');

subplot(248); imshow(f6);  title ('Chmean Q=0.2 Denoise Salt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)

subplot(241); imshow(f); title('source');

subplot(242); imshow(g2); title('salt & pepper');

subplot(243); imshow(gp); title('pepper');

subplot(244); imshow(gs); title('salt');

f7=spfilt (g2,'chmean' , 3, 3, -4.1);  %% 逆谐波滤波,窗口大小3x3,Q=-4.1

f8=spfilt (gp,'chmean' , 3, 3, -4.1);

f9=spfilt (gs,'chmean' , 3, 3, -4.1);

subplot(246); imshow(f7);  title ('Chmean Q=-4.1 Denoise Salt&Pepper');

subplot(247); imshow(f8);  title ('Chmean Q=-4.1 Denoise Pepper');

subplot(248); imshow(f9);  title ('Chmean Q=-4.1 Denoise Salt');

实验方案与步骤

1．分别产生Gaussian噪声和Rayleigh噪声，显示和观察其分布。

2．在图像中加入Gaussian噪声，然后对该噪声图像分别应用算术均值滤波、几何均值滤波、谐波均值滤波和中值滤波进行恢复，观察结果。

3．在图像中分别加入椒盐噪声、仅含椒和仅含盐的噪声，对噪声图像应用不同Q参数的逆谐波均值滤波和中值滤波进行恢复，观察结果。

4. 讨论中值滤波的窗口大小以及椒盐噪声的强度对于图像去噪效果的影响

产生Gaussian噪声 500x500像素, 均值为0, 方差为1

产生Rayleigh噪声100000×1像素，a=0，b=1。

均值为0，标准差为0.1的高斯噪声。

结果图像：

用的都是3×3的模板。

效果都不是太好。换成11×11的模板之后：

去噪效果是好了，但是图像变得太模糊了，高斯噪声还是建议使用高斯双边滤波器来滤除。

椒盐噪声参数：

噪声强度 20%--Peppers  20%--Salt

结果图像：

滤波核大小均为3×3。

仅含椒噪声参数：

噪声强度30%--Peppers

结果图像：

仅含盐噪声参数：

噪声强度30%--Salt

结果图像：

中值滤波对椒盐噪声的效果都还不错，但是一个3×3的模板太小了，效果不是非常好，模板更大效果会更好。Q>0的逆谐波均值滤波器对椒噪声非常有效，会加剧盐噪声。Q<0对盐噪声非常有效，但是也会加剧椒噪声。

中值滤波不同窗口参数以及椒盐噪声不同强度的对比图像

结论：

椒盐强度一样的时候，中值滤波窗口越大，滤除椒盐噪声的效果越好；中值滤波窗口不变的时候，椒盐强度越大，滤波效果越差。

图像处理应用

如下面左图，六点盲文是通过2列x3行排列的点进行顺序编码的，当有点出现时编码为1，否则编码为0。读出下面右图(problem3.png)，编程确定图中9个字符的0/1编码内容

代码执行过程：

首先用霍夫圆变换把圆检测出来：

然后二值化并把圆进行填充

然后进行行和列投影，最后for循环判断每个位置是否是白色，最后输出编码：

out =

3×18 char 数组

'010101110101010000'

'101111100010011011'

'011001111101010111'

也可以其它的思路：

最后是三个工业检测任务：

第一个是检测药的包装，

下面是正常的：

下面这个空了一个。

下面是包装有问题。

希望编写程序最后的结果是而已像下面框出出问题的地方：

这个说实话挺难的。

第二个：

检测这种产品的质量：

这种就是有缝了：

希望识别出缝：

1.总体思路

把原图像转成灰度图像，然后利用OSTU得到二值图像，然后取反，用区域增长来识别图像中的区域。如果有断裂的地方，那么它的面积会比正常的网格大，这样可以根据面积来找出不合格区域。

这个代码在https://gitee.com/lcl1026504480/MATLAB-imgprocess

是aaa.m，里面也有注释。

下面分别是4张图片的结果：

任务三

检测饼干的缺或者裂缝。

这个我也觉得比较难。交给你们了。

任务四

表面平整度检测。

总体思路

图片上看，本身就有轻微的凹凸不平，这会造成灰度值的变化，图片就像加了噪声一样，这会对边缘检测带来一定的困难。所以，先进行均值滤波破坏掉一些由于这些正常的轻微的凹凸不平带来的边缘，而边缘比较强的，也就是不合格的边缘，是不容易被破坏的。用sobel算子来进行边缘检测，选择合适的阈值，把图像转化为二值图像，然后再利用区域增长的方式检测目标，获得目标中心点，以一个固定半径在原图基础上标示出不合格区域。

代码还是在aaa.m：

这里面需要注意的是，在用sobel卷积核做卷积的时候，边界扩展不能是默认的，因为默认的是0来填充，而如果用0来填充，图片原来的边界的卷积结果会很大，影响后面的步骤，所以使用了镜像对称。下面分别是6张图片的结果：

由于这个二值化结果比较难以看到白色的点，所以我单独把它列出来一张图片：

总结与体会

图像处理其实很难找到比较通用的办法，因为可以干扰的因素太多了，噪声也都是随机的，现在图像识别比较常用的是利用卷积神经网络去学习的方法，不过针对这几道题目，数据集太小了，不适合用机器学习去解决，那么需要我们人工去想办法。虽然通用的办法很难找，但是对于某种比较固定的环境下，我们还是可以找到针对具体问题的处理办法的，不过这其中处理过程中的阈值需要自己事先去获得或者不断的尝试改进，直到达到我们的目标。题目2就是事先去获得，这个小矩形的面积我们是可以得到一个范围的，那么面积的阈值就容易设定了，而题目4是需要不断去尝试新的阈值。对于题目2来说，上面的方法对检测断裂通用性还是比较好的，而对于题目4，上面的方法的通用性还是有疑问的，因为这个阈值仅仅是对这6张图有用，这个阈值又是试出来的，不像题目2是一个先验知识。总之图像处理虽然理论也很多，不过很多东西还是需要在实践中去不断尝试不断修改才能确定的。总得来说还是要具体问题具体分析，在实践中不断尝试。