虹膜图像分割（一）内圆检测

        本文和下文介绍如何对一个虹膜图像进行分割提取出虹膜，在虹膜分割中，最重要的是检测两个圆，一个内圆，一个是外圆。下面是两个圆的示意图。

本文先讲解如何检测内圆，在我的方法中，内圆的检测可分为四步：

1.  用Canny边缘检测和Hough变换从经过高斯平滑处理后的图片中检测圆。在此步操作，我故意把Hough变换的阈值设得比较小，这样我们就可以检测出很多圆了。下图就是用cvHoughCircle检测出的圆。

可以看到，在此步中，我们检测到了圆，但是有很多圆，我们现在要想办法找到我们要的那个内圆。下面的步骤就是我找内圆的方法。

2.   把原图二值化，这个可以通过cvThreshold很容易实现，阈值设在五十左右就可以了。阈值化后的图像如下：

这个二值化后得到的图像是为了后面计算每个圆中包含的瞳孔点比率做准备的。

3.  计算第一步中得到的每个圆的一个比率，计算是根据第二步中二值图像进行。比率的定义如下

PupilInclusion Rate = Number of black points / Area of circle

即：用每一个circle中包含的黑点的个数除以圆的面积

4. 在计算出每一个圆的比率后，找到比率最大的那个，并把那个作为最后的内圆，这样就可以完美得检测出内圆了。

//caculate the ratio based on binary image
float CheckCirle(IplImage* img)
{
	CvPoint center;
	int width = img->width;
	int height = img->height;
	int r = width/2;
	int value; //pixel value
	int count = 0;
	float ratio = 0;
	center.x = cvRound(width/2);
	center.y = cvRound(height/2);
	for(int i=0;i<height;i++)
	{
		for(int j=0;j<width;j++)
		{
			if (sqrt(pow(float(center.x-j),2)+pow(float(center.y-i),2))< r)
			{
				value = cvGetReal2D(img,i,j);
				if(value == 0)
					count++;
			}
		}
	}
	ratio = float(count)/(3.14*r*r);
	return ratio;
}
//detect the inner circle
float* FindInnerCircle(IplImage* src)
{
	float result[3]; //r[0] r[2] are the coordinate, r[3] is radius

	float maxratio = 0;
	float ratio = 0;

	IplImage* dst = NULL;

	int threshold = FindThreshold(src);

	dst = cvCreateImage (cvGetSize(src), IPL_DEPTH_8U, 1);
	dst = cvCloneImage (src);
	
	cvSmooth (dst, dst, CV_GAUSSIAN, 5, 5);
	cvSmooth(dst,dst,CV_MEDIAN,3);


	CvMemStorage* storage = cvCreateMemStorage (0);
	//CvSeq* circles = cvHoughCircles (dst, storage, CV_HOUGH_GRADIENT, 2, dst->width / 30, 50, 80, 10, 70);
	//50 img04; 100 img009
	CvSeq* circles = cvHoughCircles (dst, storage, CV_HOUGH_GRADIENT, 2, dst->width / 30, 50, 70, 20, 90);

	//60 for 009; 60 img12;  100 img16,img19, img20
	cvThreshold(dst,dst,100,255,CV_THRESH_BINARY);
	cvNamedWindow("threshold",1);
	cvShowImage("threshold",dst);

	CvPoint center;
	int r;

	if (circles->total<0)
	{
		printf("No Circle Detected!!Please Check！！\n");
		system("pause");
	}
		
	for (int i = 0; i < circles->total; i++)
		//for (int i = 0; i < 35; i++)
	{
	 
		float* p = (float*)cvGetSeqElem (circles, i);
		CvPoint pt = cvPoint (cvRound(p[0]), cvRound(p[1]));
		center = pt;
		r = cvRound(p[2]);

		cvSetImageROI(dst,cvRect(center.x-r,center.y-r,2*r,2*r));
		IplImage* rect_img=NULL;
		rect_img = cvCreateImage(cvGetSize(dst),8,1);
		cvCopy(dst,rect_img);
		cvResetImageROI(dst);

		ratio = CheckCirle(rect_img);
		if (ratio >= maxratio)
		{
			result[0] = pt.x;
			result[1] = pt.y;
			result[2] = r;
			maxratio = ratio;
		}
		rect_img = NULL;
	}

	cvNamedWindow ("src", 1);
	cvShowImage ("src", src);

	printf("Inner circle have been found!!\n");
	cvReleaseMemStorage (&storage);
	return result;
}

Result:

完整的测试图片及代码下载：

http://download.csdn.net/detail/computerme/8158905