CN102194229A - 图像处理设备、方法及程序 - Google Patents

Abstract

本发明公开了图像处理设备、方法及程序。防止了对包含于图像中的线状结构和板状结构的错误判别。滤波单元(32)计算图像中的任意像素位置处的像素值的二阶偏导数矩阵和至少一个一阶偏导数值。估计单元(34)基于二阶偏导数矩阵的值计算像素位置上成线状结构的可能性的估计值和成板状结构的可能性的估计值，使得一阶偏导数值越大，估计值越小。

Description

图像处理设备、方法及程序

技术领域

本发明涉及一种用于判别图像中的线状结构或板状结构的图像处理设备和图像处理方法，以及涉及一种用于使计算机执行该图像处理方法的程序。

背景技术

近几年来随着医疗设备(诸如多检测器CT)的发展，高质量的三维图像被用在影像诊断中。这种三维图像由许多二维图像形成，从而具有大量的信息。因此，医生可能花费一定时间来找到期望的要被观察的部分以及对该部分进行诊断。为了解决这个问题，已经实践的方法是：提取所关注的器官并且以MIP、VR、CPR或类似方式进行显示，以增强整个器官和损伤的可见性，从而提高诊断的效率。

另一方面，作为在医学图像中提取血管和骨骼的技术，使用海赛矩阵(Hessian matrix)的海赛分析已被提出(见A.F.Frangi等人的“Multiscale vessel enhancement filtering”，MICCAI，Vol.1496/1988，PP.130-137，1998，下文中该文献将被称为非专利文献1)。海赛分析被用来分析海赛矩阵的特征值(eigenvalue)，以判别图像中的局部结构是点、线还是面，其中所述海赛矩阵包含作为其元素的使用诸如高斯核(Gaussian Kernel)之类的预定滤波器的二阶导数核计算得到的二阶偏导数。使用海赛分析使得能够将血管判别为“线状结构”，而将骨骼判别为“板状结构”。

然而，由于血管和骨骼具有各种厚度和尺寸，因此仅仅使用海赛分析可能不能实现血管和骨骼的判别。例如，在诸如脊柱之类的皮质骨作为板状结构被提取出来的情况下，如果使用了其尺寸足以检测到皮质骨的表面部分的滤波器核，则不仅骨骼100的板状结构而且厚血管102的表面都有可能被判别为板状结构，如图8中的斜线区域所示。出现该问题的原因是二阶偏导数在垂直于厚血管表面的方向上高而在平行于血管表面的两个方向上低，因此，厚血管所表现出的特征与板状结构的相同。相反，在血管作为线状结构被提取出来的情况下，如果使用了其尺寸足以检测到血管的线状部分的滤波器核，则不仅血管102的形成线状结构的部分而且还有骨骼100的角部分都有可能被判别为线状结构，如图9中的斜线区域所示。

发明内容

鉴于上述情况，本发明的目的是防止错误地判别图像中包含的线状结构和板状结构。

根据本发明的图像处理设备的第一方面包括：

导数值计算装置，用于计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；以及

估计装置，用于基于二阶偏导数矩阵的值来计算像素位置的成线状结构的可能性的估计值和/或成板状结构的可能性的估计值，其中一阶偏导数值越大，所述估计装置输出的估计值越小。

在根据本发明的图像处理设备的第一方面中，所述导数值计算装置可以使用具有不同尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值，其中用于计算一阶偏导数值的滤波器的尺寸大于用于计算二阶偏导数矩阵的滤波器的尺寸。

在根据本发明的图像处理设备的第一方面中，所述导数值计算装置可以对图像进行多分辨率变换，以获得具有不同分辨率的分辨率图像，以及在分辨率图像的每个对应像素位置上使用具有预定尺寸的滤波器计算二阶偏导数矩阵和一阶偏导数值，其中用于计算一阶偏导数值的分辨率图像的分辨率小于用来计算二阶偏导数矩阵的分辨率图像的分辨率。

在根据本发明的图像处理设备的第一方面中，导数值计算装置可以使用一维基本高斯核(basic Gaussian Kernel)、通过基本高斯核的一阶微分获得的一阶导数核和通过基本高斯核的二阶微分获得的二阶导数核来计算二阶偏导数矩阵和一阶偏导数值。

根据本发明的图像处理设备的第一方面可以进一步包括分割装置，用于通过基于估计值为图像的每个像素设置属于目标区域的可能性、属于背景区域的可能性、以及相邻像素属于同一区域的可能性来分割目标区域和背景区域。

在这种情况下，所述分割装置可以基于成线状结构的可能性的估计值来设置属于目标区域的可能性，以及可以基于成板状结构的可能性的估计值来设置属于背景区域的可能性。

而且，分割装置可以基于以下估计值来设置属于同一区域的可能性，这些估计值为：仅仅基于二阶偏导数矩阵的值计算出的像素位置上成线状结构的可能性的估计值和成板状结构的可能性的估计值，或者通过降低一阶偏导数值影响计算出来的像素位置上成线状结构的可能性的估计值和成板状结构的可能性的估计值。

根据本发明的图像处理方法的第一方面包括：

计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值，并且基于二阶偏导数矩阵的值计算像素位置上成线状结构的可能性的估计值和/或成板状结构的可能性的估计值，其中一阶偏导数值越大，估计值越小。

根据本发明的图像处理设备的第二方面包括：

导数值计算装置，用于计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；以及

估计装置，用于基于二阶偏导数矩阵的值来计算像素位置的成线状结构的可能性的估计值的可能性的估计值，其中所述估计装置基于所述一阶偏导数值的幅值(magnitude)来改变所述估计值。

在根据本发明的图像处理设备的第二方面中，所述估计装置可以根据与所述二阶偏导数矩阵的每个值的二阶偏导方向相对应的一阶偏导数值的幅值来校正所述二阶偏导数矩阵中的每个值，并且基于所述二阶偏导数矩阵的校正值来计算所述估计值。

在根据本发明的图像处理设备的第二方面中，所述导数值计算装置可以使用具有不同尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值，并且用于计算一阶偏导数值的滤波器的尺寸大于用于计算二阶偏导数矩阵的滤波器的尺寸。

在根据本发明的图像处理设备的第二方面中，所述导数值计算装置可以对图像进行多分辨率变换，以获得具有不同分辨率的分辨率图像，并且在分辨率图像的每个对应像素位置上使用具有预定尺寸的滤波器计算二阶偏导数矩阵和一阶偏导数值，其中用于计算一阶偏导数值的分辨率图像的分辨率小于用来计算二阶偏导数矩阵的分辨率图像的分辨率。

在根据本发明的图像处理设备的第二方面中，导数值计算装置可以使用一维基本高斯核、通过基本高斯核的一阶微分获得的一阶导数核和通过基本高斯核的二阶微分获得的二阶导数核来计算二阶偏导数矩阵和一阶偏导数值。

根据本发明的图像处理设备的第二方面可以进一步包括分割装置，用于通过基于估计值为图像的每个像素设置属于目标区域的可能性、属于背景区域的可能性、以及相邻像素属于同一区域的可能性来分割目标区域和背景区域。

在这种情况下，所述分割装置可以基于成线状结构的可能性的估计值来设置属于目标区域的可能性。

而且，分割装置可以基于以下估计值来设置属于同一区域的可能性，这些估计值为：仅仅基于二阶偏导数矩阵的值计算出的像素位置上成线状结构的可能性的估计值，或者通过降低一阶偏导数值的影响而计算出来的像素位置上成线状结构的可能性的估计值。

根据本发明的图像处理方法的第二方面包括：

计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；和

基于二阶偏导数矩阵的值计算像素位置上成线状结构的可能性的估计值，其中所述估计值基于所述一阶偏导数值的幅值而改变。

根据本发明的图像处理方法可以以使计算机执行图像处理方法的第一和/或第二方面的程序的形式来提供。

在使用海赛矩阵判别包含于医学图像中的诸如血管之类的线状结构期间错误地提取了骨骼结构的情况下，以及在判别包含于医学图像中的诸如皮质骨之类的板状结构期间错误地提取了血管的情况下，在错误提取的区域的亮度值图案中存在一维偏差。例如，在图8所示的错误提取的情况下，可以理解为存在这种一维偏差：血管内部区域的亮度值较高，而血管外部区域的亮度值较低。由于在亮度值图案中存在的这种一维偏差，一阶偏导数值变大。相反，理想的线状结构和理想的板状结构关于其中心对称，因此一阶偏导数值为0。根据本发明的第一方面，计算图像中的任意像素位置上的像素值的二阶偏导数矩阵和一阶偏导数值，并且基于二阶偏导数矩阵的值来计算像素位置上成线状结构的可能性的估计值和/或成板状结构的可能性的估计值，从而一阶偏导数值越大，估计装置输出的估计值越小。因此，防止了错误地判别包含于图像中的线状结构和板状结构，从而能够实现对线状结构和板状结构的准确判别。

另一方面，针对理想的线状结构，当通过对二阶偏导数矩阵的值应用特征值分解来得到三个方向上的特征值时，在垂直于线状结构形成的组织主轴的方向上的两个特征值基本上彼此相等。不过，在例如沿着心脏外围运行的冠状动脉的情况下，当线状结构出现在板状结构周围时，从线状结构到板状结构的方向(与板状结构正交的方向)上的特征值变大。因此，提供了低的成线状结构的可能性的估计值，这导致无法成功确定线状结构。根据本发明的图像处理设备和方法的第二方面，计算图像中任意像素位置处的像素值的二阶偏导数矩阵和一阶偏导数值，并且基于所述二阶偏导数矩阵的值来计算像素位置上成线状结构的可能性的估计值，其中所述估计值基于所述一阶偏导数值的幅值而变化。特别是，根据与所述二阶偏导数矩阵的每个值的二阶偏导方向相对应的一阶偏导数值的幅值来校正所述二阶偏导数矩阵中的每个值，并且基于所述二阶偏导数矩阵的校正值来计算所述估计值。这样，即使当板状结构附近存在线状结构时，也可以在沿着存在板状结构的方向以及垂直于此方向的方向上提供基本相等的二阶偏导数矩阵的值。这能够避免降低关于成线状结构的可能性的估计值，从而实现对线状结构的精确辨别。而且，在使用具有不同尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值的情况下，用来计算一阶偏导数值的滤波器的尺寸大于用来计算二阶偏导数矩阵的滤波器的尺寸，或者在通过对图像进行多分辨率变换来获得具有不同分辨率的分辨率图像以及在分辨率图像的每个对应像素位置使用具有预定尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值的情况下，用来计算一阶偏导数值的分辨率图像的分辨率小于用来计算二阶偏导数矩阵的分辨率图像的分辨率。这使得能够更容易地捕获亮度值的一维偏差，从而能够更准确地计算出成线状结构的可能性的估计值和成板状结构的可能性的估计值。

更进一步地，在使用一维基本高斯核、通过基本高斯核的一阶微分获得的一阶导数核和通过基本高斯核的二阶微分获得的二阶导数核来计算二阶偏导数矩阵和一阶偏导数值的情况下，能够通过相对简单的运算来实现二阶偏导数矩阵和一阶偏导数值的计算，从而提高运算速度。

更进一步地，在通过基于估计值为图像的每个像素设置属于目标区域的可能性、属于背景区域的可能性以及相邻像素属于同一区域的可能性来分割目标区域和背景区域的情况下，能够实现目标区域和背景区域的准确分割。

具体地讲，通过基于成线状结构的可能性的估计值来设置属于目标区域的可能性，以及基于成板状结构的可能性的估计值来设置属于背景区域的可能性，能够实现对图像中的诸如血管之类的线状结构和诸如骨骼之类的板状结构的准确分割。

附图说明

图1是示出根据本发明实施例的图像处理设备的构造的示意性框图；

图2是用于说明多分辨率变换的示图；

图3是示出高斯核的示图；

图4是用于说明线状结构的特征值的示图；

图5是用于说明板状结构的特征值的示图；

图6是用于说明Graph Cut区域分割方法的示图；

图7是示出本发明实施例中执行的处理的流程图；

图8是用于说明板状结构的错误检测的示图；以及

图9是用于说明线状结构的错误检测的示图。

图10是用于说明在板状结构附近出现线状结构的情况下对线状结构的错误检测。

具体实施方式

以下，将参照附图描述本发明的实施例。图1是示出根据本发明实施例的图像处理设备的构造的示意性框图。应当注意，图1所示的图像处理设备1的构造是通过在计算机(诸如个人计算机)上运行已装载到计算机辅助存储设备(未示出)中的程序来实现的。程序可以通过存储在诸如CD-ROM之类的信息存储介质中或者通过诸如Internet之类的网络来分布，以安装到计算机上。

图像处理设备1例如使用X射线CT设备2拍摄的大量二维图像来产生三维图像M0，并且自动分割包含于三维图像M0中的线状结构或者板状结构。图像处理设备1包括图像获得单元10、检测区域设置单元20、判别单元30、分割单元40、显示单元50和输入单元60。

图像获得单元10例如获得利用X射线CT设备2拍摄的大量CT图像(二维图像)，并且根据二维图像产生三维图像M0。应当注意，图像获得单元10不限于获得CT图像的设备，可以是获得二维图像(诸如所谓的MRI图像、RI图像、PET图像或X射线图像)的任意设备。

检测区域设置单元20首先将三维图像M0的体素大小(voxelsize)变换为各向同性体素大小。例如，如果三维图像M0的体素大小在三维图像M0的X方向、Y方向和Z方向上是0.3mm×0.3mm×0.6mm，则三维图像M0的体素大小被变换为(X，Y，Z)＝(0.5，0.5，0.5)(mm)的各向同性体素大小。将三维图像M0的体素大小以此方式变换为各向同性体素大小使得判别单元30能够在X方向、Y方向和Z方向上应用同一尺寸的核(以下将进行描述)，从而可以简化运算。

在三维图像M0变换为各向同性体素大小之后，检测区域设置单元20对三维图像M0进行多分辨率变换，以产生如图2所示的具有不同分辨率(高斯金字塔)的三维多分辨率图像Msi(i＝0至n)。应当注意，“i＝0”表示该图像的分辨率与三维图像M0的分辨率相同，而“i＝n”表示该图像具有最低分辨率。三维多分辨率图像Msi的体素大小是：按从最高分辨率开始的顺序，(X，Y，Z)＝(0.5，0.5，0.5)、(1.0，1.0，1.0)、(2.0，2.0，2.0)等等。

判别单元30包括滤波单元32和估计单元34。滤波单元32对每个三维多分辨率图像Msi执行利用高斯核的滤波，以使用海赛矩阵来执行海赛分析。即，针对具有不同分辨率的三维多分辨率图像Msi中的每一个，对同一大小的滤波器核(σ＝1.0)进行卷积。这基本上等同于对三维图像M0应用具有不同大小的滤波器核，并且允许检测具有不同大小的线状结构(例如，血管)和板状结构(例如，诸如皮质骨之类的骨骼)。

现在，描述海赛分析。海赛分析中使用的海赛矩阵是用于三维图像的3×3矩阵，如下公式(1)所示：

${\▿}^{2}I=\begin{array}{}\end{array}\left[\begin{array}{ccc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}& I_{\mathrm{xz}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}& I_{\mathrm{yz}}\\ I_{\mathrm{zx}}& I_{\mathrm{zu}}& I_{\mathrm{zz}}\end{array}\right]---\left(1\right)$ $I_{\mathrm{xx}}=\frac{{\∂}^{2}I}{{\∂x}^2},$ $I_{\mathrm{xy}}=\frac{{\∂}^{2}I}{\∂x\∂y},......$

而且，当使用高斯核函数f时，利用一维基本核、通过基本核的一阶微分获得的一阶导数核、以及通过基本核的二阶微分获得的二阶导数核来得到用于获得海赛矩阵的滤波器系数，如下公式(2)所示：

$\left.\begin{array}{c}f=\exp (-\frac{x^2}{2{\mathrm{\σ}}^2})\\ \frac{\∂f}{\∂x}=(-\frac{x}{{\mathrm{\σ}}^2})\×f\\ \frac{{\∂}^{2}f}{\∂x^2}=(\frac{x^2}{{\mathrm{\σ}}^4}-\frac{1}{{\mathrm{\σ}}^2})\×f\end{array}\right\}---\left(2\right)$

应当注意，公式(2)仅仅示出了X方向的滤波器系数。对于Y方向和Z方向中的每一个方向，能够以相同的方式得到基本核、一阶导数核和二阶导数核。如果基本核如图3中的实线所示，则一阶导数核如图3中的虚线所示，而二阶导数核如图3中的点划线所示。

例如，可以通过以下方式计算作为海赛矩阵中的二阶偏导数值的元素Ixx，该方式为：针对三维多分辨率图像Msi中要被处理的每个像素，对X方向上的二阶导数核进行卷积以及对Y方向和Z方向中的每一个方向上的基本核进行卷积，即，对三维多分辨率图像Msi中要被处理的像素进行滤波。而且，可以通过以下方式计算海赛矩阵中的元素Ixy，该方式为：针对要被处理的每个像素，对X方向和Y方向中的每一个方向上的一阶导数核进行卷积以及对Z方向上的基本核进行卷积。

已知，当对如此计算出的海赛矩阵应用特征值分解以提供特征值时，线状结构的特征值具有以下特性：如图4所示，三个特征值中的两个特征值的值较大，而一个值接近0。例如，对于线状结构形成的目标组织来说，公式(1)的特征值的关系如下公式(3)所示：

的特征值：λ₁，λ₂，λ₃

λ₁≈0

(3)

λ₂，λ₃＞＞0

λ₂≈λ₃

而且，已知板状结构的特征值具有以下特性：如图5所示，三个特征值中的一个特征值的值较大，而两个值接近0。例如，对于板状结构形成的目标组织来说，公式(1)的特征值的关系如公式(4)所示：

λ₁～0，λ₂～0    (4)

λ₃＞＞0

应当注意，在图4和图5中，e1、e2和e3代表特征值λ₁、λ₂和λ₃的特征向量的方向。

因此，可以根据特征值来确定成线状结构的可能性和成板状结构的可能性，并且确定的结果被用来在三维图像M0中分割成线状结构的血管区域和成板状结构的骨骼区域。

滤波单元32首先针对三维多分辨率图像Msi的每一个像素来计算X方向、Y方向和Z方向中的每一个方向的一阶偏导数值。针对三维多分辨率图像Msi的每一个像素，通过对X方向上的一阶导数核进行卷积以及对Y方向和Z方向中的每一个方向上的基本核进行卷积来计算X方向上的一阶偏导数值。针对三维多分辨率图像Msi的每一个像素，通过对Y方向上的一阶导数核进行卷积以及对X方向和Z方向中的每一个方向的基本核进行卷积来计算Y方向上的一阶偏导数值。针对三维多分辨率图像Msi的每一个像素，通过对Z方向上的一阶导数核进行卷积以及对X方向和Y方向中的每一个方向上的基本核进行卷积来计算Z方向上的一阶偏导数值。这样计算出的X方向、Y方向和Z方向上的一阶偏导数值分别假设为ρ_x、ρ_y和ρ_z。

而且，滤波单元32针对三维多分辨率图像Msi中的每一个像素通过计算X方向、Y方向和Z方向中的每一个方向的二阶偏导数值来计算海赛矩阵的元素。如上所述，针对三维多分辨率图像Msi的每一个像素，通过对X方向上的二阶导数核进行卷积以及对Y方向和Z方向中的每一个方向上的基本核进行卷积来计算海赛矩阵中的元素I_xx。针对三维多分辨率图像Msi中的每一个像素，通过对X方向和Y方向中的每一个方向上的一阶导数核进行卷积以及对Z方向上的基本核进行卷积来计算海赛矩阵中的元素I_xy。

估计单元34对滤波单元32计算出的海赛矩阵进行特征值分解以计算三个特征值λ₁、λ₂和λ₃。特征值λ₁、λ₂和λ₃假设为满足以下关系：|λ₁|≤|λ₂|≤|λ₃|。然后，针对三维多分辨率图像Msi的每一个像素，计算成线状结构的可能性的估计值L0(lineness)和成板状结构的可能性的估计值P0(planeness)，如下公式(5)和(6)所示：

$L0\left(\mathrm{Lineness}\right)=(1-\exp (-\frac{{R_A}^2}{2a^2}))\exp (-\frac{{R_B}^2}{2b^2})(1-\exp (-\frac{{S_{2\mathrm{nd}}}^2}{2c^2}))\exp (-\frac{{S_{1\mathrm{st}}}^2}{2d^2})---\left(5\right)$

$P0\left(\mathrm{planeness}\right)=\exp (-\frac{{R_A}^2}{2e^2})\exp (-\frac{{R_C}^2}{2f^2})(1-\exp (-\frac{{S_{2\mathrm{nd}}}^2}{2g^2}))\exp (-\frac{{S_{1\mathrm{st}}}^2}{2h^2})---\left(6\right)$

公式(5)和(6)中的符号“a”到“h”是常数。而且，R_A、R_B和R_C根据下面的公式(7)至(9)来计算：

$R_A=\left|\frac{{\mathrm{\λ}}_2}{{\mathrm{\λ}}_3}\right|---\left(7\right)$

$R_B=\frac{\left|{\mathrm{\λ}}_1\right|}{\sqrt{{\mathrm{\λ}}_2{\mathrm{\λ}}_3}}---\left(8\right)$

$R_C=\frac{\sqrt{{\mathrm{\λ}}_1{\mathrm{\λ}}_2}}{\left|{\mathrm{\λ}}_3\right|}---\left(9\right)$

S_2nd和S_1st是二阶偏导数值和一阶偏导数值的乘方(power)，如下根据公式(10)和(11)来计算：

$S_{2\mathrm{nd}}=\sqrt{{{\mathrm{\λ}}_1}^2+{{\mathrm{\λ}}_2}^2+{{\mathrm{\λ}}_3}^2},---\left(10\right)$

$S_{1\mathrm{st}}=\sqrt{{{\mathrm{\ρ}}_x}^2+{{\mathrm{\ρ}}_y}^2+{{\mathrm{\ρ}}_z}^2},---\left(11\right)$

作为用在公式(5)和(6)中的一阶偏导数值ρ_x、ρ_y和ρ_z，使用通过以下方式计算的值：使用的滤波器尺寸比用于计算二阶偏导数值的滤波器尺寸大，即，使用三维多分辨率图像Msi中分辨率较低的图像。特别地，使用三维多分辨率图像Msi中分辨率大约在高斯金字塔中低一个层次的图像来计算一阶偏导数值ρ_x、ρ_y和ρ_z。这使得能够更容易地捕获三维多分辨率图像Msi的亮度值的一维偏差，从而能够更准确地计算成线状结构的可能性和成板状结构的可能性的估计值L0和P0。

在使用海赛矩阵判别包含于医学图像中的诸如血管之类的线装结构期间错误地提取骨骼结构以及在判别包含于医学图像中的诸如皮质骨之类的板状结构期间错误地提取血管的情况下，在错误提取的区域的亮度值图案中存在一维偏差。由于亮度值图案中存在的这种一维偏差，一阶偏导数值变大。相反，理想的线状结构和理想的板状结构是关于其中心对称的，因此一阶偏导数值为0。因此，当一阶偏导数值较大，即乘方S_1st较大时，成线状结构的可能性和成板状结构的可能性小。因此，公式(5)和公式(6)中的最后一项exp(-S_1st/2d²)和exp(-S_1st/2h²)的值越小，成线状结构的可能性和成板状结构的可能性的估计值L0和P0越小。

在该实施例中，针对具有不同分辨率的多分辨率三维图像中的每一个图像，计算出成线状结构的可能性的估计值L0和成板状结构的可能性的估计值P0。这样计算出的估计值L0和P0用作原始三维图像M0的对应像素位置的估计值。针对每个多分辨率三维图像Msi的对应像素位置，计算每个估计值，以及在该实施例中，针对三维图像Msi的对应像素位置计算出的估计值中的最高估计值被用作原始三维图像M0的像素位置的估计值。

基于判别单元30计算出来的成线状结构的可能性的估计值L0和成板状结构的可能性的估计值P0，分割单元40执行三维图像M0的血管区域和血管以外的区域(包括骨骼)的区域分割。特别地，通过将血管区域设置为目标区域而将血管区域以外的区域设置为背景区域，分割单元40为三维图像M0的每个像素位置设置具有预定像素大小的判别区域，并且使用Graph Cut区域分割方法将判别区域分割成目标区域和背景区域。在Y.Y.Boykov和Marie-Pierre Jolly的“Interactive Graph Cuts for Optimal Boundary & Region Segmentationof Objects in N-D images”，“International Conference on ComputerVision”学报，Vol.I，pp.105-112，2001中对Graph Cut区域分割方法进行了描述。

在Graph Cut区域分割方法中，首先产生图6所示的图，其包括表示判别区域中的各个像素的节点N_ij、表示目标区域(血管区域)的节点S、表示背景区域的节点T、每一个均表示相邻像素的每对节点属于同一区域的可能性的N-link、将表示各个像素的各个节点N_ij连接至表示目标区域的节点S的S-link、以及将各个节点N_ij连接至表示背景区域的节点T的T-link。在图6中，为了便于说明，判别区域示为3×3的二维区域。

将表示每个像素的每个节点N_ij连接至表示目标区域的节点S的每个S-link用连线的粗度(值的大小)来表示每个像素为属于目标区域的像素的可能性。将表示每个像素的每个节点N_ij连接至表示背景区域的节点T的每个T-link用连线的粗度(值的大小)表示每个像素为属于背景区域的像素的可能性。对于判别单元30计算出的较大的成线状结构的可能性的估计值L0，设置较大的S-link值(提供较粗的连线)。类似地，对于判别单元30计算出的较大的成板状结构的可能性的估计值P0，设置较大的T-link值(提供较粗的连线)。在该实施例中，公式(5)和公式(6)中的最后一项exp(-S_1st/2d²)和exp(-S_1st/2h²)的较小值提供较小的成线状结构的可能性和成板状结构的可能性的估计值L0和P0，并且表示目标区域T-link和表示背景区域S-link不变粗。

除了估计值L0和P0之外，还可以根据如下公式(12)和(13)使用每个像素的亮度值来设置S-link和T-link的值：

S-link＝g1(L0)×g2(I)    (12)

T-link＝g3(P0)×g4(I)    (13)

在公式(12)中，“g1()”是以下函数：对于较大的成线状结构的可能性的估计值L0输出较大的值，而“g2()”是以下函数：对于三维图像M0的每个像素的在统计上较接近血管亮度值的亮度值(CT值)I输出较大的S-link值。在公式(13)中，“g3()”是如下函数：对于较大的成板状结构的可能性的估计值P0输出较大的值，而“g4()”是如下函数：对于三维图像M0的每个像素的在统计上较接近骨骼亮度值的亮度值I输出较大的T-link值

每个N-link用连线的粗度(值的大小)来表示每对相邻像素属于同一区域的可能性。对于较大的成线状结构的可能性的估计值L0’，设置较大的N-link值，其中，估计值L0′是根据类似于公式(5)的其中去除了公式(5)中的最后一项exp(-S_1st/2d²)的公式计算出来的。可以通过降低公式(5)中的最后一项exp(-S_1st/2d²)的权重来计算估计值L0′。在图6所示的示例中，表示血管区域(为该实施例的目标区域)的所有像素属于同一区域的事实由N-link反映出来。因此，根据血管区域中设置的S-link，能够将整个血管区域与三维图像M0中的其他区域分割开来。具体地讲，在有噪声的图像的情况下，即使是血管区域的像素之间也存在大的亮度差。这将导致血管区域中像素的相似度低，从而导致细的N-link。然而，通过使用根据类似于公式(5)的其中去除了公式(5)中的最后一项exp(-S_1st/2d²)的公式获得的成线状结构的可能性的估计值(Lineness)，能够提高N-link的值。

应当注意，除了估计值L0′，还可以根据如下公式(14)使用N-link连接的两个像素的亮度值来设置N-link的值：

N-link＝g5(L0′)×g6((I1+I2)/2)×g7(I1-I2)    (14)

在公式(14)中，“g5()”是如下函数：对于较大的估计值L0′输出较大的值，“g6()”是如下函数：对于在统计上较接近血管的亮度值和骨骼的亮度值的平均值的两个像素的亮度值I1和I2的平均值输出较大的N-link值，以及“g7()”是如下函数：对于两个像素的亮度值I1和I2之间的较小差值输出较大的N-link值。

在图6所示的示例中，由节点N₁₁、N₁₂、N₂₁、N₂₂和N₃₁表示的像素具有大的成线状结构的可能性的估计值L0，因此节点N₁₁、N₁₂、N₂₁、N₂₂和N₃₁通过粗的S-link连接至节点S。另外，由节点N₁₁、N₁₂、N₂₁、N₂₂和N₃₁表示的像素具有大的成线状结构的可能性的估计值L0′，其中估计值L0′是根据类似于公式(5)的其中去除了公式(5)中的最后一项exp(-S_1st/2d²)的公式计算出来的，因此节点N₁₁、N₁₂、N₂₁、N₂₂和N₃₁中的每一个通过粗的N-link连接至每个相邻节点。另一方面，由节点N₁₃、N₂₃、N₃₂和N₃₃表示的像素具有大的成板状结构的可能性的估计值P0，因此节点N₁₃、N₂₃、N₃₂和N₃₃通过粗的T-link连接至节点T。

由于血管区域和血管区域之外的区域(包括骨骼)是相互排他的区域，因此通过切断S-link、T-link和N-link中的适当连线以把节点S与节点T分开，能够将判别区域分割成目标区域和背景区域，如图6的虚线所示。在这种情况下，通过切断连线以使得将被切断的S-link、T-link和N-link的值的总数最小，能够实现理想的区域分割。

显示单元50是诸如CRT屏幕之类的显示器，用于显示二维图像或三维图像。在该实施例中，显示单元50执行作为目标区域而被分割出来的线状结构和板状结构的体积重建显示(volume-rendereddisplay)，以提供整个线状结构和板状结构的概况并且显现它们的连续性。

输入单元60例如包括键盘和鼠标。

接下来，描述在该实施例中执行的处理。图7是示出该实施例中执行的处理的流程图。首先，图像获得单元10根据X射线CT设备2拍摄的二维图像来产生三维图像M0(步骤ST1)。然后，检测区域设置单元20将三维图像M0的大小变换为各向同性体素大小并且对经过变换的三维图像M0进行多分辨率变换，以产生具有不同分辨率的三维多分辨率图像Msi(步骤ST2)。

随后，判别单元30的滤波单元32使用高斯核对每个三维多分辨率图像Msi进行滤波，以计算每个像素位置上的一阶偏导数值和二阶偏导数矩阵(步骤ST3)。然后，估计单元34计算成线状结构的可能性的估计值L0和成板状结构的可能性的估计值P0(步骤ST4)。

然后，分割单元40使用上述的Graph Cut区域分割方法将三维图像M0分割成目标区域(血管区域)和背景区域(步骤ST5)。最后，显示单元50执行被分割的目标区域和背景区域的体积重建显示(步骤ST6)，并且该处理结束。

如上所述，在该实施例中，计算三维图像M0的每个像素位置处的像素值的二阶偏导数矩阵和一阶偏导数值，并且基于二阶偏导数矩阵的特征值来计算该像素位置上的成线状结构的可能性和成板状结构的可能性的估计值L0和P0，从而一阶偏导数值越大，估计值L0和P0越小。即，形状与理想的线状结构或理想的板状结构的差别越大，估讨值L0和P0越小。因此，根据该实施例，可以防止包含于三维图像M0中的线状结构和板状结构的错误判别，以实现对包含于三维图像M0中的诸如血管之类的线状结构和诸如骨骼之类的板状结构的精确判别。

由于使用了在计算二阶偏导数矩阵时必然要计算的一阶导数核来计算一阶偏导数值，所以与执行常规的海塞分析相比，本发明不需要增加计算量。

而且，通过利用基于估计值L0和P0而为三维图像M0的每个像素设置S-link、T-link和N-link来执行将三维图像M0分割成目标区域和背景区域，其中S-link表示该像素属于目标区域的可能性，T-link表示该像素属于背景区域的可能性，以及N-link表示每对相邻像素属于同一区域的可能性，能够实现目标区域和背景区域的准确分割。

应当注意，虽然本发明在上述实施例中应用于对包含于三维图像M0中的线状结构和板状结构的判别，但是本发明还适用于对包含于二维图像中的线状结构和板状结构的判别。

进一步地，虽然在作为示例的上述实施例中对成线状结构的血管进行了判别，但是本发明还适用于对诸如支气管之类的其他线状结构进行判别。而且，本发明还适用于对除骨骼之外的诸如皮肤、叶间胸膜等之类的其他板状结构进行判别。

进一步地，虽然在上述实施例中使用Graph Cut区域分割方法来实现对包含于三维图像M0中的线状结构和板状结构的分割，但是还可以使用诸如Watershed算法之类的其他区域分割技术来实现分割。Watershed算法是以以下方式分割图像的技术：当水倾泻在将图像的像素值信息视作高度的地形上时，在不同的洼地中形成的水池之间形成边界。在这种情况下，可以通过在执行Watershed算法之前对三维图像M0的估计值L0和P0进行适当的平滑来实现线状结构和板状结构的分割。

在上述实施例中，通过对三维图像M0进行多分辨率变换以产生三维多分辨率图像Msi以及利用具有单一尺寸的滤波器计算一阶偏导数值和二阶偏导数矩阵，来判别具有各种尺寸的线状结构和板状结构。然而，可以将具有不同尺寸的滤波器应用到三维图像M0，以计算一阶偏导数值和二阶偏导数矩阵。

而且，虽然在上述实施例中不仅计算了成线状结构的可能性的估计值L0而且计算了成板状结构的可能性的估计值P0，但是可以仅仅计算成线状结构的可能性的估计值L0和成板状结构的可能性的估计值P0。

进一步地，虽然在上述实施例中通过对海塞矩阵进行特征值分解而获得的特征值被用来计算成线状结构的可能性的估计值L0和成板状结构的可能性的估计值P0，但是可以在不获得特征值的情况下使用海塞矩阵的元素来计算成线状结构的可能性和成板状结构的可能性的估计值。

针对理想的线状结构，当通过对上述海塞矩阵应用特征值分解来得到三个特征值λ₁、λ₂和λ₃时，在垂直于沿着线状结构延伸的方向的方向上的两个特征值λ₂和λ₃基本上彼此相等，如上述公式(3)所示。不过，在如图10所示沿着心脏110外围运行的冠状动脉111的情况下，当线状结构(冠状动脉111)出现在板状结构(心脏110的外表面)周围时，从线状结构到板状结构的方向(即与板状结构正交的方向)上的特征值变大。例如，如果与板状结构(心脏110的外表面)正交的方向对应于特征向量e3的方向，则这个方向上的特征值λ₃会大于与这个方向垂直的方向上的特征值λ₂。这是因为由于心脏外部的肺的存在，使得在心脏内侧和外侧的CT值上存在较大差异。因此，提供了低的关于成线状结构的可能性的估计值L0，并且这会导致无法成功地确定线状结构。

现在已经描述了作为本发明第二实施例的解决上述问题的技术。由于按照上述方式计算出的X方向、Y方向、和Z方向上的一阶偏导数值ρ_x、ρ_y和ρ_z的方向与特征值λ₁、λ₂和λ₃的特征向量e1、e2和e3的方向不相同，第二实施例中的估计单元34根据下述公式(15)来计算与特征向量e1、e2和e3的方向相对应一阶偏导数值ρ₁、ρ₂和ρ₃：

ρ₁＝ρ_x×x₁+ρ_y×y₁+ρ_z×z₁

ρ₂＝ρ_x×x₂+ρ_y×y₂+ρ_z×z₂    (15)

ρ₃＝ρ_x×x₃+ρ_y×y₃+ρ_z×z₃

假设特征值λ₃在特征值λ₂和λ₃中较大，则估计单元34根据下述公式(16)来校正特征值λ₃。也即，估计单元34根据特征向量e2和e3方向上的一阶偏导数值ρ₂和ρ₃的幅值对特征值λ₃进行校正来得到校正的特征值λ₃’。

λ3′＝λ2+(λ3-λ2)×f(|ρ2-ρ3|)    (16)

在公式(16)中，f()是输出0至1范围内的值的函数，其中针对较大的|ρ₂-ρ₃|值输出较小值。因此，|ρ₂-ρ₃|的值越大，特征值λ₃’的值越接近特征值λ₂的值，并且满足如公式(3)所示的特征值λ₂和λ₃基本彼此相等的关系。

然后，估计单元34在公式(7)至(9)中使用特征值λ₃’代替特征值λ₃来计算R_A、R_B和R_C的值，并使用计算得到的R_A、R_B和R_C的值根据公式(5)来计算关于成线状结构的可能性的估计值L0。

如上文所述，在第二实施例中，在计算关于成线状结构的可能性的估计值L0时，根据与特征值λ₂和λ₃的特征向量e2和e3的方向相对应的一阶偏导数值ρ₂和ρ₃的幅值对特征值λ₃进行校正，并基于校正的特征值λ₃’来计算估计值L0。从而，即使在板状结构附近存在线状结构时，也可以在沿着存在板状结构的方向以及垂直于此方向的方向上提供基本相等的特征值λ₂和λ₃。这能够避免降低关于成线状结构的可能性的估计值L0，从而实现对线状结构的精确辨别。

Claims (17)

1.一种图像处理设备，包括：

导数值计算装置，用于计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；以及

估计装置，用于基于二阶偏导数矩阵的值来计算像素位置上的成线状结构的可能性的估计值和/或成板状结构的可能性的估计值，其中一阶偏导数值越大，所述估计装置输出的估计值越小。

2.如权利要求1所述的图像处理设备，其中所述导数值计算装置利用具有不同尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值，以及其中用于计算一阶偏导数值的滤波器的尺寸大于用于计算二阶偏导数矩阵的滤波器的尺寸。

3.如权利要求1所述的图像处理设备，其中所述导数值计算装置对图像进行多分辨率变换，以获得具有不同分辨率的分辨率图像，以及在分辨率图像的每个对应像素位置上使用具有预定尺寸的滤波器计算二阶偏导数矩阵和一阶偏导数值，其中用于计算一阶偏导数值的分辨率图像的分辨率小于用来计算二阶偏导数矩阵的分辨率图像的分辨率。

4.如权利要求1所述的图像处理设备，其中所述导数值计算装置使用一维基本高斯核、通过基本高斯核的一阶微分获得的一阶导数核和通过基本高斯核的二阶微分获得的二阶导数核来计算二阶偏导数矩阵和一阶偏导数值。

5.如权利要求1所述的图像处理设备，还包括分割装置，用于通过基于估计值为图像的每个像素设置属于目标区域的可能性、属于背景区域的可能性、以及相邻像素属于同一区域的可能性来分割目标区域和背景区域。

6.如权利要求5所述的图像处理设备，其中所述分割装载基于成线状结构的可能性的估计值来设置属于目标区域的可能性，并且基于成板状结构的可能性的估计值来设置属于背景区域的可能性。

7.如权利要求5所述的图像处理设备，其中基于仅仅根据二阶偏导数矩阵的值计算出来的像素位置上成线状结构的可能性的估计值和成板状结构的可能性的估计值，或者基于通过降低一阶偏导数值影响计算出来的像素位置上成线状结构的可能性的估计值和成板状结构的可能性的估计值，所述分割装置设置属于同一区域的可能性。

8.一种图像处理方法，包括：

计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值，以及

基于二阶偏导数矩阵的值计算像素位置上成线状结构的可能性的估计值和/或成板状结构的可能性的估计值，其中一阶偏导数值越大，估计值越小。

9.一种图像处理设备，包括：

导数值计算装置，用于计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；以及

估计装置，用于基于二阶偏导数矩阵的值来计算像素位置的成线状结构的可能性的估计值的可能性的估计值，其中所述估计装置基于所述一阶偏导数值的幅值来改变所述估计值。

10.如权利要求9所述的图像处理设备，其中所述估计装置根据与所述二阶偏导数矩阵的每个值的二阶偏导方向相对应的一阶偏导数值的幅值来校正所述二阶偏导数矩阵中的每个值，并且基于所述二阶偏导数矩阵的校正值来计算所述估计值。

11.如权利要求9所述的图像处理设备，其中所述导数值计算装置使用具有不同尺寸的滤波器来计算二阶偏导数矩阵和一阶偏导数值，并且用于计算一阶偏导数值的滤波器的尺寸大于用于计算二阶偏导数矩阵的滤波器的尺寸。

12.如权利要求9所述的图像处理设备，其中所述导数值计算装置对图像进行多分辨率变换，以获得具有不同分辨率的分辨率图像，并且在分辨率图像的每个对应像素位置上使用具有预定尺寸的滤波器计算二阶偏导数矩阵和一阶偏导数值，其中用于计算一阶偏导数值的分辨率图像的分辨率小于用来计算二阶偏导数矩阵的分辨率图像的分辨率。

13.如权利要求9所述的图像处理设备，其中所述导数值计算装置使用一维基本高斯核、通过基本高斯核的一阶微分获得的一阶导数核和通过基本高斯核的二阶微分获得的二阶导数核来计算二阶偏导数矩阵和一阶偏导数值。

14.如权利要求9所述的图像处理设备，进一步包括分割装置，所述分割装置用于通过基于估计值为图像的每个像素设置属于目标区域的可能性、属于背景区域的可能性、以及相邻像素属于同一区域的可能性来分割目标区域和背景区域。

15.如权利要求14所述的图像处理设备，其中所述分割装置基于成线状结构的可能性的估计值来设置属于目标区域的可能性。

16.如权利要求14所述的图像处理设备，其中所述分割装置基于以下估计值来设置属于同一区域的可能性，这些估计值为：仅仅基于二阶偏导数矩阵的值计算出的像素位置上成线状结构的可能性的估计值，或者通过降低一阶偏导数值的影响而计算出来的像素位置上成线状结构的可能性的估计值。

17.一种图像处理方法，包括：

计算图像中每个像素位置上的像素值的二阶偏导数矩阵和至少一个一阶偏导数值；和

基于二阶偏导数矩阵的值计算像素位置上成线状结构的可能性的估计值，其中所述估计值基于所述一阶偏导数值的幅值而改变。

Images