OSEM 重建算法及其改进算法的研究和比较

杨 娟，王明泉，石 浪，侯慧玲

（中北大学 仪器科学与动态测试教育部重点实验室，山西 太原030051）

0 引 言

最大似然期望值最大算法 （MLEM）是一种基于统计的迭代重建算法，该算法对噪声有较强的抑制作用，但因其计算量大，收敛速度慢，在实际当中的应用较少。Hudson等将有序子集 （ordered subsets of projection data）应用到MLEM 算法中，简称OSEM。OSEM 算法将投影数据按投影角度分解成有限个有序子集，每个子集应用MLEM 算法对图像更新一次，为一次子迭代，所有的子集全部使用完，为一次完整的迭代。相对于MLEM 算法，OSEM 算法在一次迭代过程中，重建图像更新了n次，因此加快了图像的收敛速度［1，2］。

OSEM 算法子集水平的不同会对重建图像的收敛速度以及重建图像的质量产生不同的影响。如何选取子集才能使重建结果最优，许多学者对OSEM 算法做出了不同的改进。Kadrmas提出了统计自适应EM 算法 （StatREM）［3］，该算法利用投影数据的统计信息，通过检验统计量，在每次迭代时自动调整子集数目，但其统计过程比较复杂；孔慧华等将MOS－BR算法应用到OSEM 算法中，提出了子集序列EM 算法 （SSEM）［4，5］，该算法提前将子集水平设为一个单调不增序列，每次迭代顺序使用序列中的子集数，但其子集序列的选择还要视重建结果予以调整；Vaissier提出了CROSEM （Count－Regulated OSEM）算法［6］，该算法只需设定子集水平为最大值，重建过程中，根据每个像素的活度大小，对其进行更新。该算法有效地改善了OSEM 算法重建SPECT 图像时低频元素丢失的问题［7］。

本文首先研究了OSEM 算法子集的不同选取对重建结果的不同影响，然后分析了SSEM 算法和CROSEM 算法的原理，最后分别使用OSEM 算法、SSEM 算法和CROSEM算法完成对Sheep＿Logan模型和固体火箭发动机模型的重建，对比重建图像的收敛速度、重建图像的质量，分析3种算法的实用性和有效性。

1 OSEM 算法

1.1 OSEM 算法简介

设定重建图像大小为n×n，360 个投影角度，每个角度下n条投影射线。将投影数据按投影角度划分成T 个有序子集。OSEM 算法的迭代公式如下

式中：x——待重建图像，pi——第i条射线的投影数据，wij——第i条射线穿过第j个像素的长度，Sl——第l个子集，l＝1，2，...，T。

1.2 子集水平对重建结果的影响

OSEM 算法重建图像过程中，子集水平较大时，图像先恢复高频元素信息，子集水平较小时，图像先恢复低频元素信息。并且，当投影数据不含噪声时，子集水平越大，重建图像的收敛速度越快；子集水平越小，图像的收敛速度越慢。但是一般CT 设备的投影数据都含有泊松噪声，此时，使用OSEM 算法进行重建，子集水平越大，图像收敛速度依旧比较快，但是随着迭代次数的增加，图像会发散，而子集水平越小，图像收敛速度越慢，并且图像的高频元素信息会丢失［8］。

2 SSEM 算法

为了改善OSEM 算法子集水平对重建图像的收敛速度以及重建图像的质量的影响，孔慧华基于MOS－BR 算法，对OSEM 算法做出改进，提出了子集序列EM 算法 （subset sequence EM），简称SSEM。

SSEM 算法对OSEM 算法的改进之处在于：SSEM 算法的子集水平不是一个固定值，而是一个单调不增序列。在迭代重建过程中，先使用大的子集水平，恢复图像的高频元素，加快重建图像的收敛速度，然后在下一次迭代时降低子集水平，恢复图像的低频元素，避免重建图像的发散。如此循环，到子集序列结束。

这样，在重建过程中，先恢复高频元素，再恢复低频元素，既加快了重建图像的收敛速度，又避免了重建图像的发散。

SSEM 算法子集序列的不同也会对重建图像的收敛速度以及重建图像的质量产生影响。一般序列中的第一个子集水平要不大于投影角度数，大小相同的子集水平的个数还要具体视重建结果予以调整［9］。

3 CROSEM 算法

OSEM 算法重建SPECT 图像时，子集水平选择越大，在重建过程中，低活度像素的灰度值会被更新为0，导致这些像素信息的丢失。而子集水平选择越小，图像的收敛速度越慢，随着迭代次数的增加，低活度像素的灰度值仍然会被更新为0，导致信息的丢失［10－12］。

为此，Vaissier提出了CROSEM （count－regulated OSEM）算法，该算法与OSEM 算法的不同之处在于：OSEM算法每次迭代更新图像，都是对图像的每一个像素值进行修正；而CROSEM 算法在迭代更新图像时，需要判断每一个像素的活度是否大于一个设定的阈值CTV （count threshold value），如果大于，这个像素值才会被更新，否则，不予更新。并且，CROSEM 算法的子集水平是其取值范围内的较大值，如此可以加快重建图像的收敛速度。

将CROSEM 算法应用到CT 图像重建，其重建步骤是：

（1）使用MLEM 算法，迭代一次重建图像。根据重建图像估计每个像素的活度，确定阈值CTV；

（2）k＝0，并给定正的初值图像x（0），指定迭代次数Niter，指定子集水平NS等于或小于投影角度数；

（3）进入子集Sl运算，计算下列参数

（4）根据第三步的结果，计算下列参数

重复步骤 （3）、步骤 （4），直到k＝Niter。

CROSEM 算法根据每个像素的活度来更新图像，解决了子集水平过高时，低活度像素丢失的问题，并且，其子集水平为取值范围内的较大值，加快了重建图像的收敛速度［9－11］。

4 实验结果及分析

为了对比3种算法的有效性和实用性，分别使用3种算法对Sheep＿Logan头部仿真模型和实际的固体火箭发动机模型进行重建，并采用归一化均方距离对重建结果进行评价，其表达式如式 （8）所示

式中：tu，v，ru，v——原始图像和重建图像中各像素的灰度值。珋t——原始图像灰度的平均值。d 值越小表示重建图像与原始图像的偏差越小，重建图像的质量越好。

重建Sheep＿Logan头部模型，重建图像大小为256×256，投影角度是将360°均分为256份，每个投影角度下有256条投影射线。OSEM 算法选定子集为64，迭代5 次；SSEM 算法选定的子集序列为256→128→64→32→32；CROSEM 算法的子集为最大值256，迭代5次。原始图像及其第128行的灰度曲线如图1所示，3种算法的重建图像及其第128行的灰度曲线分别如图2、图3、图4所示。

图1 原始图像及灰度曲线

3种算法重建图像的归一化均方距离d 见表1。

对实际采集的投影数据进行重建，实验参数为：射线源电压为290KV，电流为1.8mA。探测器大小为1920×1536，探元大小为0.127mm，射线源到旋转中心的距离为1060mm，旋转中心到探测器的距离为140mm，工件旋转一周，采样间隔为1 度。设定重建图像大小为256×256。OSEM 算法选定子集水平为36，迭代3次，其重建结果如图5所示。

表1 3种算法重建图像的归一化均方距离

图2 OSEM 重建图像及灰度曲线

图3 SSEM 重建图像及灰度曲线

图4 CROSEM 重建图像及灰度曲线

图5 OSEM 算法重建图像

SSEM 算法，选定子集序列为90→90→60→60→30，其重建结果如图6所示。

CROSEM 算法，设定子集水平为90，迭代5次，其重建结果如图7所示。

图6 SSEM 算法重建图像

图7 CROSEM 算法重建图像

对比3种算法重建的结果以及表1所示的归一化均方距离，可知：

（1）算法的收敛速度：3种算法相同的迭代次数，所得重建图像的归一化均方距离关系为：CROSEM＜SSEM＜OSEM，因此，CROSEM 的收敛速度优于SSEM，SSEM优于OSEM；

（2）重建图像的质量：CROSEM 优于SSEM，SSEM优于OSEM；

（3）子集选择的不定性：CROSEM 算法的子集水平为其取值范围内的较大值，而OSEM 算法与SSEM 算法的子集水平，都要根据重建结果来不断的调整，因此，CROSEM 算法更有实用性。

5 结束语

OSEM 算法子集的使用加快了MLEM 算法的收敛速度，但是OSEM 算法子集水平选择较高时，重建图像收敛速度快，但随着迭代次数增加，重建图像会发散；而子集水平较低时，重建图像的收敛速度又会变慢，而且，图像的高频信息会丢失。因此，子集水平的选取对重建图像的收敛速度以及重建图像的质量有很大的影响。SSEM 算法设定子集水平为一个单调不增序列，解决了OSEM 算法子集水平带来的问题。但是SSEM 算法子集序列的选取，还要视重建结果予以调整。而CROSEM 算法只需设定子集水平为一个固定值，在重建过程中，根据每个像素的活度对图像进行更新，有效地解决了低活度像素丢失的问题。

最后，对比3种算法对Sheep＿Logan头部仿真模型和实际固体火箭发动机模型的重建图像的收敛速度以及重建图像的质量，验证了CROSEM 算法的收敛速度更快，重建质量更好，有更好的实用性。

［1］SUN Shousi，QIU Jun，GUI Zhiguo，et al.The EM iterative imaging on the discrete reconstruction points ［J］.Journal of North University of China（Natural of Science Edition），2014，35 （2）：209－217 （in Chinese）.［孙守思，邱钧，桂志国，等.重建点模型的EM 迭代成像 ［J］.中北大学学报：自然科学版，2014，35 （2）：209－217.］

［2］QUE Jiemin，WANG Yanfang，SUN Cuili，et al.Comparison of four iterative algorithm based on incomplete projection reconstruction ［J］.CT Theory and Applications，2012，21 （2）：169－178 （in Chinese）.［阙介民，王燕芳，孙翠丽，等.基于不完备投影数据重建的四种迭代算法比较研究 ［J］.CT 理论与应用研究，2012，21 （2）：169－178.］

［3］KONG Huihua，PAN Jinxiao.Statistical self－adaptive ordered subsets algorithm for image reconstruction ［J］.Journal of North University of China（National of Science Edition），2010，31 （1）：79－83 （in Chinese）.［孔慧华，潘晋孝.图像重建的统计自适应子集算法 ［J］.中北大学学报：自然科学版，2010，31 （1）：79－83.］

［4］KONG Huihua.Research on iterative algorithms to speed image reconstruction ［D］.Taiyuan：North University of China Press，2006：32－47 （in Chinese）.［孔慧华.加速图像重建的迭代算法研究 ［D］.太原：中北大学出版社，2006：32－47.］

［5］KONG Huihua，PAN Jinxiao，WU Kun.Three dimensional OSEM image reconstruction algorithm based optimizing subset order ［J］.Journal of Test and Measurement Thchnology，2011，25 （2）：169－171 （in Chinese）. ［孔慧华，潘晋孝，吴琨.基于优化子集顺序的三维OSEM 图像重建算法 ［J］.测试技术学报，2011，25 （2）：169－171.］

［6］Vaissier PEB，Goorden MC，Taylor AB，et al.Count－regulated OSEM reconstruction ［C］／／IEEE Nuclear Science Symposium and Medical Imaging Conference Record，2012：3315－3320.

［7］ZHANG Xuezhu，QI Yujin.Research on fully 3Dimage reconstruction of pinhole SPECT imaging ［J］.Chinese Science Bulletin，2010，55 （18）：1846－1855 （in Chinese）.［张雪竹，漆玉金.针孔SPECT 成像完全三维图像重建的研究 ［J］.科学通报，2010，55 （18）：1846－1855.］

［8］JI Dongjiang.Research on iterative reconstruction algorithms of 3D cone beam CT ［D］.Chongqing：Chongqing University Press，2008：18－21（in Chinese）.［冀东江.三维锥束CT 迭代重建算法研究［D］.重庆：重庆大学出版社，2008：18－21.］

［9］KONG Huihua，PAN Jinxiao.An sequence subsets simultaneous algebraic reconstruction techniques［J］.CT Theory and Applications，2008，17 （2）：42－47 （in Chinese）. ［孔慧华，潘晋孝.序列子集联合代数重建技术 ［J］.CT 理论与应用研究，2008，17 （2）：42－47.］

［10］ZHANG Quan，FU Xuejing，LI Xiaohong，et al.High quality ordered subset expectation maximization reconstruction algorithm，based on multi－resolution for PET images ［J］.Journal of Computer Applications，2013，33 （3）：648－650（in Chinese）. ［张权，付学敬，李晓红，等.基于多分辨率的PET 图像优质有序子集最大期望重建算法 ［J］.计算机应用，2013，33 （3）：648－650.］

［11］Zhao Jingwu，Su Weining.An iterative image reconstruction algorithm for SPECT ［J］.Nuclear Science and Techniques，2014，25 （3）：1－5.

［12］LIU Lu.Parallel computing for multi－pinhole SPECT reconstruction algorithm ［D］.Beijing：Tsinghua University Press，2012 （in Chinese）.［刘路.多针孔SPECT 重 建算法并 行实现 ［D］.北京：清华大学出版社，2012.］