一个简单的虹膜定位程序
2019-04-14 16:02发布
生成海报
一个简单的虹膜定位实现
%.虹膜定位程序
clear;close all;
I=imread('ip1.jpg')
f = rgb2gray(I);
imhist(f); % 求图像的直方图
[F_Size_M F_Size_N] = size(f); % 获取f的行和列
T =84;
for i = 1:F_Size_M
for j = 1:F_Size_N
if f(i,j) >= T
f(i,j) = 255;
else
f(i,j) = 0;
end
end
end
imshow(f);
f = ~f;
se = strel('disk',3);
f1 = imclose(f,se);
se = strel('disk',3);
f2 = imopen(f1,se);
f3 = imfill(f2,'holes');
f4 = im2double(f3);
f4 = medfilt2(f3,[3 3]);
% imshow(f4);
% 标记,查找,分离最大联通域
[x,num] = bwlabel(f4,4); % 标记联通域
% 标记每一个联通域的像素个数
for j = 1 : num
[r,c] = find(x == j);
rc = [r c];
long(j) = size(rc,1);
end
% 查找最大联通域
max = 0;
m=0;
for i = 1 : num
if long(i) > max
max = long(i);
m= i;
end
end
% 分离最大联通域
[row column] = size(x);
for i = 1 : row
for j = 1 : column
if x(i,j) == m
x(i,j) = 1;
else
x(i,j) = 0;
end
end
end
figure,imshow(x);
% % 利用cann算子进行边缘检测,以得到单像素的边缘图像
k = edge(x,'canny'); % 边缘检测
[m,n] = find(k == 1); % 记录边缘的坐标
figure,imshow(I);
D=[m,n];
row=length(D);%记录边缘点个数
[R,A,B]=circ(m,n,row);调用函数,圆拟合的定位显示
最小二乘圆拟合程序
function [R,A,B]=circ(x,y,N)
x1 = 0; %给定初值
x2 = 0;
x3 = 0;
y1 = 0;
y2 = 0;
y3 = 0;
x1y1 = 0;
x1y2 = 0;
x2y1 = 0;
for i = 1 : N %公式计算
x1 = x1 + x(i);
x2 = x2 + x(i)*x(i);
x3 = x3 + x(i)*x(i)*x(i);
y1 = y1 + y(i);
y2 = y2 + y(i)*y(i);
y3 = y3 + y(i)*y(i)*y(i);
x1y1 = x1y1 + x(i)*y(i);
x1y2 = x1y2 + x(i)*y(i)*y(i);
x2y1 = x2y1 + x(i)*x(i)*y(i);
end
C = N * x2 - x1 * x1;
D = N * x1y1 - x1 * y1;
E = N * x3 + N * x1y2 - (x2 + y2) * x1;
G = N * y2 - y1 * y1;
H = N * x2y1 + N * y3 - (x2 + y2) * y1;
a = (H * D - E * G)/(C * G - D * D);
b = (H * C - E * D)/(D * D - G * C);
c = -(a * x1 + b * y1 + x2 + y2)/N;
A = a/(-2); %y坐标
B = b/(-2); %x坐标
R = sqrt(a * a + b * b - 4 * c)/2;
hold on;
alpha=0:pi/20:2*pi;%角度[0,2*pi]
y=A+R*cos(alpha);
x=B+R*sin(alpha);
plot(B,A,'b-'); %标注圆心
hold on
plot(x,y,'r'); %圆显示
打开微信“扫一扫”,打开网页后点击屏幕右上角分享按钮