基于matlab指纹识别算法的实现解析
由于指纹所具有的普遍性，唯一性和不变性，以及指纹识别技术具有很高的可行性和实用性，使之成为目前最流行、也最可靠的个人身份认证技术之一。
本文主要对指纹图像进行三方面处理：图像预处理、特征提取和特征匹配。图像预处理包括四个步骤：图像分割、滤波增强、二值化、细化，对指纹图像进行预处理后，去除了原图像的冗余部分，方便后续的识别处理；特征提取主要是提取指纹图像细化后的端点和分叉点；特征匹配是利用两个指纹的图像进行特征点比较，来确定两幅图像是否来自于同一手指。
本文给出了指纹图像预处理、特征提取、特征匹配的matlab程序及处理结果。该结果证明，用matlab实现的这些算法的处理结果比较理想，满足识别的可行性和应用性。
- 指纹识别算法流程及GUI实现效果示例
指纹识别技术主要包括三大部分：指纹图像采集、指纹预处理、特征提取与匹配。如图1-4所示。
示例如下：
- 指纹识别的几个预处理过程
1、图像的分割
图像分割是从一幅图像中按一定规则将一些物体或区域加以分离，划分出我们感兴趣的部分或区域。经过分割后的图像更容易进行进一步的分类、分析和识别处理。图像分割要在指纹二值化和滤波及细化之前进行，如此可以减少计算的冗余量，提高指纹检测速度。
采集到的指纹内容分为目标区域和背景区域。在指纹图像中，由脊线和谷线组成的较清晰的部分，称之为目标区域；没有用的部分我们称之为背景区域。指纹识别中的分割就是将有用的目标区域分割出来，去掉没用的背景区域，以避免背景区域的各种干扰。指纹图像可分为四类区域：背景区、不可恢复区、可恢复区、清晰区，如下图所示。
2、 图像归一化
对指纹图像进行分割处理，消除剩下的背景区域前，首先要进行图像归一化。
对采集好的指纹图像进行归一化处理，是对指纹灰度图的灰度均值和方差做一次调整，使得不论用什么设备采集的指纹图像都可以有预期的方差和均值，从而屏蔽不必要的噪声。指纹归一化不改变指纹质量，只是方便指纹的后续处理并保证程序运行时收敛加快。
3、 图像的二值化
二值化就是将图像上的像素点的灰度值设置为0或1，也就是将整个图像呈现出明显的黑白视觉效果。指纹图像中包括目标和背景还有众多噪声，要想从原始的指纹图像中提取出目标，一般用的方法是设定一个阈值T，用T将图像中像素数据分成两部分，若输入灰度图像的函数为：
（2-4）
通过求解阈值T，从而把图像f(x,y)分成目标和背景两个区域，其中大于T的像素群为目标区域，小于等于T的像素群为背景区域，阈值的选取原则是：（1）尽可能的多保存图像信息；（2）尽可能的减少噪声。
4、指纹图像的滤波
一个优秀的指纹识别系统不仅需要高的识别准确度，还需要高的识别速度，而影响识别速度的最主要因素就是指纹图像的滤波，而滤波的好坏直接因素是增强滤波的算法，当然跟所使用的软件和硬件也有很大的关系。
5、图像细化
分割和滤波后的指纹图像再进行二值化处理后，脊线仍然有一定的宽度，指纹识别的匹配是只利用图像的点或线的特征，这些点或者特征只与脊线的走向或者纹理有关系，有一定宽度的二值化图像显得有些多余，所以需要对二值化图像进行细化处理，指纹二值化图像经过细化处理即可得到一个单一像素宽度的脊线，经过上述的细化处理，在后续的指纹特征提取和特征匹配的算法中大大的减少了计算的冗余量和出错率，使得指纹识别的速度和准确度有了很大的提高。
四、 指纹特征提取和特征匹配
特征点提取
（1）提取指纹的端点和交叉点
端点和交叉点均是指纹图像的两个细节特征，同时在指纹识别的的过程中起着重要的作用，因为识别的首要前提就是找到图像的所有端点和交叉点。先通过一函数对八个邻域的坐标位置进行定义，然后定义另一函数来找出细化后指纹图像的所有端点及交叉点。
将八邻域中的每个点依次两两相减并取其绝对值，后将所有结果加起来，因为端点处是两个点，即和为2时细化图像有端点，和为6时图像特征为交叉点。
运行完上面的和函数的程序后，能把细化图像的的端点和交叉点全部找出。在定义函数的程序中有数组txy，其中t为横坐标，x为纵坐标，y为2时为端点，y为6时为交叉点。
（2）去除图像边缘的端点
可以看出，指纹图像细化的边缘，由于采集仪器不同的关系，因此不可避免的会多出很多的端点，这些端点不仅增加了后续的工作量，还可能导致识别过程中产生错误，所以要把这些边缘的端点都去除，在matlab中这些操作都可以采用一函数来实现。
找出特征点
设置三个函数来找出图像的特征点：
（1）A函数
经过去除边缘端点的操作后进一步减少了指纹细化图像中的端点和交叉点的个数。下面就需要找出一些在细化图像中比较独特的端点来作为识别的特征点。在一幅细化的指纹图像中，如果在一个像素（该像素为端点）的周围半径为r（r为像素的个数）的圆内没有任何的端点或者交叉点，那么随着r的逐渐增大，这样的点就会越来越少，因此该点也就越来越独特。于是我们设计了A函数来找出这样独特的点。
（2）B函数
为了进一步找出特征点，我们还需定义B函数，它的主要作用就是判断某一端点在num的距离内是否还有其他的端点。
（3）C函数
A函数和B函数都是找细化图像特征点的函数，因此可以设计另一个新的last1函数，通过执行
[pxy3,error2]=C(thin,r,txy,num)
可以找出一端点以r为半径的像素内的任何端点和交叉点且沿着脊线走向的num内没有任何的其他端点和交叉点。
特征点匹配
由上文的函数可知，已经找出了指纹细化图像中的特征点，并画出了一段独特的脊线，在图像中用红色来标示。下面就是指纹匹配[12]的问题了。在此我们设置了三层匹配。
（1）脊线长度匹配
对于上面的函数即可找出细化图像中的特征点和一段脊线，沿着该段脊线走向，每隔五个像素测量一下，看到到原始端点的距离，此段距离由一juli函数得到。
函数结果会得到一数组（内有脊线的长度信息）。如果两幅指纹细化图像中的纹路是相同的，则它们就包含相同的端点和交叉点及用juli函数找出的相同的一段脊，则这两个指纹图像中的长度数组对应的位置比例会基本相等（我们选择的指纹图像大小基本相等，因此该比例选1），因此函数最终定义了一个数f=(sum(abs((d1./d2)-1)))，其中若f的值越接近于0，这两幅图像的匹配度就越高，在一定范围的阈值内我们可以认定为匹配。
（2）三角形边长匹配
找到一个指纹细化图像的特征点后，可以找出距离这个端点距离最近的两个端点或者交叉点，与这个指纹图像细化的特征点构成一个三角形，若两幅图像中的边长比例基本相等（原理同上，也选1），则说明这两幅图像匹配，越接近于1说明这两幅指纹图像越匹配。其中设置一find_point函数来找出距离最近的端点或交叉点。
函数最后定义了一个数ff=(sum(abs((dd1./dd2)-1)))，因此ff值越接近于0，这两幅指纹图像的匹配度越高，在一定范围的阈值内我们可以认定为匹配。
（3）点类型匹配
找到一个指纹细化图像的特征点后，在该端点周围找到四十个端点或者交叉点，统计在这四十个特征点中端点的个数和交叉点的个数。若有两幅指纹细化图像中的端点所占的比例近似相同，则两幅图像相匹配，越近似，则越相同。函数最终定义了一个数fff=abs(f11-f21)/(f11+f12)，所以fff值越接近于0，这两幅指纹图像的匹配度就会越高。我们也设定一阈值，在此阈值内都可以认定为匹配。
本文中取r=8，num=60，经过试验，得到f的阈值为，ff的阈值为，fff的阈值为。即两幅图像的f，ff，fff若均小于阈值，则两幅图匹配；若三个值中有至少一个值大于阈值，则不匹配。
五、部分参考源码
function img = tuxiangyuchuli(path)
M=0;var=0;
I=double(imread(path));
[m,n,p]=size(I);
for x=1:m
for y=1:n
M=M+I(x,y);
end
end
M1=M/(m*n);
for x=1:m
for y=1:n
var=var+(I(x,y)-M1).^2;
end
end
var1=var/(m*n);
for x=1:m
for y=1:n
if I(x,y)>=M1
I(x,y)=150+sqrt(2000*(I(x,y)-M1)/var1);
else
I(x,y)=150-sqrt(2000*(M1-I(x,y))/var1);
end
end
end
figure, imshow(I(:,:,3)./max(max(I(:,:,3))));title(‘归一化’)
%************************************************************************
M =3; %3*3
H = m/M; L= n/M;
aveg1=zeros(H,L);
var1=zeros(H,L); %计算每一块的平均值
for x=1:H;
for y=1:L;
aveg=0;var=0;
for i=1:M;
for j=1:M;
aveg=I(i+(x-1)*M,j+(y-1)*M)+aveg;
end
end
aveg1(x,y)=aveg/(M*M); %计算每一块的方差
for i=1:M;
for j=1:M;
var=(I(i+(x-1)*M,j+(y-1)*M)-aveg1(x,y)).^2+var;
end
end
var1(x,y)=var/(M*M);
end
end
Gmean=0;Vmean=0;
for x=1:H
for y=1:L
Gmean=Gmean+aveg1(x,y);
Vmean=Vmean+var1(x,y);
end
end
Gmean1=Gmean/(H*L); %所有块的平均值
Vmean1=Vmean/(H*L); %所有块的方差
gtemp=0;gtotle=0;vtotle=0;vtemp=0;
for x=1:H
for y=1:L
if Gmean1>aveg1(x,y)
gtemp=gtemp+1;
gtotle=gtotle+aveg1(x,y);
end
if Vmean1<var1(x,y)
vtemp=vtemp+1;
vtotle=vtotle+var1(x,y);
end
end
end
for x=2:H-1
for y=2:L-1
if e(x,y)==1
if e(x-1,y) + e(x-1,y+1) +e(x,y+1) + e(x+1,y+1) + e(x+1,y) + e(x+1,y-1) + e(x,y-1) + e(x-1,y-1) <=4
e(x,y)=0;
end
end
end
end
Icc = ones(m,n);
for x=1:H
for y=1:L
if e(x,y)==1
for i=1:M
for j=1:M
I(i+(x-1)*M,j+(y-1)*M)=G1;
Icc(i+(x-1)*M,j+(y-1)*M)=0;
end
end
end
end
end
figure, imshow(I(:,:,3)./max(max(I(:,:,3))));title('分割');
%************************************************************************
temp=(1/9)*[1 1 1;1 1 1;1 1 1]; % 模板系数、均值滤波
Im=double(I);
In=zeros(m,n);
for a=2:m-1;
for b=2:n-1;
In(a,b)=Im(a-1,b-1)*temp(1,1)+Im(a-1,b)*temp(1,2)+Im(a-1,b+1)*temp(1,3)+Im(a,b-1)*temp(2,1)+Im(a,b)*temp(2,2)+Im(a,b+1)*temp(2,3)+Im(a+1,b-1)*temp(3,1)+Im(a+1,b)*temp(3,2)+Im(a+1,b+1)*temp(3,3);
end
end
I=In;
Im=zeros(m,n);
for x=5:m-5;
for y=5:n-5;
sum1=I(x,y-4)+I(x,y-2)+I(x,y+2)+I(x,y+4);
sum2=I(x-2,y+4)+I(x-1,y+2)+I(x+1,y-2)+I(x+2,y-4);
sum3=I(x-2,y+2)+I(x-4,y+4)+I(x+2,y-2)+I(x+4,y-4);
sum4=I(x-2,y+1)+I(x-4,y+2)+I(x+2,y-1)+I(x+4,y-2);
sum5=I(x-2,y)+I(x-4,y)+I(x+2,y)+I(x+4,y);
sum6=I(x-4,y-2)+I(x-2,y-1)+I(x+2,y+1)+I(x+4,y+2);
sum7=I(x-4,y-4)+I(x-2,y-2)+I(x+2,y+2)+I(x+4,y+4);
sum8=I(x-2,y-4)+I(x-1,y-2)+I(x+1,y+2)+I(x+2,y+4);
sumi=[sum1,sum2,sum3,sum4,sum5,sum6,sum7,sum8];
summax=max(sumi);
summin=min(sumi);
summ=sum(sumi);
b=summ/8;
if (summax+summin+ 4*I(x,y))> (3*summ/8)
sumf = summin;
else
sumf =summax;
end
if sumf > b
Im(x,y)=128;
else
Im(x,y)=255;
end
end
end
for i=1:m
for j =1:n
Icc(i,j)=Icc(i,j)*Im(i,j);
end
end
for i=1:m
for j =1:n
if (Icc(i,j)==128)
Icc(i,j)=0;
else
Icc(i,j)=1;
end;
end
end
figure,imshow(double(Icc));title('二值化');
%************************************************************************
u=Icc;
[m,n]=size(u) %去空洞和毛刺
for x=2:m-1
for y=2:n-1
if u(x,y)==0
if u(x,y-1)+u(x-1,y)+u(x,y+1)+u(x+1,y)>=3
u(x,y)=1;
end
else u(x,y)=u(x,y);
end
end
end
figure,imshow(u) %title('去毛刺')
for a=2:m-1
for b=2:n-1
if u(a,b)==1
if abs(u(a,b+1)-u(a-1,b+1))+abs(u(a-1,b+1)-u(a-1,b))+abs(u(a-1,b)-u(a-1,b-1))+abs(u(a-1,b-1)-u(a,b-1))+abs(u(a,b-1)-u(a+1,b-1))+abs(u(a+1,b-1)-u(a+1,b))+abs(u(a+1,b)-u(a+1,b+1))+abs(u(a+1,b+1)-u(a,b+1))~=1 %去空洞
if(u(a,b+1)+u(a-1,b+1)+u(a-1,b))*(u(a,b-1)+u(a+1,b-1)+u(a+1,b))+(u(a-1,b)+u(a-1,b-1)+u(a,b-1))*(u(a+1,b)+u(a+1,b+1)+u(a,b+1))==0 %去毛刺
u(a,b)=0;
end
end
end
end
end
figure,imshow(u) %title('去空洞')
%*************************************************************************
v=~u;
se=strel('square',3);
fo=imopen(v,se);
v=imclose(fo,se); %对图像开操作和闭操作
img=bwmorph(v,'thin',Inf); %对图像进行细化
figure,imshow(img)
title('细化图')
function j = P (img, x, y, i)
switch (i)
case {1, 9}
j = img(x+1, y);
case 2
j = img(x + 1, y-1);
case 3
j = img(x, y - 1);
case 4
j = img(x - 1, y - 1);
case 5
j = img(x - 1, y);
case 6
j = img(x - 1, y + 1);
case 7
j = img(x, y + 1);
case 8
j = img(x + 1, y + 1);
end
%point函数
function txy=point(thin)
count = 1;
txy(count, :) = [0,0,0];
siz=min(size(thin,1),size(thin,2));
for x=40:siz - 40
for y=40:siz - 40
if (thin(y, x) )
CN = 0;
for i = 1:8
CN = CN + abs (P(thin, y, x, i) - P(thin, y, x, i + 1));
end
if (CN == 2)
txy(count, :) = [x, y,2];
count = count + 1;
end
if (CN == 6)
txy(count, :) = [x, y,6];
count = count + 1;
end
end
end
end
for i=1:count - 1
x(i) =txy(i, 1);
y(i)= txy(i, 2);
end
imshow(double(thin));
hold on;
plot(x,y,'.');
%cut函数
function txy=cut(thin,txy)
s(8,8)=0;
delta(8,8)=0;
n=size(txy,1);
txy=txy(find(txy(:,1)),:);
plot(txy(:,1),txy(:,2),'ro');
附录C 图像特征点代码
%single_point函数
function [pxy2,error]=single_point(txy,r)
error=0;
x=txy(:,1);
y=txy(:,2);
n=length(x);
d(1:n,1:n)=0;
for j=1:n
for i=1:n
if (i~=j)
d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
else
d(i,j)=2*r;
end
end
end
[a,b]=min(d);
c=find(a>r);
pxy2=txy(c,:);
pxy2=pxy2(find(pxy2(:,3)==2),:);
t=size(pxy2,1);
if t==0
error=1
else
plot(x,y,'b.');
hold on
plot(pxy2(:,1),pxy2(:,2),'r.');
end
%walk函数
function [error,a,b]=walk(thin,x0,y0,num)
error=0;
thin(y0,x0)=0;
t1=0;
for n=1:num
if error==1
break;
else
x=x0;
y=y0;
for x=x0-1:x0+1
if error==1
break;
else
for y=y0-1:y0+1
t1=sum(sum(thin(y0-1:y0+1,x0-1:x0+1)));
if (t1==0||t1>=2)
error=1;
a=x0;
b=y0;
break;
else
if (thin(y,x)==1&&(x-x0)^2+(y-y0)^2~=0)
if (t1>=2 )
error=1;
break ;
else
thin(y,x)=0;
x0=x;
y0=y;
a=x0;
b=y0;
plot(x0,y0,'r.')
end
end
end
end
end
end
end
end
%last1函数
function [pxy3,error2]=last1(thin,r,txy,num)
error=0;
[pxy2,error]=single_point(txy,r);
n=size(pxy2,1);
l=1;
error2=0;
附录D 特征点匹配代码
%distance函数程序
function d=distance(x0,y0,num,thin)
num2=fix(num/5);
for i=1:num2
[error,a,b]=walk(thin,x0,y0,5*i);
if error~=1
d(i)=sqrt((a-x0)^2+(b-y0)^2);
else
break;
end
end
%A函数
function pxy=A(x0,y0,txy,num)
x=txy(:,1);
y=txy(:,2);
n=length(x);
l(1,n)=0;
lnn=1;
pxy(num,:)=[0,0,0];
for i=1:n
l(i)=sqrt((x(i)-x0)^2+(y(i)-y0)^2);
end
ll=sort(l);
for i=1:num
xiao=ll(i+lnn);
nn=find(l==xiao);
lnn=length(nn);
pxy(i,:)=[x(nn(1)),y(nn(1)),txy(nn(1),3)];
end
plot(x0,y0,'bo');
x0;
y0;
hold on
plot(pxy(:,1),pxy(:,2),'ro');
%主函数main
close all;
tic
clear;
thin1=tuxiangyuchuli('');
thin2=tuxiangyuchuli('');
figure;
txy1=point(thin1);
txy2=point(thin2);
[w1,txy1]=guanghua(thin1,txy1);
[w2,txy2]=guanghua(thin2,txy2);
thin1=w1;
thin2=w2;
txy1=cut(thin1,txy1);
txy2=cut(thin2,txy2);
[pxy31,error2]=last1(thin1,8,txy1,60)
[pxy32,error2]=last1(thin2,8,txy2,60)
error=1;
num=20;
cxy1=pxy31;
cxy2=pxy32;
d1=distance(cxy1(1,1),cxy1(1,2),num,thin1);
d2=distance(cxy2(1,1),cxy2(1,2),num,thin2);
f=(sum(abs((d1./d2)-1)));
if ff<=
error=0;
else
error=1;
end
c11=find_point(cxy1(1,1),cxy1(1,2),txy1,1);
c12=find_point(cxy1(1,1),cxy1(1,2),txy1,2);
c21=find_point(cxy2(1,1),cxy2(1,2),txy2,1);
c22=find_point(cxy2(1,1),cxy2(1,2),txy2,2);
cxy1(2,:)=c11;
cxy1(3,:)=c12(2,:);
cxy2(2,:)=c21;
cxy2(3,:)=c22(2,:);
x11=cxy1(1,1); y11=cxy1(1,2);
x12=cxy1(2,1); y12=cxy1(2,2);
x13=cxy1(3,1); y13=cxy1(3,2);
x21=cxy2(1,1); y21=cxy2(1,2);
x22=cxy2(2,1); y22=cxy2(2,2);
x23=cxy2(3,1); y23=cxy2(3,2);
dd1(1)=juli(x11,y11,x12,y12);
dd1(2)=juli(x12,y12,x13,y13);
dd1(3)=juli(x13,y13,x11,y11);
dd2(1)=juli(x21,y21,x22,y22);
dd2(2)=juli(x22,y22,x23,y23);
dd2(3)=juli(x23,y23,x21,y21);
ff=(sum(abs((dd1./dd2)-1)))
if ff<=1
error=0;
else
error=1;
end
cxy1(2:41,:)=find_point(pxy31(1,1),pxy31(1,2),txy1,40);
cxy2(2:41,:)=find_point(pxy32(1,1),pxy32(1,2),txy2,40);
f11=length(find(cxy1(:,3)==2));
f12=length(find(cxy1(:,3)==6));
f21=length(find(cxy2(:,3)==2));
f22=length(find(cxy2(:,3)==6));
fff=abs(f11-f21)/(f11+f12)
toc