opencv 十五 红外图像中虹膜的提取
作者:mmseoamin日期:2023-12-14

一、算法需求

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

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第1张

1.1 眼睛结构说明

虹膜为圆盘状膜,中央有一黑孔称瞳孔,具体如下图所示。如果光线过强,虹膜内瞳孔括约肌收缩,则瞳孔缩小;光线变弱,虹膜开大肌收缩,瞳孔变大。

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第2张

1.2 现有方法简述

通常使用霍夫圆检测实现瞳孔定位,具体实现效果如下所示。

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第3张

其先通过二值化方法,获取瞳孔区域(包含闭运算操作,使瞳孔的圆闭合【瞳孔经常出现反光的情况】),最后在对瞳孔区域进行霍夫圆检测。

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

使用opencv的椭圆检测进行定位时发现以下情况,当瞳孔运动到眼球边缘时,其无法准确的检测到瞳孔(霍夫圆检测的黄色圆与瞳孔区域没有严格的贴合)。

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第4张

使用椭圆拟合则可以准确的圈出瞳孔区域

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第5张

二、问题分析

对现有的多个数据进行分析发现,眼球照片有以下特点:

1、在红外图像中虹膜与巩膜区域没有显著性差异性===》不可以使用现有虹膜提取方法

2、虹膜以瞳孔为中心,跟随瞳孔运动方向进行同步移动===》可以将虹膜提取转化为瞳孔提取

三、核心思路

1、读取图片为灰度图,并优化图像质量(使用滤波尽可能减少图像背景的复杂度)

2、对图像进行二值化(其可以根据调试效果设置二值化阈值,瞳孔区域与眼球其他区域存在显著的颜色差异)

3、对瞳孔区域进行优化(使用闭运算移除瞳孔中的反光区域)

4、获取图像中的轮廓,并进行椭圆拟合,并根据拟合结果排除错误的椭圆(根据拟合椭圆长轴与短轴值判定)

5、根据瞳孔与虹膜的半径比假定虹膜的椭圆轴长绘制椭圆mask,在原图中截取出虹膜区域。

四、具体实现

读取后的图片如下所示

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第6张

进行二值化后得到以下图像,可以看到瞳孔中存在黑洞,其他区域存在白色干扰点。

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第7张

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

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第8张

然后获取轮廓并进行椭圆拟合,然后将拟合的椭圆绘制在原图与mask上(画在原图上的椭圆要使用原始值,而画在mask上的椭圆需要对长轴和短轴值进行放大,使其能尽可能的盖住虹膜区域)opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第9张

使用mask与原图进行与运算可以得到以下结果

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第10张

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

opencv 十五 红外图像中虹膜的提取,在这里插入图片描述,第11张

五、完整代码

完整代码如下所示

import numpy as np
import cv2
from matplotlib import pyplot as plt
import os
    
def 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 black
    
def 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 )#生成全黑的mask
    
    for i in range(len(cnt)):
        # 椭圆拟合
        #(x, y)代表椭圆中心点的位置,(a, b)代表长短轴长度,应注意a、b为长短轴的直径,而非半径,angle 代表了中心旋转的角度
         
        ellipse= cv2.fitEllipse(cnt[i])
        (cx, cy), (a, b), angle=ellipse
        print((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)