【数字图像处理课程设计】期中、期末综合考试题目整理总结（共四个图像处理算法应用题）_试分别计算i与a和b两个模板的滤波结果

目录

一、下面两幅图像中有几处不同，编程把它们找出来、并在图中突出显示（关键步骤不能调用内置函数）。

1.算法原理

2.解题步骤

3.程序代码

4.处理结果

二、下图含有干扰条纹（moiré pattern）、并且低灰度区域的细节难以分辨。对图像进行灰度变换，再分别应用空域和频域的处理方法消除条纹干扰，并比较这两类方法去除条纹干扰的效果。

 1.算法原理

2.解题步骤 

 3.程序代码

 4.处理结果

三、下图是一张单据的扫描件，设计算法将表格横线校正为水平、并补全断裂表格线。

1.相关原理

 2.程序源码

 3.结果分析

四、有4种飞机模型，每种5架。请分别用傅里叶描述子和矩不变量提取特征，每种飞机中选3架进行训练，另外2架进行测试，用最近邻分类器进行识别。阐述所用算法的原理、关键步骤的处理结果、最终结果等，要附程序代码。

 1.所用算法原理

2.关键程序源码

 3.关键步骤的处理结果

4.最终结果

---

一、下面两幅图像中有几处不同，编程把它们找出来、并在图中突出显示（关键步骤不能调用内置函数）。

  

1.算法原理

（1）差影法：

        差影法实际是对图像进行代数运算的一种不同的叫法。

        代数运算是指对两幅输入图像进行点对点的加、减、乘、除计算得到输出图像的运算。

        图像相加的一个重要应用就是对同一场景的多幅图像求平均值,它可以有效地降低随机噪声的影响。这是因为对于一幅有噪声的图像S (x，y)，可以看成是由原始无噪声的图像F(x，y）和噪声G(x，y）叠加而成的，即

        如果叠加在图像上的噪声G(x，y）是非相关、具有零均值的随机噪声时，那么，把针对同一目标物在相同条件下，做M次重复摄取的图像相加，取平均值作为输出图像，即:

 

        图像相减:对于同一场景的两幅图S (x，y)、F(x，y），因为是同一场景，所以目标图像的背景是大致相同的，当两幅图进行相减时，相同位置的背景点因为灰度值相同，结果变成了黑点，而目标图像因为位置不同，所以相减后值不为0，处理后的图像就只在两个目标图像的位置有像素点，背景变为全黑。此时，我们计算两个目标之间的距离就非常简单了，这也就达到了检测物体运动的目的。

        乘法和除法在数字图像处理中一般应用得不多，但它们的用途也很重要。

        数字化对一幅图像各点的敏感程度可能有变化,乘和除运算有可能纠正这种影响。除运算可产生对颜色和多光谱图像分析十分重要的比率图像，而用一幅图像乘某一图像可以遮住该图像中的某些部分，仅留下感兴趣的物体。

（2）形态学图像处理：

        形态学，即数学形态学(mathematical Morphology），是图像处理中应用最为广泛的技术之一，主要用于从图像中提取对表达和描绘区域形状有意义的图像分量，使后续的识别工作能够抓住目标对象最为本质〈最具区分能力－most discriminative）的形状特征，如边界和连通区域等。同时像细化、像素化和修剪毛刺等技术也常应用于图像的预处理和后处理中，成为图像增强技术的有力补充。

        在数字图像处理中， 形态学是借助集合论的语言来描述的，其中首先较为重要的就是结构元素。

        结构元素（structure element)——设有两幅图像A, S。若A是被处理的对象， 而S是用来处理A的， 则称S为结构元素。结构元素通常都是一些比较小的图像， A与S的关系类似于滤波中图像和模板的关系。

形态学图像处理主要包括以下内容：

        a). 二值图像的基本形态学运算， 包括腐蚀、膨胀、开和闭。

        由于所有形态学运算都是针对图像中的前景物体进行的， 因而首先对图像前景和背景的认定给出必要的说明．

        注意：大多数图像，一般相对于背景而言物体的颜色（灰度）更深，二值化之后物体会成为黑色，而背景则成为白色，因此我们通常是习惯于将物体用黑色（灰度值0）表示，而背景用白色（灰度值255）表示，本章所有的算法示意图以及所有的Visual C＋＋的程序实例都遵从这种约定；但Matlab 在二位图像形态学处理中，默认情况下白色的（二位图像中灰度值为1的像素，或灰度图像中灰度值为255的像素）。

        是前景（物体），黑色的为背景，因而本次报告涉及Matlab 的所有程序实例又都遵从Matlab本身的这种前景认定习惯。

        实际上，无论以什么灰度值为前景和背景都只是一种处理上的习惯，与形态学算法本身无关。只需要在形态学处理之前先对图像反色就可以在两种认定习惯之间自由切换。

        b). 二值形态学的经典应用， 包括击中击不中变换、边界提取和跟踪、区域填充、提取连通分量、细化和像素化，以及凸壳。

        c). 灰度图像的形态学运算， 包括灰度腐蚀、灰度膨胀、灰度开和灰度闭。

（3）图像分割：

        图像分割（Segmentation）指的是将数字图像细分为多个图像子区域（像素的集合）（也被称作超像素）的过程，就是把图像分成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程。它是由图像处理到图像分析的关键步骤。图像分割的目的是简化或改变图像的表示形式，使得图像更容易理解和分析。图像分割通常用于定位图像中的物体和边界（线，曲线等）。更精确的，图像分割是对图像中的每个像素加标签的一个过程，这一过程使得具有相同标签的像素具有某种共同视觉特性。

        图像分割的结果是图像上子区域的集合（这些子区域的全体覆盖了整个图像），或是从图像中提取的轮廓线的集合（例如边缘检测）。一个子区域中的每个像素在某种特性的度量下或是由计算得出的特性都是相似的，例如颜色、亮度、纹理。邻接区域在某种特性的度量下有很大的不同。

     现有的图像分割方法主要分以下几类：基于阈值的分割、区域生长、区域分裂合并、分水岭算法、边缘分割（边缘检测）、直方图法、聚类分析、小波变换等。

        a). 大津法(OTSU)是一种确定图像分割阈值的算法，由日本学者大津于1979年提出;原理上来讲，该方法又称作最大类间方差法，有时也称之为大津算法;其按照大津法求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大；其被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。

        b). Otsu算法原理：

Otsu算法使用的是聚类的思想：

①  把图像的灰度数按灰度级分成2个部分，使得两个部分之间的灰度值差异最大，每个部分之间的灰度差异最小；

②  通过方差的计算来寻找一个合适的灰度级别来划分；

对图像I(x,y)，前景和背景的分割阈值记作T，

前景像素点数占整幅图像的比例为ω0，其平均灰度μ0；

背景像素点数占整幅图像的比例为ω1，其平均灰度μ1；

图像的总平均灰度记为μ=ω0∗μ0 +ω1∗μ1；

类间方差记为g。

假设：

背景较暗，且图像的大小为M×N

图像中，

像素灰度值小于阈值T的像素个数记作N0，

像素灰度值大于阈值T的像素个数记作N1，

则有：

ω0 = N0 / M×N ；      (1)

ω1 = N1 / M×N ；      (2)

N0 + N1 = M×N ；       (3)

ω0+ω1=1；           (4)

μ=ω0*μ0+ω1*μ1；  (5)

g =ω0(μ0-μ)^2+ω1(μ1-μ)^2；(6)

将式(5)代入式(6),得到等价公式:

g=ω0*ω1*(μ0-μ1)^2; 　(7)　这就是类间方差

采用遍历的方法得到使类间方差g最大的阈值T,即为所求。

        c). 特点：

它是按图像的灰度特性,将图像分成背景和前景两部分；

因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小；因此,使类间方差最大的分割意味着错分概率最小。

（4）SURF特征提取算法：

        SURF，英文的全称为Speed Up Robust Features,直译为:加速版的具有鲁棒特性的特征算法，是由Bay在2006年首次提出的。该算法对经典的尺度不变特征变换算法(SIFT算法)进行了改进，SIFT算法最大的缺点就是如果不借助硬件或者专门的图像处理器进行加速的话，SIFT算法很难达到实时处理的效果。SURF算法最大的特征在于采用了Haar特征以及积分图像的概念，这大大的加速了程序的运行时间。SURF算法不仅保持了SIFT算法的尺度不变和旋转不变的特性，而且对光照变化和放射变化同样具有很强的鲁棒性。SURF算法一般应用于计算机视觉中的物体识别、图像拼接、图像配准以及3D重建中。

        SIFT算法实现物体识别主要有三大工序，1、提取关键点；2、对关键点附加详细的信息（局部特征）也就是所谓的描述器；3、通过两方特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，也就建立了景物间的对应关系。

2.解题步骤

（1）读取图像，并转为灰度图。

         

 （2）由上面右图可以发现，图像并不是完全重合的，因此需要对齐矫正，这里采用了surf特征提取算法，对图像进行特征匹配，然后图像配准，最终达到对齐图像的。

（3）应用差影法对两幅图像做减法，得到两个图像的绝对差。

                   

（4）对图像差进行灰度增强，然后利用Otsu阈值分割算法和形态学处理对图像差进行进一步处理。

 右图是放大之后的细节

 5）由第4步第一个图可以看出边角有白边，所以需要边角置黑。

（6）找到图像中连通的区域，并在原图基础上绘制矩形框把不同的地方框选出来。

3.程序代码

Align.m：

function Image_Align = Align(MOVING,FIXED)  %此函数用来给图像进行配准
% 参数 moving 是要配准的图像B,参数 fixed 是固定的参考图像A, 
 
% type是使用的模式方法
% type 'no' 不使用图像注册
% type 'auto' 使用自动图像配准。
% type 'user' 使用用户定义的控制点进行注册
% 注册图像使用注册估计器应用程序自动生成的代码注册灰度图像。
%  [MOVINGREG] = registerImages(MOVING,FIXED) 寄存灰度图像
%  MOVING and FIXED 使用自动生成的代码从注册
%  Estimator App.估算器应用程序，已设置所有注册参数的值
%  在应用程序中以交互方式显示并生成存储
        % 基于特征的技术需要计算机视觉系统工具箱的许可
        checkLicense();
 
        % 默认空间引用对象
        fixedRefObj = imref2d(size(FIXED));
        movingRefObj = imref2d(size(MOVING));
 
        % 检测SURF特征并返回SURFPoints对象
        fixedPoints = detectSURFFeatures(FIXED,'MetricThreshold',750.000000,'NumOctaves',3,'NumScaleLevels',5);
        movingPoints = detectSURFFeatures(MOVING,'MetricThreshold',750.000000,'NumOctaves',3,'NumScaleLevels',5);
 
        % 特征提取
        [fixedFeatures,fixedValidPoints] = extractFeatures(FIXED,fixedPoints,'Upright',false);
        [movingFeatures,movingValidPoints] = extractFeatures(MOVING,movingPoints,'Upright',false);
        % 特征匹配
        indexPairs = matchFeatures(fixedFeatures,movingFeatures,'MatchThreshold',50.000000,'MaxRatio',0.500000);
        fixedMatchedPoints = fixedValidPoints(indexPairs(:,1));
        movingMatchedPoints = movingValidPoints(indexPairs(:,2));
        MOVINGREG.FixedMatchedFeatures = fixedMatchedPoints;
        MOVINGREG.MovingMatchedFeatures = movingMatchedPoints;
 
        % 应用转换——由于算法的随机性，运行之间的结果可能不相同
        tform = estimateGeometricTransform(movingMatchedPoints,fixedMatchedPoints,'projective');
        MOVINGREG.Transformation = tform;
        MOVINGREG.RegisteredImage = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj, 'SmoothEdges', true);
 
        Image_Align=MOVINGREG.RegisteredImage;
end
function checkLicense() %此子函数用来获取计算机视觉系统工具箱许可证
    CVSTStatus = license('test','video_and_image_blockset');
    if ~CVSTStatus
        error(message('images:imageRegistration:CVSTRequired'));
    end
end

Main.m：

clear;clc;
%读入图片
Image1=imread('hw1_painting_1.jpg');
Image2=imread('hw1_painting_2.jpg');
%变为灰度图
Image1_gary=rgb2gray(Image1);
Image2_gary=rgb2gray(Image2);
figure(1);
subplot(211);imshow(Image1_gary);title('原图像1的灰度图像');
subplot(212);imshow(Image2_gary);title('原图像2的灰度图像');
%将灰度图像对齐
Image2_Align=Align(Image2_gary,Image1_gary);
%% 差影法
Image_sub1=imabsdiff(Image2_gary,Image1_gary);%imabsdiff,两个图像的绝对差
Image_sub2=(Image2_Align-Image1_gary)+(Image1_gary-Image2_Align);%这里等同于imabsdiff
figure(2);
subplot(211);imshow(Image_sub1);title('未对齐的图像差');
subplot(212);imshow(Image_sub2);title('对齐后的图像差');
%筛选灰度值
[height,width,c]=size(Image_sub2);%height代表高度，width代表宽度，c代表通道数
A=0;B=0;
for i=1:height
    for j=1:width       
            A=A+1;%计算总灰度值
        if Image_sub2(i,j,1)<30
            B=B+1;%计算小于30的灰度值
        end
    end
end
C=B/A;
if C>0.92  %为阈值函数做准备
    Image_sub3=immultiply(Image_sub2,3);%图像的相乘，使整个图像的灰度都乘3
end
figure(3);
subplot(221);imshow(Image_sub3);title('灰度增强后的图像差');
%%  处理差影法获得的差异图
%阈值函数阈值分割
T=graythresh(Image_sub3);%为了找到合理的阈值,使用 Otsu 方法计算全局图像阈值
Image_bin=imbinarize(Image_sub3,T);%imbinarize,通过阈值化将二维灰度图像或三维体二值化
subplot(222);imshow(Image_bin);title('二值化后的图像差');
%形态学处理 -- bwmorph,针对二值图像的形态学运算
Image_bin=bwmorph(Image_bin,'clean',100);%100指的是进行100次处理，clean：去除图像中孤立的亮点，比如，一个像素点，像素值为1，其周围像素的像素值全为0，则这个孤立的亮点将被去除
Image_bin=bwmorph(Image_bin,'majority');%如果一个像素的8邻域中有等于或超过5个像素点的像素值为1， 则将该点像素值置1
Image_bin = bwareaopen(Image_bin,10,8);%bwareaopen,从二值图像中删除小对象。删除二值图像Image_bin中面积小于10的对象，默认情况下使用8邻域
subplot(223);imshow(Image_bin);title('形态学小对象处理后的二值化图像差');
% -- imdilate,膨胀图像
market_d=strel('disk',4);%创建一个4×4的圆形模板,strel,形态学结构元素
Image_bin_dilation=imdilate(Image_bin,market_d);%根据模板来膨胀
Image_bin=imfill(Image_bin_dilation,'holes');%填充，黑色背景上有白色的洞，洞便会被填满
Image_bin = bwareaopen(Image_bin,100,8);%删除二值图像Image_two中面积小于100的对象，默认情况下使用8邻域
subplot(224);imshow(Image_bin);title('形态学膨胀及大对象处理后的二值化图像差');
%边角置黑
[m,n]=size(Image_bin);%获取图像的长和宽
Image_bin(m-10:m,1:n)=0;   %下边角
Image_bin(1:m,1:7)=0;      %左边角
Image_bin(1:m,n-6:n)=0;    %右边角
Image_bin = bwareaopen(Image_bin,100,8);%处理膨胀后细小的点，删除二值图像Image_two中面积小于100的对象，默认情况下使用8邻域
figure(4);
imshow(Image_bin);title('边角置黑后的最终二值化图像差');
%% 绘图
[L,num]=bwlabel(Image_bin);%bwlabel,对二维二值图像中的连通分量进行标注。 L为和BW(Image_bin)大小相同的矩阵，包含了标记了BW中每个连通区域的类别标签，num为连通域的个数
figure(5);
mesh(L);title('连通区域分布图');
Image_T=regionprops(L,'area','boundingbox');%regionprops,测量图像连通区域的属性。area:标量，计算出在图像各个连通区域中像素总个数, boundingbox:1行 ndims(L)*2列的向量，即包含相应连通区域的最小矩形 (ndims,数组维度数目)
T_area=[Image_T.Area];%将像素个数矩阵上下拼接
rects=cat(1,Image_T.BoundingBox);%将位置矩阵上下拼接
figure(6),
subplot(1,2,1),imshow(Image1);title('原图像1');
for i=1:size(rects,1)
   rectangle('position',rects(i,:),'EdgeColor','r') %将框画在第一幅图上
end
subplot(1,2,2),imshow(Image2);title('原图像2');

4.处理结果

最后的处理结果如图所示，可以看出三处不同的地方在图中已经被框选出来了。

经过测试发现，我们的程序在处理其他图像时，也能取得很好的效果。

二、下图含有干扰条纹（moiré pattern）、并且低灰度区域的细节难以分辨。对图像进行灰度变换，再分别应用空域和频域的处理方法消除条纹干扰，并比较这两类方法去除条纹干扰的效果。

 1.算法原理：

（1）伽马变换：

        在图像处理中，常常利用伽马变换来对过曝或者曝光不足（过暗）的灰度图利用伽马变换进行对比度调节。具体年来讲：就是通过非线性变换，让图像中较暗的区域的灰度值得到增强，图像中灰度值过大的区域的灰度值得到降低。经过伽马变换，图像整体的细节表现会得到增强。

        数学公式为，其中，r为灰度图像的输入值（原来的灰度值），取值范围为[0,1]。s为经过伽马变换后的灰度输出值。c为灰度缩放系数，通常取1。γ为伽马因子大小。控制了整个变换的缩放程度。

（2）中值滤波： 

        图像处理时经常要做的一件事就是滤波，其中线性滤波器比如FIR、IIR 等类型都是研究的比较透彻的，实际使用中也有很好的效果。但是有时我们遇到的噪声比较顽固，比如说电子信号中的爆米花噪声（popcorn noise）还有图像处理中的椒盐噪声(salt-and-pepper noise)，用普通的线性滤波器只能将其压低，而无法彻底消除。这时一些非线性滤波器就体现出优势来了。

        中值滤波法是一种非线性平滑技术，它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值.中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。方法是用某种结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为，其中，分别为原始图像和处理后图像。W为二维模板，通常为3*3，5*5区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。

（3）陷波滤波器：

        条纹干扰是一种周期噪声，带阻滤波器可以用于去除周期性噪声，陷波滤波器也用于去除周期噪声，虽然带阻滤波器也能可以去除周期噪声，但是带阻滤波器对噪声以外的成分也有衰减。而陷波滤波器主要对某个点进行衰减，对其余的成分不损失。

        n阶巴特沃斯陷波滤波器的传递函数由下式给出：

其中，

        由于傅里叶的周期性，傅里叶频谱上不可能单独存在一个点的噪声，必定是关于远点对称的一个噪声对。这里的(u0,v0)和(-u0,-v0)是“陷波”的位置，即需要去除的噪声点；D0是它们的半径。

        注意，滤波器指定了频率矩形的中心，所以在使用前必须用函数fftshift进行预处理；可以同时输入多个陷波。

2.解题步骤 

（1）读取原图像，转为灰度图像，并进行归一化处理。

（2）用伽马变换进行灰度变换，增强低灰度区域的细节。

               

（3）利用中值滤波进行空域滤波处理，分别调用MATLAB的自带函数medfilt2和自己编写的函数midfilt。

（4）查看该图像的傅里叶频谱图像。

可以看出在中心上下两侧有明显的亮斑，可以推测就是条纹干扰的频谱所在。

（5）结合mesh函数绘制出的频谱分布图和在y方向的纵向截面判断噪声所在的频谱位置。

可以看出在419和553出有两个尖峰，应该就是噪声所在位置，而频谱中心位置为486，间距为67。

（6）构造n阶巴特沃斯陷波器实现频域滤波处理，通过修改半径D0和阶次n来获得更好的效果。

 3.程序代码

clear;clc;
Image=imread('hw2_rg.jpg'); %读取原始图像
figure(1);
subplot(211);imshow(Image,[0 1]); %灰度范围[0 1]
title('原始图像--带有条纹干扰');
Image_gray=rgb2gray(Image);%变为灰度图
Image_G=mat2gray(Image_gray);%归一化
%% 图像增强——灰度变换，伽马变换
r = Image_G;
c = 1;%C取值为1
gamma = 0.5;
s= c * r.^gamma; 
Image_gamma=s;
subplot(212);imshow(Image_gamma,[0 1]); %灰度范围[0 1]
title('图像增强--伽马变换后的图像');
%% 空域域处理方法 -- 中值滤波
k1=medfilt2(Image_gamma,[15,10]);%模板选择15*10，长方形的模板
k2=midfilt(Image_gamma,15);%模板选择15*10，长方形的模板
k2=histeq(k2);
figure(2);
subplot(211);imshow(k1);title('空域滤波去除条纹干扰后的图像');
subplot(212);imshow(k2);title('自己编写的中值滤波程序');
%% 傅里叶幅度谱
[M,N] = size(Image_gamma);
P=2*M;Q=2*N;
D=zeros(P,Q);
D(1:M,1:N)=Image_gamma;
u = 0:(P-1);
v = 0:(Q-1);
[V,U] = meshgrid(v,u);
j=fft2(D.*(-1).^(U+V));
figure(3);
% imshow(fftshift(abs(fft2(Image_gamma))),[0,100]);title('傅里叶幅度谱');
subplot(211);imshow(abs(j),[0,500]);title('傅里叶幅度谱');
% mesh(abs(j));title('傅里叶幅度谱的mesh分布--鼠标滑轮放大');
subplot(212);plot(abs(j(:,floor(P/2)-11:floor(Q/2)+9))); title('y方向频域纵向截面');% y方向纵向截线
%% 频域处理方法 -- 构造n阶巴特沃斯陷波器，可以修改D0和n获得更好的效果。
u = 0:(P-1);
v = 0:(Q-1);
[V,U] = meshgrid(v,u);
D0 = 50;  %陷波半径
n = 5;    %阶次  
u0 = 176; %485/2+1=243; 243-(553-486)=176
v0 = 445; %893/2+1=447; 447-2=445
u1 = 310; %485/2+1=243; 243+(553-486)=310
v1 = 449; %893/2+1=447; 447+2=449
D1 = sqrt((U-(floor(M/2)+1)-u0).^2+(V-(floor(N/2)+1)-v0).^2); %噪声点位置
D2 = sqrt((U-(floor(M/2)+1)-u1).^2+(V-(floor(N/2)+1)-v1).^2); %噪声点位置
H = 1./(1+(D0^2./(D1.*D2)).^n); %构造n阶巴特沃斯陷波器
figure(4);
subplot(211);imshow(log(1+abs(H)),[]);title('生成的n阶巴特沃斯陷波器');
% mesh(log(1+abs(H))); title('傅里叶幅度谱的mesh分布--鼠标滑轮放大');
% 频域滤波
G = j.*H;
subplot(212);imshow(abs(G),[0,500]);title('滤波后的傅里叶幅度谱');
%傅里叶逆变换恢复图像 
g = real(ifft2(G));
g=g(1:485,1:893);
figure(5);
imshow(g,[0 1]);title('频域滤波去除条纹干扰后的图像');
%% 
function d=midfilt(x,n)   
p=size(x);   %输入图像是p×q的,且p>n,q>n
x1=double(x);
x2=x1;
for i=1:p(1)-n+1
    for j=1:p(2)-n+1
        c=x1(i:i+(n-1),j:j+(n-1));  %取出x1中从(i,j)开始的n行n列元素,即模板(n×n的)
        e=c(1,:);      %是c矩阵的第一行
        for u=2:n
            e=[e,c(u,:)];     %将c矩阵变为一个行矩阵    
        end
        mm=median(e);      %mm是中值
        x2(i+(n-1)/2,j+(n-1)/2)=mm;   %将模板各元素的中值赋给模板中心位置的元素
    end
end 
%未被赋值的元素取原值
d=uint8(x2);
end

 4.处理结果：

              

        可以看出调用MATLAB系统函数的中值滤波效果很好，而自己编写的函数效果却不行，并且模板尺寸在15*15或者15*10才能滤掉条纹干扰，过大过小效果差距都很大，有时候很难准确找对模板大小。

        而频率滤波则能准确的滤掉噪声所在位置，并且对原图像细节影响很小。（使用了5阶滤波器）

三、下图是一张单据的扫描件，设计算法将表格横线校正为水平、并补全断裂表格线。

1.相关原理

（1）Radon变换

        Radon变换是一个积分变换，它将定义在二维平面上的一个函数 f(x,y) 沿着平面上的任意一条直线做线积分，相当于对函数 f(x,y) 做 CT扫描。Radon变换的本质是将原来的函数做了一个空间转换，即，将原来的XY平面内的点映射到AB平面上，那么原来在XY平面上的一条直线的所有的点在AB平面上都位于同一点。记录AB平面上的点的积累厚度，便可知XY平面上的线的存在性。公式如下:

Radon变换原理图

（2）Hough变换

        霍夫变换(Hough Transform)是图像处理中的一种特征提取技术，霍夫变换运用两个坐标空间之间的变换将在一个空间中具有相同形状的曲线或直线映射到另一个坐标空间的一个点上形成峰值，从而把检测任意形状的问题转化为统计峰值问题。

        一条直线在直角坐标系下可以用y=kx+b表示, 霍夫变换的主要思想是将该方程的参数和变量交换，即用x,y作为已知量k,b作为变量坐标，所以直角坐标系下的直线y=kx+b在参数空间表示为点(k,b)，而一个点(x1,y1)在直角坐标系下表示为一条直线y1=x1·k+b，其中(k,b)是该直线上的任意点。为了计算方便，我们将参数空间的坐标表示为极坐标下的γ和θ。因为同一条直线上的点对应的(γ,θ)是相同的，因此可以先将图片进行边缘检测，然后对图像上每一个非零像素点，在参数坐标下变换为一条直线，那么在直角坐标下属于同一条直线的点便在参数空间形成多条直线并内交于一点。因此可用该原理进行直线检测。

Hough变换原理图

 2.程序源码

（1）Radon变换简易实现：

clear;close all;
bw=imread('scanned_table.jpg');
% bw=rgb2gray(bw);
 
%% 图像校正
I=edge(bw);  %检测边缘
theta = 1:180;  
[R,xp] = radon(I,theta); %rando变换
[I,J] = find(R>=max(max(R))); %找到最亮（暗）点
angle=90-J; %确定倾斜角度
figure,imshow(bw);title('灰度图像');
bw=imrotate(bw,angle,'bilinear','crop');
 
%% 二值化
thresh=graythresh(bw);   %二值化确定阈值
pic_2b=im2bw(bw,thresh);    %图像二值化
 
%% 补全断线
[m,n]=size(pic_2b);
for i=1:m
    sum=0;
    for j=1:n %列扫描 补全行
        if (pic_2b(i,j)==0)
            sum=sum+1;
        end
    end
    if (sum/n>0.3)
        for j=200:1000  %表格长度范围
           pic_2b(i,j)=0; 
        end
    end
end
for j=1:n
    sum=0;
    for i=1:m %行扫描 补全列
        if (pic_2b(i,j)==0)
            sum=sum+1;
        end
    end
    if (sum/m>0.4)
        for i=120:550  %表格宽度范围
           pic_2b(i,j)=0; 
        end
    end
end
figure();
imshow(pic_2b);
 

（2）Hough变换均值聚类

close all;clear 
clc;clf
 
F = im2double(imread('scanned_table.jpg')); %输入归一化图像
[r,c] = size(F);
BW = edge(F,'Canny');
figure(1);
subplot(2,2,1);imshow(F);title('原图像');
subplot(2,2,2);imshow(BW);title('边缘检测');
[H,T,~] = hough(BW);   %hough变换取角度
P = houghpeaks(H,1);   %取图像横线rho、theta
theta = T(P(2));
theta = 90-abs(theta);%旋转弧度角
L = fix((pi*theta/180)*c);       %计算旋转黑边弧长
f = padarray(F,[L,L],'symmetric');    %扩大图像
N = imrotate(f,theta,'bilinear','crop'); %图像旋转
N = N(L:r+L-1,L:c+L-1);  %取原图像大小
subplot(2,2,3);imshow(N);title('校正后图像');   %校正后图像
 
%投影法
row = sum(N,1);col = sum(N,2)';     %图像矩阵 横向求和 纵向求和
[~,left] = find(row<600,1,'first');[~,right] = find(row<600,1,'last');  %找到表格边界
[~,up] = find(col<1200,1,'first');[~,down] = find(col<1200,1,'last');
sort_row = sort(row);row_data = zeros(1,100);
sort_col = sort(col);col_data = zeros(1,200);
k1 = 6;  %要画的横线数
for i = 1:200   %最大取200个 过大存在许多重复点无法找寻位置
    [~,n] = find( sort_row(i) == row );  %寻找最小值位置
    row_data(i) = n;    
end
[y1,Num_row,Mean_row] = C_mean(row_data,k1);  %C均值聚类 得到聚类数据
k2 = 16;
for i = 1:400  %最大取400个 过大存在许多重复点无法找寻位置
    [~,n] = find( sort_col(i) == col );  %寻找最小值位置
    col_data(i) = n;
end
[y2,Num_col,Mean_col] = C_mean(col_data,k2);  %C均值聚类 得到聚类数据
subplot(2,2,4);imshow(N);hold on;title('补线后图像');
for i = 1:k1  %补出竖线
    plot([y1{i}(1),y1{i}(1)],[up,down],'LineWidth',1.3,'Color','black');
end
for i = 1:k2  %补出横线
    plot([left,right],[y2{i}(1),y2{i}(1)],'LineWidth',1.3,'Color','black');
end
 function [y,Num,mean] = C_mean(data,k)  %data原数据  k为划线数  
data = data';
[m,n]  = size(data);    %获得数据维度
sum_Eirr = zeros(1,k);   %存放一到五类时的误差平方和
Mean = cell(1,k);   %聚类中心   最大M个
for i = 1:k
    Mean{1,i} = zeros(1,n);   %聚类中心初始化
end
p=[4,9,12,23,34,45,56,167,178,289,256,311,322,233,144,255];  %随机初始化K个聚类中心
for i =1:k
    Mean{1,i} = data(p(i));  
end
Eirr =0;
while 1
    Fdata = cell(1,k);     %初始化聚类
    count = zeros(1,k);     %聚类中元素的个数
    for i =1:k
        Fdata{1,i} = zeros(m,n);   %一个聚类最多有m个数据
    end
 
    for i = 1:m   %计算每一个数据跟聚类中心的距离
        for j = 1:k    
            distance(j) = sum ( ( data(i)-Mean{1,j} ).^2 );
        end
        [~,Lmin] = min( sqrt(distance) );    %寻找最小的距离 并且返回下标,Lmin表示哪一个类别
        count(Lmin) = count(Lmin) +1;
        Fdata{ 1,Lmin }(count(Lmin) , :) = data(i);   %记录下这个数据
    end
    old_Eirr = sum(Eirr);   %记录上一次误差平方和  为了和下次比较
    Eirr = zeros(1,k);   
    for i = 1:k
        for j = 1:count(i)
            Eirr(1,i) = Eirr(1,i) +(sum ( ( Fdata{1,i}(j,:)-Mean{1,i} ).^2 )); %误差平方和
        end
    end
    sum_Eirr(1,k) = sum(Eirr);     %记录每一次的误差平方和
    for i = 1:k    %计算新的聚类中心
        Mean{1,i} = sum(Fdata{1,i}) / count(i);
    end
    if old_Eirr == sum(Eirr) %当误差平方和保持不变达到最优迭代
        break;
    end
end
y = Fdata;Num = count;mean = Mean;

 3.结果分析

（1）Radon变换简易实现：

Radon变换简易实现

（2）Hough变换+均值聚类

Hough变换+均值聚类实现

        两种方案对比得知，第一种方法通过Radon变换将图像矫正后直接对二值化后的图像行列扫描进行补全，能实现要求。但是相比第二种Hough变换检测直线后均值聚类实现的方案效果差一点。存在图片存在黑边，结果是二值图像等缺陷。但是代码简单易懂。

四、有4种飞机模型，每种5架。请分别用傅里叶描述子和矩不变量提取特征，每种飞机中选3架进行训练，另外2架进行测试，用最近邻分类器进行识别。阐述所用算法的原理、关键步骤的处理结果、最终结果等，要附程序代码。

 1.所用算法原理

（1）傅里叶描绘子(Fourier Descriptors)

       傅里叶描述子的基本思想是：假定物体的形状是一条封闭的曲线，沿边界曲线上的一个动点s(k)的坐标变化x(k)+jy(k)是一个以形状边界周长为周期的函数，这个周期函数可以用傅里叶级数展开表示，傅里叶级数中的一系列系数a(u)是直接与边界曲线的形状有关的，称为傅里叶描述子。

       s(k)的离散傅里叶变换(DFT)为:

        复系数au称为边界的傅里叶描绘子，这些系数的傅里叶逆变换可以恢复s(k):

        假设我们仅使用前P个傅里叶系数，而不使用全部系数，结果是s(k)的如下近似:

        其中，k = 0,1,2,…,K-1。虽然仅使用P个傅里叶系数可得到s(k)的每一个分量，但k的范围仍是从0到K-1。换言之，近似边界中包含有相同数量的点，但在每个点的重构中却用不到如此多的系数项。高频分量决定细节部分，低频分量决定总体形状。因此，随着P的减少，边界细节的损失就会增加。

（2）矩不变量

        一幅数字图像f (x,y)的二维(p +q)阶矩定义为:

        其中，p,q=0,i,2,…，求和在跨越图像的所有空间坐标x,y的值上进行。相应的中心矩定义为:

        归一化(p+q)阶中心矩定义为:

        对平移、缩放、镜像和旋转都不敏感的7个二维不变矩的集合(ϕ1~ϕ7)可以由这些公式推导出来，因此用这7个矩就可以代表一个图像了。

（3）最近邻算法

         从方法论角度出发， KNN 认为，待分类对象的类别可以通过在它附近的训练数据的类别来确定，所以采取的策略就是找到离待分类对象最近的 K 个邻居进行分析。

        设N个样本中，来自Wc类的样本有Nc个，若K1,K1,…,Kc分别是K个近邻中属于W1,W1,…,Wc类的样本数，则我们可以定义判别函数为:

        决策规则为:     若

        则决策。这就是K近邻的基本规则。

在 KNN 的设计过程中，有四个要点需要注意:

a)用来对待分类对象所属类别进行评估的数据集合（不一定需要用到整个训练集）；

b)用来计算对象之间相似度的距离或者相似度矩阵（比如，欧式距离，马氏距离等）；

c)K 值的选取，若K=1时，即为最近邻算法；

d)用来确定待分类对象所属类别的方法（比如，距离加权与否）。

2.关键程序源码

（1）主程序代码：

clc;clear;
%每次运行都重新创建图像集的文件夹
rmdir('重构的图像集','s');
mkdir('重构的图像集');
%% 对图像数据库进行处理，获得训练集和测试集
path='.\试题3图像\';   %飞机模型数据库所在目录
for i=1:4 %4种飞机
    index=randperm(5); %随机排序数字1-5
    for j=1:3 %取前3架为训练集
        model1=char("飞机"+i+"-"+index(j)); %飞机型号
        %将型号储存到训练标签集
        TrainLable(1,(i-1)*3+j)=i;
        TrainLable(2,(i-1)*3+j)=index(j);
        imagePath1=[path model1 '.jpg']; %构成飞机图像的完整路径
        image1=imread(imagePath1); %读取飞机图像
        %第4种飞机模型是3通道的，要灰度化
        if i==4
            image1=rgb2gray(image1);
        end
        
        BW=imbinarize(image1); %将图像二值化
        Image_bin=morphology(BW);%调用morphology函数形态学处理
        %imshow(Image_bin);
        
        %以下为方法一：傅里叶描绘子
        %提取出飞机模型的边界(调用boundaries函数和bound2im函数)
        [len,wid]=size(Image_bin);
        b=boundaries(Image_bin);
        b=b{1};
        bim=bound2im(b,len,wid);
        %imshow(bim);
        
        %用傅里叶描绘子提取边界特征(调用frdescp函数和ifrdescp函数)
        z=frdescp(b);
        %z_restructure=ifrdescp(z,length(z)/2); %取一半的描述子重构（400-500）
        z_restructure=ifrdescp(z,56); %取56个描述子重构图像边界
        
        %将提取的特征储存到训练集1
        TrainSet1(:,(i-1)*3+j)=reshape(z_restructure(1:700,:),1400,1);%取1400个维度的特征
        
        %转化为可用imshow显示的类型
        img_restructure=bound2im(z_restructure,len,wid);
        %imshow(img_restructure);
        %输出到图像集文件夹 
        imwrite(img_restructure,strcat( '.\重构的图像集\',model1,'.jpg'));
        
        %以下为方法二：矩不变量
        fai=two_dim_moment(Image_bin);
        %将矩不变量储存到训练集2
        TrainSet2(:,(i-1)*3+j)=fai';
    end
    
    for j=4:5 %取后两架为测试集
        model2=char("飞机"+i+"-"+index(j)); %飞机型号
        %将型号储存到测试标签集
        TestLable(1,(i-1)*2+(j-3))=i;
        TestLable(2,(i-1)*2+(j-3))=index(j);
        imagePath2=[path model2 '.jpg']; %构成飞机图像的完整路径
        image2=imread(imagePath2); %读取飞机图像
        %第4种飞机模型是3通道的，要灰度化
        if i==4
            image2=rgb2gray(image2);
        end
        
        BW=imbinarize(image2); %将图像二值化
        Image_bin=morphology(BW);%调用morphology函数形态学处理
        %imshow(Image_bin);
        
        %以下为方法一：傅里叶描绘子
        %提取出飞机模型的边界(调用boundaries函数和bound2im函数)
        [len,wid]=size(Image_bin);
        b=boundaries(Image_bin);
        b=b{1};
        bim=bound2im(b,len,wid);
        %imshow(bim);
        
        %用傅里叶描绘子提取边界特征(调用frdescp函数和ifrdescp函数)
        z=frdescp(b);
        %z_restructure=ifrdescp(z,length(z)/2); %取一半的描述子重构（400-500）
        z_restructure=ifrdescp(z,56); %取56个描述子重构图像边界
        
        %将提取的特征储存到测试集1
        TestSet1(:,(i-1)*2+(j-3))=reshape(z_restructure(1:700,:),1400,1);%取1400个维度的特征
        
        %转化为可用imshow显示的类型
        img_restructure=bound2im(z_restructure,len,wid);
        %imshow(z_half2im);
        %输出到图像集文件夹        
        imwrite(img_restructure,strcat( '.\重构的图像集\',model2,'.jpg'));
        
        %以下为方法二：矩不变量
        fai=two_dim_moment(Image_bin);
        %将矩不变量储存到测试集2
        TestSet2(:,(i-1)*2+(j-3))=fai';
    end
end
%% 调用最近邻分类器进行识别分类
KNN(TrainSet1,TestSet1,TrainLable,TestLable,1400); %1400个维度的特征
KNN(TrainSet2,TestSet2,TrainLable,TestLable,7); %7个维度的不变矩

（2）最近邻分类器代码：

%% 最近邻分类器
function []=KNN(TrainSet,TestSet,TrainLable,TestLable,dimnum)
    if dimnum==1400
        fprintf('使用傅里叶描绘子提取特征：\n');
    elseif dimnum==7
        fprintf('使用矩不变量提取特征：\n');
    end
    k=3; %KNN，此为3近邻，取k=1时即为最近邻
    knn=zeros(8,1); %测试集一共8个飞机模型
    err_num=0;
    for x=1:8 %测试集一共8个飞机模型
        for y=1:12 %训练集一共12个飞机模型
            distance=0;
            for dim=1:dimnum
                distance=distance+(TestSet(dim,x)-TrainSet(dim,y))^2; %欧式距离算法
            end
            distance=sqrt(distance);
            knn(y,1)=distance;
        end
 
        [B,Index]=sort(knn); %升序排列，Index是原位置
        first=Index(1,1);
        
        %飞机模型数据库所在目录
        if dimnum==1400
            path='.\重构的图像集\';
        elseif dimnum==7
            path='.\试题3图像\';   %飞机模型数据库所在目录
        end     
        model_first=char("飞机"+TrainLable(1,first)+"-"+TrainLable(2,first)); %最近邻飞机型号
        imagePath_first=[path model_first '.jpg']; %构成最近邻飞机图像的完整路径
        image_first=imread(imagePath_first); %读取最近邻的飞机图像
 
        model_test=char("飞机"+TestLable(1,x)+"-"+TestLable(2,x)); %待测飞机型号
        imagePath_test=[path model_test '.jpg']; %构成待测飞机图像的完整路径
        image_test=imread(imagePath_test); %读取待测飞机图像
        
        %imshow(I,[]) 显示灰度图像 I，根据 I 中的像素值范围对显示进行转换。
        figure();
        title1=['待分类图像: ' model_test];
        title2=['最近邻图像: ' model_first];
        subplot(121);imshow(image_test,[]);title(title1);
        subplot(122);imshow(image_first,[]);title(title2);
 
        class_test=TestLable(1,x);
        class_first=TrainLable(1,first);
 
        fprintf('测试集第%d组数据原为第%d类飞机模型，最近邻识别分类后为第%d类飞机模型,',x,class_test,class_first);
        if class_test==class_first
            fprintf('分类正确。\n');
        else
            fprintf('分类错误！\n');
            err_num=err_num+1;
        end
        pause(1);
    end
    err_rate=(err_num*100)/8;
    fprintf('最近邻分类完毕！测试数据集一共%d组数据，错误率为%f。%%\n',x,err_rate);
    fprintf('\n');
end

（3）傅里叶描绘子代码： 

%% frdescp函数用于计算边界s的傅里叶描绘子
function z=frdescp(s)
    [np,nc]=size(s);
    if nc~=2
        error("S must be of size np-by-2.");
    end
    if np/2~=round(np/2)
        s(end+1,:)=s(end,:);
        np=np+1;
    end
    x=0:(np-1);
    m=((-1).^x)';
    s(:,1)=m.*s(:,1);
    s(:,2)=m.*s(:,2);
    s=s(:,1)+1i*s(:,2);
    z=fft(s);
end
%% 给定一组傅里叶描绘子s，ifrdescp函数可用给定数量nd的描绘子计算其逆变换z
function s=ifrdescp(z,nd)
    np=length(z);
    if nargin==1 || nd>np
        nd=np;
    end
    x=0:(np-1);
    m=((-1).^x)';
    d=round((np-nd)/2);
    z(1:d)=0;
    z(np-d+1:np)=0;
    zz=ifft(z);
    s(:,1)=real(zz);
    s(:,2)=imag(zz);
    s(:,1)=m.*s(:,1);
    s(:,2)=m.*s(:,2);
end
（4）矩不变量代码：
%% two_dim_moment函数用于求出矩不变量中的7个二维不变矩
function fai=two_dim_moment(image) %fai为希腊字母，moment意为矩
    [m,n]=size(image);
    image=double(image);
    
    %图像的各阶矩
    mm=zeros(4,4);   
    for y=1:m
        for x=1:n           
            for q=1:4
                for p=1:4
                    mm(q,p)=mm(q,p)+x^(p-1)*y^(q-1)*image(y,x);
                end
            end
        end
    end
    mean_x=mm(2,1)/mm(1,1);
    mean_y=mm(1,2)/mm(1,1);
    %三阶中心矩
    u00=mm(1,1);
    u11=mm(2,2)-mean_y*mm(2,1);
    u20=mm(3,1)-mean_x*mm(2,1);
    u02=mm(1,3)-mean_y*mm(1,2);
    u30=mm(4,1)-3*mean_x*mm(3,1)+2*mean_x^2*mm(2,1);
    u03=mm(1,4)-3*mean_y*mm(1,3)+2*mean_y^2*mm(1,2);
    u21=mm(3,2)-2*mean_x*mm(2,2)-mean_y*mm(3,1)+2*mean_x^2*mm(1,2);
    u12=mm(2,3)-2*mean_y*mm(2,2)-mean_x*mm(1,3)+2*mean_y^2*mm(2,1);
    %归一化中心矩
    n20=u20/u00^2;
    n02=u02/u00^2;
    n11=u11/u00^2;
    n30=u30/u00^2.5;
    n03=u03/u00^2.5;
    n12=u12/u00^2.5;
    n21=u21/u00^2.5;
    %7个不变矩
    fai(1) = n20 + n02; 
    fai(2) = (n20-n02)^2 + 4*n11^2;
    fai(3) = (n30-3*n12)^2 + (3*n21-n03)^2; 
    fai(4) = (n30+n12)^2 + (n21+n03)^2;
    fai(5) = (n30-3*n12)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n21-n03)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
    fai(6) = (n20-n02)*((n30+n12)^2-(n21+n03)^2)+4*n11*(n30+n12)*(n21+n03);
    fai(7) = (3*n21-n03)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n12-n30)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
end

 3.关键步骤的处理结果

（1）读取图像以后将图像二值化：

二值化后的飞机模型 

（2）对二值化图像进行形态学处理：

形态学处理后的飞机模型

（3）提取出飞机模型的边界：

调用boundaries函数提取的边界

（4）计算出边界的傅里叶描述子以后，分别采用一半、56、28、8个数量的傅里叶描述子进行图像重构：

不同数量描述子重构的图像结果

        通过多次运行观察，发现重构时取56个描述子，它仅保留了边界的主要特征；如果少于56将会丢失边界的主要特征(飞机的凸出部分)而显示无法接受的失真；减少到4个和2个的时候，产生的图像为一个椭圆或者圆。所以数量取56时最合适。

（5）傅里叶描述子提取特征构成的训练集：

取前1400维度构成的描述子特征训练集

（6）提取7个二维不变矩特征构成的训练集：

七个维度构成的矩不变量特征训练集

4.最终结果

（1）56个描述子重构后的图像集：

用56个傅里叶描述子重构得到的图像集

（2）运行以后的分类结果

傅里叶描述子的最近邻分类结果示例

矩不变量的最近邻分类结果示例

图识别分类完成后的最终错误率

        经过多次运行观察，傅里叶描述子的识别准确率很高，而二维不变矩偶尔会发生识别分类错误的情况。