基于3DFrangi滤波的血管强化方法（附代码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滤波，最后放上几张结果图进行对比。

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

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上面的内容)

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

原图：

结果：

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

原图：

结果：

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

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