虹膜分割算法

由上一篇https://blog.csdn.net/zqx951102/article/details/83177029
虹膜识别外圆检测 得到的画出虹膜内外圆图像：

然后需要去掉外圆之外的图像和内圆之内的瞳孔图像，然后剩下的就是我们想得到的虹膜图像了，算法思路就是遍历图像上所有的像素点，然后判断该点的距离与内圆圆心的距离是大于外圆的半径还是小于内圆的半径，然后满足这样条件的点，就修改像素值为白色，这样把图像输出就是我们需要的虹膜图像了。
部分代码如下：
我们需要新建一个图像srcImage2.然后在这个图片上进行遍历像素点。

//虹膜分割Mat srcImage2=srcImage1.clone();int width = srcImage2.cols;int height = srcImage2.rows;	for(int i=0;i<height;i++){for(int j=0;j<width;j++){if (sqrt(pow(float(center1.x-j),2)+pow(float(center1.y-i),2))>outer_r){   //value =dst.at<Vec3b>(ti,tj)[2];  RGB图像是三通道的BGR 所以0 1 2 三个通道都是255时候才是显示的白色！srcImage2.at<Vec3b>(i,j)[0]=255;srcImage2.at<Vec3b>(i,j)[1]=255; srcImage2.at<Vec3b>(i,j)[2]=255;//外圆取大于}}}for(int i=0;i<height;i++){for(int j=0;j<width;j++){if (sqrt(pow(float(center1.x-j),2)+pow(float(center1.y-i),2))<result[2]){   //value =dst.at<Vec3b>(ti,tj)[2];srcImage2.at<Vec3b>(i,j)[0]=255;srcImage2.at<Vec3b>(i,j)[1]=255; srcImage2.at<Vec3b>(i,j)[2]=255;  //内圆取小于}}}imshow("虹膜分割效果图",srcImage2);

效果图：

这个就是我们所需要的虹膜分割图像！
完整代码如下：

#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include<math.h>
using namespace cv;
using namespace std;int main()
{Mat srcImage=imread("C://1.bmp");Mat midImage,dstImage,edge;Mat srcImage1=srcImage.clone();imshow("原始图",srcImage);//外圆	 Mat kern = (Mat_<char>(5,3) << -1, 0 ,1,-1, 0 ,1,-1, 0 ,1,-1, 0 ,1,-1, 0 ,1);Mat dst;filter2D(srcImage,dst,srcImage.depth(),kern);namedWindow("卷积",WINDOW_AUTOSIZE);imshow("卷积",dst);threshold(dst,dst,10,255,CV_THRESH_BINARY);namedWindow("二值化",WINDOW_AUTOSIZE);imshow("二值化",dst);threshold(srcImage, srcImage, 30, 200.0, CV_THRESH_BINARY);//二值化imshow("二值化1",srcImage);int cn=0;//cn是圆的个数int radius=0;float ratio = 0;float maxratio = 0;float result[3];cvtColor(srcImage,midImage,COLOR_BGR2GRAY);blur(midImage,edge,Size(3,3));Canny(edge,edge,3,9,3);GaussianBlur(midImage,midImage,Size(9,9),2,2);vector<Vec3f> circles;HoughCircles(midImage,circles,CV_HOUGH_GRADIENT,2,10,200,80,0,0);for(cn=0;cn<circles.size();cn++){Point center(cvRound(circles[cn][0]),cvRound(circles[cn][1]));radius=cvRound(circles[cn][2]);//	circle(srcImage,center,3,Scalar(0,255,0),-1,8,0);//	circle(srcImage,center,radius,Scalar(155,50,255),3,8,0);//int width = srcImage.cols;int height = srcImage.rows;	int value; //pixel valueint count = 0;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))< radius){//获得某点的像素值value = srcImage.at<Vec3b>(i,j)[2];  //cvGetReal2D(img,i,j);if(value == 0)count++;}}ratio = float(count)/(3.14*radius*radius);if (ratio >= maxratio){result[0] = circles[cn][0];result[1] = circles[cn][1];result[2] = radius;maxratio = ratio;}}printf("黑色点像素的个数：%d\n",count);printf("瞳孔重合比率：%f\n",ratio);Point center1(cvRound(result[0]),cvRound(result[1]));circle(srcImage1,center1,3,Scalar(255,0,0),-1,8,0);circle(srcImage1,center1,result[2],Scalar(255,255,255),1,8,0);}//外圆int r=result[2];int value;int tmp_count = 0;	 int test=0;int outer_r = 0;int max_count = 0;int ti,tj;for(double tr = 1.8*r;tr< 3*r;tr++ ){tmp_count = 0;for(double angle = 0;angle<=60;angle++){ti = cvRound(result[1] + tr*cos(angle));tj = cvRound(result[0]+tr*sin(angle));if ((ti < dst.rows) && (ti>0)&&(dst.cols)&&(tj>0))value =dst.at<Vec3b>(ti,tj)[2];else break;if( value == 255)tmp_count++;}if (tmp_count>=max_count){max_count = tmp_count;outer_r = tr;test++;}}printf("outer radius = %d\n",outer_r); Point center1(cvRound(result[0]),cvRound(result[1]));circle(srcImage1,center1,3,Scalar(255,0,0),-1,8,0);circle(srcImage1,center1,outer_r,Scalar(255,255,255),1,8,0);imshow("画内外圆效果图",srcImage1);/*if(cn==0){printf("No Circle Detected!!Please Check！！\n");system("pause");}*///虹膜分割Mat srcImage2=srcImage1.clone();int width = srcImage2.cols;int height = srcImage2.rows;	for(int i=0;i<height;i++){for(int j=0;j<width;j++){if (sqrt(pow(float(center1.x-j),2)+pow(float(center1.y-i),2))>outer_r){   //value =dst.at<Vec3b>(ti,tj)[2];  图像是三通道的rgb 所以0 1 2 三个通道都是255时候才是显示的白色！srcImage2.at<Vec3b>(i,j)[0]=255;srcImage2.at<Vec3b>(i,j)[1]=255; srcImage2.at<Vec3b>(i,j)[2]=255;//外圆取大于}}}for(int i=0;i<height;i++){for(int j=0;j<width;j++){if (sqrt(pow(float(center1.x-j),2)+pow(float(center1.y-i),2))<result[2]){   //value =dst.at<Vec3b>(ti,tj)[2];srcImage2.at<Vec3b>(i,j)[0]=255;srcImage2.at<Vec3b>(i,j)[1]=255; srcImage2.at<Vec3b>(i,j)[2]=255;  //内圆取小于}}}imshow("虹膜分割效果图",srcImage2);waitKey(0);return 0;
}