基于3D Frangi滤波的血管强化方法（附代码python）

---

前言

Frangi滤波原文：https://www.researchgate.net/publication/2388170_Multiscale_Vessel_Enhancement_Filtering

Frangi滤波翻译讲解：
https://zhuanlan.zhihu.com/p/127951058

参考代码：https://github.com/vinhnguyen21/python_JermanEnhancementFilter

Frangi滤波原文中详细说明了3D，2D图像下的血管强化方法，但是在网上找了好久，只有找到2D滤波的代码，在做毕设的时候因为时间有限，所以对三正交平面都进行一次2D Frangi滤波的方式代替3D Frangi滤波，虽然也有效果但总是不是很舒服。

本文首先会根据参考代码中的2D Frangi滤波进行讲解，接着在此基础上按照原文的意思更改3D Frangi滤波，最后放上几张结果图进行对比。

本人水平有限，还望各位大佬批评指正。

---

一、2D Frangi滤波——原文复现

1、import

需要说明的是，我们的3D文件是.nii文件，这里使用SimpleITK进行读写

import cv2
import os
import numpy as np
from scipy import ndimage
import SimpleITK as sitk 

2、vesselness2d

vesselness2d.py(记得import上面的内容)

class vesselness2d:
    def __init__(self, image, sigma, spacing, tau):
        super(vesselness2d, self).__init__()
        # image 为numpy类型，表示n * m 的二维矩阵。
        self.image = image
        # sigma 为list 类型，表示高斯核的尺度。
        self.sigma = sigma
        # spacing 为list类型，表示.nii文件下某一切面下的体素的二维尺寸。如果输入图像本身为二维图像，则为[1,1],如果为三维图像，则为对应的space。
        self.spacing = spacing
        # tau 为float类型，表示比例系数。
        self.tau = tau
        # 图像尺寸
        self.size = image.shape
	
	# 使用特定的特定sigma尺寸下的高斯核对图像滤波
	# 这里作者并没有使用n*n的卷积核，而是分别使用n*1，1*n的卷积对图像进行x和y方向上的卷积，
	# 并且使用的是最原始的计算高斯函数得到卷积核，而不是直接用现成的高斯卷积核，
	# 通过证明可以发现在两方面的结果是等价的。
    def gaussian(self, image, sigma):
        siz = sigma * 6 # 核的尺寸
        
        # x轴方向上的滤波
        temp = round(siz / self.spacing[0] / 2)
        x = [i for i in range(-temp, temp + 1)]
        x = np.array(x)
        H = np.exp(-(x ** 2 / (2 * ((sigma / self.spacing[0]) ** 2))))
        H = H / np.sum(H)
        Hx = H.reshape(len(H), 1, 1)
        I = ndimage.filters.convolve(image, Hx, mode='nearest')
		
		# y轴方向上的滤波
        temp = round(siz / self.spacing[1] / 2)
        x = [i for i in range(-temp, temp + 1)]
        x = np.array(x)
        H = np.exp(-(x ** 2 / (2 * ((sigma / self.spacing[1]) ** 2))))
        H = H / np.sum(H[:])
        Hy = H.reshape(1, len(H), 1)
        I = ndimage.filters.convolve(I, Hy, mode='nearest')
        return I

	# 求图像的梯度
    def gradient2(self, F, option):
        k = self.size[0]
        l = self.size[1]
        D = np.zeros(F.shape)
        if option == "x":
            D[0, :] = F[1, :] - F[0, :]
            D[k - 1, :] = F[k - 1, :] - F[k - 2, :]
            # take center differences on interior points
            D[1:k - 2, :] = (F[2:k - 1, :] - F[0:k - 3, :]) / 2
        else:
            D[:, 0] = F[:, 1] - F[:, 0]
            D[:, l - 1] = F[:, l - 1] - F[:, l - 2]
            D[:, 1:l - 2] = (F[:, 2:l - 1] - F[:, 0:l - 3]) / 2
        return D
        
	# 求海森矩阵中所需要的二阶偏导数
    def Hessian2d(self, image, sigma):
        image = self.gaussian(image, sigma)
        # image = ndimage.gaussian_filter(image, sigma, mode = 'nearest')
        Dy = self.gradient2(image, "y")
        Dyy = self.gradient2(Dy, "y")

        Dx = self.gradient2(image, "x")
        Dxx = self.gradient2(Dx, "x")
        Dxy = self.gradient2(Dx, 'y')
        return Dxx, Dyy, Dxy
	
	# 求解海森矩阵的两个特征值
	# 这里作者使用求根公式，将二阶海森矩阵展开，a=1,b=-(Ixx+Iyy),c=(Ixx*Iyy-Ixy*Ixy) 
	# 首先计算 sqrt(b^2 - 4ac),通过化简得到tmp
	# 最后得到两个特征值mu1，mu2，根据大小关系，大的为mu2，小的为mu1
    def eigvalOfhessian2d(self, Dxx, Dyy, Dxy):
        tmp = np.sqrt((Dxx - Dyy) ** 2 + 4 * (Dxy ** 2))
        # compute eigenvectors of J, v1 and v2
        mu1 = 0.5 * (Dxx + Dyy + tmp)
        mu2 = 0.5 * (Dxx + Dyy - tmp)
        # Sort eigen values by absolute value abs(Lambda1) < abs(Lambda2)
        indices = (np.absolute(mu1) > np.absolute(mu2))
        Lambda1 = mu1
        Lambda1[indices] = mu2[indices]

        Lambda2 = mu2
        Lambda2[indices] = mu1[indices]
        return Lambda1, Lambda2

    def imageEigenvalues(self, I, sigma):
        hxx, hyy, hxy = self.Hessian2d(I, sigma)
        c = sigma ** 2
        hxx = -c * hxx
        hyy = -c * hyy
        hxy = -c * hxy

        # 为了降低运算量，去掉噪声项的计算
        B1 = -(hxx + hyy)
        B2 = hxx * hyy - hxy ** 2
        T = np.ones(B1.shape)
        T[(B1 < 0)] = 0
        T[(B1 == 0) & (B2 == 0)] = 0
        T = T.flatten()
        indeces = np.where(T == 1)[0]

        hxx = hxx.flatten()
        hyy = hyy.flatten()
        hxy = hxy.flatten()

        hxx = hxx[indeces]
        hyy = hyy[indeces]
        hxy = hxy[indeces]


        lambda1i, lambda2i = self.eigvalOfhessian2d(hxx, hyy, hxy)
        lambda1 = np.zeros(self.size[0] * self.size[1], )
        lambda2 = np.zeros(self.size[0] * self.size[1], )

        lambda1[indeces] = lambda1i
        lambda2[indeces] = lambda2i

        # 去掉噪声
        lambda1[(np.isinf(lambda1))] = 0
        lambda2[(np.isinf(lambda2))] = 0

        lambda1[(np.absolute(lambda1) < 1e-4)] = 0
        lambda1 = lambda1.reshape(self.size)

        lambda2[(np.absolute(lambda2) < 1e-4)] = 0
        lambda2 = lambda2.reshape(self.size)
        return lambda1, lambda2
	
	# 血管强化
    def vesselness2d(self):
        for j in range(len(self.sigma)):
            lambda1, lambda2 = self.imageEigenvalues(self.image, self.sigma[j])
            lambda3 = lambda2.copy()
            new_tau = self.tau * np.min(lambda3)
            lambda3[(lambda3 < 0) & (lambda3 >= new_tau)] = new_tau
            different = lambda3 - lambda2
            response = ((np.absolute(lambda2) ** 2) * np.absolute(different)) * 27 / (
                        (2 * np.absolute(lambda2) + np.absolute(different)) ** 3)
            response[(lambda2 < lambda3 / 2)] = 1
            response[(lambda2 >= 0)] = 0

            response[np.where(np.isinf(response))[0]] = 0
            if j == 0:
                vesselness = response
            else:
                vesselness = np.maximum(vesselness, response)
        vesselness[(vesselness < 1e-2)] = 0
        return vesselness

3、应用示例（原文）

demo.py

需要说明的是，这里使用的图像是0-255灰度图像，原文的强化针对的是背景亮，血管暗的图像，但是这里的图像是相反，所以在下面对图像进行了像素灰度值的反转。

from PIL import Image
import numpy as np
import cv2
import matplotlib.pyplot as plt
from vesselness2d import *

img_dir = 'images/test.tif' #路径写自己的

#reading image
image = Image.open(img_dir).convert("RGB")
image = np.array(image)
plt.figure(figsize=(10,10))
plt.imshow(image, cmap='gray')

#convert forgeground to background and vice-versa
image = 255-image

image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
thr = np.percentile(image[(image > 0)], 1)*0.9
image[(image <= thr)] = thr
image = image - np.min(image)
image = image / np.max(image)

sigma=[0.5,1, 1.5, 2, 2.5]
spacing = [1, 1]
tau = 2

output = vesselness2d(image, sigma, spacing, tau)
output = output.vesselness2d()

plt.figure(figsize=(10,10))
plt.imshow(output, cmap='gray')

原图：

结果：

二、3D Frangi滤波 ——三正交平面分别进行2D Frangi滤波

1、import

import cv2 as cv
import SimpleITK as sitk
from vesselness2d import *

2、main

Hessian_3D.py

# 这里使用的是MSD数据集中的肝脏血管分割数据集，并且只用已训练好的肝脏分割模型对其进行分割，
# 只保留肝脏区域，图像的灰度范围是[0,200]，血管相较于背景为白色
def edge(img,position):
    img_dt = np.zeros((len(img),len(img[0]),len(img[0][0])))
    img_dt[:] = img[:]
    origin = img_dt[0][0][0]
    img_dt[img_dt!=-origin] = 1
    img_dt[img_dt==-origin] = 0

    tmp = np.ones((len(img_dt), len(img_dt[0]), len(img_dt[0][0])))
    if position == "x":
        for i in range(len(img_dt)):
            kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))
            dst = cv.erode(img_dt[i], kernel)
            tmp[i] = dst
        img_dt[tmp == 1] = 0
    elif position == "y":
        for i in range(len(img_dt[0])):
            kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))
            dst = cv.erode(img_dt[:,i,:], kernel)
            tmp[:,i,:] = dst
        img_dt[tmp == 1] = 0
    elif position == "z":
        for i in range(len(img_dt[0][0])):
            kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))
            dst = cv.erode(img_dt[:,:,i], kernel)
            tmp[:,:,i] = dst
        img_dt[tmp == 1] = 0
    return img_dt

def frangi(img, sigma, spacing, tau,position):
    img_dt = np.zeros((len(img), len(img[0]), len(img[0][0])))
    img_dt[:] = img[:]
    result_dt = np.zeros((len(img_dt), len(img_dt[0]), len(img_dt[0][0])))
    if position == "x":
        for i in range(len(img_dt)):
            image = img_dt[i]
            output = vesselness2d(image, sigma, spacing, tau)
            output = output.vesselness2d()
            result_dt[i] = output
    elif position == "y":
        for i in range(len(img_dt[0])):
            image = img_dt[:,i,:]
            output = vesselness2d(image, sigma, spacing, tau)
            output = output.vesselness2d()
            result_dt[:,i,:] = output
    elif position == "z":
        for i in range(len(img_dt[0][0])):
            image = img_dt[:,:,i]
            output = vesselness2d(image, sigma, spacing, tau)
            output = output.vesselness2d()
            result_dt[:,:,i] = output
    return result_dt

def Hessian3D(image,sigma, tau):
    img_dt = sitk.GetArrayFromImage(image)
    stand = img_dt[0][0][0]
    img_dt[img_dt==stand] = -200
    img_dt = 200-img_dt
    img_dt[img_dt==400] = -200

    edge_x = edge(img_dt,"x")
    edge_y = edge(img_dt,"y")
    edge_z = edge(img_dt,"z")
    edge_x[edge_y == 1] = 1
    edge_x[edge_z == 1] = 1

    space = image.GetSpacing()
    spacing_x = [space[0],space[1]]
    spacing_y = [space[0],space[2]]
    spacing_z = [space[1],space[2]]

    hessian_x = frangi(img_dt, sigma, spacing_x, tau, "x")
    return

    hessian_y = frangi(img_dt, sigma, spacing_y, tau, "y")
    hessian_z = frangi(img_dt, sigma, spacing_z, tau, "z")

    result_dt = hessian_x+hessian_y+hessian_z
    result_dt[-1] = np.zeros((len(result_dt[0]), len(result_dt[0][0])))
    result_dt[edge == 1] = 0
    result_dt *= 400
    result_dt[result_dt > 200] = 200
    result_dt[img_dt == -200] = -200
    result_dt = result_dt.astype(int)

    result = sitk.GetImageFromArray(result_dt)
    result.SetSpacing(image.GetSpacing())
    result.SetOrigin(image.GetOrigin())
    result.SetDirection(image.GetDirection())
    return result


# 这里的main函数根据自己的需要改
# 这里我的直接对整个文件夹中的全部.nii文件进行处理
if __name__ == "__main__":
    sigma = [0.5, 1, 1.5, 2, 2.5]
    tau = 2
    path = "D:\\PythonProject\\Daily\\AHE"
    result_path = "F:\\3DUNet-Pytorch-master_vesselSeg\\raw_dataset\\train_seg\\hessian"
    path_list = os.listdir(path)
    for i in path_list:
        image_i_path = os.path.join(path,i)
        img = sitk.ReadImage(image_i_path)
        result = Hessian3D(img,sigma,tau)
        sitk.WriteImage(result,os.path.join(result_path,i))
        print(i + " is OK!")

原图：

结果：

---

三、3D Frangi滤波 ——原文复现

1、import

import cv2 as cv
import SimpleITK as sitk
from vesselness2d import *

2、vesselness3d

Hessian_3D.py

# 对于3D Frangi滤波，与2D Frangi不同点在于
# 1、高斯滤波考虑第三个维度
# 2、构造三阶海森矩阵，[[Ixx,Ixy,Ixz],[Ixy,Iyy,Iyz],[Ixz,Iyz,Izz]]
# 3、求解三阶海森矩阵的特征值lambda1，lambda2，lambda3，并按照绝对值的大小排序
# 4、为减小求解时间，对于Ixx+Iyy+Izz<0的情况直接将灰度置为0
# 5、使用三维滤波公式求解体素灰度值
class vesselness3d:
    def __init__(self, image, sigma, spacing):
        super(vesselness3d, self).__init__()
        self.image = image
        self.sigma = sigma
        self.spacing = spacing
        self.size = image.shape

    def gaussian(self, image, sigma):
        siz = sigma * 6
        temp = round(siz / self.spacing[0] / 2)
        # processing x-axis
        x = [i for i in range(-temp, temp + 1)]
        x = np.array(x)
        H = np.exp(-(x ** 2 / (2 * ((sigma / self.spacing[0]) ** 2))))
        H = H / np.sum(H)
        Hx = H.reshape(len(H), 1, 1)
        I = ndimage.filters.convolve(image, Hx, mode='nearest')

        # processing y-axis
        temp = round(siz / self.spacing[1] / 2)
        x = [i for i in range(-temp, temp + 1)]
        x = np.array(x)
        H = np.exp(-(x ** 2 / (2 * ((sigma / self.spacing[1]) ** 2))))
        H = H / np.sum(H[:])
        Hy = H.reshape(1, len(H), 1)
        I = ndimage.filters.convolve(I, Hy, mode='nearest')

        # processing z-axis
        temp = round(siz / self.spacing[2] / 2)
        x = [i for i in range(-temp, temp + 1)]
        x = np.array(x)
        H = np.exp(-(x ** 2 / (2 * ((sigma / self.spacing[2]) ** 2))))
        H = H / np.sum(H[:])
        Hz = H.reshape(1, 1, len(H))
        I = ndimage.filters.convolve(I, Hz, mode='nearest')
        return I

    def gradient2(self, F, option):
        k = self.size[0]
        l = self.size[1]
        h = self.size[2]
        D = np.zeros(F.shape)
        if option == "x":
            D[0, :, :] = F[1, :, :] - F[0, :, :]
            D[k - 1, :, :] = F[k - 1, :, :] - F[k - 2, :, :]
            # take center differences on interior points
            D[1:k - 2, :, :] = (F[2:k - 1, :, :] - F[0:k - 3, :, :]) / 2
        elif option == "y":
            D[:, 0, :] = F[:, 1, :] - F[:, 0, :]
            D[:, l - 1, :] = F[:, l - 1, :] - F[:, l - 2, :]
            D[:, 1:l - 2, :] = (F[:, 2:l - 1, :] - F[:, 0:l - 3, :]) / 2
        elif option == "z":
            D[:, :, 0] = F[:, :,  1] - F[:, :, 0]
            D[:, :, h- 1] = F[:, :, h - 1] - F[:, :, h - 2]
            D[:, :, 1:h - 2] = (F[:, :, 2:h - 1] - F[:, :, 0:h - 3]) / 2
        return D

    def Hessian2d(self, image, sigma):
        image = self.gaussian(image, sigma)
        self.gaus_image = image
        Dz = self.gradient2(image, "z")
        Dzz = self.gradient2(Dz, "z")

        Dy = self.gradient2(image, "y")
        Dyy = self.gradient2(Dy, "y")
        Dyz = self.gradient2(Dy, "z")

        Dx = self.gradient2(image, "x")
        Dxx = self.gradient2(Dx, "x")
        Dxy = self.gradient2(Dx, 'y')
        Dxz = self.gradient2(Dx, "z")
        return Dxx, Dyy, Dzz, Dxy, Dxz, Dyz

    def eigvalOfhessian2d(self, array):
        tmp = np.linalg.eig(array)
        lamda = sorted([(abs(tmp[0][0]),tmp[0][0]),(abs(tmp[0][1]),tmp[0][1]),(abs(tmp[0][2]),tmp[0][2])])
        Lambda1 = lamda[0][1]
        Lambda2 = lamda[1][1]
        Lambda3 = lamda[2][1]
        return Lambda1, Lambda2, Lambda3

    def imageEigenvalues(self, I, sigma):
        self.hxx, self.hyy,self.hzz, self.hxy, self.hxz, self.hyz= self.Hessian2d(I, sigma)
        hxx = self.hxx
        hyy = self.hyy
        hzz = self.hzz
        hxy = self.hxy
        hxz = self.hxz
        hyz = self.hyz

        hxx = hxx.flatten()
        hyy = hyy.flatten()
        hzz = hzz.flatten()
        hxy = hxy.flatten()
        hxz = hxz.flatten()
        hyz = hyz.flatten()

        Lambda1_list = []
        Lambda2_list = []
        Lambda3_list = []
        count = 0
        for i in range(len(hxx)):
            if hxx[i] + hyy[i] + hzz[i] <= 0:
                Lambda1, Lambda2, Lambda3 = 0, 0 ,0
            else:
                array = np.array([[hxx[i],hxy[i],hxz[i]],[hxy[i],hyy[i],hyz[i]],[hxz[i],hyz[i],hzz[i]]])
                Lambda1, Lambda2, Lambda3 = self.eigvalOfhessian2d(array)
            if Lambda1 !=0 and Lambda2!=0 and Lambda3!=0 :
                count+=1
            
            Lambda1_list.append(Lambda1)
            Lambda2_list.append(Lambda2)
            Lambda3_list.append(Lambda3)

        Lambda1_list = np.array(Lambda1_list)
        Lambda2_list = np.array(Lambda2_list)
        Lambda3_list = np.array(Lambda3_list)
        Lambda1_list[(np.isinf(Lambda1_list))] = 0
        Lambda2_list[(np.isinf(Lambda2_list))] = 0
        Lambda3_list[(np.isinf(Lambda3_list))] = 0

        # Lambda1_list[(np.absolute(Lambda1_list) < 1e-4)] = 0
        Lambda1_list = Lambda1_list.reshape(self.size)
        # Lambda2_list[(np.absolute(Lambda2_list) < 1e-4)] = 0
        Lambda2_list = Lambda2_list.reshape(self.size)
        # Lambda3_list[(np.absolute(Lambda3_list) < 1e-4)] = 0
        Lambda3_list = Lambda3_list.reshape(self.size)
        return Lambda1_list,Lambda2_list,Lambda3_list

    def vesselness3d(self):
        for k in range(len(self.sigma)):
            lambda1, lambda2, lambda3 = self.imageEigenvalues(self.image, self.sigma[k])
            c = self.gaus_image.max()/2

            item1 = (1 - np.exp(-2 * (lambda2 ** 2) / (lambda3 ** 2)))
            item2 = np.exp(-2 * (lambda1 ** 2) / np.absolute(lambda2 * lambda3))
            item3 = (1 - np.exp(-(lambda1 ** 2 + lambda2 ** 2 + lambda3 ** 2) / (2 * c ** 2)))
            
            item1[lambda3==0] = 0
            item2[lambda3==0] = 0
            
            response = item1*item2*item3

            response[np.where(np.isnan(response))] = 0

            if k == 0:
                vesselness = response
            else:
                vesselness = np.maximum(vesselness, response)

        vesselness = (vesselness / (vesselness.max()))*20000
        vesselness[vesselness>200] = 200
        return vesselness

if __name__ == "__main__":
    sigma = [0.5, 1, 1.5, 2, 2.5]
    path = "raw_hepaticvessel_001.nii.gz"
    img = sitk.ReadImage(path)
    img_data = sitk.GetArrayFromImage(img)
    space = img.GetSpacing()
    direction = img.GetDirection()
    origin = img.GetOrigin()
    img_data = 200-img_data
    v = vesselness3d(img_data,sigma,list(space))
    image_data = v.vesselness3d()
    img = sitk.GetImageFromArray(image_data)
    img.SetOrigin(origin)
    img.SetDirection(direction)
    img.SetSpacing(space)
    sitk.WriteImage(img, "Frangi_hepaticvessel_001.nii.gz")

总结

1、Frangi滤波作为经典血管强化、管状强化滤波算法，具有极好的数学证明与实验结果。
2、尽管Frangi滤波效果很好，但仍需要进行调参，如sigma的选取，2D滤波中出现过的tau，已经3D复现中使用的多个超参数。
3、从本文的结论中明显看出，实际上使用三正交平面的滤波效果优于3D Frangi，原因是2中提到的，关于超参的选取问题，而且从2D的复现中，我们能明显看出，代码原作者对Frangi原文做出了极大的改变，使其效果更优。
4、边缘会比血管更容易被增强，所以要处理边缘（在三正交平面我处理了，在3D Frangi没有处理）
5、Frangi滤波的效果告诉我们，机器学习如此发达的今天，特征工程仍必不可少。