眼底血管分割(python版本)

一、项目背景

       研究表明，各类眼科疾病以及心脑血管疾病会对视网膜血管造成形变、出血等不同程度的影响。随着生活水平的提高，这类疾病的发病率呈现逐年增长的趋势。临床上，医疗人员能够从检眼镜采集的彩色眼底图像中提取视网膜血管，然后通过对血管形态状况的分析达到诊断这类疾病的目的。但是，由于受眼底图像采集技术的限制，图像中往往存在大量噪声，再加之视网膜血管自身结构复杂多变，使得视网膜血管的分割变得困难重重。传统方法中依靠人工手动分割视网膜血管，不仅工作量巨大极为耗时，而且受主观因素影响严重。因此，利用计算机技术，找到一种能够快速、准确分割视网膜血管的算法，实现对眼底图像血管特征的实时提取，对辅助医疗人员诊断眼科疾病、心脑血管疾病等具有重要作用。

二、项目流程

0.流程图

1.视网膜眼底图像预处理

1.1 图像G通道输出
       眼底视网膜图像在获取过程中，往往受到外部条件的影响，使得采集后的图像出现灰度分布不均匀现象。再加上视网膜图像血管与背景对比度较低，使得视网膜血管难以被检测到，为了提高视网膜图像的对比度，改善视网膜图像质量，在视网膜特征提取之前需要进行预处理，在图像特征提取过程中，这一步的好坏能够直接影响到后续处理的效果。
       而图像预处理的方法有很多，比较常见的有：图像对比度增强法，图像亮度校正法以及图像归一化方法等，其中使用最广泛的是图像对比度增强法。而在图像对比度增强法中，直方图均衡又是一种常用的方法。该方法首先求解图像的全局以及局部灰度分布情况，建立图像直方图校正信息，进而达到提高视网膜图像对比度的目的。建立在这之上的，是自适应直方图均衡法。它对普通的直方图均衡算法进行了改进，即对图像的对比度进行了重新分配，所以使用该方法可以获得更多的图像细节。
       但在进行直方图均衡之前，我们首先要对原图进行处理。因为网膜血管的提取会受到众多因素的干扰，所获得的图像质量差别很大，所以，为了能够使算法很好的适用于不同的图像，在对图像预处理之前先对其进行色彩空间变换和颜色通道的提取。常用的易于分离亮度通道的色彩空间有RGB、Lab、YCb Cr和 Gaussian 色彩空间。下图为RGB和Lab通道输出的图像。

图一   原图与RGB和Lab通道输出的图像

       通过分析视网膜彩色图像，及其对应的RGB、Lab通道图像，可以发现绿色通道视网膜图像拥有最佳的目标血管与图像背景的对比度特征，而其他通道视网膜图像的质量较差，血管与图像对比度特不明显。因此，在彩色视网膜图像的预处理过程中，提取了视网膜的绿色通道。


1.2高斯滤波
       高斯滤波器是一种线性滤波器，能够有效的抑制噪声，平滑图像。其作用原理和均值滤波器类似，都是取滤波器窗口内的像素的均值作为输出。其窗口模板的系数和均值滤波器不同，均值滤波器的模板系数都是相同的为1；而高斯滤波器的模板系数，则随着距离模板中心的增大而系数减小。所以，高斯滤波器相比于均值滤波器对图像个模糊程度较小。基于二维高斯分布(如下式)，使用5*5模板进行滤波，得到滤波后的视网膜血管图。

图二   高斯滤波后的视网膜

1.3 消除不均匀光照
       消除不均匀光照的方法有很多种，比如底帽变换，对比度受限的自适应直方图均衡，同态滤波等。

1.3.1 底帽变换
       通过对f的闭操作减去f以达到消除不均匀光照的目的:

1.3.2 对比度受限的自适应直方图均衡化(CLAHE)
       直方图均衡化（HE）是一种很常用的直方图类方法，基本思想是通过图像的灰度分布直方图确定一条映射曲线，用来对图像进行灰度变换，以达到提高图像 对比度的目的。该映射曲线其实就是图像的累计分布直方图（CDF）（严格来说是呈正比例关系）。然而HE是对图像全局进行调整的方法，不能有效地提高局部 对比度，而且某些场合效果会非常差。为了提高图像的局部对比度，有人提出将图像分成若干子块，对子块进行HE处理，这便是AHE（自适应直方图均衡化）。

1.3.3 同态滤波
       同态滤波是一种在频域中进行的图像对比度增强和压缩图像亮度范围的特殊方法。同态滤波器能够减少低频并且增加高频，从而能减少光照变化并锐化边缘细节。图像的同态滤波技术的依据是图像获取过程中的照明反射成像原理。它属于频域处理,作用是对图像灰度范围进行调整,通过消除图像上照明不均的问题。非线性滤波器能够在很好地保护细节的同时,去除信号中的噪声,同态滤波器就是一种非线性滤波器,其处理是一种基于特征的对比度增强方法,主要用于减少由于光照不均匀引起的图像降质,并对感兴趣的景物进行有效地增强。

       在这里，为了进一步提高视网膜血管与背景的对比度，本实验采用了对比度受限的自适应直方图均衡法。该方法通过限制局部直方图的高度来控制局部对比度的增强幅度，不仅保证了视网膜绿色通道图像的局部对比度增强，还能够进一步的限制视网膜噪声放大。

G通道的自适应直方图均衡

2．视网膜特征提取

2.1高斯匹配滤波与图像归一化处理
       当光线通过眼底视网膜时，由于视网膜表皮组织与视网膜血管的反射系数不同，导致视网膜图像中血管呈现较暗的亮度特性。而且通过分析视网膜血管横截面信息可以知道，视网膜血管的灰度分布可以用高斯匹配滤波函数进行模拟，即利用高斯匹配滤波函数的形状以及走向，来近似视网膜图像血管。
       此外，当高斯匹配滤波器的方向与视网膜血管走向接近时，利用高斯匹配滤波器与视网膜图像进行卷积，能够获得较大的响应值。所以通过对高斯匹配滤波器进行不同方向的调整，即以相同角度间隔对高斯匹配滤波器进行旋转，能够实现对视网膜血管的追踪检测。在具体过程中，高斯匹配滤波的实现是对视网膜图像上的单个像素点，利用高斯匹配滤波器进行处理，然后通过提取高斯匹配滤波所有方向最大响应值完成的。

图三   15x15匹配滤波窗口图例

       图像经过匹配滤波处理之后，可抑制眼底图象中的背景，增强目标血管灰度，达到凸显血管结构的目的。但是为了保证眼底图象的仿射不变性，提高计算精度，还得对高斯匹配滤波后的图形进行归一化处理，如下图

图四    高斯匹配滤波+归一化处理

2.2去轮廓
       在最开始的视网膜原图基础上，利用灰度变换，获取图形的掩膜，如图四。将带轮廓的视网膜图像(图四)与其相减，再经形态学的膨胀处理，即可获得清晰的轮廓线，如图五。

然后，将图四的结果减去轮廓图，即可获得去除轮廓后的视网膜血管图。

图七   灰度拉伸后的图像

2.4 otsu二值化
       OTSU算法也称最大类间差法，有时也称之为大津算法，被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。
       它是按图像的灰度特性,将图像分成背景和前景两部分。背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。这也是前面进行灰度拉伸的原因：使得背景和前景之间的方差变大。
       实现方法是通过最大类间方差法获得阈值，再将灰度图像转换为二值图像。

图八   otsu二值化

       最后，考虑到最终分割图的噪声很少，因此不对其进行滤波，以及形态学操作。因为一旦进行滤波或者形态学操作，虽然达到了目的，但是也使得血管出现断裂，损失更多的细节。

3. 图像分割结果的评估—DICE指标

       血管分割了之后，我们要对分割的结果进行评价，虽然通过人眼可以看出结果好坏，但是这只是主观的感知。因此，我们需要对分割的结果进行量化。Dice 系数是一种评估相似度的函数，通常用于计算两个样本的相似度或者重叠度：

       一般来说，Dice指标越大，则说明两幅图像越相似。因此，我们将专家手动分割的结果与本实验分割的结果进行计算，最终算得20组样本的平均dice指标为0.735。具体指标详见附录一。

三、实验结果

       基于以上算法，我们分割了将近20组样本，效果甚佳。其中十组样本得分割结果如下。

四、结果分析

       本实验一共给出了10组视网膜图像的分割结果。第一列为视网膜原始图像，第二列为专家手动分割的视网膜标准参考图像，第三列为本实验算法分割的结果。对比分割结果，可以看出，本人的算法对视网膜细血管分割不够准确，且丢失了部分血管的连通性。原因可能是灰度拉伸，以及otsu二值化的阈值设定不准确导致的。

五、参考文献

[1]马青柯. 眼底图像检测和分析系统的设计与开发[D].暨南大学,2018.
[2]王晓红. 基于特征识别的视网膜血管分割方法研究[D].中南大学,2014.
[3]尚星宇. 视网膜眼底病变图像血管提取方法[D].湘潭大学,2012.
[4]康文炜. 冠状动脉造影图像的分割方法研究[D]. 吉林大学,2015.

python源代码

from __future__ import print_function

from multiprocessing.pool import ThreadPool

import cv2
import numpy as np
import matplotlib.pyplot as plt
import pylab as pl

def showImg(imgName, img, wsize=(400, 400)):
    cv2.namedWindow(imgName, cv2.WINDOW_NORMAL)
    cv2.resizeWindow(imgName, wsize[0], wsize[1])
    cv2.imshow(imgName, img)

def homofilter(I):
    I = np.double(I)
    m, n = I.shape
    rL = 0.5
    rH = 2
    c = 2
    d0 = 20
    I1 = np.log(I + 1)
    FI = np.fft.fft2(I1)
    n1 = np.floor(m / 2)
    n2 = np.floor(n / 2)
    D = np.zeros((m, n))
    H = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            D[i, j] = ((i - n1) ** 2 + (j - n2) ** 2)
            H[i, j] = (rH - rL) * (np.exp(c * (-D[i, j] / (d0 ** 2)))) + rL
    I2 = np.fft.ifft2(H * FI)
    I3 = np.real(np.exp(I2) - 1)
    I4 = I3 - np.min(I3)
    I4 = I4 / np.max(I4) * 255
    dstImg = np.uint8(I4)
    return dstImg

def gaborfilter(srcImg):
    dstImg = np.zeros(srcImg.shape[0:2])
    filters = []
    ksize = [5, 7, 9, 11, 13]
    j = 0
    for K in range(len(ksize)):
        for i in range(12):
            theta = i * np.pi / 12 + np.pi / 24
            gaborkernel = cv2.getGaborKernel((ksize[K], ksize[K]), sigma=2 * np.pi, theta=theta, lambd=np.pi / 2,
                                             gamma=0.5)
            gaborkernel /= 1.5 * gaborkernel.sum()
            filters.append(gaborkernel)
    for kernel in filters:
        gaborImg = cv2.filter2D(srcImg, cv2.CV_8U, kernel)
        np.maximum(dstImg, gaborImg, dstImg)
    return np.uint8(dstImg)

def process(img, filters):
    accum = np.zeros_like(img)
    for kern in filters:
        fimg = cv2.filter2D(img, cv2.CV_8U, kern, borderType=cv2.BORDER_REPLICATE)
        np.maximum(accum, fimg, accum)
    return accum

def process_threaded(img, filters, threadn=8):
    accum = np.zeros_like(img)

    def f(kern):
        return cv2.filter2D(img, cv2.CV_8U, kern)

    pool = ThreadPool(processes=threadn)
    for fimg in pool.imap_unordered(f, filters):
        np.maximum(accum, fimg, accum)
    return accum

###    Gabor特征提取
def getGabor(img, filters):
    res = []  # 滤波结果
    for i in range(len(filters)):
        res1 = process(img, filters[i])
        res.append(np.asarray(res1))

    pl.figure(2)
    for temp in range(len(res)):
        pl.subplot(4, 6, temp + 1)
        pl.imshow(res[temp], cmap='gray')
    pl.show()
    return res  # 返回滤波结果,结果为24幅图，按照gabor角度排列

def build_filters():
    filters = []
    ksize = 31
    for theta in np.arange(0, np.pi, np.pi / 16):
        # kern = cv2.getGaborKernel((ksize, ksize), 4.0, theta, 10.0, 0.5, 0, ktype=cv2.CV_32F)
        kern = cv2.getGaborKernel((ksize, ksize), 2 * np.pi, theta, 17.0, 0.5, 0, ktype=cv2.CV_32F)
        kern /= 1.5 * kern.sum()
        filters.append(kern)
    return filters

def print_gabor(filters):
    for i in range(len(filters)):
        showImg(str(i), filters[i])

def reverse_image(img):
    antiImg = np.zeros_like(img, dtype=np.uint8)
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            antiImg[i][j] = 255 - img[i][j]
    return antiImg

def pass_mask(mask, img):
    # qwe = reverse_image(img)
    qwe = img.copy()
    for i in range(mask.shape[0]):
        for j in range(mask.shape[1]):
            if mask[i][j] == 0:
                qwe[i][j] = 0
    # asd = cv2.filter2D(qwe, cv2.CV_8U, mask)
    return qwe

def showKern(filters):
    for i in list(range(16)):
        kern = filters[i]
        kern = kern - np.min(kern)
        kern = kern / np.max(kern) * 255
        kern = np.clip(kern, 0, 255)
        kern = np.uint8(kern)
        plt.suptitle('Gabor matched filter kernel')
        plt.subplot(4,4,i+1), plt.imshow(kern, 'gray'), plt.axis('off'), plt.title('theta=' + str(i) + r'/pi')
    plt.show()

def calcDice(predict_img, groundtruth_img):
    predict = predict_img.copy()
    groundtruth = groundtruth_img.copy()
    predict[predict < 128] = 0
    predict[predict >= 128] = 1
    groundtruth[groundtruth < 128] = 0
    groundtruth[groundtruth >= 128] = 1
    predict_n = 1 - predict
    groundtruth_n = 1 - groundtruth
    TP = np.sum(predict * groundtruth)
    FP = np.sum(predict * groundtruth_n)
    TN = np.sum(predict_n * groundtruth_n)
    FN = np.sum(predict_n * groundtruth)
    # print(TP, FP, TN, FN)
    dice = 2 * np.sum(predict * groundtruth) / (np.sum(predict) + np.sum(groundtruth))
    return dice

def adjust_gamma(imgs, gamma=1.0):
    # assert (len(imgs.shape)==4)  #4D arrays
    # assert (imgs.shape[1]==1)  #check the channel is 1
    # build a lookup table mapping the pixel values [0, 255] to
    # their adjusted gamma values
    invGamma = 1.0 / gamma
    table = np.array([((i / 255.0) ** invGamma) * 255 for i in np.arange(0, 256)]).astype("uint8")
    # apply gamma correction using the lookup table
    new_imgs = np.zeros_like(imgs)
    for i in range(imgs.shape[0]):
        for j in range(imgs.shape[1]):
            new_imgs[i, j] = cv2.LUT(np.array(imgs[i, j], dtype=np.uint8), table)
    return new_imgs

def build_filters2(sigma=1, YLength=10):
    filters = []
    widthOfTheKernel = np.ceil(np.sqrt((6 * np.ceil(sigma) + 1) ** 2 + YLength ** 2))
    if np.mod(widthOfTheKernel, 2) == 0:
        widthOfTheKernel = widthOfTheKernel + 1
    widthOfTheKernel = int(widthOfTheKernel)
    # print(widthOfTheKernel)
    for theta in np.arange(0, np.pi, np.pi / 16):
        # theta = np.pi/4
        matchFilterKernel = np.zeros((widthOfTheKernel, widthOfTheKernel), dtype=np.float)
        for x in range(widthOfTheKernel):
            for y in range(widthOfTheKernel):
                halfLength = (widthOfTheKernel - 1) / 2
                x_ = (x - halfLength) * np.cos(theta) + (y - halfLength) * np.sin(theta)
                y_ = -(x - halfLength) * np.sin(theta) + (y - halfLength) * np.cos(theta)
                if abs(x_) > 3 * np.ceil(sigma):
                    matchFilterKernel[x][y] = 0
                elif abs(y_) > (YLength - 1) / 2:
                    matchFilterKernel[x][y] = 0
                else:
                    matchFilterKernel[x][y] = -np.exp(-.5 * (x_ / sigma) ** 2) / (np.sqrt(2 * np.pi) * sigma)
        m = 0.0
        for i in range(matchFilterKernel.shape[0]):
            for j in range(matchFilterKernel.shape[1]):
                if matchFilterKernel[i][j] < 0:
                    m = m + 1
        mean = np.sum(matchFilterKernel) / m
        for i in range(matchFilterKernel.shape[0]):
            for j in range(matchFilterKernel.shape[1]):
                if matchFilterKernel[i][j] < 0:
                    matchFilterKernel[i][j] = matchFilterKernel[i][j] - mean
        filters.append(matchFilterKernel)

    return filters

def Z_ScoreNormalization(x, mu, sigma):
    x = (x - mu) / sigma
    return x

def sigmoid(X):
    return 1.0 / (1 + np.exp(-float(X)))

def Normalize(data):
    k = np.zeros(data.shape, np.float)
    # k = np.zeros_like(data)
    # m = np.average(data)
    mx = np.max(data)
    mn = np.min(data)
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            k[i][j] = (float(data[i][j]) - mn) / (mx - mn) * 255
    qwe = np.array(k, np.uint8)
    return qwe

def grayStretch(img, m=60.0/255, e=8.0):
    k = np.zeros(img.shape, np.float)
    ans = np.zeros(img.shape, np.float)
    mx = np.max(img)
    mn = np.min(img)
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            k[i][j] = (float(img[i][j]) - mn) / (mx - mn)
    eps = 0.01
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            ans[i][j] = 1 / (1 + (m / (k[i][j] + eps)) ** e) * 255
    ans = np.array(ans, np.uint8)
    return ans

if __name__ == '__main__':
    path = r'D:\Downloads\DRIVE\DRIVE\test\myDRIVE\\'
    numList = list(range(1, 21))
    # numList = [5]

    for num in numList:
        # 原图
        srcImg = cv2.imread(path + ('%02d' % num) + '_test.tif', cv2.IMREAD_ANYDEPTH | cv2.IMREAD_ANYCOLOR)

        # 标定图
        grountruth = cv2.imread(path + ('%02d' % num) + '_manual1.tif', cv2.IMREAD_GRAYSCALE)
        grayImg = cv2.split(srcImg)[1]

        # 提取掩膜
        ret0, th0 = cv2.threshold(grayImg, 30, 255, cv2.THRESH_BINARY)
        mask = cv2.erode(th0, np.ones((7, 7), np.uint8))
        # showImg("mask", mask)

        # 高斯滤波
        blurImg = cv2.GaussianBlur(grayImg, (5, 5), 0)
        # cv2.imwrite("blurImg.png", blurImg)

        # HE
        heImg = cv2.equalizeHist(blurImg)
        # cv2.imwrite("heImg.png", heImg)

        # CLAHE 光均衡化+对比度增强
        clahe = cv2.createCLAHE(clipLimit=2, tileGridSize=(10, 10))
        claheImg = clahe.apply(blurImg)
        # cv2.imwrite("claheImg.png", claheImg)

        # 同态滤波 光均衡化
        homoImg = homofilter(blurImg)

        preMFImg = adjust_gamma(claheImg, gamma=1.5)
        filters = build_filters2()
        # showKern(filters)
        gaussMFImg = process(preMFImg, filters)
        gaussMFImg_mask = pass_mask(mask, gaussMFImg)
        grayStretchImg = grayStretch(gaussMFImg_mask, m=30.0 / 255, e=8)

        # 二值化
        ret1, th1 = cv2.threshold(grayStretchImg, 30, 255, cv2.THRESH_OTSU)
        predictImg = th1.copy()

        # print(num)
        dice = calcDice(predictImg, grountruth)
        print(num,'',dice)
        wtf = np.hstack([srcImg, cv2.cvtColor(grountruth,cv2.COLOR_GRAY2BGR),cv2.cvtColor(predictImg,cv2.COLOR_GRAY2BGR)])
        cv2.imwrite(('m%02d' % num)+'.png', wtf)
        cv2.waitKey()

结语：这个python版的眼睛血管分割属于小组的共同成果，不过大部分都是基于本作者的matlab版本的代码实现的，代码详见matlab版本的眼睛血管分割

用于测试的图片存放在百度云网：https://pan.baidu.com/s/1y02leAk8d87A7H7sDmHP2g
提取码：b9hc