一种DSA锥束精确滤波反投影断层成像系统及成像方法与流程

本发明涉及dsa锥束精确滤波反投影断层成像领域，具体讲是一种适用于介入治疗中使用的dsa设备的x射线断层成像系统及成像方法。

背景技术：

介入治疗在dsa影像监视设备引导下，通过导管深入体内，对疾病进行诊断或直接执行治疗，它使临床疾病由创伤变为微创，减少患者的病痛和恢复时间。dsa在介入治疗中具有不可替代的重要诊疗辅助设备。它利用x射线实时显示人体某一部位的透视影像，方便医生动态观察导管、导丝和造影剂等的行程。然而，透视影像是物体三维间信息堆叠显示，很难精确定位器械及病变的空间解剖位置。在颅内介入手术中，需要及时观察颅内出血现象，这需要利用ct设备获取颅内断层图像。在大多数的医院里介入手术室不额外配置ct设备，病人需在ct室与介入手术室间往返推送，手术无菌条件难以保证，容易引起其它并发症，甚至威胁到患者的生命。

dsa设备的c型臂可作旋转运动或多轨迹运动，可实时采集x射线衰减后的信号，并且所发出的x射线在空间中形成锥束，符合ct精确重建的数据采集方式。利用dsa设备实现ct断层功能，降低介入治疗的造价，减少患者医疗费用，是一项具有广泛应用前景的重要技术。

利用dsa设备进行断层重建，最实用的扫描轨迹是单圆弧轨迹。目前有fdk类型和katsevich类型两类算法可以实现相应的重建。前者属于近似锥束重建算法，在锥角大于2°后出现了伪像。后者属于半精确的锥束重建算法，在锥角小于4°时，能取得良好的重建效果，该算法申请美国专利conebeamfilteredbackprojectionimagereconstructionmethodforshorttrajectories，其专利号为us7,203,272b2，然而，算法实现比较复杂，不利实际应用。katsevich发表论文“imagereconstructionforthecircleandlinetrajectory”是一种精确锥束断层重建方法，然而该算法需要借助螺旋扫描的pi线进行重建，使得在重建过程中，每个重建点需要单独区间判别，复杂了算法，不利重建速度的提高。

技术实现要素：

因此，为了解决上述不足，本发明在此提供一种dsa锥束精确滤波反投影断层成像系统及成像方法；是将精确锥束重建方法应用到dsa成像设备中，实现dsa成像设备具备断层成像功能。利用同一条投影射线上的重建点具有相同的构造因子以及与构造因子相关的参数在检测器的投影几何关系，将结构因子的计算映射到检测器平面上，获得构造因子在检测器上的分布规律。重建过程主要包括偏导、重排、滤波和反投影四个步骤。从而使得重建方法具有滤波反投影结果，易于加速。

本发明是这样实现的，构造一种dsa锥束精确滤波反投影断层成像系统，其特征在于：包括成像检测系统、控制系统和计算机系统；成像检测系统包括c型臂，c型臂的一端配有x射线源，用于产生锥束x射线；在c型臂另一端配置二维平面检测器，检测器由若干阵元组成，每个阵元可感光穿过躺在手术台上的病人的x射线能量衰减信号；安装在基座上的多关节机械臂可灵活运动，使c型臂作旋转运动或多轨迹运动；c型臂运动过程中，x射线源形成了多种圆弧或直线的扫描轨迹，检测器上的阵元可获取不同扫描位置的x射线衰减信号，该信号传输到控制系统的数据采集单元，转换为锥束投影数据；c型臂运动控制器，控制c型臂运动的速度、形式；x射线控制器，提供x射线的能量和时间信号；手术台控制器，控制手术台的平移和升降运动；计算机系统中的断层图像重建单元，接收来自于数据采集系统的x射线投影数据，进行断层成像重建；重建后的图像被输入到计算机中心，进一步存储到大容量存储单元；计算机中心接收操作控制台输入的参数或指令，允许操作者观察显示器重建图形和其它数据；观察者也可利用计算机中心为数据采集单元、c型臂运动控制器、x射线控制器以及手术台控制器提供控制指令和相应的参数。

一种dsa锥束精确滤波反投影断层成像方法，其特征在于；成像方法的步骤包括：成像检测系统中，机械臂可灵活运动，使安装在c型臂上的射源和检测器作旋转或平移运动，x射线源发射的x射线，检测器上阵元接受衰减信号；控制系统中，数据采集单元将阵元获取的x射线衰减信号转换为锥束投影数据；图像重建单元接收锥束投影数据，经过模块进行加权，加权后进入模块，计算关于扫描轨迹参数的偏导数；扫描轨迹若为圆弧轨迹，则进入模块计算关于检测器参数的偏导数，模块将偏导数据进行加权求和，模块对求和数据进行过端点投影重排和希尔伯特滤波；模块沿水平方向对求和数据进行希尔伯特滤波；滤波后的数据进入模块进行构造因子处理及反投影运算；若扫描轨迹为直线轨迹，模块中数据输入到模块进行加权，再转入模块，加权数据沿着过切点直线重排和希尔伯特滤波；滤波后的数据也进入模块中进行构造因子处理及反投影；反投影后，获得物体的断层图像，输入计算机中心，存到大容量存储器，并输入到显示器进行显示。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的数据采集单元是把检测器阵元接收的x射线模拟信号经过偏置校正、增益校正、坏像素校正及log运算等预处理转化成锥束投影数据。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的射源和检测器作旋转是指射源和检测器同时绕c型臂中心轴旋转形成x射线锥束投影的圆弧扫描轨迹。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的射源和检测器作平移运动是指机械臂控制射源和检测器同时沿水平方向运动形成x射线锥束投影的直线扫描轨迹。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的端点投影是指射源沿圆弧轨迹运动时，射源与圆弧扫描轨迹两个端点分别连线同检测器平面的交点。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的过端点投影重排是指射源沿圆弧轨迹运动时，在检测器平面内原水平和竖直分布的求和数据分别沿过圆弧两端点投影的直线进行线性插值重排。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的切点直线是指射源沿直线轨迹运动时，圆弧轨迹在检测器平面上的投影是一条抛物线，过抛物线上的每一点所作切线即为切点直线。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述过切点直线重排是指射源沿直线轨迹运动时，在检测器平面内原水平和竖直分布的加权数据沿过与圆弧轨迹投影的相切直线进行线性插值重排。

根据本发明所述的dsa锥束精确滤波反投影断层成像方法，其特征在于；所述的构造因子是指射源在圆弧扫描轨迹上时，其值可能为1，0.5，-0.5，射源在直线扫描轨迹上时，其值可能为1，-1；

所述的构造因子处理及反投影是指希尔伯特滤波后的数据按照两种扫描轨迹构造因子的分布特点与其值相乘，重建点按照投影关系检索检测器平面的构造因子处理后的滤波数据获得某一投影视角的重建数值。

本发明具有如下优点：本发明提供一种dsa成像设备的精确断层成像系统及成像方法，该技术方案保持了滤波反投影的框架，具备锥束精确重建的条件，计算速度较快，可以有效改善锥角伪像，提高重建精度范围；避免了fdk类型算法重建精度不高的缺点，现有半精确锥束的复杂计算过程，提高重建算法的速度；同时也避免了精确锥束重建算法的pi线检索区间的判别，为dsa设备的断层成像提供方法支撑。同时，该方法可方便适应单圆弧扫描轨迹的锥束投影数据的半精度断层重建。

附图说明

图1为dsa断层成像系统框图；

其中：1成像检测系统；

101c型臂，102x射线源，103x射线，104检测器，105阵元，106手术台，107病人，108基座，109机械臂；

2控制系统；

201数据采集单元，202c型臂运动控制器，203x射线控制器，204手术台控制器

3计算机系统；

301图像重建单元，302计算机中心，303存储单元，304操作控制台，305显示器；

图2为正交圆弧轨迹锥束扫描几何示意图；

图3锥束扫描几何参数在检测器上的投影示意图；

图4为圆弧轨迹构造因子的确定示意图；((a)与轨迹相切的临界面；(b)过圆弧端点的临界面)；

图5为圆弧轨迹检测器区域的划分示意图((a)端点投影在检测器外侧且滤波线最大夹角小于π/2；(b)端点投影在检测器外侧且滤波线最大夹角大于π/2；(c)端点投影在检测器内)；

图6为直线轨迹构造因子的确定示意图；

图7为本发明的精确滤波反投影重建方法流程图。

具体实施方式

下面将结合附图1-图7对本发明进行详细说明，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1示出了dsa断层成像系统框图，成像检测系统1包括c型臂101，c型臂101的一端配有x射线源102，用于产生锥束x射线103；在c型臂101另一端配置二维平面检测器104，检测器104由若干阵元105组成，每个阵元105可感光穿过躺在手术台106上的病人107的x射线103能量衰减信号；安装在基座108上的多关节机械臂109可灵活运动，使c型臂101作旋转运动或多轨迹运动；c型臂101运动过程中，x射线源102形成了多种圆弧或直线的扫描轨迹，检测器104上的阵元105可获取不同扫描位置的x射线衰减信号，该信号传输到控制系统2的数据采集单元201，转换为锥束投影数据；c型臂运动控制器202，控制c型臂101运动的速度、形式；x射线控制器203，提供x射线的能量和时间信号；手术台控制器204，控制手术台106的平移和升降运动；计算机系统3中的断层图像重建单元301，接收来自于数据采集系统201的x射线投影数据，进行断层成像重建；重建后的图像被输入到计算机中心302，进一步存储到大容量存储单元303；计算机中心302接收操作控制台304输入的参数或指令，允许操作者观察显示器305重建图形和其它数据；观察者也可利用计算机中心302为数据采集单元201、c型臂运动控制器202、x射线控制器203以及手术台控制器204提供控制指令和相应的参数。

图2示出了dsa成像系统简化的圆弧加直线轨迹锥束投影扫描几何。o代表旋转中心，以其为原点建立笛卡尔xyz坐标系；投影采集过程中，检测器和射源绕z轴旋转，也可沿z轴作直线运动，检测器平面始终垂直于xoy平面，射源在检测平面的投影始终垂直于检测器中心；射源到检测器中心的距离记为d，i表示圆弧c与直线l构成的复合轨迹，设且c⊥l，l位于c的起始端y0；通过重建点任意的radon平面与扫描轨迹至少有一个交点，则表达了能够精确重建的源点扫描区间。采集系统引入由标识的局部坐标系；

由此，射源和重建点确定的单位矢量方向可表达为：

则过射源沿单位矢量方向的锥束投影可表示为：

上式中的s²表示三维空间的单位球。

过源点且包含单位矢量的平面记为其单位法向矢量记为而θ表示平面在垂直于矢量局部坐标系中的极角。

katsevich在ageneralschemeforconstructinginversionalgorithmforconebeamct论文中推导出通用的具有移不变滤波反投影结构的精确锥束重建算法，其重建公式可表示为：

式中，θ表示平面上任意投影射线矢量方向，γ则为矢量方向θ与之间的夹角；称之为构造因子，其表达式为：

构造因子中，表示平面的权重值。由于过重建点的平面与轨迹的交点可能不止一个，则权重值反映了平面对重建点的冗余贡献，需要满足下面关系：

在本发明中选择不等权重策略，即重建点的平面同时与圆弧和直线轨迹相交时，直线轨迹的重建点的权重贡献值为零，圆弧点的交点个数的倒数用以确定相应的权重值。

katsevich指出重建过程仅发生构造因子不为零的投影面，该面是与轨迹相切或者通过轨迹端点的平面。对于圆弧加直线轨迹上存在三类临界面：与轨迹相切以及分别通过圆弧起始和终止端点的平面。为此需判别这三类临界面上构造因子的大小来确定投影数据的滤波方向。

图3示出了当源点为在圆弧轨迹时，锥束扫描几何参数在检测器平面上的投影，上标符号^分别表示相应参数在检测器上的投影，而πt、πs和πe分别标识了过且与轨迹相切的平面、过和圆弧起始端点s的平面和过和圆弧终止端点e的平面与检测器的交线。显然，平面法线矢量在检测器上的投影与平面同检测器的交线相垂直。此外，根据构造因子的含义可知，同一条投影射线上的重建点具有相同的构造因子，并且这些重建点在检测器上的投影地址相同，因此，可以将构造因子存储在检测器对应的阵元上。

图4(a)示出了与圆弧轨迹相切的临界面同检测器平面的交线是沿水平方向的。当该平面绕投影射线分别沿顺时针与逆时针旋转后，与圆弧轨迹的交点个数均为2，与直线轨迹的交点个数为1，因此，而的符号值发生改变，由负变为正，由此可计算出该临界面上的结构因子值为1，并且相应的滤波线是沿水平方向。

图4(b)示出了过圆弧端点的临界面相关参数在检测器上的投影。从图中可看出，重建点在检测器上的投影在检测器轴上坐标正负的不同，临界面平面穿过圆弧端点前后与轨迹的交点个数变化不同；此外，投影点位置的不同，临界面上的的符号不同。因此，计算过圆弧端点的临界面上构造因子相对要复杂些。

过圆弧端点的临界面与检测器的交线决定了投影数据的滤波方向，亦称为滤波线。端点的投影一定在轴上，根据它沿轴的坐标的位置，将滤波线在检测器上的分布划分为三种情况，如图5所示。图5(a)示出了端点投影在检测器外，最外侧两条滤波线的夹角小于等于π/2；图5(b)也示出了端点投影仍然在检测器外，但最外侧两条滤波线夹角小于等于π且大于π/2；图5(c)示出了端点投影在检测器内，滤波线通过端点投影并沿整周分布。

检测器上的投影数据需要沿这些滤波线进行重采样，而且重采样后的投影数据被权重及滤波后，还需要反采样到检测器阵列单元上以方便反投影运算。为了减少采样过程中数据的缺失，对于斜率绝对值小于等于1的滤波线，投影数据应保持检测器水平阵列线的索引不变，而沿着轴进行重采样；对于斜率绝对值大于1的滤波线，投影数据应保持检测器垂直阵列线的索引不变，而沿着轴进行重采样。

经研究得出过圆弧端点处构造因子的绝对值均为0.5，而正负取值仅与过圆弧端点在检测器上的投影的斜率为-1的直线相关。对于圆弧的起始端点，则该直线的上方对应的构造因子为0.5，而在该直线下方对应的构造因子为-0.5；而对于圆弧的终止端点，则在直线的上方对应的构造因子为-0.5，而对应在该直线下方的构造因子为0.5。

图6示出了源点在直线轨迹时确定构造因子的示例图。过圆弧终止点的临界面，绕投影射线旋转过程中，与圆弧轨迹由一个或者两个交点。根据权重约定，权重值为零，使得构造因子中为零。因此对直线轨迹而言，仅需要判别与圆弧相切时的构造因子。经推断，圆弧轨迹在检测器平面的投影为一抛物线。图6给出当切点在抛物线左侧时构造因子确定的示意图。图6(a)示出当重建点在切点左侧时，临界面绕投影方向分别沿逆时针和顺时针旋转时，扫描轨迹的导数方向与平面法矢量的夹角为锐角。逆时针微动平面仅与直线轨迹相交，权重为1；顺时针微动平面与圆弧轨迹有两个交点，权重为0；因此，构造因子为-1。图6(b)示出当重建点在切点右侧时，临界面绕投影方向分别沿逆时针和顺时针旋转时，扫描轨迹的导数方向与平面法矢量的夹角为锐角。逆时针微动平面与圆弧轨迹有两个交点，权重为0；顺时针微动平面仅与直线轨迹相交，权重为1；因此，构造因子为1。进一步分析，当切点在抛物线顶点左侧时，构造因子与切点在右侧时完全相同。从而，可以得出，源点在直线轨迹运动时，滤波线是过源点与圆弧轨迹相切平面与检测器平面的交线，切点左侧的构造因子值为-1，右侧的构造因子值为1。

确定了构造因子不为零的radon平面以及相应的数值后，还需要推导适合平面检测器上的重建公式。利用偏导数链式法则，可推导出投影数据的偏导数为：

构造因子不为零的radon平面，其上相关的运算才会对重建点有所贡献。实际上，它与检测器平面的交线确定了投影数据的滤波方向。在检测器平面上，滤波线的斜率可以是0到正无穷。为了采样精确，对于斜率绝对值小于1的滤波线，沿水平方向采样，可以推导出偏导数据的加权为：

利用hilbert变换，进一步推导出投影数据的滤波公式：

其中，hh(·)表示hilbert变换滤波核。记g^f(λ,u,v)为投影数据偏导数的滤波函数，其表达式为：

在检测器平面上，由于通过圆弧端点的滤波线存在斜率绝对值大于1的情况。投影数据需要按轴进行重排，并也需要推导相应的公式。依据上述相同的推导方式，同样获得斜率绝对值大于1的滤波线对应的滤波函数公式g^f(λ,u,v)：

由此，获得dsa断层图像的精确重建方法。其表达式如下：

重建公式中，需要沿扫描轨迹获得在检测坐标u,v处的反投影值，u,v与重建点满足的投影映射关系：

此步骤为反投影运算。

重建公式表明本发明的方法由两部分组成。第一部分针对圆弧轨迹的重建公式：

其中，c0＝1表示投影数据沿水平方向进行滤波时构造因子的取值；c1和c2分别表示投影数据沿过圆弧起始和终止端点在检测器上的投影点和的直线方向进行滤波时构造因子的取值，依据在检测器上相应的分布规律取0.5或-0.5。

另一部分是针对直线轨迹的重建公式：

其中，c0＝-1表示过临界面与圆弧轨迹的切点在检测器上的投影点左侧的投影数据构造因子的取值；c1＝1表示过临界面与圆弧轨迹的切点在检测器上的投影点右侧的投影数据构造因子的取值。

图7示出了本发明的精确滤波反投影方法的流程图。数据采集单元201将检测器阵元105采集到的x射线103衰减信号转换为锥束投影数据，锥束投影数据经过诸如偏置校正、增益校正、坏像素校正及log运算等预处理。获取适合断层重建的投影数据，并按照检测器阵元布置的二维数组形式输出给图像重建单元301。在重建单元301中，锥束投影数据经过模块3011进行加权，然后进入模块3012先进行轨迹类型判别，计算关于扫描轨迹参数的偏导数。扫描轨迹若为圆弧轨迹，则进入模块3013计算关于检测器参数的偏导数，模块3014将模块3012和模块3013的计算结果进行加权求和，从而获得投影数据的偏导数；偏导数据分别进入到模块3016和模块3017中，在模块3016对偏导数据沿着过两个圆弧端点在检测器平面投影的滤波线进行重排，再进行希尔伯特滤波；在模块3017中偏导数据直接沿水平滤波方向进行希尔伯特滤波；滤波后的数据进入模块3019进行构造因子及反投影运算。若扫描轨迹为直线轨迹，模块3012获得的数据会进入模块3015进行加权，加权后的数据输入模块3018中，沿着与圆弧轨迹在检测器投影相切的直线进行重排，重排后的数据再进行希尔伯特滤波；滤波后的数据也进入模块3019中进行构造因子处理及反投影。反投影后，获得物体的断层图像。该图像会输入计算机中心302，再存到大容量存储器303中，并输入到显示器305进行显示和相关操作。

根据上述的分析，重建过程主要可划分为投影数据的偏导数、重排、滤波以及反投影四个步骤。此外还需说明两个重建模块中，沿圆弧轨迹的重建模块可以单独运行，即在轴向重建范围小的情况下，可以采用单个沿圆弧轨迹的重建模块，可以获得较精确的重建图像，并且重建操作简单。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。