一种面向射线源直线CT扫描的源轨迹微分图像重建方法

一种面向射线源直线ct扫描的源轨迹微分图像重建方法
技术领域
1.本发明属于ct图像重建领域，更具体地说，涉及一种面向射线源直线ct扫描的源轨迹微分图像重建方法。


背景技术：

2.近年来，为了满足更高的分辨率需求，出现了多种不同的ct扫描方式。刘丰林和余海军等提出一种射线源直线ct扫描的扩大成像视场的结构，该成像几何中探测器被固定，射线源沿平行于探测器行方向进行直线平移，被测物体放于旋转台中心，并设计了基于tv最小化的迭代图像重建算法[1,2]。这种射线源直线ct扫描结构能够通过控制源平移实现扩大成像视野，易实现高精度的控制，在ct扫描技术领域具有巨大的发展潜力，随之适用于该扫描模型的一系列图像重建算法也被相继研究。
[0003]
适用于这种射线源直线ct扫描成像结构首先被提出的是基于tv最小化的迭代图像算法，其虽然可以重建无伪影的图像，但重建时间长、计算成本高，很难满足实际工程需求[1,2]。传统fbp(滤波反投影)解析类重建算法虽然可以进行高效率的重建，但在这种扩大成像视场的直线ct扫描中，采集得到的投影数据存在数据截断，进而滤波反投影得到的重建图像表现出严重的截断伪影[2-5]。为了采用fbp类重建算法避免截断伪影，现有研究采用虚拟投影的方式来实现等效全局无截断的投影，再进行图像重建[3-5]。这种通过先通过将投影重组进而转换成虚拟投影的方式虽可以避免截断伪影，但重建方法依旧不够简单直接。
[0004]
参考文献：
[0005]
[1]刘丰林,余海军,李雷,谭川东.一种新型大视场直线扫描ct系统及图像重建方法[p].重庆市：cn111839568a,2020-10-30.
[0006]
[2]h.yu,l.li,c.tan,f.liu,r.zhou,x-ray source translation based computed tomography(stct),optics express,29(2021)19743-19758.r.clackdoyle,f.noo,a large class of inversion formulae for the 2d radon transform of functions of compact support,inverse problems,20(2004)1281.
[0007]
[3]戈文杰,余海军,陈杰等.基于求导-希尔伯特变换-反投影的源直线扫描计算机断层成像解析重建[j].光学学报,2022,42(11):292-303.
[0008]
[4]李雷,余海军,谭川东,段晓礁,刘丰林.射线源平移扫描ct解析重建算法研究[j].仪器仪表学报,2022,43(02):187-195.doi:10.19650/j.cnki.cjsi.j2108157.
[0009]
[5]yu h,ni s,chen j,et al.analytical reconstruction algorithm for multiple source-translation computed tomography(mstct)[j].applied mathematical modelling,2023,117:251-266.


技术实现要素：

[0010]
为了克服上述现有技术的不足，本发明目的在于提供一种面向射线源直线ct扫描
的源轨迹微分图像重建方法，直接避免截断伪影，重建高质量无伪影的图像。
[0011]
为了实现上述目的，本发明提供二维和三维两种图像重建方法。其中，二维方法采用的技术方案为：
[0012]
s1.初始化参数i＝1，直线扫描段数t和旋转角度间隔δθ，待重建二维图像的零空间
[0013]
s2.获取第i段射线源直线扫描的投影
[0014]
s3.对采集的投影进行预加权，得到加权投影
[0015]
s4.对加权投影沿射线源平移轨迹λ方向进行微分操作；
[0016]
s5.对微分投影进行加权反投影操作，得到第i段hilbert图像
[0017]
s6.对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，重建第i段有限角图像
[0018]
s7.执行i＝i+1，并完成图像的叠加：
[0019]
s8.若i≤t，则使直线ct扫描系统旋转(i-1)
·
δθ角度，并跳转至s2，依次循环直至i＞t，完成图像的重建。
[0020]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s3中，所述投影预加权得到加权投影其计算式为
[0021][0022]
其中，θi为射线源直线扫描轨迹的方向角，θi＝(i-1)
·
δθ；λ表示射线源在直线轨迹上的坐标，λ∈[-s,s]；u表示探测器的坐标，且u∈[-d,d]，d表示探测器长度的一半；l表示沿垂直于探测器方向的射线源到旋转中心的距离；h表示旋转中心到探测器中心的距离；为冗余权重因子。
[0023]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s4中，所述加权投影沿射线源平移轨迹λ方向进行微分操作得到微分投影其计算式为
[0024][0025]
其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点的射线在射线源平移轨迹上的坐标值，
[0026][0027]
其中，l表示待重建点到射线源轨迹的垂直距离，l＝-xsinθi+ycosθi+l；h表示待重建点x到探测器的垂直距离，h＝xsinθ
i-ycosθi+h。
[0028]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在s5步骤中，所述微分投影进行加权反投影操作得到第i段hilbert图像其计算式为，
[0029][0030]
式中，表示待重建点；ηi表示对第i段hilbert图像进行有限hilbert逆变换的方向角；1/h2表示反投影加权因子；所述反投影操作在一个大于待重建图像的矩阵中完成，即在待重建图像矩阵空间每个边缘外增加p0＝0.5*i数量的零，其中i为待重建图像行或列数。
[0031]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s6中，所述对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，可以重建第i段有限角图像其计算式为
[0032][0033]
其中，y1表示一维hilbert变换方向，y1∈[ly+εy,u
y-εy]，其中[ly,uy]表示hilbert变换的有限区间，εy为一个小正数；表示重建第i段有限角图像需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为第i段hilbert图像的有限希尔伯特逆变换的具体实现包括以下步骤：
[0034]
1)将第i段hilbert图像点除-2π，并绕图像中心旋转-ηi的角度，得到旋转hilbert图像其中ηi＝θ
i-π/2(i＝1,2,...,t)，角度值为正，图像逆时针旋转；角度为负，图像顺时针旋转；
[0035]
2)计算权重矩阵w
mat
，即
[0036][0037]
3)搜索旋转hilbert图像中哪些行的数值全为0，得到元素为行号的数组ri；
[0038]
4)对旋转hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像
[0039]
5)依据s3步骤中得到的ri数值，从图像中索引出元素，沿列向计算其均值并取反，得到
[0040]
6)执行
[0041]
7)将旋转图像绕图像中心旋转角度ηi至原状态，得到
[0042]
8)提取中间的待重建图像区域，
[0043]
本发明提供的三维方法所采用的技术方案为：
[0044]
s1.初始化参数i＝1，直线扫描段数t和旋转角度间隔δθ，待重建三维图像的零空间
[0045]
s2.获取第i段射线源直线扫描的锥束投影
[0046]
s3.对采集的锥束投影进行预加权，得到加权锥束投影
[0047]
s4.对加权锥束投影沿射线源轨迹λ方向进行微分操作；
[0048]
s5.对微分锥束投影进行加权反投影操作，得到第i段三维hilbert图像
[0049]
s6.对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层地进行有限hilbert逆变换，重建第i段有限角三维图像
[0050]
s7.执行i＝i+1，并完成三维图像的叠加：
[0051]
s8.若i≤t，则使直线ct扫描系统旋转(i-1)
·
δθ角度，并跳转至s2，依次循环直至i＞t，完成三维图像的重建。
[0052]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s3中，所述锥束投影预加权得到加权锥束投影其计算式为
[0053][0054]
其中，θi为射线源直线扫描轨迹的方向角，θi＝(i-1)
·
δθ；λ表示射线源在直线轨迹上的坐标，λ∈[-s,s]；u表示面板探测器的行向坐标，且u∈[-d,d]，d表示面板探测器行向长度的一半；l表示沿垂直于面板探测器方向的射线源到旋转中心的距离；h表示旋转中心到面板探测器中心的距离；为冗余权重因子。
[0055]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s4中，所述加权锥束投影沿射线源轨迹λ方向进行微分操作得到微分锥束投影其计算式为
[0056][0057]
其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点的射线在射线源轨迹λ上的行向坐标，v表示穿过待重建点的射线在面板探测器上的列向坐标，
[0058][0059][0060]
其中，l表示待重建点到射线源轨迹的垂直距离，l＝-xsinθi+ycosθi+l；h表示待重建点到面板探测器的垂直距离，h＝xsinθ
i-ycosθi+h。
[0061]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在s5步骤中，所述微分锥束投影进行加权反投影操作得到第i段三维hilbert图像其计算式为
[0062][0063]
式中，表示待重建点；ηi表示对第i段三维hilbert图像进行有限hilbert逆变换的方向角；1/h2为反投影加权因子；所述反投影操作在一个大于待重建三维图像的矩阵中完成，即在待重建三维图像矩阵空间每个边缘外增加p0＝0.5*i数量的零，其中i为待重建图像某维度方向上的体素数。
[0064]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s6中，所述对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层进行有限hilbert逆变换，即逐zk(k＝1,2,...,i)层地对二维hilbert图像进行逆变换，进而重建出第i段有限角三维图像二维hilbert图像的逆变换也是沿xj(j＝1,2,...,i)向逐条地进行一维处理，其计算式为
[0065][0066]
其中，y1表示一维hilbert变换方向，y1∈[ly+εy,u
y-εy]，其中[ly,uy]表示hilbert变换的有限区间，εy为一个小正数；表示重建需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为c
yi
；关
于第i段三维hilbert图像的轴向逐层地进行有限hilbert逆变换，其具体实现包括以下步骤：
[0067]
1)将第i段三维hilbert图像点除-2π，并绕三维图像中心轴线旋转-ηi角度，得旋转的三维hilbert图像其中ηi＝θ
i-π/2(i＝1,2,...,t)，角度值为正，三维图像逆时针旋转；角度为负，三维图像顺时针旋转；
[0068]
2)计算一维权重序列w
mat
，即
[0069][0070]
3)for(k＝1；k≤i；k++)(即进入for循环体)：
[0071]
3.1)从旋转的三维hilbert图像索引出第k层二维hilbert图像
[0072]
3.2)搜索二维hilbert图像中哪些行的数值全为0，得到元素为行号的数组ri；
[0073]
3.3)对二维hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像
[0074]
3.4)依据s3步骤中得到的ri数值，从图像中索引出元素，沿列向计算其均值并取反，得到
[0075]
3.5)执行
[0076]
3.6)将旋转的二维图像绕图像中心旋转角度ηi至原状态，得到重建完的第k层二维
[0077]
4)提取中间的待重建三维图像区域，即
[0078][0079]
本发明的有益效果在于：
[0080]
与现有面向射线源直线ct扫描的重建算法相比，本方法能利用截断的投影数据进行精确重建，直接地避免了截断伪影，能够重建出近乎无伪影的二维或三维ct图像。
附图说明
[0081]
图1是本发明的简要流程图。
[0082]
图2是本发明面向的射线源直线ct二维扫描轨迹的示意图。
[0083]
图3是本发明以shepp-logan体模为被测对象的二维图像重建流程图。
[0084]
图4是本发明以forbild体模为被测对象的不同投影数的二维重建图像。
[0085]
图5是本发明面向的射线源直线锥束ct三维扫描轨迹的示意图。
[0086]
图6是本发明的三维hilbert图像沿射线源轨迹方向轴向逐层的有限hilbert逆变
换示意图。
[0087]
图7是本发明的三维ct图像重建流程图。
[0088]
图8是本发明以三维shepp-logan体模为被测对象且在不同投影数条件下的图像重建结果。
具体实施方式
[0089]
下面结合附图和实施例详细说明本发明的实施方式。
[0090]
一种面向射线源直线ct扫描的源轨迹微分图像重建方法，如图1所示，其二维流程包括以下步骤：
[0091]
s1.初始化参数i＝1，直线扫描段数t和旋转角度间隔δθ，待重建二维图像的零空间
[0092]
s2.获取第i段射线源直线扫描的投影
[0093]
s3.对采集的投影进行预加权，得到加权投影
[0094]
s4.对加权投影沿射线源平移轨迹λ方向进行微分操作；
[0095]
s5.对微分投影进行反投影操作，得到第i段hilbert图像
[0096]
s6.对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，重建第i段有限角图像
[0097]
s7.执行i＝i+1，并完成图像的叠加：
[0098]
s8.若i≤t，则使直线ct扫描系统旋转(i-1)
·
δθ角度，并跳转至s2，依次循环直至i＞t，完成图像的重建。
[0099]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s3中，所述投影预加权得到加权投影其计算式为
[0100][0101]
其中，θi为射线源直线扫描轨迹的方向角，θi＝(i-1)
·
δθ；λ表示射线源在直线轨迹上的坐标，λ∈[-s,s]；u表示探测器的坐标，且u∈[-d,d]，d表示探测器长度的一半；l表示沿垂直于探测器方向的射线源到旋转中心的距离；h表示旋转中心到探测器中心的距离；为冗余权重因子。
[0102]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s4中，所述加权投影沿射线源平移轨迹λ方向进行微分操作得到微分投影其计算式为
[0103][0104]
其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点的射线在射线源平移轨迹上的坐标值，
[0105][0106]
其中，l表示待重建点到射线源轨迹的垂直距离，l＝-xsinθi+ycosθi+l；h表示待重建点x到探测器的垂直距离，h＝xsinθ
i-ycosθi+h。
[0107]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在s5步骤中，所述微分投影进行加权反投影操作得到第i段hilbert图像其计算式为，
[0108][0109]
式中，表示待重建点；ηi表示对第i段hilbert图像进行有限hilbert逆变换的方向角；1/h2表示反投影加权因子；所述反投影操作在一个大于待重建图像的矩阵中完成，即在待重建图像矩阵空间每个边缘外增加p0＝0.5*i数量的零，其中i为待重建图像矩阵行或列数。
[0110]
优选地，对于面向射线源直线ct扫描的源轨迹微分二维图像重建方法，在步骤s6中，所述对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，可以重建第i段有限角图像其计算式为
[0111][0112]
其中，y1表示一维hilbert变换方向，y1∈[ly+εy,u
y-εy]，其中[ly,uy]表示hilbert变换的有限区间，εy为一个小正数；表示重建第i段有限角图像需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为第i段hilbert图像的有限希尔伯特逆变换的具体实现包括以下步骤：
[0113]
1)将第i段hilbert图像点除-2π，并绕图像中心旋转-ηi的角度，得到旋转hilbert图像其中ηi＝θ
i-π/2(i＝1,2,...,t)，角度值为正，图像逆时针旋转；角度
为负，图像顺时针旋转；
[0114]
2)计算权重矩阵w
mat
，即
[0115][0116]
3)搜索旋转hilbert图像中哪些行的数值全为0，得到元素为行号的数组ri；
[0117]
4)对旋转hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像
[0118]
5)依据s3步骤中得到的ri数值，从图像中索引出元素，沿列向计算其均值并取反，得到
[0119]
6)执行
[0120]
7)将旋转图像绕图像中心旋转角度ηi至原状态，得到
[0121]
8)提取中间的待重建图像区域，
[0122]
图2展示了射线源直线ct扫描轨迹中的前两段扫描示意图，对应于本发明的第1段射线源直线扫描和第2段射线源直线扫描。
[0123]
图3展示了以shepp-logan体模为被测对象时本发明的图像重建流程，证实了本方法有效性，能够重建无伪影的图像。
[0124]
进一步地，采用forbild体模，像素数量为512
×
512，尺寸为8.4mm
×
8.4mm。通过数值实验对本发明的图像重建效果进行评估。设定每一段射线源直线ct扫描的投影数n分别为50，100，200，400，600和800，模拟计算forbild体模的正投影并采用本发明进行重建，从图4所示的图像重建结果可知，本方法能利用射线源直线ct扫描过程中截断的投影数据重建出无伪影的高质量图像，且图像重建质量随投影数量增多而增加。
[0125]
表1数值仿真试验参数
[0126][0127]
一种面向射线源直线ct扫描的源轨迹微分图像重建方法，如图1所示，其三维方法流程包括以下步骤：
[0128]
s1.初始化参数i＝1，直线扫描段数t和旋转角度间隔δθ，待重建三维图像的零空间
[0129]
s2.获取第i段射线源直线扫描的锥束投影
[0130]
s3.对采集的锥束投影进行预加权，得到加权锥束投影
[0131]
s4.对加权锥束投影沿射线源轨迹λ方向进行微分操作；
[0132]
s5.对微分锥束投影进行加权反投影操作，得到第i段三维hilbert图像
[0133]
s6.对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层地进行有限hilbert逆变换，重建第i段有限角三维图像
[0134]
s7.执行i＝i+1，并完成三维图像的叠加：
[0135]
s8.若i≤t，则使直线ct扫描系统旋转(i-1)
·
δθ角度，并跳转至s2，依次循环直至i＞t，完成三维图像的重建。
[0136]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s3中，所述锥束投影预加权得到加权锥束投影其计算式为
[0137][0138]
其中，θi为射线源直线扫描轨迹的方向角，θi＝(i-1)
·
δθ；λ表示射线源在直线轨迹上的坐标，λ∈[-s,s]；u表示面板探测器的行向坐标，且u∈[-d,d]，d表示面板探测器行向长度的一半；l表示沿垂直于面板探测器方向的射线源到旋转中心的距离；h表示旋转中心到面板探测器中心的距离；为冗余权重因子。
[0139]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s4中，所述加权锥束投影沿射线源轨迹λ方向进行微分操作得到微分锥束投影其计算式为
[0140][0141]
其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点的射线在射线源轨迹λ上的行向坐标，v表示穿过待重建点x的射线在面板探测器上的列向坐标，
[0142][0143][0144]
其中，l表示待重建点到射线源轨迹的垂直距离，l＝-xsinθi+ycosθi+l；h表示待重建点到面板探测器的垂直距离，h＝xsinθ
i-ycosθi+h。
[0145]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在s5步骤中，所述微分锥束投影进行加权反投影操作得到第i段三维hilbert图像其计算式为
[0146][0147]
式中，表示待重建点；ηi表示对第i段三维hilbert图像进行有限hilbert逆变换的方向角；1/h2为反投影加权因子；所述反投影操作在一个大于待重建三维图像的矩阵中完成，即在待重建三维图像矩阵空间每个边缘外增加p0＝0.5*i数量的零，其中i为待重建图像某维度方向上的体素数。
[0148]
优选地，对于面向射线源直线ct扫描的源轨迹微分三维图像重建方法，在步骤s6中，所述对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层进行
有限hilbert逆变换，即逐zk(k＝1,2,...,i)层地对二维hilbert图像进行逆变换，进而重建出第i段有限角三维图像二维hilbert图像的逆变换也是沿xj(j＝1,2,...,i)向逐条地进行一维处理，其计算式为
[0149][0150]
其中，y1表示一维hilbert变换方向，y1∈[ly+εy,u
y-εy]，其中[ly,uy]表示hilbert变换的有限区间，εy为一个小正数；表示重建需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为关于第i段三维hilbert图像的轴向逐层地进行有限hilbert逆变换，其具体实现包括以下步骤：
[0151]
1)将第i段三维hilbert图像点除-2π，并绕三维图像中心轴线旋转-ηi角度，得旋转的三维hilbert图像其中ηi＝θ
i-π/2(i＝1,2,...,t)，角度值为正，三维图像逆时针旋转；角度为负，三维图像顺时针旋转；
[0152]
2)计算一维权重序列w
mat
，即
[0153][0154]
3)for(k＝1；k≤i；k++)(即进入for循环体)：
[0155]
3.1)从旋转的三维hilbert图像索引出第k层二维hilbert图像
[0156]
3.2)搜索二维hilbert图像中哪些行的数值全为0，得到元素为行号的数组ri；
[0157]
3.3)对二维hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像
[0158]
3.4)依据s3步骤中得到的ri数值，从图像中索引出元素，沿列向计算其均值并取反，得到
[0159]
3.5)执行
[0160]
3.6)将旋转的二维图像绕图像中心旋转角度ηi至原状态，得到重建完的第k层二维
[0161]
4)提取中间的待重建三维图像区域，即
[0162]
[0163]
图5展示了射线源直线锥束ct扫描轨迹的示意图，图6展示了对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层进行有限hilbert逆变换的示意图。
[0164]
图7展示了本发明的三维ct图像重建基本流程，可以证实本方法的有效性，能够重建高质量近乎无伪影的三维ct图像。
[0165]
进一步地，采用三维shepp-logan体模，像素数量设为512
×
512
×
512，尺寸为8.4mm
×
8.4mm
×
8.4mm。通过定量数值实验对本发明的图像重建效果进行评估。重建图像结果如图8所示，设定每一段射线源直线ct扫描下的不同投影数n的条件，模拟计算三维shepp-logan体模的锥束正投影并采用本发明进行重建，观察重建图像质量，其中n分别为设为200和1600，数值实验的参数设置如表2所示。从图8所示的重建结果可以可知，本发明在投影数较多时能重建高质量近乎无伪影的三维ct图像。
[0166]
表2数值仿真试验参数
[0167][0168]
本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的；本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

技术特征：
1.一种面向射线源直线ct扫描的源轨迹微分图像重建方法，其特征在于，其二维或三维方法包括以下步骤：s1.初始化参数i＝1，直线扫描段数t和旋转角度间隔δθ，待重建二维图像的零空间s2.获取第i段射线源直线扫描的投影或s3.对采集的投影或进行预加权，得到加权投影或s4.对加权投影或沿射线源平移轨迹λ方向进行微分操作；s5.对微分投影进行加权反投影操作，得到第i段hilbert图像s6.对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，重建第i段有限角图像s7.执行i＝i+1，并完成图像的叠加：s8.若i≤t，则使直线ct扫描系统旋转(i-1)
·
δθ角度，并跳转至s2，依次循环直至i＞t，完成图像的重建。2.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分图像重建方法，其特征在于，在步骤s3中，所述投影或预加权得到加权投影或其计算式为或其中，θ
i
为射线源直线扫描轨迹的方向角，θ
i
＝(i-1)
·
δθ；λ表示射线源在直线轨迹上的坐标，λ∈[-s,s]；u表示探测器的坐标，且u∈[-d,d]，d表示探测器长度的一半；l表示沿垂直于探测器方向的射线源到旋转中心的距离；h表示旋转中心到探测器中心的距离；为冗余权重因子。3.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分二维图像重建方法，其特征在于，在步骤s4中，所述加权投影沿射线源平移轨迹λ方向进行微分操作得到微分投影其计算式为
其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点x的射线在射线源平移轨迹上的坐标值，其中，l表示待重建点x到射线源轨迹的垂直距离，l＝-xsinθ
i
+ycosθ
i
+l；h表示待重建点到探测器的垂直距离，h＝xsinθ
i-ycosθ
i
+h。4.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分二维图像重建方法，其特征在于，在s5步骤中，所述微分投影或进行加权反投影操作得到第i段hilbert图像其计算式为，式中，或表示待重建点；η
i
表示对第i段hilbert图像进行有限hilbert逆变换的方向角；1/h2表示反投影加权因子；所述反投影操作在一个大于待重建图像的矩阵中完成，即在待重建图像矩阵空间每个边缘外增加p0＝0.5*i数量的零，其中i为待重建图像行或列数。5.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分二维图像重建方法，其特征在于，在步骤s6中，所述对第i段hilbert图像沿第i段射线源直线平移轨迹方向进行有限希尔伯特逆变换，可以重建第i段有限角图像其计算式为其中，y1表示一维hilbert变换方向，y1∈[l
y
+ε
y
,u
y-ε
y
]，其中[l
y
,u
y
]表示hilbert变换的有限区间，ε
y
为一个小正数；表示重建第i段有限角图像需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为第i段hilbert图像的有限希尔伯特逆变换的具体实现包括以下步骤：1)将第i段hilbert图像点除-2π，并绕图像中心旋转-η
i
的角度，得到旋转hilbert图像其中η
i
＝θ
i-π/2(i＝1,2,...,t)，角度值为正，图像逆时针旋转；角度为负，图像顺时针旋转；2)计算权重矩阵w
mat
，即
3)搜索旋转hilbert图像中哪些行的数值全为0，得到元素为行号的数组r
i
；4)对旋转hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像5)依据s3步骤中得到的r
i
数值，从图像中索引出元素，沿列向计算其均值并取反，得到c
yi
；6)执行7)将旋转图像绕图像中心旋转角度η
i
至原状态，得到8)提取中间的待重建图像区域，6.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分三维图像重建方法，其特征在于，在步骤s4中，所述加权锥束投影沿射线源轨迹λ方向进行微分操作得到微分锥束投影其计算式为其中，表示沿λ方向的微分算子，由有限差分运算实现，投影内部的数据采用中心差分，投影边界的数据采用单侧差分；λ
*
表示穿过待重建点x的射线在射线源轨迹λ上的行向坐标，v表示穿过待重建点的射线在面板探测器上的列向坐标，的射线在面板探测器上的列向坐标，其中，l表示待重建点x到射线源轨迹的垂直距离，l＝-xsinθ
i
+ycosθ
i
+l；h表示待重建点到面板探测器的垂直距离，h＝xsinθ
i-ycosθ
i
+h。7.根据权利要求1所述的面向射线源直线ct扫描的源轨迹微分三维图像重建方法，其特征在于，在步骤s6中，所述对第i段三维hilbert图像沿第i段射线源直线平移轨迹方向轴向逐层进行有限hilbert逆变换，即逐z
k
(k＝1,2,...,i)层地对二维hilbert图像进行逆变换，进而重建出第i段有限角三维图像二维hilbert图像的逆变换也是沿x
j
(j＝1,2,...,i)向逐条地进行一维处理，其计算式为
其中，y1表示一维hilbert变换方向，y1∈[l
y
+ε
y
,u
y-ε
y
]，其中[l
y
,u
y
]表示hilbert变换的有限区间，ε
y
为一个小正数；表示重建需要计算的未知常数，通过找出已知的一些位置点，将其经有限hilbert变换后的均值来作为关于第i段三维hilbert图像的轴向逐层地进行有限hilbert逆变换，其具体实现包括以下步骤：1)将第i段三维hilbert图像点除-2π，并绕三维图像中心轴线旋转-η
i
角度，得旋转的三维hilbert图像其中η
i
＝θ
i-π/2(i＝1,2,...,t)，角度值为正，三维图像逆时针旋转；角度为负，三维图像顺时针旋转；2)计算一维权重序列w
mat
，即3)for(k＝1；k≤i；k++)(即进入for循环体)：3.1)从旋转的三维hilbert图像索引出第k层二维hilbert图像3.2)搜索二维hilbert图像中哪些行的数值全为0，得到元素为行号的数组r
i
；3.3)对二维hilbert图像加权w
mat
，并沿图像矩阵列向进行hilbert变换，得到图像3.4)依据s3步骤中得到的r
i
数值，从图像中索引出元素，沿列向计算其均值并取反，得到3.5)执行3.6)将旋转的二维图像绕图像中心旋转角度η
i
至原状态，得到重建完的第k层二维4)提取中间的待重建三维图像区域，即

技术总结
一种面向射线源直线CT扫描的源轨迹微分图像重建方法，包括二维方法和三维方法。步骤如下：S1.初始化参数i＝1，直线扫描段数T和旋转角度间隔Δθ，待重建图像零空间S2.获取第i段射线源直线扫描的投影；S3.对投影数据预加权；S4.对加权投影沿源轨迹方向微分；S5.对微分投影进行加权反投影，得到第i段Hilbert图像；S6.对Hilbert图像沿第i段射线源轨迹进行有限希尔伯特逆变换，重建第i段有限角图像S7.使i＝i+1，S8.若i≤T，则使直线CT扫描系统旋转(i-1)


技术研发人员：崔俊宁 汪志胜 王顺利 边星元 赵亚敏
受保护的技术使用者：哈尔滨工业大学
技术研发日：2023.06.06
技术公布日：2023/9/14