CN110974278B - 一种dsa锥束精确滤波反投影断层成像系统及成像方法 - Google Patents

一种dsa锥束精确滤波反投影断层成像系统及成像方法 Download PDF

Info

Publication number
CN110974278B
CN110974278B CN201911331253.6A CN201911331253A CN110974278B CN 110974278 B CN110974278 B CN 110974278B CN 201911331253 A CN201911331253 A CN 201911331253A CN 110974278 B CN110974278 B CN 110974278B
Authority
CN
China
Prior art keywords
detector
data
projection
cone beam
track
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911331253.6A
Other languages
English (en)
Other versions
CN110974278A (zh
Inventor
王瑜
李迅波
闫明明
曾毅星
蔡吴缺
魏宏才
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201911331253.6A priority Critical patent/CN110974278B/zh
Publication of CN110974278A publication Critical patent/CN110974278A/zh
Application granted granted Critical
Publication of CN110974278B publication Critical patent/CN110974278B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/504—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of blood vessels, e.g. by angiography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Pathology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Public Health (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Optics & Photonics (AREA)
  • General Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Vascular Medicine (AREA)
  • Dentistry (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pulmonology (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

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

Description

一种DSA锥束精确滤波反投影断层成像系统及成像方法
技术领域
本发明涉及DSA锥束精确滤波反投影断层成像领域，具体讲是一种适用于介入治疗中使用的DSA设备的X射线断层成像系统及成像方法。
背景技术
介入治疗在DSA影像监视设备引导下，通过导管深入体内，对疾病进行诊断或直接执行治疗，它使临床疾病由创伤变为微创，减少患者的病痛和恢复时间。DSA在介入治疗中具有不可替代的重要诊疗辅助设备。它利用X射线实时显示人体某一部位的透视影像，方便医生动态观察导管、导丝和造影剂等的行程。然而，透视影像是物体三维间信息堆叠显示，很难精确定位器械及病变的空间解剖位置。在颅内介入手术中，需要及时观察颅内出血现象，这需要利用CT设备获取颅内断层图像。在大多数的医院里介入手术室不额外配置CT设备，病人需在CT室与介入手术室间往返推送，手术无菌条件难以保证，容易引起其它并发症，甚至威胁到患者的生命。
DSA设备的C型臂可作旋转运动或多轨迹运动，可实时采集X射线衰减后的信号，并且所发出的X射线在空间中形成锥束，符合CT精确重建的数据采集方式。利用DSA设备实现CT断层功能，降低介入治疗的造价，减少患者医疗费用，是一项具有广泛应用前景的重要技术。
利用DSA设备进行断层重建，最实用的扫描轨迹是单圆弧轨迹。目前有FDK类型和Katsevich类型两类算法可以实现相应的重建。前者属于近似锥束重建算法，在锥角大于2°后出现了伪像。后者属于半精确的锥束重建算法，在锥角小于4°时，能取得良好的重建效果，该算法申请美国专利cone beam filtered backprojection image reconstructionmethod for short trajectories，其专利号为US7,203,272 B2，然而，算法实现比较复杂，不利实际应用。Katsevich发表论文“Image reconstruction for the circle and linetrajectory”是一种精确锥束断层重建方法，然而该算法需要借助螺旋扫描的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构成的复合轨迹，设
Figure BDA0002329633970000061
且C⊥L，L位于C的起始端y0；通过重建点
Figure BDA0002329633970000062
任意的Radon平面与扫描轨迹至少有一个交点，则
Figure BDA0002329633970000063
表达了
Figure BDA0002329633970000064
能够精确重建的源点扫描区间。采集系统引入由
Figure BDA0002329633970000065
标识的局部坐标系；
Figure BDA0002329633970000066
由此，射源
Figure BDA0002329633970000067
和重建点
Figure BDA0002329633970000068
确定的单位矢量方向
Figure BDA0002329633970000069
可表达为：
Figure BDA00023296339700000610
则过射源
Figure BDA00023296339700000611
沿
Figure BDA00023296339700000612
单位矢量方向的锥束投影
Figure BDA00023296339700000613
可表示为：
Figure BDA00023296339700000614
上式中的S2表示三维空间的单位球。
过源点
Figure BDA00023296339700000615
且包含单位矢量
Figure BDA00023296339700000616
的平面记为
Figure BDA00023296339700000617
其单位法向矢量记为
Figure BDA00023296339700000618
而θ表示平面
Figure BDA00023296339700000619
在垂直于矢量
Figure BDA00023296339700000620
局部坐标系中的极角。
Katsevich在A general scheme for constructing inversion algorithm forcone beam CT论文中推导出通用的具有移不变滤波反投影结构的精确锥束重建算法，其重建公式可表示为：
Figure BDA00023296339700000621
式中，Θ表示平面
Figure BDA00023296339700000622
上任意投影射线矢量方向，γ则为矢量方向Θ与
Figure BDA0002329633970000071
之间的夹角；
Figure BDA0002329633970000072
称之为构造因子，其表达式为：
Figure BDA0002329633970000073
构造因子中，
Figure BDA0002329633970000074
表示
Figure BDA0002329633970000075
平面的权重值。由于过重建点的平面与轨迹的交点可能不止一个，则权重值反映了
Figure BDA0002329633970000076
平面对重建点的冗余贡献，需要满足下面关系：
Figure BDA0002329633970000077
在本发明中选择不等权重策略，即重建点的平面同时与圆弧和直线轨迹相交时，直线轨迹的重建点的权重贡献值为零，圆弧点的交点个数的倒数用以确定相应的权重值。
Katsevich指出重建过程仅发生构造因子不为零的投影面，该面是与轨迹相切或者通过轨迹端点的平面。对于圆弧加直线轨迹上存在三类临界面：与轨迹相切以及分别通过圆弧起始和终止端点的平面。为此需判别这三类临界面上构造因子的大小来确定投影数据的滤波方向。
图3示出了当源点为
Figure BDA0002329633970000078
在圆弧轨迹时，锥束扫描几何参数在检测器平面上的投影，上标符号^分别表示相应参数在检测器上的投影，而Πt、Πs和Πe分别标识了过
Figure BDA0002329633970000079
且与轨迹相切的平面、过
Figure BDA00023296339700000710
和圆弧起始端点s的平面和过
Figure BDA00023296339700000711
和圆弧终止端点e的平面与检测器的交线。显然，
Figure BDA00023296339700000712
平面法线矢量
Figure BDA00023296339700000713
在检测器上的投影
Figure BDA00023296339700000714
与
Figure BDA00023296339700000715
平面同检测器的交线相垂直。此外，根据构造因子的含义可知，同一条投影射线上的重建点具有相同的构造因子，并且这些重建点在检测器上的投影地址相同，因此，可以将构造因子存储在检测器对应的阵元上。
图4(a)示出了与圆弧轨迹相切的临界面同检测器平面的交线是沿水平方向的。当该平面绕投影射线分别沿顺时针与逆时针旋转后，与圆弧轨迹的交点个数均为2，与直线轨迹的交点个数为1，因此，
Figure BDA0002329633970000081
而
Figure BDA0002329633970000082
的符号值发生改变，由负变为正，由此可计算出该临界面上的结构因子值为1，并且相应的滤波线是沿水平方向。
图4(b)示出了过圆弧端点的临界面相关参数在检测器上的投影。从图中可看出，重建点
Figure BDA0002329633970000083
在检测器上的投影
Figure BDA0002329633970000084
在检测器
Figure BDA0002329633970000085
轴上坐标正负的不同，临界面平面穿过圆弧端点前后与轨迹的交点个数变化不同；此外，投影点
Figure BDA0002329633970000086
位置的不同，临界面上的
Figure BDA0002329633970000087
的符号不同。因此，计算过圆弧端点的临界面上构造因子相对要复杂些。
过圆弧端点的临界面与检测器的交线决定了投影数据的滤波方向，亦称为滤波线。端点的投影一定在
Figure BDA0002329633970000088
轴上，根据它沿
Figure BDA0002329633970000089
轴的坐标的位置，将滤波线在检测器上的分布划分为三种情况，如图5所示。图5(a)示出了端点投影在检测器外，最外侧两条滤波线的夹角小于等于π/2；图5(b)也示出了端点投影仍然在检测器外，但最外侧两条滤波线夹角小于等于π且大于π/2；图5(c)示出了端点投影在检测器内，滤波线通过端点投影并沿整周分布。
检测器上的投影数据需要沿这些滤波线进行重采样，而且重采样后的投影数据被权重及滤波后，还需要反采样到检测器阵列单元上以方便反投影运算。为了减少采样过程中数据的缺失，对于斜率绝对值小于等于1的滤波线，投影数据应保持检测器水平阵列线的索引不变，而沿着
Figure BDA00023296339700000810
轴进行重采样；对于斜率绝对值大于1的滤波线，投影数据应保持检测器垂直阵列线的索引不变，而沿着
Figure BDA00023296339700000811
轴进行重采样。
经研究得出过圆弧端点处构造因子的绝对值均为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平面以及相应的数值后，还需要推导适合平面检测器上的重建公式。利用偏导数链式法则，可推导出投影数据的偏导数为：
Figure BDA0002329633970000101
构造因子不为零的Radon平面，其上相关的运算才会对重建点有所贡献。实际上，它与检测器平面的交线确定了投影数据的滤波方向。在检测器平面上，滤波线的斜率可以是0到正无穷。为了采样精确，对于斜率绝对值小于1的滤波线，沿水平方向采样，可以推导出偏导数据的加权为：
Figure BDA0002329633970000102
利用Hilbert变换，进一步推导出投影数据的滤波公式：
Figure BDA0002329633970000103
其中，hH(·)表示Hilbert变换滤波核。记gF(λ,u,v)为投影数据偏导数的滤波函数，其表达式为：
Figure BDA0002329633970000104
在检测器平面上，由于通过圆弧端点的滤波线存在斜率绝对值大于1的情况。投影数据需要按
Figure BDA0002329633970000105
轴进行重排，并也需要推导相应的公式。依据上述相同的推导方式，同样获得斜率绝对值大于1的滤波线对应的滤波函数公式gF(λ,u,v)：
Figure BDA0002329633970000106
由此，获得DSA断层图像的精确重建方法。其表达式如下：
Figure BDA0002329633970000107
重建公式中，需要沿扫描轨迹获得在检测坐标u,v处的反投影值，u,v与重建点
Figure BDA0002329633970000108
满足的投影映射关系：
Figure BDA0002329633970000111
此步骤为反投影运算。
重建公式表明本发明的方法由两部分组成。第一部分针对圆弧轨迹的重建公式：
Figure BDA0002329633970000112
其中，c0＝1表示投影数据沿水平方向进行滤波时构造因子的取值；c1和c2分别表示投影数据沿过圆弧起始和终止端点在检测器上的投影点
Figure BDA0002329633970000113
和
Figure BDA0002329633970000114
的直线方向进行滤波时构造因子的取值，依据在检测器上相应的分布规律取0.5或-0.5。
另一部分是针对直线轨迹的重建公式：
Figure BDA0002329633970000115
其中，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进行显示和相关操作。
根据上述的分析，重建过程主要可划分为投影数据的偏导数、重排、滤波以及反投影四个步骤。此外还需说明两个重建模块中，沿圆弧轨迹的重建模块可以单独运行，即在轴向重建范围小的情况下，可以采用单个沿圆弧轨迹的重建模块，可以获得较精确的重建图像，并且重建操作简单。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (1)

1.一种DSA锥束精确滤波反投影断层成像系统，其特征在于：包括成像检测系统(1)、控制系统(2)和计算机系统(3)；成像检测系统(1)包括C型臂(101)，C型臂(101)的一端配有X射线源(102)，用于产生锥束X射线(103)；在C型臂(101)另一端配置二维平面检测器(104)，检测器(104)由若干阵元(105)组成，每个阵元(105)能够感光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)提供控制指令和相应的参数；
所述系统对应的成像方法的步骤包括：成像检测系统(1)中，机械臂(109)能够灵活运动，使安装在C型臂(101)上的射源(102)和检测器(104)作旋转或平移运动，X射线源(102)发射的X射线(103)，检测器(104)上阵元(105)接受衰减信号；控制系统(2)中，数据采集单元(201)将阵元(105)获取的X射线(103)衰减信号转换为锥束投影数据；图像重建单元(301)接收锥束投影数据，经过模块(3011)进行加权，加权后进入模块(3012)，计算关于扫描轨迹参数的偏导数；扫描轨迹若为圆弧轨迹，则进入模块(3013)计算关于检测器参数的偏导数，模块(3014)将偏导数据进行加权求和，模块(3016)对求和数据进行过端点投影重排和希尔伯特滤波；模块(3017)沿水平方向对求和数据进行希尔伯特滤波；滤波后的数据进入模块(3019)进行构造因子处理及反投影运算；若扫描轨迹为直线轨迹，模块(3012)中数据输入到模块(3015)进行加权，再转入模块(3018)，加权数据沿着过切点直线重排和希尔伯特滤波；滤波后的数据也进入模块(3019)中进行构造因子处理及反投影；反投影后，获得物体的断层图像，输入计算机中心(302)，存到大容量存储器(303)，并输入到显示器(305)进行显示；
所述的数据采集单元(201)是把检测器阵元接收的X射线模拟信号经过偏置校正、增益校正、坏像素校正及log运算预处理转化成锥束投影数据；
所述的射源和检测器作旋转是指射源和检测器同时绕C型臂中心轴旋转形成X射线锥束投影的圆弧扫描轨迹；
所述的射源和检测器作平移运动是指机械臂控制射源和检测器同时沿水平方向运动形成X射线锥束投影的直线扫描轨迹；
所述的端点投影是指射源沿圆弧轨迹运动时，射源与圆弧扫描轨迹两个端点分别连线同检测器平面的交点；
所述的过端点投影重排是指射源沿圆弧轨迹运动时，在检测器平面内原水平和竖直分布的求和数据分别沿过圆弧两端点投影的直线进行线性插值重排；
所述的切点直线是指射源沿直线轨迹运动时，圆弧轨迹在检测器平面上的投影是一条抛物线，过抛物线上的每一点所作切线即为切点直线；
所述过切点直线重排是指射源沿直线轨迹运动时，在检测器平面内原水平和竖直分布的加权数据沿过与圆弧轨迹投影的相切直线进行线性插值重排；
所述的构造因子是指射源在圆弧扫描轨迹上时，其值为1，0.5，-0.5，射源在直线扫描轨迹上时，其值为1，-1；
所述的构造因子处理及反投影是指希尔伯特滤波后的数据按照两种扫描轨迹构造因子的分布特点与其值相乘，重建点按照投影关系检索检测器平面的构造因子处理后的滤波数据获得某一投影视角的重建数值。
CN201911331253.6A 2019-12-21 2019-12-21 一种dsa锥束精确滤波反投影断层成像系统及成像方法 Active CN110974278B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911331253.6A CN110974278B (zh) 2019-12-21 2019-12-21 一种dsa锥束精确滤波反投影断层成像系统及成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911331253.6A CN110974278B (zh) 2019-12-21 2019-12-21 一种dsa锥束精确滤波反投影断层成像系统及成像方法

Publications (2)

Publication Number Publication Date
CN110974278A CN110974278A (zh) 2020-04-10
CN110974278B true CN110974278B (zh) 2022-02-11

Family

ID=70073842

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911331253.6A Active CN110974278B (zh) 2019-12-21 2019-12-21 一种dsa锥束精确滤波反投影断层成像系统及成像方法

Country Status (1)

Country Link
CN (1) CN110974278B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113812971B (zh) * 2021-08-27 2023-10-13 浙江大学 一种多自由度四维双能锥束ct成像系统及方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5706325A (en) * 1996-12-05 1998-01-06 General Electric Company Exact regional reconstruction of longitudinally-unbounded objects using a circle-and-line cone beam tomographic system
WO2005107598A1 (en) * 2003-12-04 2005-11-17 Research Foundation Of The University Of Central Florida, Incorporated Efficient circle and line cone beam computed tomography
CN102004111A (zh) * 2010-09-28 2011-04-06 北京航空航天大学 一种倾斜多锥束直线轨迹ct成像方法
CN102973291A (zh) * 2012-12-20 2013-03-20 电子科技大学 C型臂半精确滤波反投影断层成像方法
CN106228584A (zh) * 2016-07-20 2016-12-14 中国人民解放军信息工程大学 锥束ct圆加直线轨迹反投影滤波重建方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7305061B2 (en) * 2001-08-16 2007-12-04 Research Foundation Of The University Of Central Florida Efficient image reconstruction algorithm for the circle and arc cone beam computer tomography
US7197105B2 (en) * 2001-08-16 2007-03-27 Research Foundation Of The University Of Central Florida, Inc. Efficient image reconstruction algorithm for the circle and line cone beam computed tomography
US7359477B2 (en) * 2005-02-15 2008-04-15 Siemens Aktiengesellschaft Method for reconstructing a CT image using an algorithm for a short-scan circle combined with various lines
US7477720B2 (en) * 2005-06-28 2009-01-13 University Of Utah Research Foundation Cone-beam reconstruction using backprojection of locally filtered projections and X-ray CT apparatus
CN101011261A (zh) * 2007-02-08 2007-08-08 上海交通大学 多源螺旋ct反投影滤波精确重建系统
CN100522065C (zh) * 2007-12-06 2009-08-05 上海交通大学 三源马鞍线轨迹锥形束ct精确重建方法
CN103714578A (zh) * 2014-01-24 2014-04-09 中国人民解放军信息工程大学 针对半覆盖螺旋锥束ct的单层重排滤波反投影重建方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5706325A (en) * 1996-12-05 1998-01-06 General Electric Company Exact regional reconstruction of longitudinally-unbounded objects using a circle-and-line cone beam tomographic system
WO2005107598A1 (en) * 2003-12-04 2005-11-17 Research Foundation Of The University Of Central Florida, Incorporated Efficient circle and line cone beam computed tomography
CN102004111A (zh) * 2010-09-28 2011-04-06 北京航空航天大学 一种倾斜多锥束直线轨迹ct成像方法
CN102973291A (zh) * 2012-12-20 2013-03-20 电子科技大学 C型臂半精确滤波反投影断层成像方法
CN106228584A (zh) * 2016-07-20 2016-12-14 中国人民解放军信息工程大学 锥束ct圆加直线轨迹反投影滤波重建方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
圆弧加直线轨迹Katsevich类型的锥束CT重建算法;王瑜，陈亮，欧宗瑛，董彦磊，黄建龙;《电子科技大学学报》;20130930;第794-799页 *

Also Published As

Publication number Publication date
CN110974278A (zh) 2020-04-10

Similar Documents

Publication Publication Date Title
US7142633B2 (en) Enhanced X-ray imaging system and method
JP3373720B2 (ja) Ｘ線断層撮影装置
KR102139668B1 (ko) 단층 촬영 장치 및 그에 따른 단층 영상 복원 방법
EP1513449B1 (en) Rotational angiography based hybrid 3-d reconstruction of coronary arterial structure
RU2550542C2 (ru) Способ и устройство для формирования компьютерных томографических изображений с использованием геометрий со смещенным детектором
JP4346297B2 (ja) Ｘ線コンピュータ断層撮影装置、画像処理装置及び画像処理方法
US7251307B2 (en) Fan-beam and cone-beam image reconstruction using filtered backprojection of differentiated projection data
US6990167B2 (en) Image reconstruction method for divergent beam scanner
EP2317477B1 (en) Reconstruction of 3D image datasets from X-ray cone-beam data
EP0520778B1 (en) Tomographic image reconstruction using cross-plane rays
JPH0661327B2 (ja) 断層撮影像作成方法および装置
WO1995022115A1 (en) Reconstruction of images from cone beam data
CN102973291B (zh) C型臂半精确滤波反投影断层成像方法
US20050100124A1 (en) Methods and apparatus for artifact reduction in computed tomography imaging systems
CN1947154A (zh) 使用截短的投影和在先采集的3d ct图像的锥形束ct设备
CN101133962A (zh) 再现三维立体图像的方法和x射线设备
JP2010510833A (ja) 円状及び直線状軌道を用いた医用撮像における断層撮影法の再構成のための方法及びシステム
EP1663004A2 (en) Computer tomography method using a cone-shaped bundle of rays
US6292526B1 (en) Methods and apparatus for preprocessing volumetric computed tomography data
EP1891603A2 (en) Fast reconstruction algorithm for cone-beam ct
CN110974278B (zh) 一种dsa锥束精确滤波反投影断层成像系统及成像方法
US7272205B2 (en) Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix
JP2007198866A (ja) 広義サドルコーンビームｃｔ装置および３次元再構成法
JPH09182745A (ja) 計算機式断層撮影装置
EP1615561B1 (en) Computer tomography method for a periodically moving object

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant