基于频域分析的ct图像平移运动伪影校正方法
【专利摘要】本发明属于CT图像处理技术，尤其涉及一种基于频域分析的CT图像平移运动伪影校正方法，该方法包括：获取受平移运动影响的扇束CT图像投影数据，对投影数据进行重排得到平行束模式的投影数据；对平行束模式投影数据进行一维傅里叶变换；计算相位相关函数，利用相位相关函数求得相对位移量；利用相对位移量与平移运动参数之间的变换关系，计算获得平移运动参数；对受平移运动影响的扇束CT图像投影数据利用滤波反投影算法进行图像重建，在图像重建中加入计算得到的平移运动参数进行补偿，校正运动伪影。本发明的有益效果为：本发明解决了现有技术中对于复杂运动曲线估计效果较差的问题，能够实时准确地估计运动参数，有效校正运动伪影，且具有鲁棒性。
【专利说明】
基于频域分析的CT图像平移运动伪影校正方法
技术领域
[0001] 本发明属于CT图像处理技术，尤其涉及一种基于频域分析的CT图像平移运动伪影 校正方法。
【背景技术】
[0002] 计算机断层成像(Computed Tomography，CT)自20世纪70年代问世以来，在人类疾 病诊断和治疗方面取得巨大成就，是当今最先进的医学成像技术之一，广泛应用于临床诊 断和科学研究中，目前CT主要采用扇束扫描模式。传统的CT理论均是在假设患者静止不动 的前提下研究的，难以满足实际情况。CT成像需要一定时间，病人会发生自主性或非自主性 的运动，导致重建图像存在运动伪影，影响CT成像质量。通过大量临床实验，得出病人头部 的平均平移指数(Mean Translation Excursion，MTE)平均为2 · 5mm〇由于CT扫描需要相对 较高的辐射剂量，若因运动而重复扫描是不可取的。因此，研究运动伪影的校正方法具有重 要的理论意义和临床价值。
[0003] 常见的运动伪影校正方法大致分为三类，一类通过增加排数、X射线源和探测器的 旋转速度，加快扫描。但这种方法会极大增加硬件设计难度和硬件成本，同时具有降低图像 的信噪比（SNR)的缺点。且病人还有可能在扫描过程中运动，只能通过重复扫描修复，所以 不是一个完整的解决方案。第二类使用外部数据跟踪系统得到刚体运动数据，进行运动伪 影校正，但需要额外增加设备。第三类是基于运动参数估计的，采用高阶多项式曲线模型描 述病人的平移运动，寻找运动参数与投影数据间的关系，实现运动参数估计，然后在图像重 建过程中进行运动补偿，最经典的方法是利用数据一致性原理(HLCC)的估计方法挖掘投影 数据内在关联估计出运动参数，对于单纯的平移运动和简单的运动曲线估计精度高，但对 于复杂运动曲线估计效果较差。
[0004] 滤波反投影算法是一种常用的图像重建方法，它基于傅里叶切片定理，该定理能 够将投影数据的一维傅里叶变换与重建图像的二维傅里叶变换联系起来，对于每一个投影 数据，首先获得滤波后的投影数据，然后通过求解逆傅里叶变换以获得重建图像。
[0005] 中国专利CN102521801中提出一种CT图像的弧形伪影校正方法，以及中国专利 CN102768759中提出一种术中CT图像射束硬化伪影校正方法，以及中国专利CN103310432中 提出一种基于四阶全变分流的CT图像归一化的金属伪影校正的方法，分别是对弧形伪影、 射束硬化伪影和金属伪影的校正方法，不是对运动伪影进行的校正方法，目前国内专利中 对运动伪影校正方法的研究很少，存在很大的研究空间。
【发明内容】
[0006] 为了解决上述技术问题，本发明克服已有基于运动参数估计的伪影校正技术中存 在的不足，针对复杂运动曲线估计效果较差的情况，提供了一种基于频域分析的CT图像平 移运动伪影校正方法，在投影数据含噪且运动曲线复杂情况下，本发明仍可以准确估计运 动参数，有效校正运动伪影。
[0007] 一种基于频域分析的CT图像平移运动伪影校正方法，其特征在于按照下述步骤进 行：
[0008] 步骤1:获取受平移运动影响的扇束CT图像投影数据gF(i3，u)，其中下标F表示扇束 模式，β为扇束模式的扫描视角βε[0,2π)， u为探测器测得任意光束与射线源中线的距离；
[0009] 步骤2:对投影数据gF(i3，u)进行重排得到平行束模式的投影数据gp(0，u)，其中下 标P表示平行束模式，θ为重排后平行束模式的扫描视角；
[0010] 步骤3:对平行束模式投影数据gP(0，u)进行一维傅里叶变换，得到FdgpO，!!）]，将 Fi[gp(9，u)]分为两部分Fi[gp(0 lk，u)]和Fi[gp(02k，u)]，其中，0lk和02k均为重排后平行束模 式的扫描视角，只是大小范围不同，
和S2k的长度，γ为任意光束 与射线源中线的夹角，Ymax为γ的最大张角；
[0011] 步骤4:利用Fi[gp(0lk，u)]和?1|^(021{，11)]之间的变换关系计算获得相位相关函 数，进而求得相对位移量；
[0012] 步骤5:利用相对位移量与平移运动参数之间的变换关系，计算获得平移运动参 数；
[0013] 步骤6:利用步骤1获取的投影数据gF(i3，u)由滤波反投影算法进行图像重建，在图 像重建中加入步骤5计算得到的平移运动参数进行运动补偿，校正运动伪影。
[0014] 进一步地，所述步骤2中：
[0015] 利用线性插值方法将扇束模式投影数据gF(i3，u)重排转换为平行束模式的投影数 据gp(9，U)，具体重排过程为：gp(0，u) =ClgF(01，U)+C2gF(02，U)，其中，C1，C2为插值系数， /?，/t e /? (= [0,2瓦)，.Θ e [ γ max，2π- γ max]。
[0016] 进一步地，所述步骤3中：
[0017]扫描过程中被须_发生平移运动为d=(dx，dy)时，平行束模式的投影数据g P(0，u) 的一维傅里叶变换为：
[0018]
[0019] 其中，ω为u变换到频域对应的角频率，5(xcos9+ysin9-u)为变量xcos9+ysin9-u 的狄拉克函数，f(X，y)表示被测物的图像函数，FU，r〇是f(x，y)的象函数，ξ，?!分别为x，y变 换到频域对应的角频率；
[0020:
[0021] 其中，dx(tlk)，dy(tlk)分别表示扫描视角为0 lk时被测物在X、y方向的位移量，
Ω表示CT射线源旋转角速度，dx(t 2k)，dy(t2k)分别表示扫描视角为02k时被测物在 X、y方向的位移量，&=·^。
[0022] 进一步地，所述步骤4中：
[0023] Fitgp(Q1I^u)]和Fi [gp( 02k，u)]有如下关系
，贝丨jFi[gp(9lk，u)]^PIFi[gp (921{，11)]的互功率谱为：
[0024]
[0025] 计算互功率谱一维傅里叶逆变换，得到相位相关函数为 ω )]=3(u_u〇k)，其中，u〇k为相对位移量
UokG [u0i，UQ2，K，U()m]，k = 1，2，1(，111，5(11-11他)为变量11-11他在任意视角911{时的狄拉克函数，3 (u-UQk)对任意视角0lk在UQk处均不为零，在除UQk以外其它位置均为零，则根据相位相关函数 的峰值位置求取相对位移量
[0026] 进一步地，由于相对位移量受采样率影响，通过搜索得到的最大值可能并不能定 位正确的峰值，从而引起估计误差，相对位移量u 〇 k的精确估计u 〇 7k为
[0027] 进一步地，所述步骤5中：
[0028] 平移运动函数dx(t)和dy(t)通过基本函数奶(/).炉:(/),!< ,炉、(/)展开表示，即：
[0029]
[0030] 其中，Pij和P2j( j = l，2，K，N)为平移运动参数，选取= ^ 同时考虑初始时刻 被测物是静止的，所以dx(0)=0，dy(0)=0,则平移运动函数为：
[0031]
[0032] 将平移运动函数代. 方程：
[0034] 方程组包含2N个未知数Plj和?2』〇 = 1，2，1(少），111个独立的方程，实际中111要远大于 2N，所以平移运动参数可以由下式确定：
[0035] P=(HtH)^1(HtU);
[0036] 其中，11=(1/01,11(/2, Λ，U(/m)T，P=(pii，pi2, Λ，Ρ1Ν，Ρ21，Ρ22, Λ，P2N)T，
[0038] 进一步地，所述步骤6中：
[0039] 将视角β下的投影数据gF(i3，u)先滤波，后利用步骤5估计得到的平移运动参数进 行坐标系平彳
累加到图像函数上，直至最后一个扫描视角。
[0040] 本发明的有益效果是：（1)首次提出基于频域分析的CT图像平移运动伪影校正方 法，利用插值重排将扇束投影变为平行束投影，可以将平行束投影数据频域分析方法拓展 到扇束模式，然后根据运动参数与相对位移量的关系进行估计，利用修正的相对位移量公 式可以提高参数估计精度;最后在重建过程中改进滤波反投影算法，加入估计出的运动参 数进行补偿，有效校正运动伪影。（2)通过实验证明了本发明的实时准确性和鲁棒性。在运 动曲线复杂情况下，仍可以准确估计运动参数，有效校正运动伪影。（3)本发明的方法简单 可行，可以实时准确地估计运动参数，有效校正运动伪影。
【附图说明】
[0041] 图1是CT成像坐标系；其中，图1(a)是扇形光束成像坐标系，图1(b)是平行光束成 像坐标系；
[0042] 图2是投影数据重排示意图；其中，图2(a)是平行光束投影数据网格图，图2(b)是 扇形光束投影数据网格图；
[0043]图3是模拟头部模型的示意图；
[0044] 图4为本发明实施例中4组实验运动参数估计结果；
[0045] 其中，图4(a)是第I组实验真实运动曲线和估计运动曲线，图4(b)是第II组实验真 实运动曲线和估计运动曲线，图4(c)是第III组实验真实运动曲线和估计运动曲线，图4(d) 是第IV组实验真实运动曲线和估计运动曲线；
[0046] 图5为本发明实施例中4组实验伪影校正结果；
[0047] 其中，图5(a)是第I组实验原始图像、伪影图像和本发明校正后图像，图5(b)是第 II组实验原始图像、伪影图像和本发明校正后图像，图5(c)是第III组实验原始图像、伪影 图像和本发明校正后图像，图5(d)是第IV组实验原始图像、伪影图像和本发明校正后图像。
【具体实施方式】
[0048] 下面结合附图和实施例对发明的【具体实施方式】做出说明。
[0049] 本发明提供的基于频域分析的CT图像平移运动伪影校正方法包括下列步骤：
[0050] 步骤1、获取受平移运动影响的扇束CT图像投影数据gF(i3，u)，其中下标F表示扇束 模式，β为扇束模式的扫描视角βε[0,2π)， u为探测器测得任意光束与射线源中线的距离， 如图1(a)所示；
[0051]步骤2:对投影数据gF(i3，u)进行重排得到平行束模式的投影数据gP(0， u)，其中下 标P表示平行束模式，Θ为重排后平行束模式的扫描视角，如图I (b)所示；
[0052] 其中，步骤2中重排方法具体包括：
[0053] 利用线性插值方法将扇束模式投影数据gF(i3，u)重排转换为平行束模式的投影数 据gp(9，U)，具体重排过程为：gp(0，u) =ClgF(01，U)+C2gF(02，U)，其中，C1，C2为插值系数， /?丨，/；?2￡/?(=[0,2/了) ;为了防止插值数据缺失，0￡[丫繼，231-丫111狀]，丫为任意光束与射线源 中线的夹角，Ymx为γ的最大张角，如图2(a) (b)所示。
[0054] 需要说明的是，重排得到的投影数据gP(0，u)是相邻原投影数据的加权和，其运动 信息近似表示为：
[0055] d[t(9，u) ] = cid[t(0i，u)]+C2d[t(&，u)]
[005?]其中，?(βι,ιι)和t(fe，u)分别表不扇束模式扫描视角βι和(62的扫描时间，t(0，u)表 示新的平行束视角Θ的扫描时间，/= $，Ω表示CT射线源旋转角速度。
[0057]进一步说明，因物体在扫描过程中发生运动，每个邻近点均含有不同的运动信息， 插值后可能会使其受到破坏。但由于运动的连续性，相邻的两次投影运动位移接近，线性插 值方式使得运动信息损坏较少，且仅与采样点时间有关。因此，较短的时间差内，可以认为 平移运动变化不大，从而忽略gp(9，u)中各Θ视角下投影数据运动信息受到的破坏。
[0058]步骤3:对平行束模式投影数据gP(0，u)进行一维傅里叶变换，得到FjgpO，!!)]，将 Fi[gp(9，u)]分为两部分Fi[gp(0lk，u)]和Fi[gp(02k，u)]，其中，0lk和02k均为重排后平行束模
与射线源中线的夹角，Ymax为γ的最大张角；
[0059]步骤3中，扫描过程中被测物发生平移运动为d=(dx，dy)时，平行束模式的投影数 据gP(9，u)的一维傅里叶变换为：
[0060]
[0061J 其中，ω为U变换到频域对应的角频率，6(xcose+ysine-u)为变量xcose+ysine-u 的狄拉克函数，f(X，y)表示被测物的图像函数，FU，r0是f(x，y)的象函数，ξ，?!分别为x，y变 换到频域对应的角频率；
[0063] 其中，dx(tlk)，dy(tlk)分别表示扫描视角为G lk时被测物在x、y方向的位移量，
Ω表示CT射线源旋转角速度，dx(t 2k)，dy(t2k)分别表示扫描视角为02k时被测物在 x、y方向的位移量：
[0064] 步骤4:利用？1匕4011{，11)]和?1匕[ 5(021{，11)]之间的变换关系计算获得相位相关函 数，进而求得相对位移量；
[0065] 其中，步骤4中计算相对位移量的方法具体包括：
[0066] MgpOm)]和 MgP(02k，u)]有如下关系
贝 IjF1 [ gP (Θ lk，u)]和?丄[gP (Θ 2k，u)]的互功率谱为：
[0067]
[0068]计算互功率谱一维傅里叶逆变换，得到相位相关函数为ω )] = S(U-UQk)，其中，UQk为相对位移量
[u0i，UQ2，K，u0m]，k = l，2，K，m，5( u-UQk)为变量U-Uok在任意视角 Θ ik 时的狄拉克函数，δ (u-u〇k) 对任意视角9lk在uQk处均不为零，在除uQk以外其它位置均为零，则根据相位相关函数的峰值 位置求取相对位移量.
[0069]需要说明的是，相对位移量受采样率影响，通过搜索得到的最大值可能并不能定 位正确的峰值，从而引起估计误差，利用公；E :青确估 计相对位移量UOk。
[0070]步骤5、利用相对位移量与平移运动参数之间的变换关系，计算获得平移运动参 数；
[0071 ]其中，步骤5中计算平移运动参数的方法具体包括：
[0072] 平移运动函数dx(t)和dy(t)通过基本函数於別,釣(?),Κ,%(〇展开表示，即：
[0073]
[0074] 其中，Plj和p2j (j = 1，2，K，Ν)为平移运动参数，选取a (/) = / ν ,同时考虑初始时刻 被测物是静止的，所以dx(0)=0，dy(0)=0,则平移运动函数为：
[0075]
[0076] 将平移运动函数代Λ 得到m个 方程：
[0077]
[0078] 方程组包含2N个未知数Plj和p2j(j = 1，2，K，N)，m个独立的方程，实际中m要远大于 2N，所以平移运动参数可以由下式确定：
[0079] P=(HtH)^1(HtU);
[0080] 其中，U= (U7 01，U(/ 2, Λ，U(/ m)T，P= (Ρ11，Ρ12, Λ，Ρ1Ν，Ρ21，Ρ22, Λ，P2N)T，
[0081]
[0082]步骤6、利用步骤I获取的投影数据gF(i3，u)由滤波反投影算法进行图像重建，在图 像重建中加入步骤5计算得到的平移运动参数进行运动补偿，校正运动伪影。。
[0083]其中，步骤6中图像重建方法具体包括:将视角β下的投影数据gF(i3，u)先滤波，后 利用步骤5估计得到的平移运动参数进行坐标系平彳3累加到图像函数 上，直至最后一个扫描视角。
[0084] 为了验证本发明的有效性，做了充分的实验验证，实验结果表明在投影数据含噪 且运动曲线复杂情况下，本发明仍可以准确估计运动参数，有效校正运动伪影。
[0085] 实施例：为了验证本发明CT图像平移运动伪影校正效果，本实施例采用模拟头部 模型进行仿真实验，如图3所示。模拟头部模型采用shepp-logan头部模型模拟人体活动，在 模型上增加多个更小的圆，用于模拟内耳等细微结构。
[0086] 仿真参数设置:有效视野半径为128mm，扫描起始时间t = 0，扫描时间tm= 1。射线 源沿着半径R = 500mm的圆周以Ω = 2Jirad/S的角速度旋转，在2JT范围内获取720组投影，探 测器单元为1000个，探测器孔径Um=0.3mm，Ymax = O. 184rad/s，投影数据加入Poisson分布 噪声，射线源射向每个探测器的光子数1〇6。
[0087] 评价指标:平均平移指数(Mean Translation Excursion，MTE)、平均补偿平移偏 离指数(Compensative Mean Translation Excursion，CMTE)和平均相对平移偏离指数 (Relative Mean Translation Excursion，RMTE)，艮P
[0088]
[0089]
[0090] RMTE = CMTE/MTE；
[0091] 其中，(1(〇和|@分别为真实平移运动和估计平移运动位移大小。MTE越大说明病 人运动幅度越大，图像伪影越严重;MTE为零时，CMTE越小估计结果越精确;MTE不为零时， CMTE越小即RMTE越小，估计结果越精确。
[0092] 本发明以4组不同的运动曲线为例，第I组无运动，第11组为简单运动，第III、IV组 为复杂运动。计算平移运动参数估计结果的准确性，并记录运行时间，每组运动进行30次仿 真实验取平均值，如表1所示。
[0093] 衷1参数估计结果和运行时间
[0095] 由表1可知，第I组实验表明，在扫描过程不存在运动的情况下，本发明方法不会引 入误差;第II组实验表明，在运动曲线简单的情况下，本发明方法可以准确估计运动参数； 第m、 IV组实验表明，在运动曲线复杂情况下，本发明方法估计结果稳定，且仍可以估计运 动参数。4组实验的运行时间都很短，表明本发明方法实时性好。
[0096] 将表1中4组实验的真实运动曲线和估计运动曲线进行比较，如图4所示。由图4(a) 可知，在第I组无运动实验中，本发明方法估计曲线与原曲线完全重合；由图4(b)可知，在第 II组简单运动曲线实验中，本发明方法估计得到曲线几乎与原曲线重合；由图4(c)、4(d)可 知，在第III、IV组复杂运动曲线实验中，本发明方法估计结果与原曲线十分接近;综合4组 实验结果，本发明在投影数据含噪且运动曲线复杂情况下，仍能准确估计运动参数，表明本 发明方法具有鲁棒性。
[0097] 将表1中4组实验伪影校正后的结果和原始图像进行比较，如图5所示。由图5(a)、5 (b)可知，在第I、II组无运动或简单运动曲线实验中，本发明方法均有很好的校正效果，校 正后图像十分接近原始图像；由图5(c)、5(d)可知，在第III、IV组复杂运动曲线实验中，本 发明方法运动参数估计准确，因此本发明方法仍有很好的校正效果，内部细节清晰；综合4 组实验伪影校正结果，表明本发明可以有效校正CT图像平移运动伪影。
[0098] 以上对本发明的实例进行了详细说明，但所述内容仅为本发明的较佳实施例，不 能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均 应仍归属于本发明的专利涵盖范围之内。
【主权项】
1. 一种基于频域分析的CT图像平移运动伪影校正方法，其特征在于包括以下步骤： 步骤1:获取受平移运动影响的扇束CT图像投影数据gF(i3，u)，其中下标F表示扇束模式， β为扇束模式的扫描视角βε[〇,2π)，u为探测器测得任意光束与射线源中线的距离； 步骤2:对投影数据gF(i3，u)进行重排得到平行束模式的投影数据gP(0，u)，其中下标P表 示平行束模式，Θ为重排后平行束模式的扫描视角； 步骤3:对平行束模式投影数据gP(0，u)进行一维傅里叶变换，得到 (Θ，u)]分为两部分Fl[gp(0lk，u)]和Fl[gP(02k，U)]，其中，0lk和02k均为重排后平行束模式的 扫描视角，只是大小范围不同，化e [义,6*n, K , 3 ;" ] (= ，心e低爲,K :,02m ] .G +，知_ - ，k = 1，2，K，m，m为θ 1 k和θ 2 k的长度，丫为任意光束与射线源中线的夹 角，Tmax为γ的最大张角； 步骤4:利用？1|^(011<，11)]和?1|^(021{，11)]之间的变换关系计算获得相位相关函数，进 而求得相对位移量； 步骤5:利用相对位移量与平移运动参数之间的变换关系，计算获得平移运动参数； 步骤6:利用步骤1获取的投影数据gF(i3，u)由滤波反投影算法进行图像重建，在图像重 建中加入步骤5计算得到的平移运动参数进行运动补偿，校正运动伪影。2. 根据权利要求1所述的CT图像平移运动伪影校正方法，其特征在于，所述步骤2中：利 用线性插值方法将扇束模式投影数据gF(i3，u)重排转换为平行束模式的投影数据g P(0，u)， 具体重排过程为：gP(θ，U)=ClgF(βl，U)+C2gF(β2，U)，其中，Cl，C2为插值系数， 々? e e [(Χ2;τ), Θ e [ γ max，231- γ max]。3. 根据权利要求1所述的CT图像平移运动伪影校正方法，其特征在于，所述步骤3中：扫 描过程中被测物发生平移运动为d=(d x，dy)时，平行束模式的投影数据gP(0，u)的一维傅里 叶变换为：其中，ω为u变换到频域对应的角频率，5(xcos9+ysin9-u)为变量xcos9+ysin9-u的狄 拉克函数，f(x，y)表示被测物的图像函数，FU，r〇是f(x，y)的象函数，ξ，τ?分别为x，y变换到 频域对应的角频率；其中，dx(tlk)，dy(tlk)分别表示扫描视角为9 lk时被测物在x、y方向的位移量，Ω 表不CT射线源旋转角速度，dx(t2k)，dy(t2k)分别表不扫描视角为92k时被测物在x、y方向的位 移量4. 根据权利要求1所述的方法，其特征在于，所述步骤4中：Fi [gp(Θlk，u)]和Fi [gp(Θ2k， u )]有如下关系 5 [办(九，《)] = e-(~丨乂"+|i/' (0lk，u)]和MgpO^u)]的互功率谱为：计算互功率谱一维傅里叶逆变换，得到相位相关函数为 (U-U〇k)，其中，UOk为相对位移量，"w =[<(/丨.t )-c/,(i1;.)]cos 九 +[<(/4)^ [U〇l，UQ2，K，UQm]，k=l，2，K，m，5( U-UQk )为变量U-UQk在任意视角 Θ lk 时的狄拉克函数，δ (u-u〇k) 对任意视角0lk在UQk处均不为零，在除UQk以外其它位置均为零，则根据相位相关函数的峰值 位置求取相对位移量为= u5. 如权利要求4所述的方法，其特征在于，相对位移量uQk的精确估计V (^为6. 如权利要求1所述的方法，其特征在于，所述步骤5中，平移运动函数dx(t)和dy(t)通 过基本函数仍(?),朽(i),L.(；).展开表示，即：其中，pij和P2j(j = 1，2，Κ，N)为平移运动参数，选取炉? = ，.同时考虑初始时刻被测物 是静止的，所以dx(0)=0，dy(0)=0,则平移运动函数为：方程组包含2N个未知数Plj和P2j(j = 1，2，K，N)，m个独立的方程，实际中m要远大于2N，所 以平移运动参数可以由下式确定： P=(HTH)-WU); 其中，11=(1^ 01,1/02，Λ Λ，Ρ1Ν，Ρ21，Ρ22, Λ，P2N)T，7.如权利要求1所述的方法，其特征在于，所述步骤6中：将视角β下的投影数据gF(e，u) 先滤波，后利用步骤5估计得到的平移运动参数进行坐标系平移'累加到 图像函数上，直至最后一个扫描视角。
【文档编号】G06T11/00GK106056645SQ201610352005
【公开日】2016年10月26日
【申请日】2016年5月25日
【发明人】孙云山, 张立毅, 张媛, 费腾, 张勇, 陈雷
【申请人】天津商业大学