基于3D Frangi滤波的血管强化方法（附代码python）
Posted HitStuHan
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了基于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(im以上是关于基于3D Frangi滤波的血管强化方法（附代码python）的主要内容，如果未能解决你的问题，请参考以下文章
图像增强基于matlab Frangi滤波器血管图像增强含Matlab源码 2108期
深度强化学习 Policy Gradients 模型解析，附Pytorch完整代码