CN102144927B - 基于运动补偿的ct设备和方法 - Google Patents

Abstract

本发明公开了一种基于运动补偿的CT设备和方法。本发明利用基于立体视觉的运动测量系统获得成像目标的运动参数，然后通过基于重建图像矩阵变换的技术实现运动补偿，消除运动伪影，获得清晰地二维/三维CT图像。本发明能够有效消除CT扫描过程中由于被扫描物体自身运动带来的运动伪影，并且能够方便地嵌入到现有的CT扫描装置中，实现无缝衔接。本发明能够有效提高现有CT的影像质量，特别是对某些运动无法自控的特殊目标CT成像，具有尤为重要意义，例如帕金森患者、幼儿、活体小老鼠等，以及超高分辨率人体成像等。

Description

基于运动补偿的CT设备和方法

技术领域

本发明的实施例涉及辐射成像，具体涉及一种基于运动补偿的CT设备和方法。

背景技术

自从1972年Hounsfield发明了第一台CT机，CT技术给医学诊断和工业无损检测带来了革命性的影响，CT已经成为医疗、生物、航空航天、国防等行业重要的检测手段之一。随着技术的进步，CT扫描模式和成像方法也在不断地改进，三维锥束CT已经成为研究和应用的主流。X射线锥束CT已经在医学临床、安全检查、无损检测等领域得到了广泛的应用，特别是在医学临床诊断中，螺旋CT已经成为不可或缺的检查手段之一。

1989年，螺旋CT开始投入医学临床应用，由于螺旋CT的巨大优势，使得它逐步替代了以前的断层CT，螺旋CT相对于断层CT的优势在于：螺旋CT可以连续不间断地采集投影数据，并通过专门设计的重建算法得到物体的三维体数据，使得CT扫描的时间大大缩短，提供了重建图像的Z轴分辨率，减少了运动伪迹。1991年，Elscint公司在单层螺旋CT基础上，首先推出了双层螺旋CT，从此揭开了多层螺旋CT(Multi-slice CT，MSCT)飞速发展的序幕。

MSCT与单层螺旋CT的主要区别在于单层螺旋CT的检测器是单排的，每次只能采集一层扇束投影数据，而MSCT的检测器是多排的，可以同时采集多层锥束投影数据。因此，MSCT和单层螺旋CT相比在性能上有了很大的提升，大大增加了X射线束的覆盖范围，有效地提高X射线的利用率，缩短了扫描时间，能够得到更高质量的三维重建图像。1998年，GE、Siemens、Toshiba、Philips公司推出了4层螺旋CT，揭开了各大医疗设备公司激烈竞争的序幕。2001年，GE公司率先推出了8层。2002年，GE、Siemens、Toshiba、Philips公司分别推出了16层螺旋CT。2005年，Toshiba公司推出了256层螺旋CT。2007年，在美国芝加哥的第93届北美放射学会议上Toshiba公布了其最新推出的320层螺旋CT产品。而Siemens公司则在2005年推出了首台双能量螺旋CT。目前最新的MSCT的扫描速度达到了每秒4周，已经广泛地被应用于人体三维成像、血管造影成像、心脏成像、脑灌注成像等领域。在MSCT技术上还发展起来了计算机辅助手术、虚拟内窥镜技术和辅助放射治疗等新技术。

尽管目前CT技术已经工业、安检、医疗等领域取得了巨大的成功，但由于工程应用条件的复杂性和多样性，对CT技术的进一步发展提出了更高的要求。在CT扫描过程中经常会遇到被扫描物体自身运动的情况，特别是在生物活体CT成像中。例如在医疗CT扫描过程，由于人体器官的固有运动，CT影像经常存在运动伪影，导致图像模糊、重影等，影响医生对疾病的准确诊断。人体器官在正常的新陈代谢过程中会存在多种生理性的运动，包括：心脏的跳动、肺部的呼吸运动、脉搏的跳动、肠胃的蠕动等等。这些器官的运动，除去呼吸可以通过摒息来控制，其他多数运动都很难使之暂时停止。另外，一些无法自我控制身体运动的特殊群体，例如帕金森患者、幼儿、没有自我控制能力的精神病人等，都可能在CT扫描过程中存在较大幅度的肢体运动。所有上述运动会不可避免地在MSCT图像中形成运动伪影，导致图像质量变差，甚至扫描失败。因此，研究抑制和消除运动伪影的技术对于CT成像具有重要现实意义。

在摒住呼吸的情况下，健康人体在10秒钟内其头部的平均运动位移约为0.35mm，在目前MSCT图像只有0.4mm分辨率的情况下，不会在图像中产生严重的运动伪影。而人体躯干部分的器官，由于受心跳、肠胃蠕动等运动的影响更大，运动会比头部更明显。特别是对病人而言，由于其身体的控制能力变弱，其运动幅度可能会超过1mm甚至更多，可能会导致严重的运动伪影，影响CT图像的质量。

减少运动伪影的根本办法是提高成像的速度，以减少扫描过程中器官自身运动带来的影响。但是，现有MSCT通过滑环技术实现X光机和探测器围绕人体快速旋转，以完成CT扫描。目前最快的扫描速度是0.25秒/圈，受离心力、材料强度等多方面限制，这个速度已经是业内的极限，在未来几年内很难有大的突破。除了MSCT，其他医疗CT、工业CT扫描速度更是远远低于0.25秒/圈。因此，单纯地依靠提高扫描速度是无法做到完全消除运动伪影的。

目前抑制运动伪影的另外一个技术是门控技术(Gating)，该技术在心脏成像中已经有相当应用，通过门控技术(Gating)把MSCT的扫描数据按照心脏相位进行分组，分别进行图像重建，然后利用图像融合技术得到最终的CT图像。这种门控技术带来的问题就是一些X射线投影数据没有得到有效地使用，大大增加了心脏成像的辐射剂量。而医疗照射已经成为全民最大的人为电离辐射来源，减少CT检查的X射线剂量是关系到全体社会公众及其后代健康的重大课题。

发明内容

本发明的目的是提出一种基于运动补偿的CT设备和方法，它利用立体视觉系统测量被扫描物体运动坐标，进而利用基于重建图像矩阵变换的运动补偿技术来消除CT图像中的运动伪影，获得清晰的CT图像。

在本发明的一个方面，提出了一种基于运动补偿的CT设备，包括：CT扫描设备，用于对被扫描物体进行CT扫描；立体视觉测量系统，用于与所述CT扫描设备的扫描同步地测量被扫描物体的运动参数；成像装置，基于所述运动参数，通过对重建图像矩阵进行变换来在CT成像过程中对被扫描物体的运动进行补偿。

根据本发明的实施例，所述立体视觉测量系统包括第一摄像机和第二摄像机，与所述CT扫描同步地采集被扫描物体的2D/3D图像，其中所述成像装置基于采集的2D/3D图像来计算被扫描物体的运动参数。

根据本发明的实施例，所述第一摄像机和第二摄相机与CT扫描设备之间的同步是基于CT扫描设备的探测器的采集时序的。

根据本发明的实施例，当CT扫描设备使用的探测器的采集速度高于所述第一摄像机和第二摄像机的采集速度时，将同步采集脉冲触发信号周期调整到探测器采集周期的整数倍，作为所述第一摄像机和第二摄像机的采集触发。

根据本发明的实施例，在进行CT重建时，通过插值得到和CT投影一一对应的被扫描物体的三维坐标以及运动参数。

根据本发明的实施例，当CT扫描设备使用的探测器的采集速度低于或等于所述第一摄像机和第二摄像机时，将探测器的采集信号当作所述第一摄像机和第二摄像机同步触发信号进行图像采集。

根据本发明的实施例，所述成像装置在每次反投影之前将重建图像矩阵按照测量到的运动参数做相应的矩阵变换，使得图像矩阵和该时刻物体形态一致。

根据本发明的实施例，所述运动参数包括被扫描物体的至少一部分的旋转和平移。

根据本发明的实施例，所述运动补偿是在FBP和ART算法中实现的。

在本发明的另一方面，提出了一种用于CT设备的方法，包括步骤：利用CT扫描设备对被扫描物体进行CT扫描；与所述CT扫描设备的扫描同步地测量被扫描物体的运动参数；基于所述运动参数，通过对重建图像矩阵进行变换来在CT成像过程中对被扫描物体的运动进行补偿。

本发明实施例的设备可以作为独立模块嵌入到已有或新开发的CT设备中，利用本发明实施例的设备和方法可获得消除运动伪影的CT图像，能够有效消除在CT扫描过程中由于物体运动而产生的运动伪影。

附图说明

从下面结合附图的详细描述中，本发明的上述特征和优点将更明显，其中：

图1示出了根据本发明实施例的基于运动补偿的CT设备的示意性框图；

图2示出了双目立体视觉系统中两台摄像机的相对位置关系和各坐标系示意图；

图3是描述使用立体视觉系统进行三维空间坐标测量的流程的示意性流程图；

图4示出了控制X射线探测器和摄像机同步采集的时序图，其中：

(a)表示当CT系统使用的探测器采集速度高于摄像机时的时序图；

(b)表示当CT系统使用的探测器采集速度摄像机相等时的时序图；

(c)表示当CT系统使用的探测器采集速度低于摄像机时的时序图；

图5示出了根据本发明实施例的图像重建方法的流程图；

图6示出了测试本发明实施例的设备和方法所用的圆盘模型，其中：

(a)表示初始时刻模型；

(b)表示一段时间以后变化的模型；以及

图7示出了利用根据本发明实施例的设备对运动物体进行CT重建的仿真结果，其中：

(a)传统FBP重建结果；

(c)使用本专利的运动补偿FBP重建结果；

(b)和(d)分别是重建结果中的水平和竖直中心剖面线；

(e)传统ART重建结果；

(g)使用本专利的运动补偿ART重建结果；

(f)和(h)分别是重建结果中的水平和竖直中心剖面线。

具体实施方式

下面，参考附图详细说明本发明的优选实施方式。在附图中，虽然示于不同的附图中，但相同的附图标记用于表示相同的或相似的组件。为了清楚和简明，包含在这里的已知的功能和结构的详细描述将被省略，否则它们将使本发明的主题不清楚。

根据本发明的实施例，提出了一种基于运动补偿的CT成像设备和方法。图1示出了根据本发明实施例的基于运动补偿的CT设备的示意性框图。

如图1所示，根据本发明实施例的CT设备配备了立体视觉系统。在CT扫描控制模块30的控制下，在CT扫描过程中由图像采集卡11同步采集被扫描物体的运动，例如通过摄像机A 20和摄像机B 21并配合存储在工控机10的存储器(未示出)的算法来计算物体的运动参数。进而，基于计算的运动参数来在图像重建过程中对运动进行补偿，消除运定伪影。

如图1所示，根据本发明实施例的设备还具备进行CT扫描所需的X光机50、探测器40和滑环60。在CT扫描过程中，设置在X光机50和探测器40之间的被扫描物体会发生一些部位例如头部的运动。

根据本发明的实施例，通过一套立体视觉系统，在CT扫描过程中同步采集被扫描物体的2D光学图像，通过基于双目/三目的快速三维坐标标定算法，记录每个CT投影对应的被扫描物体的三维坐标，然后利用测量的被扫描物体运动参数，在CT重建过程中通过重建图像矩阵变换技术进行运动补偿，进而消除最终重建2D/3D图像中的运动伪影，获得高精度的CT影像。

立体视觉是利用计算机通过多幅2D图像重建三维场景的方法，是计算机视觉研究领域的重要研究目标，尤其是在近二十年来取得了大量的研究成果。立体视觉主要分为被动和主动两种立体测量技术，前者中景物的照明由物体周围光线提供，后者则使用专门的光源装置来提供目标物体周围的照明。应用较多的是基于摄像机的立体视觉测量技术，根据使用到的摄像机数目的不同，又分为单目法、双目法、多目法。后面具体讲到的是利用两台摄像机进行立体视觉测量的双目技术，是一种基于双视角的三维坐标快速定位方法。但利用立体视觉技术进行运动三维坐标的测量，并不仅仅限于双目立体视觉。本领域的普通技术人员应该意识到，也可以采用其他的立体视觉系统来与CT扫描过程同步地采集被扫描物体的三维坐标和运动参数，而不是必须要使用上述的双目/三目三维坐标标定算法。

双目立体视觉是计算机视觉的一个重要分支，即由不同位置的两台摄像机拍摄同一幅运动场景，通过计算空间点在两幅图像中的视差，获得该点的三维坐标值。80年代美国麻省理工大学的Marr提出了一种视觉计算理论并应用在双目匹配上，使两张有视差的平面图产生了有深度的立体图形，奠定了双目立体视觉发展的理论基础。相比其他体视方法，双目视觉直接模拟人类双眼处理景物的方式，可靠简便，在许多领域具有重要的应用价值，如工业非接触式测量、机器人导航与航测、三维测量学、虚拟现实等。

如图3所示，双目立体视觉技术的实现可分为以下步骤：图像获取步骤S10、摄像机标定步骤S11、特征提取步骤S12、图像匹配步骤S13和三维坐标计算步骤S14。

根据本发明的实施例，针对CT扫描的环境和目标物体的运动特点，提出了用于CT成像的双目立体视觉测量系统和方法，然后利用后面提出的基于重建图像矩阵变换的技术来实现运动物体的准确CT成像，消除运动伪影。

根据本发明的实施例，使用一套立体视觉测量系统，在X射线源进行CT扫描的同时，对被扫描物体进行三维坐标的测量，获得该物体的运动参数，然后利用提出的运动补偿技术消除运动伪影，获得清晰的CT图像。本发明实施例的立体视觉测量系统是基于双目立体视觉的原理，从两个不同位置的视点观察同一景物，以获取CT扫描物体在不同视角下的感知图像，通过几何成像原理，计算图像像素间的位置偏差(视差)，来获取物体的三维坐标信息。如图1所示，立体视觉测量系统主要包括以下软硬件组成部分：

两台参数性能相同的摄像机20和21，图像采集卡11、数据传输线等辅助设备。在两个不同视角同时获取被扫描物体的两幅图像，然后计算出被扫描物体的三维坐标。如图1所示，两台摄像机A、B分开一定距离固定布置，一般可采用CCD摄像机，分别通过视频采集卡11与计算机连接，将数据传输到计算机10，在计算机10中通过执行软件程序来进行基于运动补偿的图像重建过程。一般立体视觉系统只要求两台摄像机同步采集，而本发明实施例中，不仅仅要求这两台摄像机20和21要同步采集，而且其采集必须保证能够与CT扫描系统中的探测器40采集同步，以便后面计算得到的被扫描物体三维坐标能够和采集到的CT投影数据一一对应，方便进行CT运动补偿重建。因此，为了实现两台摄像机和CT扫描系统的同步采集，可以从CT扫描控制模块30中引出与CT探测器40同步的采集脉冲信号，作为外触发源实现两台摄像机20和21的同步采集。

为方便物体的三维坐标计算，一般固定两台摄像机A、B的相互位置关系如附图2。首先定义被扫描物体所在的三维空间内的直角坐标系为xyz，其次分别定义x_Ay_A和x_By_B为摄像机A和B的成像平面，由于两台摄像机性能参数完全相同，可以使得它们的光轴z_A、z_B平行，同时使得成像平面x_Ay_A、x_By_B在同一个平面内，并且它们的水平轴x_A、x_B与三维空间坐标系的x轴完全重合。P_A、P_B分别表示物体点P在摄像机A、B成像平面内的成像点，而O_A、O_B分别表示摄像机A、B的焦点，在x轴上相对坐标原点对称，焦点到成像平面的距离，即焦距f由像机自身的硬件决定，为已知条件。两台像机的距离为D，可通过测量得到。

空间物体点P的三维坐标记为(x_p，y_p，z_p)，在某一时刻，两台摄像机A、B同时拍摄到点P的图像，在其各自成像平面内分别记为P_A＝(x_pA，y_pA)、P_B＝(x_pB，y_pB)。由于A、B的图像在同一平面内，且x轴重合，因此，点P在A、B上的图像点的y坐标相同，即y_pA＝y_pB。由三角几何关系可以得到：

$\left\{\begin{array}{c}{x}_{\mathrm{pA}}=f\·\frac{x_p+D/2}{z_p}\\ x_{\mathrm{pB}}=f\·\frac{x_p-D/2}{z_p}\\ y_{\mathrm{pA}}=y_{\mathrm{pB}}=f\·\frac{y_p}{z_p}\end{array}\right.\·\·\·\·\·\·\left(1\right)$

求解上式可以计算出物体点P的三维空间坐标为：

$\left\{\begin{array}{c}{x}_p=\frac{D}{2}\·\frac{x_{\mathrm{pA}}+x_{\mathrm{pB}}}{x_{\mathrm{pA}}-x_{\mathrm{pB}}}\\ y_p=\frac{y_{\mathrm{pA}}\·D}{x_{\mathrm{pA}}-x_{\mathrm{pB}}}\\ z_p=\frac{f\·D}{x_{\mathrm{pA}}-x_{\mathrm{pB}}}\end{array}\right.\·\·\·\·\·\·\left(2\right)$

由此，通过左右两台同步摄像机20和21上拍摄的两幅图像就能够计算出物体任意点的三维空间坐标，并且上述方法是完全点对点的运算，只要在两台摄像机拍摄的图像中找到P点的匹配点，就可以利用上述算法计算出其对应的三维坐标。基于双目的立体视觉测量流程和三维坐标计算由下面方法完成。

如图3所示，在步骤S10，图像同步获取是双摄像机进行三维坐标测量的基础，当摄像机按照图1布置完毕，通过CT扫描控制模块30输出的同步采集脉冲信号，两台摄像机20和21同时采集目标物体的2D图像，通过图像采集卡11传输到计算机10，同时记录图像数据采集的时间。

计算机立体视觉成像一般采用双目图像，有时也采用三目或多目成像技术，考虑到CT扫描的物体一般比较固定，例如人体、包裹行李等，因此，基于双摄像机的立体视觉系统一般能够满足实际应用需求。

在步骤S11，执行摄像机标定。摄像机标定的目的是为了准确确定摄像机的位置、属性参数，并建立成像模型，以便确定空间坐标系中物体点同它所在像平面上像点的对应关系，即确定变换矩阵。一般情况下，该对应关系由摄像机几何模型决定，虽然有些参数是可以通过测量和阅读摄像机参数说明得到的，例如两台摄像机的距离、摄像机焦距等，但为了更准确的确定这些参数，需要我们通过实验对这些参数信息进行测量和计算，即摄像机的定标。因此，摄像机定标的过程就是确定摄像机内部几何和光学特性(内部参数)，以及和相对世界坐标系的摄像机坐标系三维位置和方向(外部参数)。为了更准确地确定这些参数，一般地，在实验中设置10个左右标记点，这些点在世界坐标系中的三维坐标已知，通过测量这些10个标记点在像机成像坐标系中的坐标，来确定上述的全部参数，完成摄像机标定。

在步骤S12，执行图像特征提取。图像特征提取是为了得到匹配赖以进行的图像特征，主要的匹配特征有点状特征、线状特征和区域特征等。提取何种特征需要根据实际的被扫描物体来确定，例如：金属工件、精密器件等一般选择线特征，行李包裹等选择区域特征。有些物体不具备明显的特征，比如医疗CT对人体腹部成像，一般需要在人体表面粘贴辅助标记，例如黑白相间的“田”字图案等。需要说明的是：由于摄像机采集的2D图像往往包含了各种各样的随机噪声和畸变，因此在进行特征提取前需要对图像进行预处理，突出有用信息、抑制无用信息，提高图像清晰度，改善图像质量；使图像便于计算机处理，进行各种特征分析和提取。

在步骤S13，执行立体匹配。立体匹配就是要将两台摄像机同步拍摄的两幅图像中的图像特征对应起来，同一物体在不同视点下的图像会有较大不同，根据前面确定的匹配特征，在两幅图像间建立正确的匹配关系。

在步骤S14，计算三维坐标。在确定好匹配点后，就可以利用公式(2)确定特征点的三维坐标信息，进而获得物体各特征点的三维运动参数。

如上所述，摄像机20和21的图像采集与CT扫描系统的扫描要同步。同步采集保证两台摄像机能够和CT系统的探测器实现同步采集，由于摄像机和X射线探测器的采集速度往往是不同的，因此，必须专门进行时序控制设计。不同CT系统使用的探测器采集速度往往存在较大差异，比如目前医院使用的主流MSCT产品，其探测器采集速度能够达到1000Hz左右，而某些锥束CT系统使用的平板探测器的采集速度只有1～30Hz左右。而摄像机的图像拍摄速度一般在25或者30Hz，高速摄像机能够达到超过100Hz甚至更高的速度。在进行同步采集时需要根据采集频率的不同，设计不同的同步时序。

如图4所示，有下面三种选择方式：(1)当CT系统使用的探测器采集速度高于摄像机时可采用图4中的(a)时序，将同步采集脉冲触发信号周期调整到X射线探测器40采集周期的整数倍，作为摄像机20和21的采集触发，此时摄像机20和21拍摄到的物体图像要少于CT扫描的投影个数，在后面进行运动补偿CT重建时，需要进行插值，得到和CT投影一一对应的被三维物体三维坐标。(2)当CT系统使用的探测器40采集速度摄像机相等时，可采用图4的(b)时序，此时CT投影和摄像机图像一一对应。(3)当CT系统使用的探测器40采集速度低于摄像机20和21时可采用图4的(c)时序，将X射线探测器40的采集信号当作摄像机20和21同步触发信号进行图像采集，此时摄像机图像和CT投影也是一一对应的。

在通过立体视觉测量计算出物体在CT扫描过程中的三维坐标，即可利用后面描述的基于图像矩阵变换的CT重建算法对运动物体进行重建，获得消除了运动伪影的清晰CT图像。

在知道了物体在CT扫描过程中的运动情况就可以使用本发明实施例的重建方法进行CT重建。需要说明的是：本发明实施例的方法对于物体做何种运动是有限制的，要求被扫描物体所做的运动可以分解为任意小单元的刚性运动，即被扫描物体的运动可以看作是组成该物体的小单元做彼此无关的刚性运动。实际上现实生活中能够遇到的绝大多数物体的运动都可以看作是该类运动，例如金属工件的震动、心脏的跳动、人体头部的轻微摆动等等。因此，本实施例提出的CT图像重建方法能够满足实际CT应用所能够碰到的绝大多数的运动伪影消除需求。当被扫描物体作为一个整体做刚性运动时，可以看作是物体所有小单元一起运动的特殊形式，因此，后面以物体所有小单元各自做刚性运动的一般情况对本发明的实施例的方法进行描述。

另外需要说明的是：由于该重建方法与CT扫描的具体形式无关，即该方法可以用于平行束CT图像重建，也可用于扇束、锥束CT图像重建，因此，为了清楚说明，下面以平行束为例详细说明该方法。

首先，考虑二维平面内平行束投影可以用下面公式表示：

$p(\widehat{\mathrm{\θ}},s)={\∫}_{R^2}f\left(\stackrel{\→}{x}\right)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)d\stackrel{\→}{x}\·\·\·\·\·\·\left(3\right)$

其中，

是被扫描物体的函数表示，

是垂直于投影方向的单位方向向量，s是某条X射线投影到探测器上的投影位置。在CT图像重建中应用最广泛的FBP(滤波反投影)算法需要首先要对投影数据进行滤波，如下：

$g(\widehat{\mathrm{\θ}},s)={\∫}_{-s_0}^{s_0}p(\widehat{\mathrm{\θ}},s^{\′})\·h(s-s^{\′}){\mathrm{ds}}^{\′}\·\·\·\·\·\·\left(4\right)$

其中，h(s)是滤波函数，一般采用Ramp滤波器。然后将滤波后的投影数据进行反投影即可重建出

反投影公式如下：

$f\left(\stackrel{\→}{x}\right)={\∫}_{\mathrm{\θ}}g(\widehat{\mathrm{\θ}},s)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)d\widehat{\mathrm{\θ}}\·\·\·\·\·\·\left(5\right)$

公式(5)就是传统CT中使用最多的FBP重建算法。当物体在CT扫描过程中存在运动时，X射线探测器所记录的是运动中的物体投影数据，此时利用上述FBP算法重建得到的混合了CT扫描过程中物体的不同运动形态。因此，重建图像中混叠了运动伪影，导致图像质量下降。为了在重建过程中消除运动伪影，首先定义多体刚性运动模型如下：被扫描物体可以分成若干个小单元，每个单元各自做刚性运动，单个单元的刚性运动可以用如下公式表示：

刚性运动在上式中被分解为两个运动：旋转和平移。其中，

表示旋转运动中物体单元在t时刻转过的角度，表示平移运动中物体单元在平移运动向量。定义：

将上式代入式(6)中，可以得到：

$f_t\left(\stackrel{\→}{x}\right)=f\left({\mathrm{\Γ}}_t\left(\stackrel{\→}{x}\right)\right)\·\·\·\·\·\·\left(8\right)$

将所有的物体单元累加即可得到完整物体的运动表达式：

$f_t\left(\stackrel{\→}{x}\right)={\mathrm{\Σ}}_i{f}_{t,i}\left(\stackrel{\→}{x}\right)\·\·\·\·\·\·\left(9\right)$

其中，

表示第i个物体单元在t时刻的分布函数。

当物体存在多体运动时，其平行束投影相应的变为：

$p_t(\widehat{\mathrm{\θ}},s)={\∫}_{R^2}{\mathrm{\Σ}}_i{f}_{t,i}\left(\stackrel{\→}{x}\right)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)d\stackrel{\→}{x}$

$={\mathrm{\Σ}}_i{\∫}_{R^2}{f}_{t,i}\left(\stackrel{\→}{x}\right)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)d\stackrel{\→}{x}$

$={\mathrm{\Σ}}_i{p}_{t,i}(\widehat{\mathrm{\θ}},s)\·\·\·\·\·\·\left(10\right)$

即：多体运动情况下的物体投影等于所有物体单元分别投影之和。由于在CT扫描过程中物体存在运动，因此每个角度下的CT投影对应的物体形态都是不一样的，这样就无法直接利用传统的CT重建方法进行重建，否则重建的结果将是各个不同形态的物体的简单叠加。

本发明实施例提出了一种CT重建方法，在将各个不同时刻采集的投影数据进行反投影之前，需要将被扫描物体的重建图像矩阵(初始值为预先设置值)按照测量到的运动情况对应不同时刻做相应的矩阵变换，使得图像矩阵和该相应时刻物体形态一致，然后再进行反投影，该时刻数据反投影完成后，再将此时的图像矩阵变换到下一个时刻对应的物体形态，进行下一组投影数据的反投影操作，直到全部投影数据反投影完毕，再将此时的图像矩阵变换到t＝0时刻形态，即得到了最终消除了运动伪影的t＝0时刻物体形态的CT重建结果。由于在CT扫描过程中物体不同时刻的运动参数已经测量得到，因此可以利用重建得到的t＝0时刻重建结果通过图像矩阵变换得到任意时刻的物体重建图像。

根据本发明的实施例，也可以在每次反投影操作之前将更新后的重建图像矩阵变换到与初始时刻的形态。

图5示出了根据本发明实施例的图像重建过程的流程图。如图5所示，在步骤S20，将被扫描物体的重建图像矩阵的元素设置为全零或1，作为图像重建矩阵的初始值。

在步骤S21，取某一时刻(角度)的投影数据和该时刻的被扫描物体的三维运动参数。

在步骤S22，在计算机10中利用三维运动参数对上述重建图像矩阵进行矩阵变换，使得重建图像矩阵和该时刻的物体的形态一致；

在步骤S23，利用投影数据，通过执行反投影算法或者其他算法来对变换后的重建图像矩阵进行更新；

在步骤S24，将更新后的重建图像矩阵变换到初始时刻的形态；

在步骤S25，判断是否重建结束；

如果结束，则在步骤S26输出重建图像，例如显示在显示器上。如上所述，得到了初始时刻的重建图像，通过矩阵变换可以得到其他任意时刻的重建图像；

如果没有结束，则在步骤S27取另一时刻(角度)的投影数据和相应的三维运动参数，然后流程转入到步骤S22，重复执行上述的操作过程。

下面以FBP算法和ART(代数迭代)算法两类算法为例，详细说明根据本发明实施例重建过程的例子。

1.FBP类型的运动补偿算法：

考虑FBP重建公式(5)，由于各个角度投影数据对应的物体形态不同，因此必须将其统一才能对各个角度的反投影结果进行累加。公式(5)可以重新写为：

$f\left(\stackrel{\→}{x}\right)={\∫}_{\mathrm{\θ}}{f}_{\mathrm{FBP},\widehat{\mathrm{\θ}}}\left(\stackrel{\→}{x}\right)d\widehat{\mathrm{\θ}}\·\·\·\·\·\·\left(11\right)$

其中：

$f_{\mathrm{FBP},\widehat{\mathrm{\θ}}}\left(\stackrel{\→}{x}\right)=g(\widehat{\mathrm{\θ}},s)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)\·\·\·\·\·\·\left(12\right)$

对于多体运动物体而言，每个角度下的投影数据反投影对应的是不同时刻下物体的形态

$f_{\mathrm{FBP},t,\widehat{\mathrm{\θ}}}\left(\stackrel{\→}{x}\right)=g_t(\widehat{\mathrm{\θ}},s)\mathrm{\δ}(\stackrel{\→}{x}\·\widehat{\mathrm{\θ}}-s)\·\·\·\·\·\·\left(13\right)$

因此，(13)式结果不能直接代入到(11)中，考虑到通过前面的立体视觉系统已经计算得到每个投影角度时刻对应的物体三维坐标，因此，可以得到t时刻物体变换到t＝0时刻的刚性变换公式：

$f_{\mathrm{FBP},\widehat{\mathrm{\θ}}}\left(\stackrel{\→}{x}\right)=f_{\mathrm{FBP},t,\widehat{\mathrm{\θ}}}\left({\mathrm{\Γ}}_t^{-1}\left(\stackrel{\→}{x}\right)\right)\·\·\·\·\·\·\left(14\right)$

其中，Γ_t ^-1是刚性变换式(7)的逆变换。代入公式(11)就得到了最终的经过运动补偿之后的FBP重建公式：

$f\left(\stackrel{\→}{x}\right)={\∫}_{\mathrm{\θ}}{f}_{\mathrm{FBP},t=0,\widehat{\mathrm{\θ}}}\left(\stackrel{\→}{x}\right)d\widehat{\mathrm{\θ}}\·\·\·\·\·\·\left(15\right)$

2.ART类型的运动补偿算法：

ART算法是除了FBP算法之外经常在实际CT系统中应用的一种迭代算法，下面推导如何在ART算法中加入重建图像矩阵变换思想，得到一种新的能够用于运动补偿的ART算法。ART算法的基本思想是根据每一条X射线的投影值通过一定规则修正迭代图像，最终得到符合收敛条件的重建图像。此处，记第j条X射线的投影为：

            p_j＝H_jf        ……(16)

其中，H_j是这条X射线的投影矩阵，记录了该X射线对其经过的所有重建矩阵像素的投影贡献。该新型ART重建算法计算过程如下：

(1)设置迭代图像初值f₀ ^k＝0，一般为全零或全1，也就是将图像变换矩阵的元素设置为零或1；

(2)进行第k步迭代，对应X射线投影为第j个投影值，此过程分为3步：

(2.1)将f₀ ^k变换到j投影对应的t时刻形态：

$f_t^k\left(\stackrel{\→}{x}\right)=f_0^k\left({\mathrm{\Γ}}_t^{-1}\left(\stackrel{\→}{x}\right)\right)\·\·\·\·\·\·\left(17\right)$

(2.2)将f_t ^k代入下面公式计算得到f_t ^k+1

$\left\{\begin{array}{cc}{f}_t^{k+1}=f_t^k+\frac{p_j-H_j{f}_t^k}{{\left|\right|H_j\left|\right|}^2}{H}_j^T& \left|\right|H_j\left|\right|\≠0\\ f_t^{k+1}=f_t^k& \left|\right|H_j\left|\right|=0\end{array}\right.\·\·\·\·\·\·\left(18\right)$

(2.3)将f_t ^k+1变换到t＝0时刻图像矩阵：

$f_0^k\left(\stackrel{\→}{x}\right)=f_t^k\left({\mathrm{\Γ}}_t\left(\stackrel{\→}{x}\right)\right)\·\·\·\·\·\·\left(19\right)$

(3)k＝k+1，进行下一步迭代，直到重建图像f₀ ^k符合收敛条件。

需要特别说明的是：上述重建算法并不仅仅限于平行束CT，可以方便地扩展到扇束和锥束CT图像重建。下面是利用计算机仿真实验对上述CT图像重建算法进行的实验验证结果。

如附图6所示，构建了一个包含了3层共36个圆盘的模型，每层圆盘的旋转速度是不一样的，圆盘内线性衰减稀疏设为1，其他部分设为0。图6的(a)为模型初始时刻的图像，(b)为经过一段时间运动后的模型图像。

在实验中，三层圆盘的旋转速度从外至内分别为2ω₀、3ω₀和5ω₀。ω₀为X射线源和探测器的旋转扫描速度，在360度内采集了260度的投影，探测器单元个数为256，重建图像矩阵像素为256×256。图7示出了重建结果：

(a)是使用传统FBP算法，直接利用运动物体的投影数据进行重建的结果，可以发现各个圆盘完全被运动伪影掩盖了，无法进行分辨；

(c)是使用本专利的基于图像矩阵旋转的新型FBP重建结果，可以发现图像清晰地重建出了所有的圆盘，运动伪影被很好地补偿掉了；

(b，d)分别是(a，c)重建结果中的水平和竖直中心剖面线，可以发现新型运动补偿的FBP算法结果重建数值也是准确的；

(e)是使用传统ART算法，直接利用运动物体的投影数据进行重建的结果，所有圆盘也完全被运动伪影掩盖了，无法进行分辨；

(g)是使用本专利的基于图像矩阵旋转的新型ART重建结果，图像清晰地重建出了所有的圆盘，运动伪影被同样被很好地补偿掉了；

(f，h)分别是(e，g)重建结果中的水平和竖直中心剖面线，该新型运动补偿的ART结果重建数值同样也是准确的。

如上所述，本发明实施例提出的采用立体视觉系统的CT设备和方法是针对CT图像中的运动伪影而提出的一种解决方法。

当前，CT技术已经在很多技术领域得到了应用，CT应用经常面临的一个问题就是在CT扫描过程中物体存在或多或少的运动(特别是在生物活体成像中)，导致CT图像中存在运动伪影而降低图像质量。

本发明实施例的硬件部分能够方便地嵌入到现有的CT扫描设备中，进行系统功能升级。本发明实施例的设备和方法通过立体视觉系统计算得到的物体三维运动参数代入到提出的新型运动补偿FBP或ART算法后，能够有效抑制甚至消除运动伪影，提高CT影像的质量。特别是在某些运动无法自控的特殊目标群体，本发明实施例的设备和方法尤为重要，例如帕金森患者、幼儿、活体小老鼠等。另外，对于超高分辨率CT成像，即时轻微的运动都可能会造成严重的运动伪影，此时本发明实施例的设备也是一种非常好的选择。

上面的描述仅用于实现本发明的实施方式，本领域的技术人员应该理解，在不脱离本发明的范围的任何修改或局部替换，均应该属于本发明的权利要求来限定的范围，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (18)

1.一种基于运动补偿的CT设备，包括：

CT扫描设备，用于对被扫描物体进行CT扫描；

立体视觉测量系统，用于与所述CT扫描设备的扫描同步地测量被扫描物体的运动参数；

成像装置，基于所述运动参数，通过对重建图像矩阵进行变换来在CT成像过程中对被扫描物体的运动进行补偿。

2.如权利要求1所述的CT设备，其中所述立体视觉测量系统包括第一摄像机和第二摄像机，与所述CT扫描同步地采集被扫描物体的2D/3D图像，其中所述成像装置基于采集的2D/3D图像来计算被扫描物体的运动参数。

3.如权利要求2所述的CT设备，其中所述第一摄像机和第二摄相机与CT扫描设备之间的同步是基于CT扫描设备的探测器的采集时序的。

4.如权利要求3所述的CT设备，其中，当CT扫描设备使用的探测器的采集速度高于所述第一摄像机和第二摄像机的采集速度时，将同步采集脉冲触发信号周期调整到探测器采集周期的整数倍，作为所述第一摄像机和第二摄像机的采集触发。

5.如权利要求4所述的CT设备，其中，在进行CT重建时，通过插值得到和CT投影一一对应的被扫描物体的三维坐标以及运动参数。

6.如权利要求3所述的CT设备，其中当CT扫描设备使用的探测器的采集速度低于或等于所述第一摄像机和第二摄像机时，将探测器的采集信号当作所述第一摄像机和第二摄像机同步触发信号进行图像采集。

7.如权利要求1所述的CT设备，其中所述成像装置在每次反投影之前将重建图像矩阵按照测量到的运动参数做相应的矩阵变换，使得图像矩阵和该时刻物体形态一致。

8.如权利要求7所述的CT设备，其中所述运动参数包括被扫描物体的至少一部分的旋转和平移。

9.如权利要求7所述的CT设备，其中所述运动补偿是在FBP和ART算法中实现的。

10.一种用于CT设备的方法，包括步骤：

利用CT扫描设备对被扫描物体进行CT扫描；

利用立体视觉测量系统，与所述CT扫描设备的扫描同步地测量被扫描物体的运动参数；

基于所述运动参数，通过对重建图像矩阵进行变换来在CT成像过程中对被扫描物体的运动进行补偿。

11.如权利要求10所述的方法，其中所述立体视觉测量系统包括第一摄像机和第二摄像机，通过所述第一摄像机和第二摄像机，与所述CT扫描同步地采集被扫物体的2D/3D图像，其中基于采集的2D/3D图像来计算被扫描物体的运动参数。

12.如权利要求11所述的方法，其中所述第一摄像机和第二摄相机与CT扫描设备之间的同步是基于CT扫描设备的探测器的采集时序的。

13.如权利要求12所述的方法，其中，当CT扫描设备使用的探测器的采集速度高于所述第一摄像机和第二摄像机的采集速度时，将同步采集脉冲触发信号周期调整到探测器采集周期的整数倍，作为所述第一摄像机和第二摄像机的采集触发。

14.如权利要求13所述的方法，其中，在进行运动补偿CT重建时，通过插值得到和CT投影一一对应的被扫描物体的三维坐标以及运动参数。

15.如权利要求14所述的方法，其中当CT扫描设备使用的探测器的采集速度低于或等于所述第一摄像机和第二摄像机时，将探测器的采集信号当作所述第一摄像机和第二摄像机同步触发信号进行图像采集。

16.如权利要求10所述的方法，其中在每次反投影之前将重建图像矩阵按照测量到的运动参数做相应的矩阵变换，使得图像矩阵和该时刻物体形态一致。

17.如权利要求16所述的方法，其中所述运动参数包括被扫描物体的至少一部分的旋转和平移。

18.如权利要求16所述的方法，其中所述运动补偿是在FBP和ART算法中实现的。

Images