CN101822545B - 一种数字减影血管造影运动伪影消除方法及其系统 - Google Patents

Abstract

一种数字减影血管造影运动伪影消除方法，步骤如下：步骤1、读取；步骤2、选点；步骤3、构建DSA时空体；步骤4、时空切片；步骤5、连轨迹；步骤6，分析DSA像素的时空运动特性；步骤7、三角剖分；步骤8、仿射变换；步骤9、时空卷绕；步骤10、优化；步骤11、配准；步骤12、灰度校正；步骤13，对数减影。本发明属于图像处理技术。采用时空分析的方法，使DSA图像配准更加精确，有效去除运动伪影，得到清晰的血管造影图像。进而提高医生的诊断准确性，提高工作效率。

Description

一种数字减影血管造影运动伪影消除方法及其系统

技术领域

本发明涉及医学图像处理技术领域。特别是一种造影图像中运动伪影的消除方法和伪影消除系统。

背景技术

随着现代科学技术的发展，一些先进的技术手段和计算机科学技术正在不断地应用于医学领域。特别是在放射诊断和治疗学科中，计算机图像处理技术正发挥着越来越重要的作用。                                            

目前，数字减影血管造影（Digital Subtraction Angiography, DSA）是血管可视化的关键技术，国内最早于1984年引进了DSA设备，并在全国迅速推广。在普通的X光图像中，由于血管与其周围组织的对比度非常小，若要单纯对血管进行细微观察就较为困难。为了获得更清晰的血管图像，通常将造影剂（不透光的碘溶液）注入需要诊断的血管中以提高血管的对比度。病人在检查台上接受多次X射线的照射, 第一次是未注射造影剂时拍摄的X光片，称为蒙片(Mask Image)。经导管向病变区域或检查区域血管内快速注入有机碘水造影剂，注入血管造影剂之后拍摄的X光片称为活片或盈片(Live Image)。从注入造影剂开始到血管内造影剂浓度达到峰值这段时间内，使检查部位连续成像。在这一系列图像中，取蒙片图像与造影剂浓度达到峰值的盈片图像执行减法操作，可以去除非感兴趣组织（如肌肉、骨骼）对X射线吸收形成的背景图像，从而得到清晰的血管造影图像。

假设在血管周围组织、成像条件完全不变的情况下，将蒙片与活片相减，如果活片与蒙片的背景（如肌肉、骨骼）完全对应，且有相同的灰度分布，则这些背景将被完全消除，即理论上可以得到仅包含血管的清晰数字减影图像。但在临床应用中，这种理论假设总是不成立的。由于蒙片和活片序列在不同时刻拍摄，在成像的同时, 病人会存在各种运动，其中有些运动是不可避免的，比如呼吸、吞咽、眼球转动、肌肉运动、脏器蠕动、病人的躁动不配合以及在造影剂注人体时由于病人的自然反应而引起的局部运动等。这些运动会使活片和蒙片之间存在并非由造影剂产生的差异，从而使减影后的图像出现难以预料的背景，称为运动伪影（Motion Artifact）。

在早期DSA中，减少运动伪影的技术主要是在如何去避免病人的运动，以及改进减影方法（如：双能量减影技术）。近年来的研究中，运动伪影主要通过图像配准（Image Registration）的方法被校正。图像配准就是将不同条件下获取的两幅或多幅图像进行匹配、叠加的过程，实质上是两个图像坐标系之间的矩阵变换的计算过程，可以为刚性配准或非刚性配准。图像配准也可以看作是一个运动估计问题，为了获得两幅图像之间的对应关系，必须进行运动估计。运动估计就是自动计算图像中局部运动或某个结构位移的技术。

根据DSA原理，造影剂浓度达到峰值的盈片图像与蒙片图像相减，从而得到清晰的血管造影图像。从注入造影剂开始到血管内造影剂浓度达到峰值这段时间内（约有数秒），检查部位连续成像，将产生一系列血管逐渐清晰的盈片图像，即在一个DSA的DICOM文件中包含了几十张盈片图像数据。在DSA图像序列中，运动像素不仅与同一幅盈片图像中周围的像素相关（空间相关性），而且也与前后几幅盈片图像中周围的像素相关（时间相关性）。基于时空信息的分析方法近年来在数字视频处理中已经得到了广泛的应用，但是由于DSA图像序列中大部分是利用价值不高的、造影剂浓度没有达到峰值的盈片图像（包括到达峰值前以及部分消退图像），所以导致消除运动伪影的效果不佳。

发明内容

本发明的目的是提供一种数字减影血管造影运动伪影消除方法及其系统，要解决数字减影血管造影图像中的运动伪影的消除技术问题，并解决盈片图像与蒙片图像精细配准的技术问题。

为实现上述目的，本发明采用如下技术方案：一种数字减影血管造影运动伪影消除方法，其特征在于步骤如下：步骤1、读取；将数字减影血管造影仪器与计算机连接，在计算机中安装操作系统和图像处理系统，用图像处理软件读入一组DSA图像，并将DSA图像进行预处理。

步骤2、选点；采用边缘检测算法找到每幅DSA图像中血管影像的边缘和边缘点，计算出边缘点的梯度值，根据阈值来判断该点是否作为备选控制点，并选取至少10个备选控制点。

步骤3、构建DSA时空体；将该组DSA二维图像以时间T为轴按拍摄顺序间隔排列成平行六面体空间。

步骤4、时空切片；用切片指令对DSA时空体做垂直方向的时空切片。

步骤5、连轨迹；用图像处理软件将各时空切片对应的各控制点连线，得到各个控制点的运动轨迹。

步骤6，分析DSA像素的时空运动特性；选择运动轨迹变化最大的10个反映DSA图像背景结构运动规律的控制点。

步骤7、三角剖分；用三角剖分的方法将所有的控制点作为三角形的顶点，把图像剖分成不均匀的三角形网格。

步骤8、仿射变换(affine transfomation)；采用仿射变换法进行初步几何校正。

步骤9、时空卷绕；用Warping非线性变换法参考前一幅图像中控制点的变形位置，逐幅确定图像非控制点的位移。

步骤10、优化；采用Powell优化算法来优化相似度准则函数求极值的过程。

步骤11、配准；选择三次样条插值算法进行亚像素级的配准。

步骤12、灰度校正；逐幅对图像的灰度失真进行回朔校正。

步骤13，对数减影；用对数减影法进行减影，获得消除伪影的减影图像。

所述边缘检测算法是Canny算子或高斯-拉普拉斯算子（Gauss-Laplacian operator）。

一种数字减影血管造影运动伪影消除系统，包括DSA设备、实现运动伪影消除的DSA软件工作站、DSA胶片打印机、图像显示系统、外部存储设备、网络传输系统，它的DSA软件工作站包括以下部分：输入接口：用于接收DSA设备输出的DICOM图像。

存储器：用于存储操作系统、图像处理软件、图像分析数据和图像资料；

数字减影部分：采用本发明的算法进行数字减影。

DSA软件工作站：用于读入DSA图像，并将DSA图像进行预处理；采用边缘检测算法找到每幅DSA图像中血管影像的边缘和边缘点，计算出边缘点的梯度值，根据阈值来判断该点是否作为备选控制点，并选取至少10个备选控制点；用于构建DSA时空体；将该组DSA二维图像以时间T为轴按拍摄顺序间隔排列成平行六面体空间；用于进行时空切片；用切片指令对DSA时空体做垂直方向的时空切片；用于将各时空切片对应的各控制点连线，得到各个控制点的运动轨迹；用于分析DSA像素的时空运动特性；选择运动轨迹变化最大的10个反映DSA图像背景结构运动规律的控制点；用三角剖分的方法将所有的控制点作为三角形的顶点，把图像剖分成不均匀的三角形网格；用于采用仿射变换法进行初步几何校正；用于作时空卷绕；用Warping非线性变换法参考前一幅图像中控制点的变形位置，逐幅确定图像非控制点的位移；采用Powell优化算法来优化相似度准则函数求极值的过程；选择三次样条插值算法进行亚像素级的配准；用于逐幅对图像的灰度失真进行回朔校正；用对数减影法进行减影，获得消除伪影的减影图像。

胶片打印控制部分：控制DSA图像的胶片打印。

图像处理部分：对DICOM图像进行适当的处理。

网络传输控制：对DSA图像在PACS系统中的网络传输进行控制。

报告系统：自动生成DSA图文检查报告。

显示器：用于显示DSA图像。

输出接口：用于接打印机。

还包括用于与远程计算机或互联网连接的通讯接口。

与现有技术相比本发明具有以下特点和有益效果：这种数字减影血管造影运动伪影消除方法，采用时空体（Spatio-Temporal Volume）结构来描述DSA图像序列；在采用常规方法（如Canny算子）得到的备选控制点的基础上，对DSA时空体进行时空切片（Spatio-Temporal Slice）得到的各个控制点的运动轨迹，选择运动轨迹变化最大的若干个反映DSA图像背景结构运动规律的控制点；控制点不是仅仅分布在在蒙片或减影盈片上，而是贯穿在所有的蒙盈片序列上；然后基于这些控制点进行图像配准和数字减影以消除运动伪影。

我们选择的控制点不是仅仅分布在在蒙片或减影盈片上，而是贯穿在所有的蒙盈片序列上，更便于分析DSA图像的时空运动特性，也为下一步进行时空卷绕（Warping）打下基础。在DSA图像序列中的每个图像上，以所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格。

在采用常规方法得到的备选控制点的基础上，在DSA图像序列中的每个图像上，以所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格；从第一幅盈片图像开始，结合前一幅图像中控制点的位置来确定本幅图像的Warping参数，Warping采用多项式Warping方式；以此类推，直到完成对最后一幅图像的Warping，即以时空卷绕的方式逐幅图像地对盈片图像进行非线性校正，逐步达到与蒙片的精细配准；然后进行数字减影以消除运动伪影。

本发明克服了传统不能充分利用DSA图像序列中大部分造影剂浓度没有达到峰值的盈片图像的缺点，可以进行更加精细的图像配准。

由于采用的控制点不是仅仅分布在在蒙片或减影盈片上,而是贯穿在所有的蒙盈片序列上，更便于分析DSA图像的时空运动特性，也为下一步进行时空卷绕打下基础。在DSA图像序列中的每个图像上，以所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格，本发明具有充分利用时空信息，通过时空卷绕处理，逐幅图像地对盈片图像进行非线性校正，逐步提高与蒙片的精细配准，数字减影后有效地去除运动伪影的技术进步。

本发明的应用能够有效地去除运动伪影，改善DSA图像质量，进而提高医生的诊断准确性，提高工作效率。

附图说明

下面结合附图对本发明做进一步详细的说明。

图1：减影图像中出现的伪影示意图。

图2：DSA蒙片与盈片的部分序列示意图。

图3：采用边缘检测算法选取控制点示意图。

图4：时空卷绕的示意图。

图5：时空体（Spatio-temporal Volume）的结构示意图。

图6：DSA图像的Delaunay三角剖分示意图。

图7：图像卷绕的示意图。 

图8：基于时空切片的DSA运动伪影消除方法的流程图。

图9：本发明的伪影消除与DSA设备减影效果对比图。

图10：仿射变换示意图。

图11：本发明的系统示意图。

具体实施方式

消除伪影的主要手段就是要找到背景部分像素的运动规律，实现精确的图像配准并减影。虽然DSA图像序列中大部分是造影剂浓度没有达到峰值的盈片图像，如图1所示。但这些图像中只是血管部分没有被完全突出出来（变化中），其它部分如骨骼、肌肉的灰度值并没有明显变化（部分像素的位置会发生变化）。我们完全可以利用DSA图像序列中的时空相关信息，研究位于背景部分的控制点（Control Points）的运动规律，再通过图像卷绕实现图像配准。

参见图8所示，本发明的具体实施步骤如下：DSA图像资料由河南大学附属医院淮河医院（三级甲等医院）提供，采用的DSA设备型号为美国GE公司LCV Plus DSA系统。

步骤1、读入DSA图像序列（DICOM格式）并进行预处理。预处理包括灰度拉伸等简单的处理，如果原始图像质量比较好也可不做预处理。如图2所示。

步骤2、采用Canny算子找到图像的边缘点。由于没必要把所有用Canny 算子检测出的点都作为控制点，所以在Canny 边缘点计算出该点的梯度值后，根据阈值来判断该点是否作为备选控制点。本实例选取约100个备选控制点。如图3所示。

    Canny 边缘检测算子是John F. Canny于 1986 年开发出来的一个多级边缘检测算法，Canny边缘检测算法：step1:用高斯滤波器平滑图象；step2:用一阶偏导的有限差分来计算梯度的幅值和方向； step3:对梯度幅值进行非极大值抑制； step4:用双阈值算法检测和连接边缘。为了满足这些要求 Canny 使用了变分法，这是一种寻找满足特定功能的函数的方法。最优检测使用四个指数函数项的和表示，但是它非常近似于高斯函数的一阶导数。

当噪声点对边缘检测有较大的影响，效果更好的边缘检测器是高斯-拉普拉斯(Lo G)算子。它把高斯平滑滤波器和拉普拉斯锐化滤波器结合起来，先平滑掉噪声，再进行边缘检测，所以效果更好。

步骤3、利用DSA图像序列构建DSA时空体，将该组DSA二维图像以时间T为轴按拍摄顺序间隔排列成平行六面体空间。

步骤4、对时空体做若干个垂直方向的时空切片（Spatio-Temporal-Slice）。DSA时空体中运动点目标下一时刻的运动方向可归纳为：与时空切片平行；与时空切片垂直；与时空切片成夹角三种可能性。如图5所示。

当目标控制点沿时空切片方向运动时，控制点将被多个时空切片依次捕获，呈现在时空切片图像中为一条斜线。当目标控制点沿水平方向运动时，如果控制点运动速度小于时空切片间隔，则控制点将被同一个时空切片多次捕获，在时空切片图中表现为一段有限长线段。当目标控制点运动方向与时空切片成夹角时，它可以分解成平行方向和垂直方向，它的运动轨迹为上述两个方向轨迹的合成。

步骤5、将各时空切片对应的各控制点连线，得到各个控制点的运动轨迹。

步骤6，选择运动轨迹变化最大的10个反映DSA图像背景结构运动规律的控制点；本发明选择的控制点不是仅仅分布在在蒙片或减影盈片上，而是贯穿在所有的蒙盈片序列上，更好地记录了DSA图像的时空运动特性。如图7所示。

步骤7、用三角剖分的方法（Delaunay Triangulation Method）。将所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格；如图6所示。在DSA图像序列中的每个图像上，以所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格。三角剖分是代数拓扑学里最基本的研究方法，是极为重要的一项预处理技术。通过Delaunay三角剖分，可以将复杂的多边形转化相对比较简单的三角形集合。在DSA伪影区域会出现大量的三角形，这些区域正是要重点控制的区域。

步骤8、剖分完成后，采用仿射变换进行初步几何校正。如图10所示。仿射变换是空间直角坐标变换的一种，它是一种二维坐标到二维坐标之间的线性变换，保持二维图形的“平直线”和“平行性”，其可以通过一系列的原子变换的复合来实现，包括平移(Translation)、缩放（Scale）、翻转（Flip）、旋转（Rotation）和剪切(Shear)。几何上，两个向量空间之间的一个仿射变换或者仿射映射由一个线性变换接上一个平移组成。在有限维的情况，每个仿射变换可以由一个矩阵A和一个向量b给出，它可以写作A和一个附加的列b。一个仿射变换对应于一个矩阵和一个向量的乘法，而仿射变换的复合对应于普通的矩阵乘法，只要加入一个额外的行到矩阵的底下，这一行全部是0除了最右边是一个1，而列向量的底下要加上一个1。

步骤9、用Warping非线性变换法参考前一幅图像中控制点的变形位置，逐幅确定图像非控制点的位移；如图4所示。从第一幅盈片图像开始，结合前一幅图像中控制点的位置来确定本幅图像的Warping参数。Warping采用多项式Warping方式，多项式的项数与控制点数相同。本实例中控制点数为10，多项式如下：

解线性方程组，得系数后矩阵求逆，然后按控制点逐点Warping。以此类推，直到完成对最后一幅图像的Warping。由于综合考虑了DSA时空体的整体和局部时空信息，以时空Warping的方式能逐幅图像地对盈片图像进行非线性校正，逐步达到与蒙片的精细配准。

步骤10、本实例中配准中所选取相似性测度为互信息（Mutual Information）准则。本发明利用时空信息进行高精度DSA配准，计算量较大，所以采用Powell优化算法来优化相似度准则函数求极值的过程，使得有效地克服互信息函数的局部极值,大大地提高了配准精度,达到亚像素级。 Powell优化算法的配准速度快、精度高，鲁棒性较好。

步骤11、在DSA图像配准中，即使亚像素（或子像素）的错误配准也会在减影图像上产生大量的伪影。为配合精确的时空卷绕，选择较为复杂的三次样条插值算法进行亚像素级的配准。三次样条插值（简称Spline插值）是通过一系列形值点的一条光滑曲线，数学上通过求解三弯矩方程组得出曲线函数组的过程。

步骤12、逐幅对图像的灰度失真进行回朔校正；采用对数减影方法进行减影，即在进行减影运算前先将盈片和蒙片所有点的灰度值取反，分别做对数运算，获得的值作为当前点的灰度值，然后盈片减蒙片获得减影结果图像。

步骤13，对数减影；用对数减影法进行减影，获得消除伪影的减影图像。参见图9所示，左侧的两幅图像是DSA设备减影的效果，右侧的两幅图像是采用本发明方法的进行DSA运动伪影的消除后的效果。

参见图11所示，一种数字减影血管造影运动伪影消除系统，包括DSA设备、实现运动伪影消除的DSA软件工作站、DSA胶片打印机、图像显示系统、外部存储设备、网络传输系统，它的DSA软件工作站包括以下部分：输入接口：用于接收DSA设备输出的DICOM图像。

存储器：用于存储操作系统、图像处理软件、图像分析数据和图像资料。

数字减影部分：采用本发明的方法进行数字减影。

中央处理单元：用于读入DSA图像，并将DSA图像进行预处理；采用边缘检测算法找到每幅DSA图像中血管影像的边缘和边缘点，计算出边缘点的梯度值，根据阈值来判断该点是否作为备选控制点，并选取至少10个备选控制点；用于构建DSA时空体；将该组DSA二维图像以时间T为轴按拍摄顺序间隔排列成平行六面体空间；用于进行时空切片；用切片指令对DSA时空体做垂直方向的时空切片；用于将各时空切片对应的各控制点连线，得到各个控制点的运动轨迹；用于分析DSA像素的时空运动特性；选择运动轨迹变化最大的10个反映DSA图像背景结构运动规律的控制点；用三角剖分的方法将所有的控制点作为三角形的顶点，把图像剖分成不均匀的三角形网格；用于采用仿射变换法进行初步几何校正；用于作时空卷绕；用Warping非线性变换法参考前一幅图像中控制点的变形位置，逐幅确定图像非控制点的位移；采用Powell优化算法来优化相似度准则函数求极值的过程；选择三次样条插值算法进行亚像素级的配准；用于逐幅对图像的灰度失真进行回朔校正；用对数减影法进行减影，获得消除伪影的减影图像。

胶片打印控制部分：控制DSA图像的胶片打印。

图像处理部分：对DICOM图像进行适当的处理。

网络传输控制：对DSA图像在PACS系统中的网络传输进行控制。

报告系统：自动生成DSA图文检查报告。

显示器：用于显示DSA图像。

输出接口：用于接打印机。

还包括用于与远程计算机或互联网连接的通讯接口。

参见图4所示，图像卷绕是一种非线性的方法，是通过指定一系列控制点的位移来定义空间变换的图像变形处理。非控制点的位移根据控制点进行插值来确定。图像卷绕一般是在两维图像上进行，但是在DSA时空体上的控制点具有空间和时间上的相关性，在对某一盈片图像进行图像卷绕时，如果我们还考虑到其它图像对本次图像卷绕的影响，就可以使得图像卷绕的效果更加精细。因此本发明中把卷绕推广为时空卷绕（Spatio-Temporal Warping），即在对一幅盈片进行图像卷绕时，要结合前一幅图像中相关点的位置关系。

参见图5所示，本发明采用了时空体（Spatio-Temporal Volume）结构来描述DSA图像序列。时空体定义了一个包括X，Y和T（时间维度）三维坐标轴的平行六面体空间，把空间信息（XY坐标）和时间信息（T坐标）整合在一个数据结构中。采用时空体的好处是分析图像序列特征时，把整个图像序列看作一个整体，而不是离散的帧，可以充分保证空间和时间的连续性，能更好地表现图像序列中结构对象的运动特征。

在DSA图像时空体的基础上，本发明进一步对时空体进行时空切片（Spatio-Temporal Slice）来分析DSA背景部分像素的运动规律。时空切片是由图像序列中相同位置的像素条带按时序组成的一幅图像，或者说是从图像序列的同一个位置提取出的一些相关像素组合而成的一幅二维图像。再形象一些来说，假如有一个（X, Y, T）三维DSA时空体，其中（X, Y）是图像维，T是时间维。我们沿平行于T轴的方向对DSA时空体进行切割，这样将得到一个剖面，称为时空切片，时空切片可以看作是包括时间维的一幅二维图像。时空切片平行于T轴，但相对于X、Y轴的位置，即切割的方向、位置，可以根据具体情况来确定。引入时空切片可以方便地获取DSA图像序列中背景部分运动的动态信息，并能减少计算复杂性。

为了精确配准图像，需要选择控制点，这一部分采用常规方法。DSA图像序列中的运动目标是在空间域和时间域上做连续的运动，因此它在DSA时空体（Spatio-Temporal-Volume）上的轨迹点具有空间和时间上的相关性，本发明利用时空切片对用常规方法得到的控制点进行运动估计，即分析运动像素的时空相关特征。通过对控制点运动轨迹的研究，选择一批最有代表性的控制点，既可以保证配准精度，也提高了DSA 图像配准的鲁棒性。

参见图6所示，本发明提出的基于时空切片的控制点的选择方法与其它方法的一个显著的区别是，我们选择的控制点不是仅仅分布在在蒙片或减影盈片上，而是贯穿在所有的蒙盈片序列上，更便于分析DSA图像的时空运动特性，也为下一步进行时空卷绕（Warping）打下基础。在DSA图像序列中的每个图像上，以所有的控制点作为三角形的顶点，把图像剖分成不均匀的Delaunay三角形网格。

参见图7所示，通过时空卷绕处理，逐幅图像地对盈片图像进行非线性校正，逐步达到与蒙片的精细配准。由于充分利用时空信息，提高了配准精度，数字减影后可以有效地去除运动伪影。

Claims (2)

1.一种数字减影血管造影运动伪影消除方法，其特征在于步骤如下：

步骤1、读取：将数字减影血管造影仪器与计算机连接，在计算机中安装操作系统和图像处理系统，用图像处理软件读入一组DSA图像，并将DSA图像进行预处理；

步骤2、选点：采用边缘检测算法找到每幅DSA图像中血管影像的边缘和边缘点，计算出边缘点的梯度值，根据阈值来判断该点是否作为备选控制点，并选取至少10个备选控制点；

步骤3、构建DSA时空体：将该组DSA二维图像以时间T为轴按拍摄顺序间隔排列成平行六面体空间；

步骤4、时空切片：用切片指令对DSA时空体做垂直方向的时空切片；

步骤5、连轨迹：用图像处理软件将各时空切片对应的各控制点连线，得到各个控制点的运动轨迹；

步骤6、分析DSA像素的时空运动特性；选择运动轨迹变化最大的10个反映DSA图像背景结构运动规律的控制点；

步骤7、三角剖分：用三角剖分的方法将所有的控制点作为三角形的顶点，把图像剖分成不均匀的三角形网格；

步骤8、仿射变换：采用仿射变换法进行初步几何校正；

步骤9、时空卷绕：用Warping非线性变换法参考前一幅图像中控制点的变形位置，逐幅确定图像非控制点的位移；

步骤10、优化：采用Powell优化算法来优化相似度准则函数求极值的过程；

步骤11、配准：选择三次样条插值算法进行亚像素级的配准；

步骤12、灰度校正：逐幅对图像的灰度失真进行回朔校正；

步骤13、对数减影：用对数减影法进行减影，获得消除伪影的减影图像。

2.根据权利要求1所述的一种数字减影血管造影运动伪影消除方法，其特征在于：所述边缘检测算法采用Canny算子或高斯-拉普拉斯算子。

Images