我下载了daugman的虹膜定位算法，我运行源文件的图可以出结果，但是运行自己的图就出不了结果。
救救孩子吧。
function [ci,cp,out]=thresh(I,rmin,rmax)
%proceed by scaling down images
scale=1;
rmin=rmin*scale;
rmax=rmax*scale;
%convert image to double
I=im2double(I);
pimage=I;
I=imresize(I,scale);
%help remove reflection.
I=imcomplement(imfill(imcomplement(I),'holes'));
rows=size(I,1);
cols=size(I,2);
[X,Y]=find(I<0.5);
s=size(X,1);
for k=1:s
if (X(k)>rmin)&&(Y(k)>rmin)&&(X(k)<=(rows-rmin))&&(Y(k)<(cols-rmin))
A=I((X(k)-1):(X(k)+1),(Y(k)-1):(Y(k)+1));
M=min(min(A));
if I(X(k),Y(k))~=M
X(k)=NaN;
Y(k)=NaN;
end
end
end
v=find(isnan(X));
X(v)=[];
Y(v)=[];
index=find((X<=rmin)|(Y<=rmin)|(X>(rows-rmin))|(Y>(cols-rmin)));
X(index)=[];
Y(index)=[];
N=size(X,1);
%defines two arrays to store the maximum value of blur the corresponding radius
maxb=zeros(rows,cols);
maxrad=zeros(rows,cols);
for j=1:N
[b,r,~]=partiald(I,[X(j),Y(j)],rmin,rmax,'inf',600,'iris');
maxb(X(j),Y(j))=b;
maxrad(X(j),Y(j))=r;
end
[x,y]=find(maxb==max(max(maxb)));
%the function search searches for the centre of the pupil and its radius by scanning a 10*10 window around the iris centre for establishing the pupil's centre and hence its radius
ci=search(I,rmin,rmax,x,y,'iris');
ci=ci/scale;
cp=search(I,round(0.1*r),round(0.8*r),ci(1)*scale,ci(2)*scale,'pupil');%Ref:Daugman's paper that sets biological limits on the relative sizes of the iris and pupil
cp=cp/scale;
out=drawcircle(pimage,[ci(1),ci(2)],ci(3),600);
out=drawcircle(out,[cp(1),cp(2)],cp(3),600);
```function [cp]=search(im,rmin,rmax,x,y,option)
rows=size(im,1);
cols=size(im,2);
sigma=0.5;%(standard deviation of Gaussian)
R=rmin:rmax;
maxrad=zeros(rows,cols);
maxb=zeros(rows,cols);
for i=(x-5):(x+5)
for j=(y-5):(y+5)
[b,r,~]=partiald(im,[i,j],rmin,rmax,0.5,600,option);
maxrad(i,j)=r;
maxb(i,j)=b;
end
end
B=max(max(maxb));
[X,Y]=find(maxb==B);
radius=maxrad(X,Y);
cp=[X,Y,radius];
function [b,r,blur]=partiald(I,C,rmin,rmax,sigma,n,part)
R=rmin:rmax;
count=size(R,2);
for k=1:count
[L(k)]=lineint(I,C,R(k),n,part);%computing the normalized line integral for each radius
if L(k)==0%if L(k)==0(this case occurs iff the radius takes the circle out of the image)
%In this case,L is deleted as shown below and no more radii are taken for computation
%(for that particular centre point).This is accomplished using the break statement
L(k)=[];
break;
end
end
D=diff(L);
D=[0 D];
%append one element at the beginning to make it an n vector
%Partial derivative at rmin is assumed to be zero
if strcmp(sigma,'inf')==1%the limiting case of the gaussian with sigma infinity(pls remember to change the code)strcmp syntax is different
f=ones(1,7)/7;
else
f=fspecial('gaussian',[1,5],sigma);%generates a 5 member 1-D gaussian
end
blur=convn(D,f,'same');%Smooths the D vecor by 1-D convolution
%'same' indicates that size(blur) equals size(D)
blur=abs(blur);
[b,i]=max(blur);
r=R(i);
b=blur(i);
%calculates the blurred partial derivative
function [L]=lineint(I,C,r,n,part)
theta=(2*pi)/n;% angle subtended at the centre by the sides
%orient one of the radii to lie along the y axis
%positive angle is ccw
%Author:Anirudh S.K.
%Department of Computer Science and Engineering
%Indian Institute of Techology,Madras
rows=size(I,1);
cols=size(I,2);
angle=theta:theta:2*pi;
x=C(1)-r*sin(angle);
y=C(2)+r*cos(angle);
if (any(x>=rows)||any(y>=cols)||any(x<=1)||any(y<=1))
L=0;
return
%This process returns L=0 for any circle that does not fit inside the image
end
%lines 34 to 42 compute the whole line integral
if (strcmp(part,'pupil')==1)
s=0;
for i=1:n
val=I(round(x(i)),round(y(i)));
s=s+val;
end
L=s/n;
end
%lines 44 onwards compute the lateral line integral(to prevent occlusion affecting the results,the pixel average is taken only along the lateral portions)
if(strcmp(part,'iris')==1)
s=0;
for i=1:round((n/8))
val=I(round(x(i)),round(y(i)));
s=s+val;
end
for i=(round(3*n/8))+1:round((5*n/8))
val=I(round(x(i)),round(y(i)));
s=s+val;
end
for i=round((7*n/8))+1:(n)
val=I(round(x(i)),round(y(i)));
s=s+val;
end
L=(2*s)/n;
end
function [O]=drawcircle(I,C,r,n)
if nargin==3
n=600;
end
theta=(2*pi)/n;% angle subtended at the centre by the sides
%orient one of the radii to lie along the y axis
%positive angle ic ccw
rows=size(I,1);
cols=size(I,2);
angle=theta:theta:2*pi;
%to improve contrast and help in detection
x=C(1)-r*sin(angle);%the negative sign occurs because of the particular choice of coordinate system
y=C(2)+r*cos(angle);
if any(x>=rows)||any(y>=cols)||any(x<=1)||any(y<=1)%if circle is out of bounds return image itself
O=I;
return
end
for i=1:n
I(round(x(i)),round(y(i)))=1;
end
O=I;
I=imread('F:\iris2.jpg');
%设置最小半径值和最大半径值
rmin=35;
rmax=400;
%瞳孔巩膜虹膜分割的调用函数
[ci,cp,out]=thresh(I,rmin,rmax);
%显示结果
figure;
imshow(out);