opencv 十五红外图像中虹膜的提取

一、算法需求

在医疗检测中，需要使用红外相机拍摄眼睛照片，然后提取出虹膜的区域。在拍摄过程瞳孔需要进行运动，其通常不在正前方，无法形成圆形，不能使用常规的霍夫圆检测进行提取定位。且在在红外图像中，虹膜区域与巩膜区域差别不明显(具体如下图所示)，故需要设计出算法提取红外图像中的虹膜区域。
在这里插入图片描述

1.1 眼睛结构说明

虹膜为圆盘状膜，中央有一黑孔称瞳孔，具体如下图所示。如果光线过强，虹膜内瞳孔括约肌收缩，则瞳孔缩小；光线变弱，虹膜开大肌收缩，瞳孔变大。
在这里插入图片描述

1.2 现有方法简述

通常使用霍夫圆检测实现瞳孔定位，具体实现效果如下所示。
在这里插入图片描述
其先通过二值化方法，获取瞳孔区域（包含闭运算操作，使瞳孔的圆闭合【瞳孔经常出现反光的情况】），最后在对瞳孔区域进行霍夫圆检测。

参考链接：https://blog.csdn.net/cungudafa/article/details/119726505

使用opencv的椭圆检测进行定位时发现以下情况，当瞳孔运动到眼球边缘时，其无法准确的检测到瞳孔（霍夫圆检测的黄色圆与瞳孔区域没有严格的贴合）。
在这里插入图片描述
使用椭圆拟合则可以准确的圈出瞳孔区域

二、问题分析

对现有的多个数据进行分析发现，眼球照片有以下特点：
1、在红外图像中虹膜与巩膜区域没有显著性差异性===》不可以使用现有虹膜提取方法
2、虹膜以瞳孔为中心，跟随瞳孔运动方向进行同步移动===》可以将虹膜提取转化为瞳孔提取

三、核心思路

1、读取图片为灰度图，并优化图像质量（使用滤波尽可能减少图像背景的复杂度）
2、对图像进行二值化（其可以根据调试效果设置二值化阈值，瞳孔区域与眼球其他区域存在显著的颜色差异）
3、对瞳孔区域进行优化（使用闭运算移除瞳孔中的反光区域）
4、获取图像中的轮廓，并进行椭圆拟合，并根据拟合结果排除错误的椭圆（根据拟合椭圆长轴与短轴值判定）
5、根据瞳孔与虹膜的半径比假定虹膜的椭圆轴长绘制椭圆mask，在原图中截取出虹膜区域。

四、具体实现

读取后的图片如下所示
在这里插入图片描述
进行二值化后得到以下图像，可以看到瞳孔中存在黑洞，其他区域存在白色干扰点。

先找到图像中最大面积的连通域，然后进行闭运算，最终得到的结果如下所示

然后获取轮廓并进行椭圆拟合，然后将拟合的椭圆绘制在原图与mask上（画在原图上的椭圆要使用原始值，而画在mask上的椭圆需要对长轴和短轴值进行放大，使其能尽可能的盖住虹膜区域）在这里插入图片描述
使用mask与原图进行与运算可以得到以下结果

根据连通域获取外接矩形，将虹膜区域裁剪出来得到以下图片

五、完整代码

完整代码如下所示

import numpy as np
import cv2
from matplotlib import pyplot as plt
import osdef find_topK_areo(img,k=1):ret,result=cv2.threshold(img,128,255,cv2.THRESH_BINARY)contours, hierarchy = cv2.findContours(result,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)#找到最大面积连通域areos=[]for i in range(len(contours)):area = cv2.contourArea(contours[i])areos.append({'area':area,'id':i})areos.sort(key=lambda x:x['area'],reverse=True)topk_areo=areos[:k]black=np.zeros(result.shape,np.uint8)for f in topk_areo:cv2.drawContours(black,contours,f['id'],(255,255,255),-1)return blackdef getHoughCircle(img):blur = cv2.GaussianBlur(img, (3, 3), 5) # 高斯模糊，给出高斯模糊矩阵和标准差gray = cv2.cvtColor(blur, cv2.COLOR_BGR2GRAY)# 灰度化# 图像二值化，全局自适应阈值:对输入的单通道矩阵逐像素进行阈值分割#ret, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_TRIANGLE)ret,binary=cv2.threshold(gray, 50, 255, cv2.THRESH_BINARY_INV)dst=find_topK_areo(binary,1)kernel=np.ones((3,3),np.uint8)dst_close=cv2.morphologyEx(dst, cv2.MORPH_CLOSE, kernel)cnt, hierarchy = cv2.findContours(dst_close, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)print(len(cnt))  #   得到该图中总的轮廓数量mask= np.zeros(dst.shape, np.uint8 )#生成全黑的maskfor i in range(len(cnt)):# 椭圆拟合#（x, y）代表椭圆中心点的位置,（a, b）代表长短轴长度，应注意a、b为长短轴的直径，而非半径,angle 代表了中心旋转的角度ellipse= cv2.fitEllipse(cnt[i])(cx, cy), (a, b), angle=ellipseprint((cx, cy), (a, b), angle) #椭圆拟合结果有一些非瞳孔区域，需要跳过。经过观察，其a与b的值特别小if a+b<40:continue# 绘制椭圆，使用ellipse(img, ellipse,color, 2)方法，不要使用另外一种多参数的用法cv2.ellipse(img, ellipse,(0,0,255), 2)cv2.drawContours(img,cnt,i,(0,0,255),1)#将椭圆区域进行放大，使其转换虹膜的mask圆ellipse=((cx, cy), (a*2.5, b*2.5), angle)cv2.ellipse(mask, ellipse,(255,255,255), -1)res=cv2.bitwise_and(gray,mask) #与灰度图进行与运算，提取目标区域（虹膜)contours, _ = cv2.findContours(mask,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)for cont in contours:# 外接矩形x, y, w, h = cv2.boundingRect(cont)#裁剪出椭圆区域crop=res[y:y+h,x:x+w]return crop if __name__=="__main__":path = 'vedio/tor3.avi/'crop1=getHoughCircle(cv2.imread(path+'tor3.avi1.jpg',1))cv2.imshow('crop',crop1)cv2.waitKey(0)