一个简单的虹膜定位程序

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'); %圆显示