基于先验知识的自适应阈值DSA图像分割

论文简介

Knowledge-based adaptive thresholding segmentationof digital subtraction angiography images
由华中科技大学提出的基于先验知识的自适应阈值DSA图像分割方法。
论文有三个关键步骤：

1. 将图像分割成小的区域，针对每个区域单独进行大津二值化
2. 对二值化的结果进行评估，判断是否接受
3. 对接受的二值化阈值，进行计算来选取最优阈值

代码

实现了论文中提出的方法二：BusinessMeasure

def BussinessMeasureSegment(img:np.ndarray, a: int = 16, b: int = 32, k: float = 0.2) -> np.ndarray:
    '''Business Measure

    args:
        img:输入图像
        a: El的大小为(a,a),要求d<a<2d,d为血管最大宽度
        b: bkl的大小为(b,b),要求a<b<3a
        k: 计算最优阈值时的参数,按论文中的经验值取0.2
    returns:
        output:分割结果
    '''
    
    output = np.zeros_like(img)
    # 论文中的El
    E = np.zeros([b, b, a, a], dtype=np.uint8)
    for j in range(img.shape[0]//a):
        for i in range(img.shape[1]//a):
            E[i, j] = img[i*a:i*a+a, j*a:j*a+a]

    for i in range(img.shape[0]//a):
        for j in range(img.shape[1]//a):
            bkl = np.zeros([4, b, b], dtype=np.uint8)  # 论文中的bkl
            fkl = np.zeros([4, b, b], dtype=np.uint8)  # 论文中的fkl,为bkl的二值化结果
            Tkl = np.zeros([4])  # 论文中的Tkl,为bkl的二值化阈值
            tkl = np.zeros([4, b, b], np.uint8)  # 论文中的tkl

            # 按边界条件计算E[i,j]的8邻域
            El = np.zeros([9, a, a], dtype=np.uint8)
            El[4] = E[i, j].copy()
            if(i > 0 and j > 0):
                El[0] = E[i-1, j-1].copy()
            if(i > 0):
                El[1] = E[i-1, j].copy()
            if(i > 0 and j < 31):
                El[2] = E[i-1, j+1].copy()
            if(j > 0):
                El[3] = E[i, j-1].copy()
            if(j < 31):
                El[5] = E[i, j+1].copy()
            if(i < 31 and j > 0):
                El[6] = E[i+1, j-1].copy()
            if(i < 31):
                El[7] = E[i+1, j].copy()
            if(i < 31 and j < 31):
                El[8] = E[i+1, j+1].copy()

            # 计算bkl
            bkl[0] = np.vstack([np.hstack([El[0], El[1]]),
                              np.hstack([El[3], El[4]])])
            bkl[1] = np.vstack([np.hstack([El[1], El[2]]),
                              np.hstack([El[4], El[5]])])
            bkl[2] = np.vstack([np.hstack([El[3], El[4]]),
                              np.hstack([El[6], El[7]])])
            bkl[3] = np.vstack([np.hstack([El[4], El[5]]),
                              np.hstack([El[7], El[8]])])

            # 计算fkl和Tkl
            for n in range(4):
                Tkl[n], fkl[n] = cv2.threshold(bkl[n], 0, 255, cv2.THRESH_OTSU)
            fkl = 255 - fkl
            # 将fkl改为float型,避免uint8计算0-1时溢出
            fkl = fkl.astype(float)
            fkl[fkl != 0] = 1

            # 计算tkl,为二值图像中相邻像素值变化的频率
            # if |fkl(i,j)-fkl(i,j+1) = 1| or |fkl(i,j)-fkl(i+1,j)| = 1, tkl(i,j) = 1
            # if |fkl(i,j)-fkl(i,j+1) = 0| and |fkl(i,j)-fkl(i+1,j)| = 0, tkl(i,j) = 0
            # 论文中没有提到边界的处理策略,此处默认边界为0,故只计算为1的情况
            for n in range(4):
                for x in range(b-1):
                    for y in range(b-1):
                        if(fkl[n, x, y]-fkl[n, x, y+1]) or (fkl[n, x, y]-fkl[n, x+1, y]):
                            tkl[n, x, y] = 1

            # 论文中的business measure,频率越高说明是背景的可能性越高
            busy = np.sum(tkl, axis=(1, 2))
            # 论文中的Hσ,通过Hσ和busy来确认二值化的结果是否被接受,k为经验值
            # busy(bkl) <= Hσ时,bkl二值化结果被接受,认为bkl含有分割目标,否则认为是背景
            H = k*busy.max() + (1-k)*busy.min()
            num = np.sum(busy <= H)
            if num == 0:
                output[i*a:i*a+a, j*a:j*a+a] = 0
            elif num == 1:
                idx = np.where(busy <= H)[0][0]
                _, binary = cv2.threshold(
                    E[i, j].copy(), Tkl[idx], 255, cv2.THRESH_BINARY)
                output[i*a:i*a+a, j*a:j*a+a] = 255 - binary
            else:
                # 计算最优阈值
                threshold = 0
                Vkl_max = 0
                # 被接受的bkl索引
                idx_list = np.where(busy <= H)[0]

                for idx in idx_list:
                    # 计算Ekl: E[i,j]按Tkl的阈值分割的结果
                    Ekl = E[i, j].copy()
                    Ekl[Ekl < Tkl[idx]] = 0
                    Ekl[Ekl >= Tkl[idx]] = 1
                    # 计算Vkl: 选择使Vkl最大的Tkl作为最优阈值
                    Vkl = np.sum(np.power(Ekl-np.sum(Ekl)/(a*a), 2)/(a*a))
                    if Vkl >= Vkl_max:
                        threshold = Tkl[idx]
                _, binary = cv2.threshold(
                    E[i, j].copy(), threshold, 255, cv2.THRESH_BINARY)
                output[i*a:i*a+a, j*a:j*a+a] = 255 - binary
    return output

分割结果

血管部分的分割效果还是很好的，但是图像背景区域的噪声还是无法有效的分割，没有论文里的效果那么好。