基于混合高斯/泊松最大似然函数的CBCT三维重建方法与流程

本发明属于图像重建
技术领域：
，主要涉及到混合高斯/泊松噪声下的cbct三维重建方法。
背景技术：
：计算机断层成像(computedtomography，ct)技术通过无损方式获取物体内部结构信息、材质、缺陷等信息，广泛用于医学辅助诊断、工业无损检测、安全检查等领域。随着技术的快速发展出现了平板探测器，锥束ct(conebeamcomputedtomographyct，cbct)系统通过它得到的投影数据是二维的，可以直接利用三维重建算法得到三维图像，是一种真正的三维重建。且cbct具有体积小，重量轻，移动灵活，扫描速度快，射线利用率高、重建分辨率高等显著优点，广泛用于介入手术治疗、无损检测等。cbct重建图像各向同性，提高了图像准确性和重建的实时性，能够更精确快速的分析病灶，提高患者的生活质量。cbct价格低廉，可紧凑安装在口腔外科诊所和私人诊所，更方便为患者提供服务。为了降低cbct成像中射线的辐射问题，低辐射剂量的cbct成像方法应运而生，但是这种情况下信号的信噪比将会随之降低。因此，有效地去除低剂量条件下的噪声问题，对于cbct重建方法而言是非常重要的。在低剂量条件下，系统中最常见的噪声包括高斯噪声和泊松噪声。高斯噪声主要来源于图像采集期间，往往由于不良照明或者高温引起的传感器噪声。泊松噪声是是由图像传感器中的光电转换过程引起的，在低剂量的cbct重建过程中，影响更为严重。1984年，feldkamp、davis和kress针对圆形扫描轨迹提出了最著名的锥束三维图像重建算法--fdk算法。它基于圆轨迹扫描，是一种易于实现、重建速度快的算法，并且它也是近似算法中最为经典的算法。至今fdk算法仍然是商业ct应用的主流，在ct重建算法中有着无可替代的地位。该算法是一种近似的三维图像重建算法，fdk算法在小锥角时具有良好的重建效果，但是当锥角增大时，与中心层距离越远重建效果越差。针对fdk算法的缺点，出现了多种fdk算法的衍生算法，改善了fdk算法重建图像质量，然而，这些重建算法不足之处在于低剂量条件下复杂噪声模型对重建图像质量的影响，且随着噪声增加重建图像质量退化严重。以fdk为主的解析算法相当于是对一个积分方程的求解，而后出现迭代算法相当于是对一个大型线性方程组的求解。由于线性方程组对成像几何结构和成像物理效应的刻画比积分方程更加的真实，还可以加入先验知识，约束条件等，所以迭代算法比解析算法对扫描模式的依赖性小，更适合于实际成像问题，且迭代算法得到的结果更加准确，对噪声的抑制效果更好。由于迭代算法的重建速度慢是其致命缺点，很多学者对迭代算法的加速进行了研究，通过设计提高计算机的并行处理能力，加速或者利用专门的图像处理器进行加速或是对算法进行改进加速。随着硬件技术的不断更新和各种专用服务器的开发使得计算机的运行速度越来越快，迭代算法的优势也将更好的发挥。迭代重建算法分为代数重建和统计重建，统计迭代重建算法(sir，statisticaliterativereconstruction)是通过对投影数据的统计特性进行分析建立目标函数，利用一定的估计方法重建出高质量的图像。统计迭代算法的研究，国内主要是针对放射型ct进行的。本发明将利用统计迭代重建算法进行锥形束ct的重建，提出一种基于混合高斯/泊松最大似然函数的cbct三维重建算法，从而重建出高质量的图像。技术实现要素：本发明提出一种基于混合高斯/泊松最大似然函数的cbct三维重建算法，属于统计迭代重建算法，是通过对投影数据的统计特性进行分析建立目标函数，利用一定的估计方法重建出高质量的图像。主要包含以下步骤：步骤1，设置系统参数和重建图像估计值，并加载真实投影数据；步骤2，对真实投影数据进行anscombe变换；步骤3，根据预设的迭代次数t，若当前迭代次数在预设的迭代次数内且未达到迭代终止条件e，则进入迭代循环内具体包括步骤4和步骤5，对重建结果进行迭代计算逼近更优值，否则退出；步骤4，对估计值的投影数据进行anscombe变换，并以真实投影数据的anscombe变换结果和估计值投影数据的anscombe变换结果为基础，以求解重建图像x的最大似然函数作为目标函数，求解最大似然函数更新图像的重建结果；步骤5，对重建图像进行3dtv正则化平滑去噪，优化重建结果，更新迭代步长进行新一轮的迭代计算。进一步的，步骤1中设置的系统参数包括，锥形束ct仪器的设置，射源到旋转中心和探测器的距离，旋转角及采样情况，探测器相关参数。进一步的，步骤1中，针对真实投影数据，通过调用cbct传统重建算法fdk获得重建图像作为重建图像估计值。进一步的，步骤4的具体实现方式如下，s41，设估计值为xk，对估计值进行投影得proj_xk，对估计值的投影进行anscombe变换得proj_traxk；s42，以重建图像x的最大似然函数作为生成重建图像的目标函数，由于第i次实测投影yi是服从混合泊松-高斯分布，因此其最大似然函数可以表示为：基于gat得到的概率密度函数可以近似为：这里x代表三维重建数据，aij代表第i个角度的投影矩阵的第j行，[p]j代表向量p中的第j个元素；在此基础上，负的log最大似然函数可以表示为其中，m代表投影个数，n2代表单个投影上元素的总个数，也代表第i个投影角度下投影矩阵a的行数，是指经过anscombe变换后第i个投影的第j个元素；对x求偏导数可得依据上式，利用真实投影数据和估计值投影数据及其ancombe变换结果求解最大似然函数，具体为：求解对其解的无穷大值和空值进行置零处理，并采用截断函数对梯度过大的值进行截断，截断后进行滤波反投影得最大似然函数的解，即图像的梯度信息gk；s43，判断gk是否满足迭代终止的条件，即梯度的模长是否小于迭代终止条件e，满足则终止迭代，否则利用估计值xk和反gk更新图像，得新估计值xkj1，xkj1＝xk-dk*gk(8)初始dk为迭代循环前赋予初值，初次循环dk为初值，非初次循环为上一次迭代后期更新的dk。进一步的，步骤5的具体实现方式如下，s51，将xkj1作为重建结果对其进行误差评估，判断误差是否随着迭代次数的增加而减小，否则通过松弛因子调整迭代步长dk＝s*dk，s为初始设置的一个常数值；重复步骤4中的s43，直至评估误差满足减小的条件；s52，如果评估误差满足减小的条件，则对重建图像进行3dtv正则化平滑去噪，正则项λtv(xkj1)其中k为图像的像素点索引，x,y,z分别代表三维图像三个方向，λ为正则化参数，表示图像xkj1的全变分，表示图像的梯度信息；作为全变分约束正则项，保证解具有一定的正则性，实现对图像噪声的平滑，去除噪声；s53，更新迭代步长，xk＝xkj1；gk_1＝gk，以进行新一轮的迭代，迭代步长dk通过barzilar-borwein方法确定，当迭代次数k>1时，gk_1为上一轮迭代的gk；其中，yk_1为本次迭代和上次迭代过程中梯度的差值，将其作为中间变量运用于迭代步长的更新中。本发明首先通过anscombe变换将服从高斯-泊松混合分布的实测投影数据变换为近似服从高斯分布，达到稳定噪声方差的目的，再通过梯度下降法求解混合高斯/泊松最大似然函数，结合三维全变分正则化方法(3d-totalvariation,3d-tv)进行迭代重建，从而达到改善重建图像质量的目的。附图说明图1为本发明实施例重建方法流程图；图2为fdk算法在不同大小高斯泊松混合噪声模型下重建结果；图3为投影数据bm3d去噪后fdk算法在不同大小高斯泊松混合噪声模型下重建结果；图4为sart算法在不同大小高斯泊松混合噪声模型下重建结果；图5为本发明方法在不同大小高斯泊松混合噪声模型下重建结果；图6为fdk算法、bm3d_fdk、sart、本发明方法在不同大小高斯泊松混合噪声模型下重建结果对比。图7为步骤3中迭代方法示意图。具体实施方式下面结合附图和实施例对本发明的技术方案作进一步说明。本发明基于混合高斯/泊松最大似然函数的cbct三维图像重建，利用anscombe变换将服从高斯-泊松混合分布的实测投影数据变换为近似服从高斯分布，达到稳定噪声方差的目的，通过梯度下降法求解混合高斯/泊松最大似然函数，结合三维全变分正则化方法(3d-totalvariation,3d-tv)进行迭代重建，本方法充分考虑了每次迭代估计值的不同，通过barzilar-borwein方法调整迭代步长，迭代多次来不断逼近最佳重建结果，改善重建图像质量。步骤1、设置系统参数，包括锥形束ct仪器的设置，射源到旋转中心和探测器的距离，旋转角及采样情况等，探测器相关参数。采集投影数据proj，初始化作为重建图像估计值的xk，anscombe变换等相关参数。步骤2、对投影数据proj进行anscombe变换。设置初始迭代步长，添加松弛因子来减小步长使得算法更加稳定。设置3dtv正则化率λ。设置初始误差参考值。设置迭代次数t和迭代终止条件e，在步骤4和步骤5的迭代次数范围内若达到了迭代终止的条件则终止迭代循环。anscombe变换方法如下，由于i代表了第i个投影数据，yi代表第i次实测投影，代表受到泊松噪声干扰的信号，服从泊松分布，ni代表系统的高斯噪声，ni～ν(0,δ2)。基于gat(anscombetransform)公式(1)，可以得到投影数据的anscombe变换的结果。步骤3，根据预设的迭代次数t，若当前迭代次数在预设的迭代次数内且未达到迭代终止条件e，则进入迭代循环内具体包括步骤4和步骤5，对重建结果进行迭代计算逼近更优值，否则退出；具体迭代方法如下，以求解重建图像x的最大似然函数作为目标函数，基于该目标函数形成外部迭代循环，每次迭代使用梯度下降法求解目标函数得图像梯度信息，计算重建图像，重建图像结合3dtv正则化项形成内部迭代循环，对重建的含噪图像平滑去噪，提高重建图像质量，如图7所示。步骤4、对估计值的投影数据进行anscombe变换，并以真实投影数据的anscombe变换结果和估计值xk投影数据的anscombe变换结果为基础，以求解重建图像x的最大似然函数作为目标函数，求解最大似然函数更新图像得重建结果。步骤4具体为：s41:对估计值进行投影得proj_xk，对估计值的投影进行anscombe变换得proj_traxk。s42:以重建图像x的最大似然函数作为生成重建图像的目标函数，由于第i次实测投影yi是服从混合泊松-高斯分布，因此其最大似然函数可以表示为：基于gat得到的概率密度函数可以近似为[1]：[1]marnissiy,zhengy,pesqueij.fastvariationalbayesiansignalrecoveryinthepresenceofpoisson-gaussiannoise[j].icassp,ieeeinternationalconferenceonacoustics,speechandsignalprocessing-proceedings,2016,2016-may:3964-3968.这里x代表三维重建数据，aij代表第i个角度的投影矩阵的第j行，[p]j代表向量p中的第j个元素。在此基础上，负的log最大似然函数可以表示为其中，m代表投影个数，n2代表单个投影上元素的总个数，也代表第i个投影角度下投影矩阵a的行数，是指经过anscombe变换后第i个投影的第j个元素。对x求偏导数可得依据上式，利用真实和估计值投影数据的及其ancombe变换结果求解最大似然函数，具体为：求解在实际的程序运算中即为对其解的无穷大值和空值进行置零处理，并采用截断函数对梯度过大的值进行截断，截断后进行滤波反投影得最大似然函数的解即图像的梯度信息gk。s43:判断gk是否满足迭代终止的条件，即梯度的模长是否小于迭代终止条件e。满足则终止迭代，程序运行结束，否则利用估计值xk和反gk更新图像，得新估计值xkj1：xkj1＝xk-dk*gk(8)初始dk为迭代循环前赋予初值，初次循环dk为初值，非初次循环为上一次迭代后期更新的dk。步骤5、对重建图像进行3dtv正则化平滑去噪，优化重建结果，更新迭代步长进行新一轮的迭代计算。步骤5具体为：s51:xkj1作为重建结果对其进行误差评估，判断误差是否随着迭代次数的增加而减小，否则通过松弛因子调整迭代步长dk＝s*dk，s为程序初始设置的一个常数值。重复步骤4的s3，直至评估误差满足减小的条件。s52:是则对重建图像进行3dtv正则化平滑去噪，正则项为λtv(xkj1)其中k为图像的像素点索引，x,y,z分别代表三维图像三个方向。λ为正则化参数，表示图像xkj1的全变分，表示图像的梯度信息。作为全变分约束正则项，保证解具有一定的正则性，实现对图像噪声的平滑，去除噪声。正则化参数λ可以用来调节正则项的程度，λ取值比较大时，正则项程度更高，平滑效果更加明显，λ取值比较小时，数据保真程度更高，边缘的细节区域更容易得到保持。当图像中噪声比较大时取较大的λ值，当图像中噪声比较小时取较小的λ值。s53:更新迭代步长，xk＝xkj1；gk_1＝gk，以进行新一轮的迭代。迭代步长dk可以通过barzilar-borwein方法确定，当迭代次数k>1时：gk_1为上一轮迭代的gk，其中，yk_1为本次迭代和上次迭代过程中梯度的差值，将其作为中间变量运用于迭代步长的更新中。本发明提供的方法能够用计算机软件技术实现流程。实施例以修定的shepp-logan体模模拟生成低剂量ct投影数据proj，采用基于圆轨迹的近似重建算法fdk获得重建图像reconimg为初始重建图像估计值，以下对本发明的流程进行一个具体的阐述，实施例具体实施方案为如下：通过修定的shepp-logan128×128×128体模模拟生成低剂量ct投影数据，设xk为重建图像估计值，其初值设置为fdk获得重建图像reconimg，dk为迭代步长，s为快速更新代步长的松弛因子，t为迭代次数，e为迭代终止条件。sigma为所添加高斯噪声大小的参数，alpha为所添加泊松噪声的大小，λ为正则化率。步骤1、加载shepp-logan模型，生成211张128×128cbct投影数据proj，调用cbct传统重建算法fdk获得重建图像为reconimg作为初始重建图像估计值xk。步骤2、设置对投影数据添加噪声的参数，分别设置sigma为2,6,10,14,18的高斯噪声，lambda为1e5的高斯泊松噪声混合模型。设置初始迭代步长dk为1，松弛因子s＝0.5。设置3dtv正则化率λ为0.01。设置初始误差参考值qm11为0.2。设置迭代次数t＝10和迭代终止条件e＝0.2。根据添加噪声情况设置进行anscombe变换的相关参数。对投影数据proj进行anscombe变换变换得proj_tra。步骤3、根据预设的迭代次数t，若当前迭代次数在预设的迭代次数内则进入迭代循环内，否则退出程序。即迭代次数小于等于10时进入迭代循环内，对重建结果进行不断迭代计算逼近更优值。步骤4s41:进入迭代循环，对重建图像估计值xk投影得proj_xk进行anscombe变换得proj_traxk。s42:利用proj_tra、proj_traxk，及投影数据proj求取目标函数偏导数的因子fx1,即：对偏导数的因子的无穷大值和空值进行置零处理，并采用截断函数进行截断，截断值此处设置为0.05，对截断后的偏导因子进行滤波反投影得目标函数的解gk。s43:gk若满足迭代终止的条件则终止迭代，程序运行结束，否则更新重建图像xkj1＝xk-dk*gk。步骤5s51:对重建图像xkj1进行误差评估判断，判断误差是否随着迭代次数的增加而减小，否则通过松弛因子调整迭代步长dk＝s*dk，重复步骤xkj1＝xk-dk*gk，直至评估误差满足减小的条件。s52:对重建图像xkj1进行3dtv正则化平滑去噪。s53:更新迭代步长dk，xk＝xkj1；gk_1＝gk，以进行新一轮的迭代。实施例过程涉及数据如下表：1:在泊松噪声分布大小为lambda＝1e5,高斯噪声大小分别为sigma＝2，sigma＝6，sigma＝10，sigma＝14，sigma＝18的不同混合噪声模型下使用fdk算法进行cbct三维重建。结果数据如表1，实验结果如图2.2:在泊松噪声分布大小为lambda＝1e5,高斯噪声大小分别为sigma＝2，sigma＝6，sigma＝10，sigma＝14，sigma＝18的不同混合噪声模型下使用bm3d对投影数据进行预处理然后通过fdk算法进行cbct三维重建。结果数据如表2，实验结果如图3.3:在泊松噪声分布大小为lambda＝1e5,高斯噪声大小分别为sigma＝2，sigma＝6，sigma＝10，sigma＝14，sigma＝18的不同混合噪声模型下使用sart算法进行cbct三维重建。结果数据如表3，实验结果如图4.4:在泊松噪声分布大小为lambda＝1e5,高斯噪声大小分别为sigma＝2，sigma＝6，sigma＝10，sigma＝14，sigma＝18的不同混合噪声模型下使用本文算法:基于混合高斯/泊松最大似然函数的cbct三维重建算法(ml_tv)进行cbct三维重建。结果数据如表4，实验结果如图5.fdk、投影数据bm3d(三维块匹配)去噪fdk，sart(同时迭代重建迭代算法)和本文算法(ml_tv)实验结果比较：表1利用fdk算法进行cbct三维重建实验结果sigma-lambda2-1e56-1e510-1e514-1e518-1e5rmse0.0477740.0525040.0568280.0608430.064607表2利用bm3d_fdk算法进行cbct三维重建实验结果sigma-lambda2-1e56-1e510-1e514-1e518-1e5rmse0.0477610.0486540.0495150.0502110.050598表3利用sart算法进行cbct三维重建实验结果sigma-lambda2-1e56-1e510-1e514-1e518-1e5rmse0.0400930.0419420.0435830.045330.046899表4利用本发明方法(ml_tv)进行cbct三维重建实验结果sigma-lambda2-1e56-1e510-1e514-1e518-1e5rmse0.0242420.0269430.0305910.0356020.035886lambda0.010.050.050.050.05sigma4040404245实施例过程对比了在不同大小高斯泊松混合噪声模型下不同重建方法的重建结果，涉及重建结果图如下：fdk算法在不同大小高斯泊松混合噪声模型下重建结果如图2；投影数据bm3d去噪后fdk算法在不同大小高斯泊松混合噪声模型下重建结果如图3；sart算法在不同大小高斯泊松混合噪声模型下重建结果如图4；本发明方法在不同大小高斯泊松混合噪声模型下重建结果如图5；fdk算法,bm3d_fdk、sart、本发明方法在不同大小高斯泊松混合噪声模型下重建结果对比如图6；cbct三维重建方法目前分为解析法和迭代重建方法，fdk算法为最经典的解析重建方法，而bm3d为目前去除高斯噪声效果非常好的去噪方法，用bm3d对投影数据进行预处理再通过fdk重建，由表1、表2可以看出运用bm3d去噪方法预处理后，fdk重建图像质量有所提升，由图2图3可以看出运用bm3d后重建结果优化的效果。由表3可以看出，迭代算法sart得到的结果更加准确，对噪声的抑制效果更好，重建图像质量更好些，即使与使用很好的去噪方法bm3d的fdk重建结果相比也要好一些，由表4可以看出本文的方法重建结果比sart还要好很多。通过图4图5可以明显看出两种迭代重建算法重建结果的差异。通过gpu的加速克服迭代重建速度慢的问题后，迭代重建算法比解析算法更适合于实际成像问题，且迭代算法得到的结果更加准确，对噪声的抑制效果更好。本文提出的方法属于统计迭代重建方法，使得重建图像的性能有了很大的提升，重建图像质量更佳，若运用于实际更能方便医护人员区分病灶所在。本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属
技术领域：
的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。当前第1页12