一种圆轨道下不完全CT投影数据的恢复方法与流程

本发明涉及一种医学图像处理，具体涉及图像的增强或复原，特别是涉及一种圆轨道下不完全ct投影数据的恢复方法。
背景技术：
：锥形束计算机断层扫描是一种可以快速且低辐射剂量获得人体的三维解剖图像的医学成像方法。数字平面探测器的出现使得真正的锥形束计算机断层扫描(cone-beamcomputedtomography,cbct)在临床应用成为可能。由于cbct采用锥形束扫描结构，所采集的投影数据具有高度冗余性，想要利用数据的冗余性来达到高效率的目的就必须满足一致性条件(consistencyconditions)。约翰于1938年从数学的角度推导出了一个超双曲线偏微分方程(ultrahyperbolicpartialdifferentialequation,pde)作为线积分的一致性条件，故此方程称为约翰方程。在70年代，随着ct的出现和兴起，约翰方程(john'sequation)被广泛地应用在ct领域并进行了更深入地推导和延伸。时至今日，已经有许多学者提出了多种基于约翰方程的一致性条件，如patch提出一种基于第三代螺旋ct的一致性条件(sarahk.patch,computationofunmeasuredthird-generationvctviewsfrommeasuredviews,2002,ieeetrans.med.imaging)，在此文献中，作者通过使用第三代螺旋ct参数来替换基于笛卡尔坐标系的约翰方程中的参数，并在此基础上进行了一系列地深入推导，将二阶偏微分方程转换成了一组一阶常微分方程(first-orderordinarydifferentialequation)，对此方程进行傅里叶变换，然后应用其频率域的形式可以从获得的投影中计算出未测量的投影数据。根据此理论，文献中使用“旋转+螺旋+旋转”的射线源轨道，螺距为0.128m，每转一圈采集984个投影，然后针对获得的投影数据使用推导方程对应的傅立叶变换形式进行外插，从而获得未测量的投影数据。但是patch推导出的一致性条件是基于螺旋ct的几何结构推导而成，螺旋结构要求变量z沿着纵轴在不停地变化，这会在公式中引入一个变化量，使得公式计算复杂，解决问题困难。文献中由于采用“旋转+螺旋+旋转”扫描方式，需要在“螺旋”这个扫描阶段前后分别进行工作台移动速度的增加和减小，以保证在这三个扫描阶段数据采集的快速转换和平滑过渡，但是在实际应用中，一般螺旋ct的工作台移动速度都是恒定的，因此，文献中的一致性条件很难具体应用。此外，haoyan等人在使用移动的波束停止阵列(beamstoparray,bsa)来进行散射校正时采用了patch基于三代螺旋ct的一致性条件(haoyan,xuanqinmou,shaojietang,qiongxuandmariazankl,projectionmovingbeamstoparray,2010,phys.med.biol.)。波束停止阵列是一个遮挡板，上面有规律地排列着遮挡点。x射线经过遮挡板照射到物体上，其中经过遮挡板上遮挡点的x射线会被屏蔽掉。如果在探测器遮挡点对应位置出现数值，则这个数值是散射值。由于使用了bsa，得到的投影数据是不完整的，因此需要使用一致性条件对其进行恢复。此文献采用patch基于三代螺旋ct的推导公式[patch(2002)]，并结合空间插值方法，提出一种pc_si方法用于散射校正领域中。但如前面提到的那样，由于螺旋ct的z轴在不停地变化，导致此公式很难具体应用，所以此文献参照defrise等人的做法，用探测器的纵坐标来近似代替z坐标，但是这种近似做法，会导致精度的下降。从几何结构上来看，只需要令z＝0，螺旋ct就可以退化成圆轨道cbct，但是文献中的公式是对z求导，当z为常数时导数没有意义，所以此文献中所用的方法并不适用于圆轨道cbct。因此，针对现有技术不足，提供一种圆轨道下不完全ct投影数据的恢复方法以克服现有技术不足甚为必要。技术实现要素：本发明所要解决的技术问题是提供一种圆轨道下不完全ct投影数据的恢复方法，和传统的空间插值方法相比，本发明的方法能够很好地恢复丢失信息中的高频成分，进而大幅度减少图像中的条形伪影，显著提高图像质量。在圆轨道下使用patch的方法需要先使用探测器的纵坐标来近似代替z坐标，而本方法采用的一致性条件是基于圆轨道cbct几何结构推导得到的，能够避免因为近似z坐标而导致的精度下降，并且patch的方法中z轴在不停变化，使得此方法具体应用困难。与之相反，本发明方法只是一个简单的一阶求导公式，实现简单，具有可应用性。本发明采用如下技术方案：一种圆轨道下不完全ct投影数据的恢复方法，通过如下步骤进行：第一步，获得坏像素点位置信息；第二步，对每个投影中的坏像素点使用水平一维三次样条插值方法进行线性插值，以恢复坏像素点的低频信息，经过此步骤得到每个角度下的初始化投影；第三步，令n＝1，进入第四步，第四步，令m＝n，m为当前迭代次数；将每个角度下的初始化投影进行二维傅里叶变换，获得每个角度下的初始化投影在频率域中的频域数据；第五步，分别将前后两个相邻投影的频率域形式带入到一致性条件对应的公式中计算投影结果；第六步，对加权后的投影进行傅里叶反变换，得到第m次迭代终值；第七步，判断当前迭代次数m是否大于等于迭代终止次数s，如果是，则进入第九步；否则进入第八步；第八步，将第m次迭代终值作为初始化投影，令n＝n+1，进入第三步；第九步；以第m次迭代终值为最终恢复的数据结果。优选的，上述的圆轨道下不完全ct投影数据的恢复方法，第五步中使用的一致性条件，具体形式如下：其中，其中ρ是射线源到旋转中心的距离，d是旋转中心到探测器的距离，(α1,α2)是平板探测器坐标，θ是方位角。优选的，上述第五步具体是使用相同的权重因子ω对前后两个投影结果进行加权，得到加权更新后的投影结果。优选的，上述的圆轨道下不完全ct投影数据的恢复方法，第五步，具体使用一致性条件即公式(28)对每一个初始化投影的频域数据进行恢复，公式(28)的具体形式如下：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0……公式(28)；其中，u*(θ+dθ)和u*(θ)分别表示θ+dθ和θ角度下初始化投影的频域数据，u*(θ)为已知量，u*(θ+dθ)是要求的未知量，dθ代表相邻两个角度的间隔，k1、k2是平板探测器的频率域坐标，i代表复数的虚部；公式(28)的具体使用方法如下：首先将u*(θ+dθ)代入到公式(28)中计算u*(θ)：u*(θ)＝u*(θ+dθ)-dθ·u；然后用u*(θ-dθ)代替公式(28)中的u*(θ)，u*(θ)代替公式(28)中的u*(θ+dθ)计算得到：u*(θ)＝u*(θ-dθ)+dθ·u；由于选取前后两个相邻投影频域数据的任一个都可以实现投影数据恢复，为了平衡恢复效果，使用权重因子ω，ω∈[0,1]，对上述两个公式的计算结果进行加权平均，得到加权后的结果ur*(θ)，即:ur*(θ)＝ω(u*(θ-dθ)+dθ·u)+(1-ω)(u*(θ+dθ)-dθ·u),k2≠0；由于公式(28)只有在k2≠0时才有效，所以对频域数据u*(θ)中k2＝0的部分，直接使用θ角度下初始化投影频域数据中的相应部分来补偿；第六步，具体是对ur*(θ)进行傅里叶反变换，得到第m次迭代终值ur(θ)；ur(θ)＝f-1(ur*(θ))……(1)；其中公式(1)中ur(θ)表示θ角度下加权后的频域数据ur*(θ)的反傅里叶变换结果，f-1代表傅里叶反变换；第七步，具体是判断当前迭代次数m是否大于等于迭代终止次数s，如果是，则进入第九步；否则进入第八步；第八步，具体是将第m次迭代终值ur(θ)作为初始化投影，令n＝n+1，进入第三步；第九步，具体是以第m次迭代终值ur(θ)为最终恢复的数据结果。优选的，上述的圆轨道下不完全ct投影数据的恢复方法，第一步，获得坏像素点位置信息具体是将实际采集的投影数据转换成对应的弦图，即：弦图的x轴为投影的像素点位置，y轴为该像素位置对应的投影角度位置；对某个探测器单元存在损坏的条件下，则在弦图中对应的像素位置上会出现一条竖线，根据此理论在弦图中获得坏像素点的位置信息，并记录坏像素点位置。其中本方法中应用的公式(28)推导过程如下：由于本发明提出的一致性条件是基于约翰方程，所以先对约翰方程进行简单介绍和分析。fjohn于1938年在文献《theultrahyperbolicdifferentialequationwithfourindependentvariables》中提出一个超双曲线二阶偏微分方程，作为线积分的一致性条件，具体形式如下：公式中i,j为下标；u是函数f穿过ε和η两点的线积分：u(ε；η)＝∫f(ε+t(η-ε))dt(3)；以下为圆轨道下基于约翰方程的一致性条件推导步骤：首先采用圆轨道cbct坐标参数对ε和η进行参数化：其中ρ是射线源到旋转中心的距离，d是旋转中心到探测器的距离，(α1,α2)是平板探测器坐标，θ是方位角。分别求ε和η的每个分量求偏导：为书写方便，定义如下算子：把公式(6)-(11)代入到公式(2)中得到：由公式(13)-(16)对h1，h2,的定义，代入到公式(17)-(19)中，得到：从上述公式中，得出：对上面公式(21)移项合并同类项：求解公式(23)，最终推导结果如下：公式(24)即为我们推导的圆轨道下基于约翰方程的一致性条件。对公式(24)进行傅里叶变换，得到对应的频率域的表达形式如下：在上述公式中，uθ*表示u对θ求导后的傅里叶变换，(k1,k2)是(α1,α2)的对应频率域形式。为了方便表达和应用，将公式(25)右边式子用u代替：考虑到求导公式的一般定义，可以将公式(25)改下成如下形式：即：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0(28)；公式(28)是本发明改进后的一致性条件，本发明使用公式(28)进行投影数据的恢复。传统的空间插值技术使用同一个角度的投影的其他完整像素点数据来恢复丢失的数据，所使用重建信息只局限于丢失数据所在角度的投影，所以插值方法只能恢复丢失数据中的低频信息，对丢失的高频信息基本没有恢复。本发明方法应用的一致性条件，即不仅基于同一个投影中的数据来恢复丢失数据，还利用相邻角度投影中的对应位置数据信息。因此本发采用一致性条件可以很好地恢复丢失数据中的高频信息，能基本消除使用插值技术获得的重建图像中的条形伪影，显著提高图像质量。附图说明图1为实例中圆轨道cbct扫描几何示意图及参数示意图，图中，ρ是射线源到旋转中心的距离，d是旋转中心到探测器的距离，(α1,α2)是平板探测器坐标，θ是方位角。图2是投影数据对应的弦图，横坐标x轴为像素点位置，纵坐标y轴为投影角度。图3为本发明的数据恢复方法流程图。图4为三种方法的仿真结果示意图，图4(a)为理想图像，图4(b)为未经任何处理的原始图像，图4(c)为线性差值获得的结果，图4(d)、(e)、(f)为本发明方法一次迭代、二次迭代、三次迭代后的结果图，图(g)、(h)、(i)为patch方法一次迭代、二次迭代、三次迭代后的结果图。具体实施方式实施例1。一种圆轨道下不完全ct投影数据的恢复方法，通过如下步骤进行：第一步，获得坏像素点位置信息，将实际采集的投影数据转换成对应的弦图，即：弦图的x轴为投影的像素点位置，y轴为该像素位置对应的投影角度位置；对某个探测器单元存在损坏的条件下，则在弦图中对应的像素位置上会出现一条竖线，根据此理论在弦图中获得坏像素点的位置信息，并记录坏像素点位置。第二步，对每个投影中的坏像素点使用水平一维三次样条插值方法进行线性插值，以恢复坏像素点的低频信息，经过此步骤得到每个角度下的初始化投影。第三步，令n＝1，进入第四步，第四步，令m＝n，m为当前迭代次数；将每个角度下的初始化投影进行二维傅里叶变换，获得每个角度下的初始化投影在频率域中的频域数据；第五步，分别将前后两个相邻投影的频率域形式带入到一致性条件对应的公式中计算投影结果；具体是使用相同的权重因子ω对前后两个投影结果进行加权，得到加权更新后的投影结果。第五步中使用的一致性条件，形式如下：其中，其中ρ是射线源到旋转中心的距离，d是旋转中心到探测器的距离，(α1,α2)是平板探测器坐标，θ是方位角，如图1所示。第五步，可具体使用一致性条件即公式(28)对每一个初始化投影的频域数据进行恢复，公式(28)的具体形式如下：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0……公式(28)；其中，u*(θ+dθ)和u*(θ)分别表示θ+dθ和θ角度下初始化投影的频域数据，u*(θ)为已知量，u*(θ+dθ)是要求的未知量，dθ代表相邻两个角度的间隔，k1、k2是平板探测器的频率域坐标，i代表复数的虚部；公式(28)的具体使用方法如下：首先将u*(θ+dθ)代入到公式(28)中计算u*(θ)：u*(θ)＝u*(θ+dθ)-dθ·u；然后用u*(θ-dθ)代替公式(28)中的u*(θ)，u*(θ)代替公式(28)中的u*(θ+dθ)计算得到：u*(θ)＝u*(θ-dθ)+dθ·u；由于选取前后两个相邻投影频域数据的任一个都可以实现投影数据恢复，为了平衡恢复效果，使用权重因子ω，ω∈[0,1]，对上述两个公式的计算结果进行加权平均，得到加权后的结果ur*(θ)，即:ur*(θ)＝ω(u*(θ-dθ)+dθ·u)+(1-ω)(u*(θ+dθ)-dθ·u),k2≠0；由于公式(28)只有在k2≠0时才有效，所以对频域数据u*(θ)中k2＝0的部分，直接使用θ角度下初始化投影频域数据中的相应部分来补偿；第六步，具体是对ur*(θ)进行傅里叶反变换，得到第m次迭代终值ur(θ)；ur(θ)＝f-1(ur*(θ))……(1)；其中公式(1)中ur(θ)表示θ角度下加权后的频域数据ur*(θ)的反傅里叶变换结果，f-1代表傅里叶反变换。第七步，具体是判断当前迭代次数m是否大于等于迭代终止次数s，如果是，则进入第九步；否则进入第八步。第八步，具体是将第m次迭代终值ur(θ)作为初始化投影，令n＝n+1，进入第三步。第九步，具体是以第m次迭代终值ur(θ)为最终恢复的数据结果。其中本方法中应用的公式(28)推导过程如下：由于本发明提出的一致性条件是基于约翰方程，所以先对约翰方程进行简单介绍和分析。fjohn于1938年在文献《theultrahyperbolicdifferentialequationwithfourindependentvariables》中提出一个超双曲线二阶偏微分方程，作为线积分的一致性条件，具体形式如下：u是函数f穿过ε和η两点的线积分：u(ε；η)＝∫f(ε+t(η-ε))dt(3)；以下为圆轨道下基于约翰方程的一致性条件推导步骤：首先采用圆轨道cbct坐标参数对ε和η进行参数化：其中ρ是射线源到旋转中心的距离，d是旋转中心到探测器的距离，(α1,α2)是平板探测器坐标，θ是方位角。分别求ε和η的每个分量求偏导：为书写方便，定义如下算子：把公式(6)-(11)代入到公式(2)中得到：由公式(13)-(16)对h1，h2,的定义，代入到公式(17)-(19)中，得到：从上述公式中，得出：对上面公式(21)移项合并同类项：求解公式(23)，最终推导结果如下：公式(24)即为我们推导的圆轨道下基于约翰方程的一致性条件。对公式(24)进行傅里叶变换，得到对应的频率域的表达形式如下：在上述公式中，uθ*表示u对θ求导后的傅里叶变换，(k1,k2)是(α1,α2)的对应频率域形式。为了方便表达和应用，将公式(25)右边式子用u代替：考虑到求导公式的一般定义，可以将公式(25)改下成如下形式：即：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0(28)；公式(28)是本发明改进后的一致性条件，本发明使用公式(28)进行投影数据的恢复。传统的空间插值技术使用同一个角度的投影的其他完整像素点数据来恢复丢失的数据，所使用重建信息只局限于丢失数据所在角度的投影，所以插值方法只能恢复丢失数据中的低频信息，对丢失的高频信息基本没有恢复。本发明方法应用的一致性条件，即不仅基于同一个投影中的数据来恢复丢失数据，还利用相邻角度投影中的对应位置数据信息。因此本发采用一致性条件可以很好地恢复丢失数据中的高频信息，能基本消除使用插值技术获得的重建图像中的条形伪影，显著提高图像质量。和传统的空间插值方法相比，本发明的方法能够很好地恢复丢失信息中的高频成分，进而大幅度减少图像中的条形伪影，显著提高图像质量。在圆轨道下使用patch的方法需要先使用探测器的纵坐标来近似代替z坐标，而本方法采用的一致性条件是基于圆轨道cbct几何结构推导得到的，能够避免因为近似z坐标而导致的精度下降，并且patch的方法中z轴在不停变化，使得此方法具体应用困难。与之相反，本发明方法只是一个简单的一阶求导公式，实现简单，具有可应用性。实施例2。本实例分别使用传统空间插值方法、patch的方法(
背景技术：
中有详细说明)、以及本方法对一组不完整的投影数据进行恢复，不完整的投影数据由一台平板探测器具有坏像素点的cbct采集得到，其ct扫描参数具体如下：射线源到旋转中心的距离ρ为500mm，旋转中心到探测器的距离d为500mm，平板探测器大小为850×200mm，像素点尺寸为1×1mm，围绕物体旋转360度采集1080个投影。本发明技术的方法具体实施方式如下：第一步，获得坏像素点位置信息。将实际采集的投影数据转换成对应的弦图，如图2所示。图中黑色竖线所对应的横坐标即为坏像素点位置，通过此方式可以找到所有坏像素点的位置信息，并记录。第二步，对每一个投影中的坏像素点使用水平一维三次样条插值方法线性插值，并将插值后的结果作为该坏像素点的初始值。第三步，将插值后的投影以及前后两个相邻投影进行傅里叶变换，得到对应的频率域形式，以方便应用本发明推导出的一致性条件。第四步，将上一步得到的前后两个相邻投影的频率域形式代入到公式(28)中，恢复中间投影的丢失信息(尤其是高频成分信息)，并使用同一个权重因子ω对两个结果进行加权，由于没有更多的约束条件表明前后两个相邻投影哪个对中间投影的贡献更大，所以本实例令ω＝0.5。然后用计算加权后的结果修正该投影初始值。第五步，对修正后的值进行傅里叶反变换，得到第一次迭代最终值。第六步，将上一次迭代的最终结果作为初值，重复上面第三步到第五步，本实例共迭代三次。具体步骤流程图如图3所示。仿真结果如图4所示，图4(a)为理想图像，图4(b)为未经任何处理的原始图像，图4(c)为线性差值获得的结果，图4(d)、(e)、(f)为本发明方法一次迭代、二次迭代、三次迭代后的结果图，图(g)、(h)、(i)为patch方法一次迭代、二次迭代、三次迭代后的结果图。使用绝对误差(absoluteerror)对重建图像进行分析:其中n为像素点个数，ri是重建图像第i个点对应的灰度值，ti是理想图像第i个点对应的灰度值。结果如下：方法第101层的mae线性空间插值1.3368×10+6patch方法一次迭代3.0641×10+5patch方法二次迭代3.1477×10+5patch方法三次迭代2.6800×10+5本方法一次迭代2.9724×10+5本方法二次迭代2.9658×10+5本方法三次迭代2.6187×10+5从重建图像效果来看，使用空间插值方法获得的重建图像条形伪影严重，严重影响了图像诊断信息的提取，而patch方法和本方法则可以显著地改善图像质量。使用绝对误差对重建图像进行定量分析可以发现：和patch方法相比，本发明方法获得的图像质量更接近理想图像质量，每一次的迭代效果都优于patch方法，而且后面的迭代效果优于前一次的迭代效果。最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。当前第1页12