1.Introduction

通常对CTA和MRA的血管成像方法中，可视化3维血管结构的方法是MIP。但是MIP的缺点是会将非血管结构展现出来并且对于对比度比较低的小血管展现不太好。这也正是TOF-MRA的缺点。

这篇论文的目的是为了增强血管结构，最终达到血管分割的目的。血管增强作为MIP和血管分割的预处理步骤，将会提升对于小血管的展示，并且减少器官的过度投影。

2.Method

这篇论文的方法中，将血管增强设想成为一个过滤过程，搜索管状几何结构。由于血管尺寸不同，因此需要引入一个在一定范围内变化的尺度。

在分析图像L信息时，通常是考虑图像中一个像素点 $x0$ 附件的泰勒展开式。

上述公式表示了图像 $x0$ 点处的二阶导数结构。 $\triangledown o,s$ 和 $Ho,s$ 是s尺度下在图像中 $x0$ 处计算得到的梯度向量和Hessian矩阵。

1.Hessian矩阵的由来和定义

由高等数学知识可知，若一元函数 $f(x)$ 在 $x=x^{(0)}$ 点的某个领域内具有任意阶导数，则 $f(x)$ z在 $x^{^{(0)}}$ 点处的泰勒展开式为：

$f(x)=f(x^{(0)})+f^{'}(x^{(0)})\Delta x+1/2f^{''}(x^{(0)})(\Delta x)^{2}+....$

其中， $\Delta x=x-x^{(0)},\Delta x^{2}=(x-x^{(0)})^{2}$ 。

二元函数 $f(x1,x2)$ 在 $\left ( x_{1}^{(0)}, x_{2}^{(0)}\right )$ 点处的泰勒展开式为：

其中， $\Delta x_{1}=x_{_{1}}-x^{(0)}_{1}$ , $\Delta x_{2}=x_{2}-x_{2}^{^{(0)}}$ ,将上述泰勒展开式写成矩阵形式，则有：

也就是：

其中：

$G\left ( X^{(0)} \right )$ 是 $f(x_{1},x_{2})$ 在 $\left ( x_{1} ^{(0)},x_{2}^{^{(0)}}\right )$ 点处的Hessian矩阵。可以将Hessian矩阵简写为：

为了以适当的方式计算图像 $L$ 的这些微分算子，作者使用了线性尺度空间理论的概念。（尺度空间理论的基本思想是：在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓提取，并且以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取。尺度空间理论的特点是：将传统的单尺度图像信息处理技术纳入尺度不断变化的动态分析框架中，更容易获取图像的本质特征。尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程）在该理论下，对图像做偏导就等于对原图以特定的高斯核做卷积操作。这样对于图的Hessian运算将极大降低计算量如下：

D维的高斯被定义为：

2.基于尺度理论的Hessian简化算法

对于二维图像，Hessian矩阵描述每个像素点在主方向上的二维导数：

根据尺度空间理论，二阶导数可以通过图像和高斯函数的卷积获得，例如，在点（x,y）处有

其中， $G\left ( x,y;\sigma \right )$ 为尺度为 $\sigma$ 的高斯函数。

参数 $\gamma$ 定义了一组归一化导数。这种标准化对于在多个尺度上比较微分算子的响应特别重要。当没有首选的尺度时， $\gamma$ 应该设置为1。

分析二阶导数信息 (Hessian ) 在检测血管的情景中有一个直观的依据。尺度 s 下，高斯核的二阶导数可以作为一个探针核在导数方向上测量区域 (−s,s) 内部区域和外部区域之间的对比度，如图 1。

公式（1）的第三部分给出了二阶有向导数，如下：

Hessian矩阵特征值分析背后隐含的想法是在图像中可以分解出局部二阶导数结构的位置提取主方向。由于该方法可以直接给出最小曲率的方向（沿血管方向），因此避免在多个方向上使用多个过滤器。

$\lambda s,k$ 是Hessian矩阵 $Ho,s$ 中第k个归一化特征向量 $\mu s,k$ 的特征值，所有的计算均是在尺度s下进行的。

通过分析公式4-6可以得到一个很好的几何解释。特征值分解提取了三个正交的方向，当这三个方向被Hessian矩阵映射时，在同一个尺度因子下是不变的。特别的，以 $x0$ 为球心，半径为1的球形领域 $Nx0$ 将被 $Hx0$ 映射到一个椭球体，该椭球体各个轴的方向与Hessian矩阵的特征向量的方向一致，并且对应轴的半长与对应特征值的大小一致。这个椭球局部描述了图像的二阶结构，可以作为几何相似性度量设计的直观工具。

3.Hessian矩阵特征值的求解方法

根据定义求解二阶矩阵的特征值：对于矩阵A，它的特征值满足 $\left | \lambda E-A \right |=0$ ,其中E是二阶对角阵。现在我们表示A为： $\begin{bmatrix} a&b \\ c& d\end{bmatrix}$ ，那么 $\left | \lambda E-A \right |=\left ( \lambda -a \right )\left ( \lambda -d \right )-bc=\lambda ^{2}-\left ( a+d \right )\lambda +ad-bc=0$

这个是一元二次方程，可以计算得到有两个解，分别为：

$\lambda 1=\left ( a+d+\sqrt{\left ( a-d \right )^{2}+4bc} \right )/2$

$\lambda 2=(a+d-\sqrt{(a-d)^{^{2}}+4bc})/2$

前面已经介绍过，Hessian矩阵的定义：

并且根据图像的特性，可以得到

$I_{xy}=I_{yx}$

带入以上方程得到Hessian矩阵的特征值的解：

4.Hessian矩阵特征值的图像性质

一个Hessian矩阵可以分解为两个特征值以及定义的特征向量。 $\lambda _{1}$ 和 $\lambda _{2}$ （ $\left | \lambda _{1} \right |\geqslant \left | \lambda _{2} \right |$ ），其中最大的绝对特征值 $\left | \lambda _{1} \right |$ 表示最大的局部灰度变化。其特征向量代表它方向，可以认为是切线方向；而较小的那个代表垂直方向，也就是法线方向。

在本文的剩下部分， $\lambda k$ 表示第k个小的特征值（ $\left | \lambda 1 \right |\leq \left | \lambda 2 \right |\leq \left | \lambda 3 \right |$ ）。在这个假设下，表1总结了为了检测不同的结构，Hessian矩阵各个特征值之间必须遵循的关系。特别的，当一个像素属于血管区域时，各个特征值符合以下关系： $\lambda 1$ 小，（理想情况下为0），而 $\lambda 2$ 和 $\lambda 3$ 大且符号相等（此处符号表示亮/暗）对应的特征向量指出了奇异方向： $\mu 1$ 表示血管方向（强度变化最小的方向）， $\mu 2\mu 3$ 构成了正交平面的基础。在MRA和CTA，血管以明亮的管状结构出现在暗色的环境中。该先验信息与图像模态相关，可以作为一致性检测，以此抛弃数据库中存在的极性与所寻找的极性不同的结构。因此，我们要寻找 $\lambda 1\lambda 2$ 同时为负的结构。

总的来说，在3D结构中，理想的管状结构如下：

原文强调了这三个特征值在局部方向模式的识别中都起着重要的作用。作者提出的差异度量考虑了基于二阶椭球体的两个几何比率。第一个比率给出了血管与斑点状结构的偏差，但是无法区分直线和平板状的图案：

对于斑点状结构，该比率将达到最大值，为0。当 $\lambda 1$ $\approx$ 0或者 $\lambda 2\lambda 3$ 趋于无穷大时，该比率为0.

第二个比率涉及到二阶椭球体的最大面积界面，该比率给出了两个最大的二阶导数的长宽比。该比率对于区分板状结构和线状结构至关重要，因为只有在线状图案中，该比率为0，如下：

当 $\lambda 2> \lambda 3$ 时为板状结构，当 $\lambda 2< \lambda 3$ 时为线状结构。

这两个几何比率是灰度不变的。这就保证了作者提出的度量只获取了图像的几何信息。但是，在MRA和CTA图像中还有额外的可用信息：血管结构比背景亮，占据整个体素中相对较小的体积。如果不合并这些信息，由于随机噪声波动，背景像素将产生不可预测的滤波器响应。然而，背景像素的一个显著特征是导数的大小很小。为了量化这一点，作者建议通过计算Hessian 范式来度量。作者使用Frobenius矩阵范式来计算Hessian矩阵的范式。因为当矩阵是实数且对称时，可以很容易根据特征值来计算。作者定义了以下二阶导数结构的度量规则：

其中D表示图像的维度。这个公式的度量在没有显示结构的背景区域中的值会小，并且缺乏对比度，导致特征值也是小的。与背景有高对比度的区域，其范式会更大。因此，至少其中一个特征值会很大。因此，作者提出以下的度量组合来定义血管函数：

其中， $\alpha \beta c$ 是控制线性滤波对度量 $R_{A}$ , $R_{B}$ , $S$ 的灵敏度的阈值。公式13隐含的思想是将公式10-12中的特征根据根据不同标准映射至血管的似然估计。作者使用乘法将各个不同的标准结合在一起。确保滤波器的响应只在3个标准均满足时2才能达到最大值。在作者的实验中， $\alpha$ 和 $\beta$ 的值为0.5。 $c$ 的值取决于图像的灰度范围，在多数情况下，Hessian范式最大值的一般是有效的。但是后续研究将是自动化选择阈值。