课题介绍
随着生物识别技术的不断发展，人们发现每个人的指纹具有唯一性和不变性。因此指纹识别技术逐步发展为一种新的身份识别方式，并且凭借其良好的安全可靠性，大有取代传统身份识别方式的趋势。
本文简要介绍了指纹识别的基本步骤，分别是指纹图像预处理、指纹特征提取、指纹匹配。在图像预处理中，依次介绍了规格化处理、图像增强、二值化处理和细化处理的方法。预处理后将得到一幅宽度为一个像素的细化二值图像，然后通过特定的端点和交叉点的特征进行指纹匹配。实验表明，该方法效果良好。
二、算法流程
传统的身份识别技术包括基于密码的身份识别技术和基于令牌的身份识别技术。但是这两种方式都存在着一定的安全风险：令牌容易丢失，密码容易忘记，且它们都容易被窃取或冒充。近年来，依赖人体特征的新型身份识别技术正在兴起，凭借其唯一且不变的特点，生物识别技术逐步成为人们研究的热点。
而在生物识别技术中，指纹识别技术已经被广泛的研究和应用。它的基本原理框图如图一所示。
指纹图像具体实现及源码
3.1图像规格化
受采集设备参数和环境的影响，采集到的指纹图像可能总体对比度较差。图像归一化的作用就是使所有指纹图像具有相同的灰度均值和方差，从而将每一幅图像的灰度调整到统一的范围，方便后续处理。归一化的算法是：
（1） 先计算图像的平均值和方差。这部分工作主要通过统计图像中各点像素值得到该图像的直方图，然后利用直方图来计算指纹图像的相关指标。
（2） 指定期望经过处理后的图像均值和方差，计算归一化后的图像G。
% 部分参考源码
% 归一化
M=0;var=0;
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(uint8(I)) ;
带入指纹图像’zhiwen1.bmp’后得到的原始图像和归一化图像见图2和图3：
3.2图像分割
指纹图像通常包括纹线区域和周边的无用区域。如果保留这些无用区域，会致使计算冗余度增加。为此要先将其从待处理区域中去除。通常用的分离方法有像素领域特征的方法和基于像素分布概率的方法等。本文采取对各像素邻域特征进行分析，然后分割图像算法。
吧图像分成T*T的非重叠块，计算每一块的灰度均值M和方差Var，通过下列条件区分前景与背景：
（1） 如果M>M1且Var<Var1，则认为背景
（2） 如果M<=M1且Var>=Var1，则认为前景，否则待定；
（3） 对所有分块划分完毕后，在对待定块进行判决。如果在8邻域中，背景景小于等于4则认为是前景，否则是背景
在编程过程中通过几次调试后对图像进行3*3分块处理
Matlab程序如下：
% 分割
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
G1=gtotle/gtemp;V1=vtotle/vtemp;
gtemp1=0;gtotle1=0;vtotle1=0;vtemp1=0;
for x=1:H
for y=1:L
if G1<aveg1(x,y)
gtemp1=gtemp1-1;
gtotle1=gtotle1+aveg1(x,y);
end
if 0<var1(x,y)<V1
vtemp1=vtemp1+1;
vtotle1=vtotle1+var1(x,y);
end
end
end
G2=gtotle1/gtemp1;V2=vtotle1/vtemp1;
e=zeros(H,L);
for x=1:H
for y=1:L
if aveg1(x,y)>G2 && var1(x,y)<V2
e(x,y)=1;
end
if aveg1(x,y)< G1-100 && var1(x,y)< V2
e(x,y)=1;
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(uint8(I));title('分割');
得到的图像如图4：
图4经过分割处理后的图像
3.3图像二值化
一般的指纹图像都有比较清晰的方向场，方向场估计得准确性直接决定了图像增强算法的效果。
为估计方向场，我们把指纹脊线的走向分为如下8个方向，如下图：
图5 在一个像素处的8个指纹脊线方向
我们先对分割后的图像进行了平均滤波，然后对图像的每一个像素，为确定在该像素出的脊线方向，在以该像素为中心的9*9窗口内，分别计算8个方向上的经过处理后的灰度值，即将图5中标了i（i=0,1,…7分别代表8个方向）的位置的像素灰度值去他们中最大summax和最小值summin，若满足(summax+summin+ 4*I(x,y))> (3*summ/8),则该像素点的脊线方向为summin,否则为summax.
确定完方向后就根据该向场对图像进行二值化。
Matlab程序如下：
%二值化¯
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('二值化');
得到的二值化图像如图6：
3.4图像增强
在当前的指纹采集条件下，不可避免的会受到环境，皮肤上的油脂、水分、污渍的影响，使采集到的指纹图像出现纹线粘连、纹线断裂等缺陷，对后续的指纹特征提取带来很大困难。所以我们必须对指纹图像进行图像增强处理，例如分离粘连的纹线，连接断裂的纹线，平滑纹线的边缘等，以保证指纹特征提取的可靠性。
实验中进行了初步去除空洞和毛刺的处理，程序如下：
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('去除空洞')
得到的去除毛刺后的图像和再去除空洞后的图像分别如图7和图8：
3.5图像细化
二值化后的纹线仍然具有一定的宽度，因为指纹识别只与纹线的走向有关，所以需要对二值图像进行细化，以减少冗余信息，突出纹线的有效特征，便于后续的特征提取。
程序如下：
v=~u;
se=strel('square',3);
fo=imopen(v,se);
v=imclose(fo,se); %对图像进行开操作和闭操作
w=bwmorph(v,'thin',Inf);%对图像进行细化
figure,imshow(w)
title('细化图')
得到的图像如图9：
4．特征点提取
4.1 找出所有的端点和交叉点
因为特征点必然是从端点和交叉点里找出，所以在得到细化的图像后，我们首先要找出所有的端点和交叉点。
程序如下：
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
然后找出细化后图像的所有端点。
将一个点的8个邻域依次两两相减并取绝对值，并将所有结果相加，从细化图像的特征来说，和为2时为端点，和为6时为交叉点。
程序如下：
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,'.');
运行完程序后，将所有的端点和交叉点全部找出。定义的数组txy第一项为横坐标，第二项为纵坐标，第三项为2或6（2为端点、6为交叉点）。
得到的标出端点和交叉点的图像如图10：
4.2 纹线光滑处理
在指纹图像预处理中，已经对指纹图像进行过去除毛刺和空洞的处理。进一步对细化后的图像进行光滑处理。基本原理为：找到每个端点，使其沿着纹线的方向移动5个像素，如果在5个像素之内遇到交叉点，则认为此端点为毛刺，去除此点。
for j=1:5
txy=point(thin);
pxy=txy(find(txy(:,3)==2),:);
n=size(pxy,1);
for i=1:n
error=0;
error=walk(thin,pxy(i,1),pxy(i,2),5);
if error==1
thin(pxy(i,2),pxy(i,1))=0;
end
end
end
w=thin;
imshow(w);
光滑后的图像如图11所示：
光滑完后需再画出新的端点，如下图所示：
图12 光滑后标出端点和特征点的图像
4.3 去除图像边缘的端点
可以看出，在指纹图像的边缘，由于采集仪器的关系，不可避免的多出很多端点，一方面增加了后续工作量，另一方面还可能产生错误，所以有必要将这些边缘的端点去除。
函数如下：
s(8,8)=0;
delta(8,8)=0;
n=size(txy,1);
for i=1:8
for j=1:8
mp{i,j}=thin(1+31*(i-1):31+31*(i-1),1+31*(j-1):31+31*(j-1));
s(i,j)=sum(sum(mp{i,j}))/(31*31);
mp{i,j}=(mp{i,j}-s(i,j)).^2;
delta(i,j)=sum(sum(mp{i,j}));
if delta(i,j)<=70
for k=1:n
if (txy(k,1)>=1+31*(i-1)&&txy(k,1)<=31+31*(i-1)&&txy(k,2)>=1+31*(j-1)&&txy(k,2)<=31+31*(j-1)&&txy(k,3)==2)
txy(k,:)=[0,0,0];
end
end
end
end
end
txy=txy(find(txy(:,1)),:);
plot(txy(:,1),txy(:,2),'ro');
得到的图像如下：
图13 去除图像边缘端点后的图像
4.4 找出特征点
4.4.1 经过光滑处理和去除边缘端点后进一步减少了端点和交叉点的个数。下面就要找出一些独特的端点来作为特征点。在一幅细化的指纹图上，如果一个端点的周围半径为r个像素的圆内没有任何端点或交叉点，那么随着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
比如在执行后，会显示出下图：
由运行结果可知，有3个端点满足周围半径为20个像素的圆内没有任何端点和交叉点的条件。
这三个点在下图中用红点标出。
图15 满足条件的特征点
4.4.2 用于进一步找出特征点，它的作用是判断离某一端点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
4.4.3 综合以上两个找特征点的函数，可以设计一个新的last1函数。通过找出周围半径为r个像素的圆内没有任何端点或交叉点，并且沿纹线走num个像素内没有另一个端点或交叉点的端点。
程序如下所示：
error=0;
[pxy2,error]=single_point(txy,r);
n=size(pxy2,1);
l=1;
error2=0;
for i=1:n
[error,a,b]=walk(thin,pxy2(i,1),pxy2(i,2),num);
if error~=1
pxy3(l,1)= pxy2(i,1);
pxy3(l,2)= pxy2(i,2);
pxy3(l,3)= pxy2(i,3);
l=l+1;
error2=0;
plot(pxy2(i,1) ,pxy2(i,2),'r+');
end
end
得到的Matlab运行结果显示为：
图16 last1函数的运行结果
得到的图像如下：
图17 指纹图像的显示
其中红色的纹线为沿上述端点走一定距离后的纹线，一部分可以一直走到规定距离（这些即是要找的特征点），还有一部分走到端点或叉点停止。
至此，每幅指纹图像的特征点均可找出。
5．特征点匹配
5.1 纹线长度匹配
上面，我们已经可以找出每幅指纹图像的特征点，并画出一段纹线。下面就是匹配的问题了。在此，我们设置了三层匹配。
首先是纹线长度匹配。对于上面找出的特征点和纹线，每沿着纹线走5个像素测量一下到原始端点的距离。
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
最后会得到一个装有长度信息的数组。试想如果两幅指纹图中的指纹是一样的，则它们会包含相同的特征点和从这个特征点出发画出的纹线，则这两个长度数组对应位置的比例应基本相等（考虑到老师所给的24幅指纹图像大小相同，这个比例应近似为1）；最终的函数中定义了一个数f=(sum(abs((d1./d2)-1)))，所以f值越接近于0，这两幅图像的匹配度越高。
5.2 三角形边长匹配
找到一个特征点后，可以找出距离其最近的2个端点或交叉点，与原特征点构成三角形，若两幅图像的三角形边长比例相等（本实验中都为1），则说明这两幅图像匹配。
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');
在最终程序中定义了一个数ff=(sum(abs((dd1./dd2)-1)))，所以ff值越接近于0，这两幅图像的匹配度越高。
5.3 点类型匹配
找到一个特征点后，在其周围找到40个端点或交叉点，统计这40个点中端点和交叉点的个数。若两幅图中端点占的比例近似相同，则两幅图像相匹配。在最终函数中定义了一个数fff=abs(f11-f21)/(f11+f12)，所以fff值越接近于0，这两幅图像的匹配度越高。
最终程序如下（里面调用到前面一些函数）：
close all;
tic
clear;
thin1=tuxiangyuchuli('zhiwen8.bmp');
thin2=tuxiangyuchuli('zhiwen9.bmp');
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 f<=0.5
error=0;
else
error=1;
end
f
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);
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
6．实验结果
本实验中，取r=8，num=60，经过反复试验，可以得到f的阈值为0.5，ff的阈值为1.5，fff的阈值为0.2。即三幅图像的f，ff，fff若均小于阈值，则两幅图匹配，显示出MATCH。若三个值中有至少一个值大于阈值，则不匹配，显示出NOT MATCH。
实验中我们进行了数次匹配（每次匹配大概需要70秒），下面列出两组结果：
zhiwen8和zhiwen7：f=0.3525；ff=0.6803；fff=0.1463。所以匹配。
Zhiwen8和zhiwen1：f=0.1389；ff=1.5441；fff=0.2927。所以不匹配。
