matlab实现人眼定位程序

%-------------------------------------------------------------------------
clear all
x=imread('C:\wuzun.jpg'); %读取原始图像
figure(1);subplot(1,4,1);imshow(x);title('原图像1');
y=rgb2gray(x);
subplot(1,4,2);imshow(y);title('图像1的灰度图'); % 图1灰度图
u1=imnoise(y,'salt & pepper',0.07);
subplot(1,4,3);imshow(u1);title('图像1加噪声图'); %给图1加椒盐噪声
zz=medfilt2(u1,[3 3]);z=medfilt2(zz,[5 5]);%(2次中值滤波)
subplot(1,4,4);imshow(z);title('图像1中值滤波'); %图1中值滤波（3*3窗口）
% figure(6);subplot(1,2,1);imshow(zz);title('一次滤波');
% subplot(1,2,2);imshow(z);title('二次滤波');
x1=imread('C:\wuzun1.jpg');
figure(2);subplot(1,4,1);imshow(x1);title('原图像2');
y1=rgb2gray(x1);
subplot(1,4,2);imshow(y1);title('图像2灰度图'); %图2灰度图
u2=imnoise(y1,'salt & pepper',0.13);
subplot(1,4,3);imshow(u2);title('图像2叫噪声图'); %图2加椒盐噪声
zz1=medfilt2(u2,[3 3]);z1=medfilt2(zz1,[5 5]);%(2次中值滤波)
subplot(1,4,4);imshow(z1);title('图像2滤波图'); %图2中值滤波
% figure(5);subplot(1,2,1);imshow(zz1),title('一次滤波');
% subplot(1,2,2);imshow(z1);title('二次滤波');
% figure(3);subplot(1,2,1);imshow(z);title('处理后图1');
% figure(3);subplot(1,2,2);imshow(z1);title('处理后图2');
f=imsubtract(z1,z); %图1图2差分
figure(4);subplot(1,1,1);imshow(f);title('差分后图像');
% agin=medfilt2(f,[3 3]);figure(10);imshow(agin);f=agin; %差分图中值滤波
%-------------------------------------------------------------------------




%-------------------------------------------------------------------------
%眼睛的粗定位
%标注一只眼睛
[h w]=size(f);
Amax=0;m=0;n=0;
for i=10:h
for j=10:w
if Amax<=f(i,j)
Amax=f(i,j);m=i;n=j;
end
end
end
%标另外一只眼睛
Bmax=0;m1=0;n1=0;
for i1=10:m-10
for j1=10:w-10
if Bmax<=f(i1,j1)
Bmax=f(i1,j1);m1=i1;n1=j1;
end
end
end
for i2=m+10:h-10
for j2=10:w-10
if Bmax<=f(i2,j2)
Bmax=f(i2,j2);m1=i2;n1=j2;
end
end
end
%判断两眼的参数是否符合要求
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++






%----------------------------------------------------------------------

figure(7);imshow(z);hold on;plot(n,m,'+',n1,m1,'+');title('标定眼睛图');
%给眼睛画上矩形框
figure(8);imshow(z),title('眼睛粗定位');hold on;
Bvertex=[n-25 m-15;n+25 m+15]; %矩形的左上顶点坐标和右下顶点坐标
plot([Bvertex(1,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(1,2)],'r')
plot([Bvertex(2,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(2,2)],'r')
plot([Bvertex(2,1),Bvertex(1,1)],[Bvertex(2,2),Bvertex(2,2)],'r')
plot([Bvertex(1,1),Bvertex(1,1)],[Bvertex(1,2),Bvertex(2,2)],'r');
ho

ld on;
Bvertex=[n1-25 m1-15;n1+25 m1+15]; %矩形的左上顶点坐标和右下顶点坐标
plot([Bvertex(1,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(1,2)],'r')
plot([Bvertex(2,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(2,2)],'r')
plot([Bvertex(2,1),Bvertex(1,1)],[Bvertex(2,2),Bvertex(2,2)],'r')
plot([Bvertex(1,1),Bvertex(1,1)],[Bvertex(1,2),Bvertex(2,2)],'r');
%----------------------------------------------------------------------





%------------------------------------------------------------------------
%最大类间方差定二值化阈值
[height,width]=size(z);
hist=zeros(256,1);%直方图
%threshold=128; %初始阈值
%计算直方图
for i=1:height
for j=1:width
m=z(i,j)+1; %zz(i,j)灰度值从0-255
hist(m)=hist(m)+1;%某级灰度的像素数
end
end
hist=hist/(height*width);%落在每一灰度级上的概率
avg=0;
for m=1:256
avg=avg+(m-1)*hist(m);
end
temp=0;
for i=1:256
p1=0;
avg1=0;
avg2=0;
T_current=i-1;%当前分割阈值
for m=1:T_current-1
p1=hist(m)+ p1;%低灰度级概率总和
end
p2=1-p1;%高灰度级概率总和
for m=1:256
if mavg1=avg1+(m-1)*hist(m);
else
avg2=avg2+(m-1)*hist(m);
end
end
avg1=avg1/p1;
avg2=avg2/p2;

D=p1*(avg1-avg)^2+p2*(avg2-avg)^2;
if D>=temp
finalT=T_current;
temp=D;
end
end


% 滤波后差分图二值化
f1=z;
for i=1:h
for j=1:w
if f1(i,j)<=finalT
f1(i,j)=0;
else
f1(i,j)=255;
end
end
end
figure(13);imshow(f1);title('差分后图二值化');
%------------------------------------------------------------------------






%------------------------------------------------------------------------
%提取眼睛轮廓
%精确定位眼睛矩形，最高点、最低点、最左边、左右边四点确定的矩形。
%一只眼轮廓
[h1,w1]=size(f1);
zg=0;zd=0;zzb=0;zyb=0;
%i3=0;j3=0;
xj=50;yj=30;%垂直和水平方向中心点矩形大小一半
x3=zeros(1,100);y3=zeros(1,60);
x4=zeros(1,100);y4=zeros(1,60);
x5=zeros(1,100);y5=zeros(1,60);
x6=zeros(1,100);y6=zeros(1,60);
zgx=0;zdx=0;zzbx=0;zybx=0;
zgy=0;zdy=0;zzby=0;zyby=0;

%最高的点y
for i3=n-50:n+50
for j3=m-30:m+30
if f1(i3,j3)==255&f1(i3,j3+1)==0
x3(i3-n+51)=i3;y3(j3-m+31)=j3;
end
end
end

for k=1:1:100
if zgyzgy=y3;
end
end

%最低点y
for j3=m-30:m+30
for i3=n-50:n+50
if f1(i3,j3)==255&f1(i3,j3-1)==0
x4(i3-n+51)=i3;y4(j3-m+31)=j3;
end
end
end

for k=1:1:100
if zdyzdy=y4-1;
end
end

%最左点x
for j3=m-30:m+30
for i3=n+50:1:n-50
if f1(i3,j3)==255&f1(i3+1,j3)==0
x5(i3-n+51)=i3;y5(j3-m+31)=j3;
end
end
end

for k=1:1:100
if zzbxzzbx=x5+1;
end
end

%最右点x
for j3=m-30:m+30
for i3=n+50:1:n-50
if f1(i3,j3)==255&f1(i3-1,j3)==0
x6(i3-n+51)=i3;y6(j3-m+31)=j3;
end
end
end

for k=1:1:100
if zybxzybx=x6-1;
end
end

%精确画出眼睛矩形区域
figure(9);imshow(f1),title('眼睛轮廓');hold on;
Bvertex=[215 163;241 175]; %矩形的左上顶点坐标和右下顶点坐标
plot([Bvertex(1,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(1,2)],'r')
plot([Bvertex(2,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(2,2)],'r')
plot([Bvertex(2,1),Bvertex(1,1)],[Bvertex(2,2),Bvertex(2,2)],'r')
plot([Bvertex(1,1),Bvertex(1,1)],[Bvertex(1,2),Bvertex(2,2)],'r');
%---------------------------------------------------------------------------

%另一只眼轮廓
hold on;
Bvertex=[139 154;170 167]; %矩形的左上顶点坐标和右下顶点坐标
plot([Bvertex(1,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(1,2)],'r')
plot([Bvertex(2,1),Bvertex(2,1)],[Bvertex(1,2),Bvertex(2,2)],'r')
plot([Bvertex(2,1),Bvertex(1,1)],[Bvertex(2,2),Bvertex(2,2)],'r')
plot([Bvertex(1,1),Bvertex(1,1)],[Bvertex(1,2),Bvertex(2,2)],'r');
%---------------------------------------------------------------------------