CN107016709A - 多源摆动动态ct成像方法 - Google Patents

多源摆动动态ct成像方法 Download PDF

Info

Publication number
CN107016709A
CN107016709A CN201710233458.5A CN201710233458A CN107016709A CN 107016709 A CN107016709 A CN 107016709A CN 201710233458 A CN201710233458 A CN 201710233458A CN 107016709 A CN107016709 A CN 107016709A
Authority
CN
China
Prior art keywords
image
data
dynamic
reconstruction
projection
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.)
Granted
Application number
CN201710233458.5A
Other languages
English (en)
Other versions
CN107016709B (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.)
Chongqing University
Original Assignee
Chongqing University
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 Chongqing University filed Critical Chongqing University
Priority to CN201710233458.5A priority Critical patent/CN107016709B/zh
Publication of CN107016709A publication Critical patent/CN107016709A/zh
Application granted granted Critical
Publication of CN107016709B publication Critical patent/CN107016709B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种多源摆动动态CT成像方法,首先把采集的多对X射线源-探测器一次摆动扫描的所有投影数据进行重排;接着对重排后的投影数据采用FBP算法重建出CT图像;把上述CT图像作为先验图像,通过带TV约束的ART算法把相应的时间帧的投影数据进行图像重建获得对应时间帧的CT图像;把所有时间帧的CT图像按时间顺序组合获得该次X射线源-探测器摆动时间段内检测对象的动态CT图像;重复上述过程,最后获得整个检测对象变化过程的动态CT图像。本发明所提出MS‑PICCS算法不仅能够从高度欠采样的数据恢复真实图像而且还具有良好的抗噪性能。

Description

多源摆动动态CT成像方法
技术领域
本发明涉及CT扫描领域,具体涉及一种多源摆动动态CT成像方法。
背景技术
实际应用中存在过程可视化检测需求,如增材制造缺陷形成、铸件冷却缺陷形成、材料(新材料、矿石)拉压破坏、空隙材料(石油岩心、高速路面)液体渗透、多相流体动态流动、固体火箭发动机喷管喉颈材料烧蚀等过程,传统无损检测方法都有一定的局限性。近年来,工业CT技术已成为工业无损检测的重要手段,但由于上述可视化检测过程中检测对象内部发生变化,而CT扫描过程中要求检测对象稳定,传统工业CT技术不能直接用于过程可视化检测。然而,上述需要可视化检测的过程有以下共同特点:(1)变化过程有一定持续时间(非瞬时);(2)同时其过程变化为检测对象内部的局部变化。这些特点使动态CT检测成为可能。
动态过程的可视化检测研究主要包括:RB Arthur和MA Cheverton通过融合多个高速相机采集的不同数据类型从而对增材制造过程的变化进行成像;光滑粒子流体动力学方法和大涡流模拟方法用于预测和可视化铸造时流体瞬时流动以及氧化物含量的动态变化;通过融合X射线温度梯度阶段(XTGS),X射线显微成像技术和尺度模拟方法观察铝-铜合金的固化冷却过程中缩松以及锁孔的生长率;Xuanhao Sun等人通过高放大倍率光学显微镜和原子力显微镜(AFM)研究皮质骨的屈服后的变形过程等。另外,微观粒子图像速度测定技术也被用于研究多孔微型模型中动态溶液的相互作用。然而,已有的这些技术都有一定局限性,如AFM只能观察变化过程的表面结构,对其内部结构的变化显得无能为力。又如平滑粒子流体动力学方法空间分辨率较低。
发明内容
鉴于此,本发明提供一种多源摆动动态CT成像方法。
本发明的目的是通过以下技术方案实现的,一种多源摆动动态CT成像方法,该方法以基于先验图像约束的压缩感知图像重建算法为基础,具体包括以下步骤为:步骤一、把采集的多对X射线源-探测器单次摆动的所有投影数据进行重排;步骤二、对重排后的投影数据采用FBP算法重建出CT图像;步骤三、把上述CT图像作为先验图像,通过带TV约束的ART算法把相应的时间帧的投影数据进行图像重建获得对应时间帧的CT图像;步骤四、把所有时间帧的CT图像按时间顺序组合获得该次X射线源-探测器摆动时间段内检测对象的动态CT图像;步骤五、重复步骤三、四,最后获得整个检测对象变化过程的动态CT图像。
进一步,基于先验图像约束压缩感知图像重建算法被设计为以下无约束最小化问题:
minκ||TV(I)||1+(1-κ)||TV(I-Ip)||1such that AI=P (1)
其中,||·||1表示矩阵L1范数,A代表系统矩阵,I,P分别表示待重建图像和测量的投影数据,TV(I),TV(I-Ip)分别表示待重建图像以及待重建图像与先验图像的总变分,κ是由先验图像Ip中的伪影程度决定的经验参数,取值范围为[0,1]。
进一步,通过采用梯度下降搜索方法对公式()进行求解,图像中(i,j)点的梯度方向hi,j可由分量ai,j,bi,j,ci,j,di,j和ei,j的计算得到:
其中,I(k) i,j,I(k) i-1,j,I(k) i,j-1分别表示在第k次迭代获得的图像中位置为(i,j),(i-1,j)以及(i,j-1)的灰度值,ε为很小的正数,防止分母等于0的情况。
其中,I(k) i+1,j,I(k) i-1,j-1分别表示在第k次迭代获得的图像中位置为(i+1,j)和(i-1,j-1)的灰度值。
其中,I(k) i,j+1,I(k) i-1,j+1分别表示在第k次迭代获得的图像中位置为(i,j+1)和(i-1,j+1)的灰度值
其中,I(k)表示在第k次迭代获得的图像。
其中,Ip i,j表示先验图像中位置为(i,j)对应的图像值。
总梯度表示为
β=max(|κI(k) i,j+(1-κ)(I(k) i,j-Ip i,j)|)/max(|hi,j|)
Ii,j (k+1)=Ii,j (k)-σ×β×hi,j
其中,β表示当前图像梯度的最大值以及先验图像与重建图像靠近的程度。
σ=σ×σs
σs表示梯度下降的调整步长。
由于采用了上述技术方案,本发明具有如下的优点:
本发明所提出MS-PICCS算法不仅能够从高度欠采样的数据恢复真实图像而且还具有良好的抗噪性能。
附图说明
为了使本发明的目的、技术方案和优点更加清楚,下面将结合附图对本发明作进一步的详细描述,其中:
图1为MSSCT系统配置平面图;
图2为MSSCT系统的投影数据模型,不同颜色点对应不同的时间帧;
图3为原始图像与重建图像,其中图a、b、c分别为时间帧1、9、20的原始图像,图d、e、f分别为采用FBP算法在系统配备7个X射线源时采用该时间帧全数据的重建图像,图g、h、i分别为采用FBP算法在系统配备9个X射线源时采用该时间帧全数据的重建图像;
图4为中心圆盘面积动态变化曲线;
图5为对于9个X射线源的欠采样因子为40的先验图像;
图6为重建图像,其中图a、d、g为ART算法的重建图像,欠采样因子分别为50,25,10;图b、e、h为TVM-SD算法的重建图像,欠采样因子分别为50,25,10;图c、f、i为MS-PICCS算法的重建图像,欠采样因子分别为50,25,10;
图7为重建图像,其中图a、d、g为ART算法的重建图像,欠采样因子分别为40,20,8;图b、e、h为TVM-SD算法的重建图像,欠采样因子分别为40,20,8;图c、f、i为MS-PICCS算法的重建图像,欠采样因子分别为40,20,8;
图8为重建图像,与图6相同,但是投影数据添加噪声;
图9为重建图像,与图7相同,但是投影数据添加噪声;
图10为左栏(a):采用FDK重建的荞麦三维结构图;(b)和(c)(即中间和右列)表示采用FBP算法使用全数据重建完整的切片,每个重建的切片图像由320×320像素组成,每个覆盖31.44×31.44μm2,矩形表示均方根误差RMSE测量的子区域;
图11为对于9个X射线源的欠采样因子50的先验重建图像(左),先验图像和重建图像之间的差异图像(右);
图12为对于显示窗口中的7个X射线源/检测器对,利用欠采样因子26(第一行),13(中间行)和10(第三行)在时间帧0处的重建图像,第一列至第三列分别表示在相同的模拟实验条件下使用ART,TVM-SD和MS-PICCS算法在200次迭代后的重建图像,同时,每个重建图像由320×320像素组成;
图13为采用9个X射线源/检测器对的MSSCT系统,且使用欠采样因子50(第一行),25(中间行)和10(第三行)在时间帧0处的重建图像,第一列到第三列分别表示在相同的模拟实验条件下使用ART,TVM-SD和MS-PICCS算法在200次迭代之后的图像重建,同时,每个重建图像由320×320像素组成;
图14为图12时间帧9处的重建图像;
图15为图13时间帧9处的重建图像。
具体实施方式
以下将结合附图,对本发明的优选实施例进行详细的描述;应当理解,优选实施例仅为了说明本发明,而不是为了限制本发明的保护范围。
传统的X射线CT受扫描速度的限制,其时间分辨能力较低。为此,提出一种“多X射线源摆动动态CT(multi-source swing CT,MSSCT)”检测方法用于过程可视化检测。一方面,使用多个X射线源,一次采样就可以获得多个角度的投影数据,从而提高CT检测的时间分辨能力;另一方面,摆动扫描方法可以避免由于检测对象系统复杂(存在附属设备),难以将检测对象系统放置到传统工业CT转台(检测对象360°旋转)上进行检测,或采用传统CT系统滑环(X射线源\探测器放置在滑环上旋转)成本高、难度大、系统复杂等问题,易于工程实现。
MSSCT系统由N对圆周均匀分布的X射线源/检测器组成(N取奇数)。每个X射线源以及对应探测器固定在旋转机构上,同心且共面。在CT一次扫描时,所有X射线源/探测器对同时同向摆动相应角度(例如N=9,每个X射线源的摆动角度为40o)。由于该系统中每个射线X射线源/探测器的扫描运动方式为往复摆动,所以系统不需要配置传统的CT滑环;另外,CT扫描时,系统一次采样可同时获得N个圆周均布的投影数据用于图像重建,因此该系统可以实现高时间分辨率的动态CT成像。根据CT常识,如果对象的状态在扫描期间保持不变,当系统沿着一个方向旋转2π/N的角度,则可获得图像重建所需要的完备数据。如图1,假设从X射线源到检测视场中心的距离是h。检测视场中心到探测器的距离为d,则视场半径R为
为了进一步分析该系统,可以进一步假设对象内存在连续变化的动态区域,如图1所示。首先,X射线源被固定在给定的初始位置并且记为时间帧0,所有射线X射线源曝光则获得时间帧0的投影数据,其对应于图2中的点1。然后,X射线源旋转到相邻位置并获得时间帧1(图2中的点2)的投影数据。而当X射线源到达点3时,假设X射线源已经到达其单向旋转极限位置,就可获得时间帧T处的投影数据。可以这样考虑,当物体保持静止时,组合上述所有时间帧数据即可获得完整的投影数据,此时被称为运动上半周期。最后X射线源沿着相反方向按照同样方式回到初始位置,此时段称之为下半周期。
在本研究中,为了获得连续动态对象的超高时间分辨率重建图像,以经典图像重建算法:基于先验图像约束压缩感知(PICCS)图像重建算法为基础,提出适合MSSCT系统的MS-PICCS算法。
经典PICCS算法通过引入与真实图像相似的先验图像可以从高度欠采样的数据中重建出稀疏图像,其中先验图像一般通过完备数据下FBP算法获得。
假设先验图像Ip,元素J×K,PICCS算法被设计以下无约束最小化问题:
minκ||TV(I)||1+(1-κ)||TV(I-Ip)||1such that AI=P (1)
其中||·||1被表示为矩阵L1范数,A代表系统矩阵,I,P分别表示待重建的图像和测量到的投影数据,TV(I-Ip)表示先验图像与待重建图像的差图像的总变分,κ(0≤κ≤1)是由先验图像Ip中的伪影程度决定的经验参数。也就是说当先前图像包含伪影较少,κ应该接近1。相反,如果先前图像被高度污染,则κ将趋向于0。事实上,当使用FBP算法获得先验图像时,κ=0.5适应绝大部分情况。TV是公式(1)中的离散梯度变换,可以写为以下表达式:
其中J,K分别表示图像矩阵的像素边长,求解公式(1)的方法很多。在本发明中通过采用梯度下降搜索方法来实现。MS-PICCS算法伪代码如下:
S1重排单次CT扫描(摆动)的投影数据;
S2采用FBP算法从重排数据重建CT图像作为先验图像。
S3初始化迭代次数k=0,平衡因子κ;
S4k=k+1和I(k)=I(k-1);k为当前迭代次数S5ART算法重建图像;
S6初始化最陡下降的最大步长σ和减小σ的步长σs=0.997;
S7重复TV最小化循环;
计算梯度的方向,
其中,I(k) i,j,I(k) i-1,j,I(k) i,j-1分别表示在第k次迭代获得的图像中位置为(i,j),(i-1,j)以及(i,j-1)的灰度值,ε为很小的正数,防止分母等于0的情况。
其中,I(k) i+1,j,I(k) i-1,j-1分别表示在第k次迭代获得的图像中位置为(i+1,j)和(i-1,j-1)的灰度值。
其中,I(k) i,j+1,I(k) i-1,j+1分别表示在第k次迭代获得的图像中位置为(i,j+1)和(i-1,j+1)的灰度值
其中,I(k)表示在第k次迭代获得的图像。
其中,Ip i,j表示先验图像中位置为(i,j)对应的图像值。
总梯度表示为
β=max(|κI(k) i,j+(1-κ)(I(k) i,j-Ip i,j)|)/max(|hi,j|)
Ii,j (k+1)=Ii,j (k)-σ×β×hi,j
其中,β表示当前图像梯度的最大值以及先验图像与重建图像靠近的程度,
σ=σ×σs
σs表示梯度下降的调整步长。
S8返回(4),满足给定条件停止迭代。
TVM-SD(total variation minimum steep discent)算法已被证明在多X射线源CT有限角内部成像系统的图像重建中有着巨大优势。一方面,与TVM-SD算法不同的是,MS-PICCS算法需要计算当前图像I(k)和差分图像I(k)-Ip的混合梯度h,而前者只需要计算当前图像的梯度。其次,参数β的更新不仅取决于当前图像和差分图像,而且还与参数κ有关。另一方面,与传统的PICCS算法相比,MS-PICCS算法的优点是一系列参数的自适应性,即MS-PICCS能够搜索梯度的最佳步长。
为了验证MS-PICCS算法在多X射线源摆动X射线成像系统动态图像重建方面的优势,在实验中给出了ART和TVM-SD算法重建图像,并且将TVM-SD和ART的初始猜测图像设置为零。将MS-PICCS算法的初始猜测图像设置为采用FBP算法从一次单方向获得数据的重建图像。采用数值模拟和荞麦实际数据进行研究,并在PC上采用matlab编程实现。
动态图像由17个圆组成,如图3所示。在该实验中,通过在模体中引入9个不同半径的圆来评估空间分辨率,同时引入绕中心圆分布的五个大小相同圆测试算法密度分辨率。假设中心圆盘的面积连续变化,如图4所示,在本实验中,假定射线X射线源/探测器对的数量为7或9,扫描的详细参数如表1。从表1可以推断得出,当射线X射线源/探测器数目为7个时,在X射线源的圆形台架上均匀分布350个投影,当X射线源数目为9个时候,均匀分布360个投影。
图3第一行包含来自不同时间帧1,9和20的原始图像。第二行和第三行分别是采FBP算法在系统配备7个和9个X射线源时采用该时间帧全数据的重建图。
表1数值模拟参数
为了验证9射线X射线源系统的MS-PICCS图像重建算法,选择三个欠采样因子w=8,20,40,其分别对应于每个时间帧中采集45,18和9个视角。同样的,当只有7个X射线源时,如果选择w=10,25,50,相应地每一时间帧分别具有35,14和7个视角。根据时间不同,组合半个时间周期投影数据产生360°数据集,重排投影数据,采用FBP算法重建先验图像,如图5所示。显然,动态变化中心圆附近存在由于数据不一致引起的条形伪影。在本发明中只给出不同的算法在时间帧6处迭代500次迭重建图像。且实验中MS-PICCS重建参数κ设置为0.52。TV次数和最陡下降最大步长分别设置为5和0.015,重建图像如图6和图7所示。可以观察到,通过利用ART和TVM-SD算法,重建图像中出现严重的有限角伪影。相比之下,当采用MS-PICCS算法时,即使只有7或9个角度,依然可以获得高质量的重建图像。
为了测试MS-PICCS在抗噪性方面的优势,通过加入均匀分布的高斯噪声使投影数据更符合实际情况,并且选择无噪声投影数据最大值的1.8%作为高斯噪声标准方差。相关重建结果显示如图8和图9。从这些结果,可以推断MS-PICCS算法具有良好的抑制噪声的性能。
图6对于7个X射线源的欠采样因子50(第一行),25(中间行)和10(第三行)的图像重建。第一列至第三列分别表示使用ART,TVM-SD和MS-PICCS算法的图像重建,每个重建图像有256×256个像素。
图7针对9个X射线源的欠采样因子40(第一行),20(中间行)和8(第三行)的图像重建。第一列到第三列分别表示使用ART,TVM-SD和MS-PICCS算法的图像重建,每个重建图像由256×256个像素组成。
为了量化图像重建的质量,通过计算所选择的子区域(如图5中所示的矩形区域)的均方根误差(RMSE)来定量评估重建图像质量。
其中I*和Ir分别表示参考图像和重建图像,J1和K1分别表示所选区域大小的像素边长。由于数值模拟研究中,原始图像不包含任何噪声和伪像,可以选择作为参考图像。然而,由于理想图像在实际数据研究中是未知的,所以通过使用FBP或ART算法从完整的数据重建的图像可以近似作为参考图像。注意,来自实际数据的这个参考图像可以包含较少的噪声,甚至误差。相关结果展示在表2中,从表2中可以看出,配备9个X射线源/检测器对的系统优于7个X射线源/检测器对结构。因为9个X射线源/探测器对系统可以一次获得比7个X射线源/探测器对系统更均匀数据,然而如果X射线源的数量越大,则系统的视场将越小。对于给定数量的X射线源的系统,每个X射线源在单独时间帧处扫描角度越大,重建的图像将更接近真实的图像,但较大的扫描角度会大大降低系统的时间分辨率。因此,需要在X射线源/检测器对数量,视场和时间分辨率之间做出折中。
表2不同采样因子以及X射线源的RMSE数值对比
为了证明扫描结构和相关算法在实际应用中的有用性和可行性,利用单个X射线源和平面探测器的微焦CT系统(μCT),采用锥形束扫描了萌芽荞麦。图10(a)显示通过使用FDK算法的3D重建图像。在本发明中,从微焦点X射线源到平面检测器的距离为578.40mm,X射线源与旋转轴之间的距离为50.30mm,旋转一圈获得采集1000个投影即每个相邻采样角度为0.36°。此外,平板探测器由1536×876个元素组成,每个元素覆盖0.075×0.075平方毫米的面积。检测器在旋转轴方向上的偏移量为0.525mm,在横向偏移为0.01275mm。则最大锥形角τmax为:
τmax=2arctan((876-1)*0.075/578.40)=6.4938
因此,当只使用接近中间平面的检测器的一部分时,锥角将非常小。也就是说,如果在中心平面附近选择50行检测器,则最大锥角可以减小到0.3715°。实际上,可以使用由平板探测器的每个切片采集的数据来执行近似模拟动态物体图像重建。由于荞麦的结构相似性,即除了小部分的萌芽变化之外,大部分结构是相似的。因此,可以从中间平面附近的50行线阵列提取投影数据,完成相关算法的试验研究。更准确地说,如果X射线源的数量是9,将提取50个1536×900的正弦图矩阵(每10个投影中丢弃1个投影),欠采样因子设置为10,25,50。当X射线源的数量降到7时,可以提取50个1536×910正弦图矩阵(每100个投影中丢弃9个投影),并且相应的欠采样因子为10,13,26。事实上,每个正弦图矩阵可以被认为是不同时间帧的投影数据。然后,根据X射线源的数量N和给定的欠采样因子w,可以进一步提取每个数据矩阵的一部分以生成新的正弦图矩阵。注意新的正弦图的数据一致性损坏,这导致来自动态部分的重建图像模糊(图11)。在这个实验中,只分析时间帧0和时间帧9的结果。图10(b)和(c)显示了使用经典FBP方法从原始全正弦图重建的图像,其被认为是参考图像。参考图像被噪声干扰,它可能影响定量化的结果。选择图10(b)和(c)中由矩形表示的子区域以量化图像重建质量。与上述数值模拟实验相同,使用ART算法,TVM-SD算法和MS-PICCS算法重建时间帧0的图像,结果示于图12和13。同样,时间帧9的结果是如图14和15所示。
图10左栏(a):采用FDK重建的荞麦三维结构图。(b)和(c)(即中间和右列)表示采用FBP算法使用全数据重建完整的切片。每个重建的切片图像由320×320像素组成,每个覆盖31.44×31.44μm2。矩形表示RMSE测量的子区域。
图11对于9个X射线源的欠采样因子50的先验重建图像(左)。先验图像和重建图像之间的差异图像(右)。
图12对于显示窗口中的7个X射线源/检测器对,利用欠采样因子26(第一行),13(中间行)和10(第三行)在时间帧0处的重建图像。第一列至第三列分别表示在相同的模拟实验条件下使用ART,TVM-SD和MS-PICCS算法在200次迭代后的重建图像。同时,每个重建图像由320×320像素组成。
图13采用9个X射线源/检测器对的MSSCT系统,且使用欠采样因子50(第一行),25(中间行)和10(第三行)在时间帧0处的重建图像。第一列到第三列分别表示在相同的模拟实验条件下使用ART,TVM-SD和MS-PICCS算法在200次迭代之后的图像重建。同时,每个重建图像由320×320像素组成。
通过与ART和TVM-SD算法重建图像进行比较,可以推断出MS-PICCS方法即使在高度欠采样的数据中也具有恢复出高质量图像的能力。例如,配备9个X射线源/检测器对系统中18个视图或9个X射线源/检测器对的21个视图的重建图像质量依然能够满足要求。为了量化重建图像误差,选择图像中矩形区域,如图10所示,并计算参考图像之间与重建图像之间的均方根误差,相关结果列在表3中,
表3荞麦选定的子区域均方根误差(单位:10-3)
本发明针对增材制造缺陷形成过程等可视化无损检测需求,提出一种“多X射线源摆动动态CT”方法及系统。该系统使用多个X射线源,一次采样就可以获得多个角度的投影数据,可以提高CT检测的时间分辨能力;同时,摆动扫描方法可以避免难以将检测对象系统放置到传统工业CT转台(检测对象360度旋转)上进行检测,或采用传统CT系统滑环(X射线源\探测器放置在滑环上旋转)成本高、难度大、系统复杂等问题,易于工程实现。
本发明提出多X射线源摆动动态CT模型,提出基于先验图像的MS-PICCS图像重建算法。具体做法:重排一次摆动的搜游投影数据,采用FBP算法对重排数据重建出先验图像,运用先验图像采用带TV约束优化的ART迭代重建算法分别重建出各时段图像。
具体通过设计气泡铸件冷却过程中生长与消失过程来模拟铸件冷却过程中缺陷的形成过程,该过程主要通过模体中圆半径的大小来实现。此外,为了进一步凸显MS-PICCS算法在多X射线源动态CT图像重建中的优势,在试验中提供ART、TVM-SD算法的图像重建结果,通过对实验结果的定性与定量分析表明(如表2),本发明所提出MS-PICCS算法不仅能够从高度欠采样的数据恢复真实图像而且还具有良好的抗噪性能。
为表明该系统的实际可行性以及算法对实际数据的适应性,本发明采用荞麦数据实验,具体做法即把物体三维结构当做随时间变化的二维图像,并把单X射线源采集实际数据用于该系统算法的研究。该实验进一步验证MSSCT系统的可行性。此外,通过与相同条件采用ART、TVM-SD算法的重建结果定量与定性分析表明(表3),本发明所提出MS-PICCS算法也能很好的适应实际数据并且高精度的重建出图像。通过以上研究分析表明本发明提出多X射线源摆动动态CT系统即MSSCT系统以及针对该系统的MS-PICCS算法在理论和实际上都是完全可行的,能够满足增材制造成形等过程可视化无损检测需求,有望进一步推动工业过程工艺发展。
以上所述仅为本发明的优选实施例,并不用于限制本发明,显然,本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样,倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内,则本发明也意图包含这些改动和变型在内。

Claims (3)

1.一种多源摆动动态CT成像方法,其特征在于:该方法以基于先验图像约束的压缩感知图像重建算法为基础,具体包括以下步骤为:
步骤一、把采集的多对X射线源-探测器单次摆动的所有投影数据进行重排;
步骤二、对重排后的投影数据采用FBP算法重建出CT图像;
步骤三、把上述CT图像作为先验图像,通过带TV约束的ART算法把相应的时间帧的投影数据进行图像重建获得对应时间帧的CT图像;
步骤四、把所有时间帧的CT图像按时间顺序组合获得该次X射线源-探测器摆动时间段内检测对象的动态CT图像;
步骤五、重复步骤三、四,最后获得整个检测对象变化过程的动态CT图像。
2.根据权利要求1所述的多源摆动动态CT成像方法,其特征在于:基于先验图像约束压缩感知图像重建算法被设计为以下无约束最小化问题:
min κ||TV(I)||1+(1-κ)||TV(I-Ip)||1such that AI=P (1)
其中,||·||1表示矩阵L1范数,A代表系统矩阵,I,P分别表示待重建图像和测量的投影数据,TV(I),TV(I-Ip)分别表示待重建图像以及待重建图像与先验图像的总变分,κ是由先验图像Ip中的伪影程度决定的经验参数,取值范围为[0,1]。
3.根据权利要求2所述的多源摆动动态CT成像方法,其特征在于:通过采用梯度下降搜索方法对公式(1)进行求解,图像中(i,j)点的梯度方向hi,j可由分量ai,j,bi,j,ci,j,di,j和ei,j的计算得到:
其中,I(k) i,j,I(k) i-1,j,I(k) i,j-1分别表示在第k次迭代获得的图像中位置为(i,j),(i-1,j)以及(i,j-1)的灰度值,
其中,I(k) i+1,j,I(k) i-1,j-1分别表示在第k次迭代获得的图像中位置为(i+1,j)和(i-1,j-1)的灰度值,
其中,I(k) i,j+1,I(k) i-1,j+1分别表示在第k次迭代获得的图像中位置为(i,j+1)和(i-1,j+1)的灰度值,
其中,I(k)表示在第k次迭代获得的图像
其中,Ip i,j表示先验图像中位置为(i,j)对应的图像值,
总梯度表示为
β=max(|κI(k) i,j+(1-κ)(I(k) i,j-Ip i,j)|)/max(|hi,j|)
Ii,j (k+1)=Ii,j (k)-σ×β×hi,j
其中,β表示当前图像梯度的最大值以及先验图像与重建图像靠近的程度,
σ=σ×σs
σs表示梯度下降的调整步长。
CN201710233458.5A 2017-04-11 2017-04-11 多源摆动动态ct成像方法 Active CN107016709B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710233458.5A CN107016709B (zh) 2017-04-11 2017-04-11 多源摆动动态ct成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710233458.5A CN107016709B (zh) 2017-04-11 2017-04-11 多源摆动动态ct成像方法

Publications (2)

Publication Number Publication Date
CN107016709A true CN107016709A (zh) 2017-08-04
CN107016709B CN107016709B (zh) 2020-06-19

Family

ID=59445261

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710233458.5A Active CN107016709B (zh) 2017-04-11 2017-04-11 多源摆动动态ct成像方法

Country Status (1)

Country Link
CN (1) CN107016709B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110827370A (zh) * 2019-11-09 2020-02-21 中北大学 一种非等厚构件的多能ct循环迭代重建方法
CN111552002A (zh) * 2020-05-19 2020-08-18 重庆大学 一种三源摆动螺旋安检ct成像装置及方法
CN115097535A (zh) * 2021-07-07 2022-09-23 清华大学 检查系统和方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103136773A (zh) * 2013-02-05 2013-06-05 南方医科大学 一种稀疏角度x射线ct成像方法
CN104167007A (zh) * 2013-05-17 2014-11-26 上海联影医疗科技有限公司 基于部分扫描的ct图像重建方法、装置及ct设备
US8989469B2 (en) * 2010-12-20 2015-03-24 The Board Of Trustees Of The Leland Stanford Junior University Systems and methods for simultaneous acquisition of scatter and image projection data in computed tomography

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8989469B2 (en) * 2010-12-20 2015-03-24 The Board Of Trustees Of The Leland Stanford Junior University Systems and methods for simultaneous acquisition of scatter and image projection data in computed tomography
CN103136773A (zh) * 2013-02-05 2013-06-05 南方医科大学 一种稀疏角度x射线ct成像方法
CN104167007A (zh) * 2013-05-17 2014-11-26 上海联影医疗科技有限公司 基于部分扫描的ct图像重建方法、装置及ct设备

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
GUANG-HONG CHEN 等: "Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets", 《NIH PUBLIC ACCESS AUTHOR MANUSCRIPT》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110827370A (zh) * 2019-11-09 2020-02-21 中北大学 一种非等厚构件的多能ct循环迭代重建方法
CN110827370B (zh) * 2019-11-09 2023-06-06 中北大学 一种非等厚构件的多能ct循环迭代重建方法
CN111552002A (zh) * 2020-05-19 2020-08-18 重庆大学 一种三源摆动螺旋安检ct成像装置及方法
CN115097535A (zh) * 2021-07-07 2022-09-23 清华大学 检查系统和方法

Also Published As

Publication number Publication date
CN107016709B (zh) 2020-06-19

Similar Documents

Publication Publication Date Title
JP4675869B2 (ja) 結像システム
JP4668161B2 (ja) 直線軌跡走査を採用した画像復元システム及び方法
US6028910A (en) High resolution areal tomosynthesis
US7215731B1 (en) Fast backprojection/reprojection with hexagonal segmentation of image
CN100464707C (zh) 三维锥束ct图像重建的处理系统
US20050078861A1 (en) Tomographic system and method for iteratively processing two-dimensional image data for reconstructing three-dimensional image data
JPH08117220A (ja) X線透視撮影方法およびx線装置
JP2016002467A (ja) X線ct装置及びその制御方法
JPH0728862B2 (ja) Ct装置
CN107016709A (zh) 多源摆动动态ct成像方法
CN102004111A (zh) 一种倾斜多锥束直线轨迹ct成像方法
Misawa et al. Ultra-fast x-ray tomography for multi-phase flow interface dynamic studies
CN100581471C (zh) 用于检查周期性运动的对象的ct方法
US7170966B2 (en) Practical implementation of a CT cone beam algorithm for 3-D image reconstruction as applied to nondestructive inspection of baggage, live laboratory animal and any solid materials
WO2005010824A1 (en) A fourier tomographic image reconstruction method for fan-beam data
CN106228601A (zh) 基于小波变换的多尺度锥束ct图像快速三维重建方法
CN107796835A (zh) 一种x射线柱面三维锥束计算机层析成像方法及装置
Lau et al. Ultrafast X-ray tomographic imaging of multiphase flow in bubble columns-Part 1: Image processing and reconstruction comparison
JPH0943354A (ja) 物体の三次元画像の再構成方法
JP2007529258A (ja) 複数の焦点取得方法及び装置
Zhou et al. Computerized tomography
Wang et al. Sparse-view cone-beam CT reconstruction by bar-by-bar neural FDK algorithm
Kaestner et al. Recent developments in neutron imaging with applications for porous media research
Zou et al. Dual helical cone-beam CT for inspecting large object
JP5053389B2 (ja) ラドンデータに基づくイメージ関数を処理する方法並びに撮像方法

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