前言
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")
文章来源：https://www.toymoban.com/news/detail-403885.html
总结
1、Frangi滤波作为经典血管强化、管状强化滤波算法，具有极好的数学证明与实验结果。
2、尽管Frangi滤波效果很好，但仍需要进行调参，如sigma的选取，2D滤波中出现过的tau，已经3D复现中使用的多个超参数。
3、从本文的结论中明显看出，实际上使用三正交平面的滤波效果优于3D Frangi，原因是2中提到的，关于超参的选取问题，而且从2D的复现中，我们能明显看出，代码原作者对Frangi原文做出了极大的改变，使其效果更优。
4、边缘会比血管更容易被增强，所以要处理边缘（在三正交平面我处理了，在3D Frangi没有处理）
5、Frangi滤波的效果告诉我们，机器学习如此发达的今天，特征工程仍必不可少。文章来源地址https://www.toymoban.com/news/detail-403885.html
到了这里，关于基于3D Frangi滤波的血管强化方法（附代码python）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
