扇形模式下SPECT图像重建算法之比较研究

【摘要】  　　目的： 比较研究扇形几何模式下单光子发射断层成像（SPECT）中三种典型重建算法的衰减补偿性能. 方法： 描述并分析扇行投影方式下FBP, OS?EM和Novikov逆变换三种算法的重建公式，对Shepp?Logan模型进行重建，并对重建时间及图像质量进行比较. 结果： 基于Novikov逆变换的定量解析重建算法得到的图像质量与OSEM迭代算法近似，而重建时间大大缩短. 结论： 定量解析重建算法可快速有效补偿非均匀衰减因素影响，具有广泛应用前景.

【关键词】  体层摄影术 发射型机 单光子 有序子集最大期望 滤波反投影 定量重建

     0引言

    单光子发射断层成像（single photon emission computer tomography, SPECT）技术广泛应用在核医学临床诊断中. 由于人体组织对发射的光子具有吸收衰减作用，如果在重建过程中不考虑该因素，将导致出现假阳性结果. 以往对非均匀衰减的补偿主要是通过迭代算法来实现［1-3］. 近年来Novikov［4］给出了平行投影下、非均匀衰减Radon逆变换的求解公式，才使得具有任意真实衰减分布的SPECT解析重建算法成为可能. Kunyansky［5］利用Novikov的逆Radon变换求解公式，提出一种可校正非均匀衰减的SPECT解析重建算法，对类似人的胸腔这样复杂的非均匀衰减分布，也能进行准确的补偿.  本研究在将Novikov逆变换公式扩展至扇行投影方式的基础上［6］，对三种典型重建算法，即经典滤波反投影法（FBP），迭代算法的代表有序子集最大期望值法（OS?EM），及基于Novikov逆变换的解析重建算法的衰减补偿性能、重建图像质量及重建时间进行比较评价，为今后该领域的研究工作提供参照.

    1对象和方法

    1.1对象对扇行投影方式下FBP，OS?EM以及Novikov求逆变换公式三种重建算法进行描述与推导，并对重建结果进行定量分析.

    1.2方法

    1.2.1FBP滤波反投影重建算法的基本思想是：在某一投影角度下取得投影函数（一维函数），对此一维投影函数做滤波处理，得到一个经过修正的投影函数. 然后再将此修正后的投影函数做反投影，得到重建后的图像［7］：

    f(x,y)=∫π〖〗0dθ ∫+∞〖〗-∞gθ(R)δ(xcosθ+ysinθ-R)dR(1)

    其中，gθ(R)=∫+∞〖〗-∞Fθ(ρ)|ρ|e2πjρRdρ，Fθ(ρ)为在θ角度下投影函数（极坐标）的一维傅立叶变换. δ(xcosθ+ysinθ-R)代表直线xcosθ+ysinθ=R

　　滤波反投影在实现图像重建时，只需做一维傅立叶变换，且可并行进行，图像重建速度快，因而在临床中得到广泛应用. 但是，由于该算法不能补偿衰减等因素对图像的影响，重建结果只能提供定性分析，图像中存在伪轮廓现象. 扇形投影方式下的滤波反投影研究可参阅［7］.

    1.2.2有序子集最大期望值（OS?EM）算法作为迭代算法的代表，OS?EM 以最大似然?最大期望值方法(ML?EM)为基础. 由于在ML?EM算法中，每一次对所有投影数据计算的结果只能更新重建图像一次，而在OS?EM算法中，投影数据被划分为G个有序子集Sg∶g=1,2,Λ,G，对每个子集的计算结果都将重新更新一次图像. 这样，对所有的子集都计算一遍后，就相当于对初始图像更新了G次，从而大大提高收敛的速度.

    对有序子集g=1,2,Λ,G

    投影：p(I, g)lmn=∑〖〗i?j?k?f(I, g)i?j?k?hi?j?k?, lmn (l, m, n)∈Sg（2）

    反投影：

    f(I, g+1)ijk=f(I, g)ijk〖〗∑〖〗lmn∈SGhijk,lmn∑〖〗lmn∈SGhijk,lmnplmn〖〗p(I, g)lmn（3）

    其中，hijk,lmn是点(i,,j,k)在探测头(l,m,n)上的投影. I为迭代次数. 每次迭代完成后，f(I, G+1)ijk作为新的f(I+1,1)ijk用于下次迭代计算.

   在利用OS?EM算法重建SPECT图像的过程中，子集的选取极为关键. 子集数量过少，将影响收敛速度；过多，则可能导致不收敛或收敛到局部收敛点. 在实际操作中，子集数目通常取8的倍数. 迭代算法能够处理复杂的真实成像模型，重图像质量好，但由于计算量较大，且存在正则及收敛问题，目前还未在临床广泛应用.

    1.2.3Novikov逆变换公式平行投影方式下Novikov的逆变换公式可参阅文献［4］. 在扇形投影方式下，投影线、探测器及重建图像间有对应几何关系（图1）.  其中O，S分别为坐标原点与扇形束焦点，P为重建图像中任意点. σ为OS与投影线夹角，σ?为OS与PS夹角. β为OS与y轴夹角，D为OS间距离.

    基于推导卡迪尔坐标与极坐标间的偏微分方程，我们将Novikov逆变换公式进一步推广到扇形模式下［6］：

    f(r,φ)=1〖〗4πRe∫2π0Wβ(r,φ,σ)〖〗πKdβ∫π/2-π/2D2cos2σ〖〗sin(σ?-σ)?g(σ,β)〖〗?σdσ+1〖〗4πRe∫2π0W1β(r,φ,σ)〖〗πKdβ∫π/2-π/2g(σ,β)Dcosσ〖〗sin(σ?-σ)dσ (4)

    其中，K=r2+D2+2rDsin(β-φ)(5)

    σ?=arcsinrcos(β-φ)〖〗K(6)

    g(σ,β)=e1〖〗2［(I+iΓ)k］(σ,β) p(Dsinσ,σ+β)(7)

    Wβ(r,φ,σ?)=eaθ(s,t)-h(s,θ)|θ=σ?+β,s=rcos(θ-φ),t=rsin(θ-φ)(8)

    W1β(r,φ,σ?)=?eaθ (s,t)-h(s,θ)〖〗?s|θ=σ?+β,s=rcos(θ-φ),t=rsin(θ-φ)(9)

    Novikov逆变换公式解决了一直以来通过解析算法无法解决的、非均匀衰减情况下逆Radon变换的精确求解问题. 基于该公式，近年来提出了大量SPECT解析重建算法.

    2结果

    比较上述三种重建算法的性能，采用Shepp?Logan模型进行数字仿真实验，得出所用模型及其衰减分布（图2）.

    图2Shepp?Logan模型(A)及其衰减系数分布(B)

    衰减图中包含了三种不同的灰度值，0.25，0.75和1，分别代表肺，软组织和骨的衰减系数. D=2，σ=60°. Shepp?Logan模型和衰减分布图完全包含在扫描区域内.

    FBP, OS?EM及Novikov逆公式对上述模型的扇形投影数据进行重建，得到图像（图3）. 衰减因子在FBP重建结果中造成了明显的伪轮廓（图3A），而在其余两种算法中，则得到了较好的补偿（图3B，C）. A: FBP; B:OS?EM; C: Novikow.

    图3投影数据重建

    为定量评价不同算法的重建性能，实验中采用感兴趣区域（region of interest, ROI）的偏差(bias)?方差(variance)特性对算法进行评估. ROI的偏差与方差分别定义为：

    Bias(%)=100(Bz(R)/Z(R))(10)

    Variance(%)=100(Vz(R)/Z2(R))(11）

    其中，Bz(R)=1〖〗N∑N〖〗k=1(fk(x,y)-f(x,y))(12)

    Vz(R)=1〖〗N-1∑N〖〗k=1(fk(x,y)-f(x,y))2(13)

    Z(R)=1〖〗NR∑〖〗x,y∈Rf(x,y)(14)

　　R为感兴趣区域，NR为感兴趣区域的面积（像素数）. N为加噪声后的实现次数. fk(x,y)为第次重建结果. fk(x,y)为N次实现的均值.

　　仿真实验中，选取了具有不同灰阶的三个ROI区域（图2A,面积均为12个像素）. 对模型理想投影数据进行了200次Poisson加噪实现，并用三种算法分别重建，得到的ROI偏差?方差结果（表1）. 表1扇形投影方式下不同ROI的偏差?方差相同条件下3种算法所需的重建时间，FBP所需时间最短，Novikov次之，OS?EM算法所需时间最长（表2）. 表2扇行投影模式下三种算法重建图像所用时间

　　3讨论

　　FBP, OS?EM及基于Novikov逆变换的重建算法代表了目前SPECT重建领域中的方向. 其中，FBP算法因运算速度快，得到了最广泛的应用. 但由于该算法不能对各种退化因素特别是衰减因子进行补偿，因此一直以来，对SPECT的定量分析还是以迭代算法为主. 以OS?EM为代表的迭代算法能够对重建过程中的各种影响因素进行建模及补偿，从而最大限度的逼近统计意义下的最优解. 但其迭代特性决定了运算量较大，重建图像需时较长，限制了在临床的广泛应用. 近年来很多研究提出了多个迭代加速算法，但考虑到图像分辨率的提高将导致迭代算法的量呈几何指数增长，因此对于高分辨率及动态图像，迭代算法仍很难满足实时性要求. Novikov逆变换公式的给出，使得理论上能够在逆Radon变换中对衰减因子进行精确的补偿，是今后SPECT重建技术的重要发展方向. 如果能进一步与噪声、探测器响应、散射等退化因素的解析补偿算法相结合，这种“FBP类型”的定量解析重建将非常可能成为临床新的重建标准. 由于本研究并未对Novikov的算法进行性能与速度上的优化. 因此，所得到的数据并不能完全体现该算法的优越性. 另外，由于受到具体实现因素（如插值计算、Hilbert变换方法实现等）的影响，重建结果中部分区域存在偏离真实值较大的现象，这些都有待进一步研究及完善.

【】
  ［1］ Hudson HM, Larkin SR. Accelerated image reconstruction using ordered subsets of projection data［J］.IEEE Trans Med Imag， 1994，13:601-609.

［2］ Byrne L. Block?Iterative Methods for Image Reconstruction from Projections ［J］. IEEE Trans Image Proc， 1994，5:792-794.

［3］ Browne J， De Pierro AR. A row?action alternative to the EM algorithm for maximization likelihood in emission tomography［J］. IEEE Trans Med Imag，1996， 15（5）: 687-699.

［4］ Novikov R. An inversion formula for the attenuated X?ray transformation［J］. Ark Math, 2002,40:145-167.

［5］ Kunyansky L. A new SPECT reconstruction algorithm based on the Novikov’s explicit inversion formula［J］. Inverse Problems, 2001, 17: 293-306.

［6］ You J. FBP algorithms for attenuated fan?beam projections［J］. Inverse Problems, 2005, 21: 1179-1192.

［7］ Horn BKP. Fan?beam reconstruction methods［J］. Proc IEEE, 1979, 67:1616-1623.