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 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)