文章目录

前言
一、2D Frangi滤波——原文复现
- 1、import
- 2、vesselness2d
- 3、应用示例（原文）
二、3D Frangi滤波 ——三正交平面分别进行2D Frangi滤波
- 1、import
- 2、main
三、3D Frangi滤波 ——原文复现
- 1、import
- 2、vesselness3d
总结

前言

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 pointsD[1:k - 2, :] = (F[2:k - 1, :] - F[0:k - 3, :]) / 2else: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]) / 2return 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，小的为mu1def eigvalOfhessian2d(self, Dxx, Dyy, Dxy):tmp = np.sqrt((Dxx - Dyy) ** 2 + 4 * (Dxy ** 2))# compute eigenvectors of J, v1 and v2mu1 = 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 = mu1Lambda1[indices] = mu2[indices]Lambda2 = mu2Lambda2[indices] = mu1[indices]return Lambda1, Lambda2def imageEigenvalues(self, I, sigma):hxx, hyy, hxy = self.Hessian2d(I, sigma)c = sigma ** 2hxx = -c * hxxhyy = -c * hyyhxy = -c * hxy# 为了降低运算量，去掉噪声项的计算B1 = -(hxx + hyy)B2 = hxx * hyy - hxy ** 2T = np.ones(B1.shape)T[(B1 < 0)] = 0T[(B1 == 0) & (B2 == 0)] = 0T = 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] = lambda1ilambda2[indeces] = lambda2i# 去掉噪声lambda1[(np.isinf(lambda1))] = 0lambda2[(np.isinf(lambda2))] = 0lambda1[(np.absolute(lambda1) < 1e-4)] = 0lambda1 = lambda1.reshape(self.size)lambda2[(np.absolute(lambda2) < 1e-4)] = 0lambda2 = 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_taudifferent = lambda3 - lambda2response = ((np.absolute(lambda2) ** 2) * np.absolute(different)) * 27 / ((2 * np.absolute(lambda2) + np.absolute(different)) ** 3)response[(lambda2 < lambda3 / 2)] = 1response[(lambda2 >= 0)] = 0response[np.where(np.isinf(response))[0]] = 0if j == 0:vesselness = responseelse:vesselness = np.maximum(vesselness, response)vesselness[(vesselness < 1e-2)] = 0return 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-imageimage = 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 = 2output = 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] = 1img_dt[img_dt==-origin] = 0tmp = 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] = dstimg_dt[tmp == 1] = 0elif 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,:] = dstimg_dt[tmp == 1] = 0elif 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] = dstimg_dt[tmp == 1] = 0return img_dtdef 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] = outputelif 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,:] = outputelif 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] = outputreturn result_dtdef Hessian3D(image,sigma, tau):img_dt = sitk.GetArrayFromImage(image)stand = img_dt[0][0][0]img_dt[img_dt==stand] = -200img_dt = 200-img_dtimg_dt[img_dt==400] = -200edge_x = edge(img_dt,"x")edge_y = edge(img_dt,"y")edge_z = edge(img_dt,"z")edge_x[edge_y == 1] = 1edge_x[edge_z == 1] = 1space = 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")returnhessian_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_zresult_dt[-1] = np.zeros((len(result_dt[0]), len(result_dt[0][0])))result_dt[edge == 1] = 0result_dt *= 400result_dt[result_dt > 200] = 200result_dt[img_dt == -200] = -200result_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 = 2path = "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 = imageself.sigma = sigmaself.spacing = spacingself.size = image.shapedef gaussian(self, image, sigma):siz = sigma * 6temp = round(siz / self.spacing[0] / 2)# processing x-axisx = [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-axistemp = 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-axistemp = 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 Idef 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 pointsD[1:k - 2, :, :] = (F[2:k - 1, :, :] - F[0:k - 3, :, :]) / 2elif 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, :]) / 2elif 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]) / 2return Ddef Hessian2d(self, image, sigma):image = self.gaussian(image, sigma)self.gaus_image = imageDz = 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, Dyzdef 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, Lambda3def imageEigenvalues(self, I, sigma):self.hxx, self.hyy,self.hzz, self.hxy, self.hxz, self.hyz= self.Hessian2d(I, sigma)hxx = self.hxxhyy = self.hyyhzz = self.hzzhxy = self.hxyhxz = self.hxzhyz = self.hyzhxx = hxx.flatten()hyy = hyy.flatten()hzz = hzz.flatten()hxy = hxy.flatten()hxz = hxz.flatten()hyz = hyz.flatten()Lambda1_list = []Lambda2_list = []Lambda3_list = []count = 0for i in range(len(hxx)):if hxx[i] + hyy[i] + hzz[i] <= 0:Lambda1, Lambda2, Lambda3 = 0, 0 ,0else: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+=1Lambda1_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))] = 0Lambda2_list[(np.isinf(Lambda2_list))] = 0Lambda3_list[(np.isinf(Lambda3_list))] = 0# Lambda1_list[(np.absolute(Lambda1_list) < 1e-4)] = 0Lambda1_list = Lambda1_list.reshape(self.size)# Lambda2_list[(np.absolute(Lambda2_list) < 1e-4)] = 0Lambda2_list = Lambda2_list.reshape(self.size)# Lambda3_list[(np.absolute(Lambda3_list) < 1e-4)] = 0Lambda3_list = Lambda3_list.reshape(self.size)return Lambda1_list,Lambda2_list,Lambda3_listdef vesselness3d(self):for k in range(len(self.sigma)):lambda1, lambda2, lambda3 = self.imageEigenvalues(self.image, self.sigma[k])c = self.gaus_image.max()/2item1 = (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] = 0item2[lambda3==0] = 0response = item1*item2*item3response[np.where(np.isnan(response))] = 0if k == 0:vesselness = responseelse:vesselness = np.maximum(vesselness, response)vesselness = (vesselness / (vesselness.max()))*20000vesselness[vesselness>200] = 200return vesselnessif __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_datav = 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滤波的效果告诉我们，机器学习如此发达的今天，特征工程仍必不可少。

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

深度学习篇之数据集划分方法-附代码python详细注释
在深度学习训练模型过程中,我们第一步就是要收集相应的数据集,之后我们就是要将数据划分为训练集train和验证集val,但是有时间我们时常面临数据量庞大的问题,手动划分显然是不现实的,因为太麻烦了,而且 ...
2020华为杯E题——基于灰色预测的大雾能见度预测模型（附代码）
2020华为杯E题--基于灰色预测的大雾能见度预测模型(附代码) 一.E题赛题二.赛题分析三.基于灰色预测的大雾能见度预测模型 3.1灰色预测模型的定义 3.2 灰色预测模型的目的 GM(1,1) ...
路径规划算法：基于灰狼优化(GWO)的路径规划算法- 附代码
路径规划算法:基于灰狼优化(GWO)的路径规划算法- 附代码文章目录路径规划算法:基于灰狼优化(GWO)的路径规划算法- 附代码 1.算法原理 1.1 环境设定 1.2 约束条件 1.3 适应度函 ...
FPGA零基础学习：基于FPGA的二进制转BCD设计（附代码）
FPGA零基础学习:基于FPGA的二进制转BCD设计(附代码) 本系列将带来FPGA的系统性学习,从最基本的数字电路基础开始,最详细操作步骤,最直白的言语描述,手把手的"傻瓜式"讲 ...
路径规划算法：基于入侵杂草优化的路径规划算法- 附代码
路径规划算法:基于入侵杂草优化的路径规划算法- 附代码文章目录路径规划算法:基于入侵杂草优化的路径规划算法- 附代码 1.算法原理 1.1 环境设定 1.2 约束条件 1.3 适应度函数 2.算法 ...
基于遗传算法优化的Elman神经网络数据预测-附代码
基于遗传算法优化的Elman神经网络数据预测 - 附代码文章目录基于遗传算法优化的Elman神经网络数据预测 - 附代码 1.Elman 神经网络结构 2.Elman 神经用络学习过程 3.电力负 ...
java实现生日提醒_asp实现的可以提醒生日的几种方法附代码
asp实现的可以提醒生日的几种方法附代码更新时间:2008年06月20日 08:44:52 作者: asp提醒生日的原理一般就是把用户的出生日期和当前日期对比,如果少于几天则提醒显示,下面是几种 ...
轨迹系列2——一种基于中值滤波的轨迹纠偏方法和几点思考
文章版权由作者李晓晖和博客园共有,若转载请于明显处标明出处:http://www.cnblogs.com/naaoveGIS/ 1.背景在无路网的情况下,如何进行轨迹纠偏也是一个很多人在研究的内容, ...
图像处理/计算机视觉/ python环境下如何用滤波器（/逆滤波/均值滤波/低通滤波/高通滤波）处理图片【附代码】
计算机视觉滤波器实操基础知识一. 计算机视觉技术中常见的几种滤波器二.滤波器相关知识应用一:算术均值.几何均值.谐波逆谐波一.问题分析二.结果图三.代码附录应用二:维纳滤波,逆滤波一 ...

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

文章目录

前言

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

1、import

2、vesselness2d

3、应用示例（原文）

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

1、import

2、main

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

1、import

2、vesselness3d

总结

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

最新文章

热门文章