1. 清华大学 工程物理系, 北京100084;
2. 北京永新医疗设备有限公司, 北京 102206
收稿日期：2022-06-10
基金项目：核能开发科研项目(〔2019〕1342号); 国家重大科研仪器研制项目(81727807); 国家重点研发计划项目(2016YFC0105405, 2019YFF0302503); 核退役科研项目(〔2018〕1521号)
作者简介：程李(1992—), 男, 博士
通讯作者：刘亚强, 研究员, E-mail：liuyaqiang@tsinghua.edu.cn
摘要：碘-131单光子发射计算机断层扫描(single photon emission computed tomography, SPECT)平片成像已广泛应用于甲状腺癌的临床诊断和治疗评估中。由于碘-131发射的Gamma光子能量较高, 穿透SPECT准直器隔壁的概率较大, 易导致最终图像上出现放射状伪影。为解决该问题, 该文提出了一种点扩展函数解卷积结合先验正则化的方法, 并基于Monte Carlo模拟和临床数据对该方法进行了验证。模拟数据表明该算法具有很好的单调收敛性。当采用127×127的点扩展函数(像素大小2.2 mm×2.2 mm), 设置正则化参数为0.01时, 模拟数据、临床点源数据和人体数据的放射状伪影区域计数值可分别降低至原始值的4%、35%和28%, 同时还可以维持背景区域较低的噪声水平。研究结果表明：该文算法可以有效解决碘-131 SPECT平片成像中的放射状伪影问题, 从而提高临床诊断精度。
关键词：碘-131单光子发射计算机断层扫描放射状伪影解卷积先验正则化
Artifacts correction algorithm for iodine-131 SPECT planar imaging
CHENG Li1, LIU Fan2, GAO Lilei2, LIU Hui1, LIU Yaqiang1
1. Department of Engineering Physics, Tsinghua University, Beijing 100084, China;
2. Beijing Novel Medical Equipment Ltd., Beijing 102206, China
Abstract: Objective Iodine-131 SPECT (single photon emission computed tomography) planar imaging has been widely used in the clinical diagnosis and treatment evaluation of thyroid cancer. Because of the high-energy emissions of iodine-131, the photons have a high probability of penetrating the collimator septa of SPECT, causing "spoke" artifacts in the final result. The "spoke" artifacts make it difficult to distinguish accurate concentration regions of iodine-131, and they may o b scure lower uptake regions nearby, such as metastatic spread to lymph nodes. In this paper, a deconvolution method based on the point spread function was com b ined with a priori regularization to suppress the spoke artifacts and to improve the diagnostic accuracy in clinical studies. Methods This study is based on the NET632 SPECT system with the corresponding high-energy general purpose collimator. The collected data follow the Poisson distribution, and they consist of two parts, a forward projection of the true activity distribution and the scatter data. The forward projection progress can be well modeled using a shift-invariant PSF (point spread function). An objective function is built based on the aforementioned approximation, and a priori function is introduced to regularize the reconstruction. A monotonic and convergent algorithm is derived to iteratively solve the o b jective function. In contrast, the conventional deconvolution method regularizes the solution using a total variation term, and the ob jective function is optimized based on the "one-step-late" algorithm; thus, nonnegativity and convergence can not be guaranteed. The triple-energy window method is employed to estimate scattering data, and PSFs of different sizes are generated based on Monte Carlo simulations. Simulated NEMA torso phantom data are reconstructed with different parameters to validate the monotonicity and convergence of the proposed method. Moreover, the dataset is also used to evaluate the effects of PSF size and regularization strength on reconstruction images. Normalized spoke counts and background noise are calculated for quantitative comparison. Simulation data are also used to compare the reconstruction performance of the proposed method and the conventional deconvolution method. With the optimized parameters determined by simulation data, the proposed method is further validated by clinical point data and volunteer data. Results With different reconstruction parameters, the objective function value increased monotonically, and the image differences between two adjacent iterations rapidly reduced to a value close to zero. The simulation study also demonstrated that a 127×127 PSF size could provide performance similar to a 255×255 PSF size, which was significantly better than 63×63 and 31×31 PSF sizes. A study on different regularization strengths suggested an optimized regularization parameter, 0.01. When the PSF size and regularization parameter were set as 127×127 and 0.01, the mean spoke counts could be reduced to 4% of the original value with a low background noise level. A comparison study based on simulation data showed the superiority of the proposed method over the conventional deconvolution method. The clinical point data and volunteer data also validated the performance of the proposed method, and the mean spoke counts could be reduced to 35% and 28% of the original values, respectively. Conclusions The proposed method suppresses the spoke artifacts in iodine-131 imaging using a PSF that models the physical response, and it also introduces a priori regularization to suppress noise amplified by deconvolution. The derived algorithm can guarantee the monotonicity and convergence of the iterative reconstruction. Studies on simulation data and clinical data have demonstrated that the proposed method can achieve the desired performance and is expected to improve the diagnostic accuracy in clinical studies.
Key words: iodine-131single photon emission computed tomography"spoke" artifactsdeconvolutiona priori regularization
单光子发射计算机断层扫描(single photon emission computed tomography，SPECT)[1-2]是一种重要的分子影像技术，可以通过追踪放射性药物在人体内的分布反映生理代谢情况。目前临床主要使用双探头平行孔SPECT系统。准直器作为SPECT系统的重要部件，可以阻止非特定方向的光子入射到探测器表面，从而对人体内放射性药物进行定位。
对于甲状腺癌的诊断和治疗评估，碘-131 SPECT平片成像是一种常用且有效的技术手段[3-4]。碘-131衰变纲图如图 1所示，临床SPECT显像主要利用364.5 keV光子。对于采用铅合金和六边形设计的平行孔准直器(见图 2a)，较高能量的Gamma光子容易穿透准直器隔壁到达探测器，从而形成放射状伪影[5](见图 2b)。该伪影一方面会导致难以区分真实的碘-131富集区域，另一方面可能会掩盖周围较低摄取区域，例如淋巴结转移。因此，消除放射状伪影对于辨别真实的碘-131分布以及开展进一步的临床诊断、治疗评估都具有十分重要的意义。
 |
| 图 1 碘-131衰变纲图 |
| 图选项 |
 |
| 图 2 平行孔准直器对应的放射状伪影示意图 |
| 图选项 |
针对放射状伪影问题，一种常用的解决方法是采用点扩展函数(point spread function，PSF)对准直器穿透效应进行建模，并结合散射校正和先验函数实现解卷积算法，从而消除放射状伪影[5-6]。先验函数的选择和相应迭代算法的设计都会对最终效果产生影响。传统方法采用全变分先验改善重建结果，但是在迭代算法方面存在一定问题，迭代算法的收敛性无法保证。本文提出了一种新的解卷积算法用于解决碘-131成像中的放射状伪影问题，选取适当的先验函数并推导了单调收敛的迭代算法；对提出的算法，考察了不同大小PSF和不同正则化参数对结果的影响，并基于AMIDE软件[7]对感兴趣区(region of interest, ROI)进行定量分析；对于优化后的参数，采用临床数据对算法效果进行了进一步评估。
1 SPECT系统简介本研究主要针对北京永新医疗设备有限公司的NET632 SPECT系统[8]进行探索。该系统的探测器采用碘化钠(NaI)晶体，固有空间分辨率和能量分辨率分别为3.4 mm和9.8%@140 keV。该系统配备的用于碘-131成像的高能通用准直器参数如表 1所示。
表 1 高能通用准直器参数
| 参数名称 | 参数描述 |
| 芯体孔形状 | 六边形 |
| 芯体孔内径/mm | 3.40 |
| 芯体孔长度/mm | 48.00 |
| 芯体隔壁厚/mm | 1.65 |
| 芯体材料 | 铅合金 |
表选项
2 伪影消除算法考虑碘-131成像时的准直器穿透效应和噪声特性，平片图像y和真实图像x之间的关系可表示为
| $\boldsymbol{y} \sim \operatorname{Poisson}(\boldsymbol{x} \otimes \boldsymbol{h}+\boldsymbol{s}).$ | (1) |
目标函数L(x)可表示为
| $L(\boldsymbol{x})=\sum\limits_i y_i \ln \left[(\boldsymbol{x} \otimes \boldsymbol{h})_i+s_i\right]-\left[(\boldsymbol{x} \otimes \boldsymbol{h})_i+s_i\right] .$ | (2) |
此外，为了压低噪声并保持图像边界清晰，改善重建图像质量，引入先验函数
| $U(\boldsymbol{x})=\frac{1}{4} \sum\limits_j \sum\limits_{k \in N_j} \omega_{j k} \sqrt{\left(x_j-x_k\right)^2+\delta^2} .$ | (3) |
综上所述，最终基于目标函数(4)对原始图像进行重建，具体可表示为
| $\psi(\boldsymbol{x})=L(\boldsymbol{x})-\beta U(\boldsymbol{x}) .$ | (4) |
2.1 迭代算法对于目标函数(4)采用迭代算法求极大值。为保证迭代过程单调收敛，采用优化转移的方法[9-10]，寻找目标函数的替代函数并推导迭代公式。L(x)的替代函数可表示为：
| $Q_L\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=\sum\limits_j p_j \left( {\hat x_{j, {\rm{EM}}}^{n + 1}{\rm{ln}}{x_j} - {x_j}} \right), $ | (5) |
| $\hat{x}_{j, \mathrm{EM}}^{n+1}=\frac{x_j^n}{p_j}\left[\boldsymbol{h}^{\mathrm{T}} \otimes \bf{b} \bf{p}\right]_j, $ | (6) |
| $\bf{b p}_i=\frac{y_i}{[\boldsymbol{h} \otimes \boldsymbol{x}+\boldsymbol{s}]_i}, $ | (7) |
| $\boldsymbol{h}^{\mathrm{T}}(a, b)=\boldsymbol{h}(N-1-a, N-1-b), $ | (8) |
| $\boldsymbol{p}=\boldsymbol{h}^{\mathrm{T}} \otimes \boldsymbol{I} .$ | (9) |
U(x)的替代函数可表示为：
| $Q_U\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=\frac{1}{2} \sum\limits_j w_j^n\left(x_j-\hat{x}_{j, \text { reg }}^{n+1}\right)^2, $ | (10) |
| $w_j^n=\sum\limits_{k \in N_j} w_{j, k}^n, $ | (11) |
| $\hat{x}_{j, \text { reg }}^{n+1}=\frac{1}{2 w_j^n} \sum\limits_{k \in N_j} w_{j, k}^n\left(x_j^n+x_k^n\right), $ | (12) |
| $w_{j, k}^n=\frac{\omega_{j, k}}{\sqrt{\left(x_j^n-x_k^n\right)^2+\delta^2}} .$ | (13) |
综上可得目标函数ψ(x)的替代函数，
| $Q_\psi\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=Q_L\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)-\beta Q_U\left(\boldsymbol{x} ; \boldsymbol{x}^n\right) .$ | (14) |
| $Q_\psi\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)-Q_\psi\left(\boldsymbol{x}^n ; \boldsymbol{x}^n\right) \leqslant \psi(\boldsymbol{x})-\psi\left(\boldsymbol{x}^n\right), $ | (15) |
| $\nabla Q_\psi\left(\boldsymbol{x}^n ; \boldsymbol{x}^n\right)=\nabla \psi\left(\boldsymbol{x}^n\right).$ | (16) |
| $x_j^{n+1}=\frac{\sqrt{\left(1-\beta_j^n \hat{x}_{j, \mathrm{reg}}^{n+1}\right)^2+4 \beta_j^n \hat{x}_{j, \mathrm{EM}}^{n+1}}-\left(1-\beta_j^n \hat{x}_{j, \mathrm{reg}}^{n+1}\right)}{2 \beta_j^n}, $ | (17) |
| $\beta_j^n=\frac{\beta w_j^n}{p_j} .$ | (18) |
2.2 相关工作对比一种传统的PSF解卷积算法是基于采集数据的Poisson分布模型并结合全变分先验构建目标函数，具体可表示为：
| $J_{\mathrm{reg}}(\boldsymbol{x})=\int|\nabla \boldsymbol{x}(\boldsymbol{r})| \mathrm{d} \boldsymbol{r}, $ | (19) |
| $ \frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}^n=-\operatorname{div}\left(\frac{\nabla \boldsymbol{x}^n(\boldsymbol{r})}{\left|\nabla \boldsymbol{x}^n(\boldsymbol{r})\right|}\right). $ | (20) |
采用一步延迟法(one-step-late, OSL)解决目标函数难以优化的问题[6]。在具体考虑散射数据的影响后[5]，可将最终迭代公式表示为：
| $x_j^{n+1}=\frac{x_j^n}{1+\beta\left[\frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}\right]_j^n}\left[\boldsymbol{h}^{\mathrm{T}} \otimes \bf{b}^{\prime}\right]_j, $ | (21) |
| ${\bf{b}}_i^{\prime}=\frac{\left[y_i-s_i\right]_{+}}{[\boldsymbol{h} \otimes \boldsymbol{x}]_i} .$ | (22) |
对比式(6)和(21)可以发现，该算法与传统EM算法相比，主要有3种不同做法：1) 将不同像素位置处的灵敏度因子pj全部视为1；2) 引入全变分先验的导数对pj进行修正，以体现先验函数对重建的影响；3) 基于扣除散射数据后的投影进行重建。
由于实际考察图像的尺寸有限，导致不同像素位置的灵敏度因子pj存在一定差异。第1种做法会导致图像边缘处出现一定差异；第2种做法导致选择参数β时要十分谨慎，并且无法保证迭代过程的非负性和单调收敛性。直接扣除散射数据的投影，理论上不再服从Poisson分布；而且由于噪声影响，直接扣除散射数据的做法可能导致某些像素位置出现负值，式(22)中的[*]+操作会强制将这些像素值置为0，这些都会对最终结果造成一定影响。
2.3 散射校正对于散射数据可基于多能窗法[11-12]进行估计，
| $\boldsymbol{s}=\left(\frac{C_{\text {lower }}}{W_{\text {lower }}}+\frac{C_{\text {upper }}}{W_{\text {upper }}}\right) \times \frac{W_{\text {main }}}{2}.$ | (23) |
对于本研究的SPECT系统，主能窗设置为328~400 keV, 相邻的低能窗和高能窗分别设置为304~ 328 keV和400~424 keV。
2.4 PSF获取获取准确的PSF是伪影消除算法的关键。本文采用Monte Carlo模拟方法，基于GATE v8.0平台[13]对SPECT探测器和高能通用准直器进行模拟，并采集了碘-131点源的投影数据。考虑到SPECT平片成像时甲状腺到准直器前表面的典型距离，模拟时点源与准直器前表面的典型距离设置为10 cm。
对于模拟采集的数据，提取主能窗中未经散射的数据用于生成PSF。为了降低PSF图像的噪声，采用未经散射的数据约976×104计数，获得的碘-131点源PSF图像如图 3所示，PSF图像大小为255×255，像素为2.2 mm×2.2 mm。
 |
| 图 3 碘-131点源PSF图 |
| 图选项 |
此外，应用时还分别截取了大小为127×127、63×63、31×31的PSF，用于评估不同大小的PSF对于伪影消除效果的影响。这3种大小的PSF用到的事件数(指用于生成PSF的数据)分别占总计数的96.9%、84.0%和69.8%，占比越小表明舍弃的有效信息越多。
3 模拟数据评估3.1 Monte Carlo模拟为了评估算法性能并优化参数，本文基于GATE v8.0平台模拟了美国国家电气制造商协会(National Electrical Manufacturers Association，NEMA)躯干模型，如图 4所示, 在直径分别为13、17和22 mm的3个球体内充满碘-131放射源。图 5展示了模拟采集的数据能谱和各能窗投影图。主能窗、高能窗和低能窗的光子计数分别为22.8×104、2.1×104和4.3×104。
 |
| 图 4 NEMA躯干模型横断面示意图(单位：mm) |
| 图选项 |
 |
| 图 5 Monte Carlo模拟数据 |
| 图选项 |
3.2 收敛性验证为了验证算法收敛性，分别选用4组不同参数进行重建，并计算相对函数值增量和相对图像差异随迭代次数的变化：
| $\text { 相对函数值增量 }=\frac{\left|\psi\left(\boldsymbol{x}^n\right)-\psi\left(\boldsymbol{x}^0\right)\right|}{\left|\psi\left(\boldsymbol{x}^*\right)-\psi\left(\boldsymbol{x}^0\right)\right|} \text {, }$ | (24) |
| $\text { 相对图像差异 }=\frac{\left|\boldsymbol{x}^{n+1}-\boldsymbol{x}^n\right|}{\left|\boldsymbol{x}^{n+1}\right|} \text {. }$ | (25) |
 |
| 图 6 基于不同正则化参数重建的算法收敛性验证 |
| 图选项 |
3.3 结果分析和参数优化图 7展示了不同大小PSF的重建结果。由图可知，直接扣除散射数据可一定程度降低放射状伪影，由于准直器穿透效应，热球附近区域仍有明显计数，而使用PSF解卷积可以明显改善放射状伪影，降低热球附近区域计数。但是PSF过小时，仍会存在一定的放射状伪影，且显著放大背景区域噪声。PSF大小为127×127和255×255时的结果相对更好，背景噪声较低，无明显放射状伪影。
 |
| 注：1为对应原始数据的结果；2为对应原始数据直接扣除散射数据的结果；3—6分别为对应大小为31×31、63×63、127×127、255×255的PSF的重建结果。 图 7 不同大小PSF的重建图像 |
| 图选项 |
为了对结果进行定量评估，在重建图像的放射状条纹、条纹间隙和背景区域分别勾画了ROI，如图 8所示。通过计算大小为160 mm×50 mm的矩形ROI(蓝色)的标准差，对背景噪声进行评估。此外，为了定量描述放射状伪影的改善情况，在放射状条纹区域放置4个直径为11 mm的圆形ROI(黄色)，在条纹间隙放置3个直径为11 mm的圆形ROI(绿色)，并根据式(26)计算归一化伪影计数，具体可表示为
| $\text { 归一化伪影计数 }=\frac{(\text { 条纹区均值 }- \text { 条纹间隙均值 })_{\text {校正后 }}}{(\text { 条纹区均值 } - \text { 条纹间隙均值 })_{\text {校正前 }}}.$ | (26) |
 |
| 图 8 定量分析ROI勾画图 |
| 图选项 |
校正图像定量分析结果如表 2所示。由表 2可知，随着PSF增大，归一化伪影计数逐渐降低。但是相对于仅做散射校正的结果，大小为31×31的PSF处理结果的归一化伪影计数更大。对比图 7中编号2和3可知，基于大小为31×31的PSF解卷积虽然可以压低条纹间隙计数，但是对于条纹区域计数的压低效果并不理想，导致归一化伪影计数增大。另外，虽然更大的PSF可以降低背景噪声，但是相比于原始数据，过小的PSF反而会导致背景噪声更大。综合来看，PSF大小为127×127和255×255时效果十分接近, 明显优于另外2组选择。因此下一步将针对大小为127×127的PSF引入先验函数进行进一步优化。
表 2 校正图像定量分析结果
| 类别 | 归一化伪影计数 | 背景噪声 |
| 原始数据 | 1.00 | 2.45 |
| 散射校正 | 0.61 | 2.48 |
| 31×31 PSF | 0.81 | 3.46 |
| 63×63 PSF | 0.03 | 3.13 |
| 127×127 PSF | 0.01 | 2.06 |
| 255×255 PSF | 0.01 | 2.09 |
表选项
基于大小为127×127的PSF，图 9展示了不同β的重建图像，编号1的2条线表示剖线位置，用于后续分析。可以发现，随着β增大，背景噪声得到一定程度的压制，但是β过大时会导致过平滑。图 10展示了β为2.5×10-3、1.0×10-2、4.0×10-2共3组重建图像剖线图，剖线位置如图 9中编号1所示，图 10a对应水平剖线，图 10b对应竖直剖线。由图 10a可知，随着β增大，相邻的充满放射源的玻璃管中间背景区域的像素值会明显增大，导致对比度降低。由图 10b可知，随着β增大，条纹伪影区域的值会逐步增大。这是由于式(3)中的先验函数在降低背景噪声的同时还要保持边界，因此放射状伪影区域的条纹特征也会被当作边界进行一定程度的强化。图 11给出了不同β(0, 1.25×10-3，2.5×10-3，5.0×10-3，1.0×10-2，2.0×10-2，4.0×10-2)重建的定量分析结果。前期随着β增大，归一化伪影计数并无明显变化，背景噪声明显减小。但是随着β进一步增大，归一化伪影计数也会随之增大，这与图 10展示的条纹伪影被强化的现象一致。
 |
| 注：1—6分别对应β为1.25×10-3、2.5×10-3、5.0×10-3、1.0×10-2、2.0×10-2、4.0×10-2。 图 9 不同β的重建图像 |
| 图选项 |
 |
| 图 10 重建图像剖线图 |
| 图选项 |
 |
| 图 11 不同β重建的定量分析结果 |
| 图选项 |
综上所述，最终选择PSF为127×127、β为1.0×10-2用于临床评估。此时基于Monte Carlo模拟数据的归一化伪影计数为0.04，背景噪声为0.96。
3.4 传统算法与本文算法对比对于2.2节所述传统算法和本文算法的差异，基于Monte Carlo模拟数据对2种算法的重建结果进行对比评估。选择同样的初始值、PSF(大小为127×127)和迭代次数(50次)，考察不同正则化强度、不同算法的重建结果，不同解卷积算法效果对比如图 12所示，β相同时传统算法背景噪声相对更大，且传统算法对β值较为敏感，β显著增大时重建结果会明显出现问题。相比之下，本文算法综合表现更好，β增大后图像变得更加平滑，符合预期，且没有出现传统算法结果中的点阵伪影。
 |
| 注：1、2、3对应传统解卷积算法结果，4、5、6对应本文算法结果，从左往右分别对应β为2.0×10-1、2.0×10-2和2.0×10-3。 图 12 不同解卷积算法效果对比 |
| 图选项 |
4 临床数据评估对于临床采集的点源数据和人体数据，分别采用本文算法进行处理。其中点源数据基于连续床位运动采集。表 3展示了临床数据各能窗计数的情况。图 13展示了临床数据的处理结果。采集的真实投影数据大小为1 024×256，图 13仅展示关注部分，大小为256×256。不论是点源数据还是人体数据，引入解卷积算法后放射状伪影均得到了有效压制。对于点源数据，图像分辨率得到了明显改善，并且由于计数较多，点源数据解卷积后并未明显导致背景噪声放大。相对而言，由于人体数据计数低、背景噪声大，单纯的解卷积算法会导致背景噪声显著增大，而引入先验函数可使背景噪声变小，在图像上表现为背景区域更加“干净”。
表 3 临床数据各能窗计数
| 临床数据 | 低能窗 | 主能窗 | 高能窗 |
| 点源数据 | 340 478 | 2 004 681 | 265 882 |
| 人体数据 | 32 745 | 154 612 | 20 885 |
表选项
 |
| 注：1—原始数据；2—仅做散射校正；3—解卷积；4—解卷积+先验正则化的结果。 图 13 临床数据的处理结果 |
| 图选项 |
图 14和表 4分别展示了临床数据ROI勾画示意图和临床数据定量分析结果。条纹伪影和条纹间隙区域ROI直径为11 mm, 背景区域ROI大小为160 mm×50 mm。由表 4可知，相对于仅做散射校正，解卷积校正结果的归一化伪影计数明显更小。对于低计数的人体数据，先验正则化可以显著改善解卷积导致的高背景噪声，与图 13观测结果一致，代价是略微提升归一化伪影计数。此外，点源数据本身放射状伪影区域较小，ROI区域条纹特征相对已不是特别明显，因此引入先验函数后反而会被进一步平滑，导致归一化伪影计数进一步降低。
 |
| 图 14 临床数据ROI勾画示意图 |
| 图选项 |
表 4 临床数据定量分析结果
| 类别 | 归一化伪影计数 | 背景噪声 | |||
| 点源数据 | 人体数据 | 点源数据 | 人体数据 | ||
| 原始数据 | 1.00 | 1.00 | 1.57 | 1.03 | |
| 散射校正 | 0.79 | 0.72 | 0.96 | 0.95 | |
| 散射估计+解卷积 | 0.40 | 0.22 | 0.83 | 1.74 | |
| 散射估计+解卷积+先验正则化 | 0.35 | 0.28 | 0.62 | 0.92 | |
表选项
5 结论针对碘-131 SPECT平片成像中的放射状伪影问题，本文提出了一种PSF解卷积结合先验正则化的方法，构建目标函数并推导了迭代算法。在Monte Carlo模拟实验中，对算法的单调收敛性进行了验证，研究了不同大小PSF、不同正则化参数β对重建结果的影响，并验证了本文算法相对于传统算法的优势。结果表明：虽然PSF可以有效降低放射状伪影，但也会引入背景噪声，PSF越小背景噪声越明显。而先验约束在降低背景噪声的同时也会一定程度保持放射状伪影特征。综合考虑选择PSF大小为127×127、正则化参数β为1.0×10-2时，Monte Carlo模拟数据处理结果的归一化伪影计数仅为0.04。临床数据处理结果也证明了算法的有效性，点源数据和人体数据处理后的归一化伪影计数分别为0.35、0.28。
本研究基于Monte Carlo模拟方法获取临床典型距离(10 cm)的PSF，而临床采集时实际位置会有一定偏差，这也可能是临床数据结果相对模拟数据较差的原因。此外，先验函数中的超参数δ用于区分边界和噪声，像素差异远大于该参数时，会被认为是图像边界进行维持，而像素差异远小于该参数时，会被当作背景噪声进行平滑。2组实验数据测试中，由于伪影特征强度不同，导致引入先验函数对伪影区域计数的影响也存在一定差异。下一步工作一方面要考察不同位置的PSF对算法结果的影响，另一方面也要基于更多临床数据对算法进行进一步的参数优化和有效性验证。
参考文献
| [1] | AHMADZADEHFAR H, BIERSACK H J, HERRMANN K. Clinical applications of SPECT-CT[M]. Berlin: Springer, 2014. |
| [2] | 金永杰, 马天予. 核医学仪器与方法[M]. 哈尔滨: 哈尔滨工程大学出版社, 2010. JIN Y J, MA T Y. Nuclear medicine instruments and methods[M]. Harbin: Harbin Engineering University Press, 2010. (in Chinese) |
| [3] | AVRAM A M. Radioiodine scintigraphy with SPECT/CT: An important diagnostic tool for thyroid cancer staging and risk stratification[J]. Journal of Nuclear Medicine Technology, 2014, 42(3): 170-180. |
| [4] | AHMED N, NIYAZ K, BORAKATI A, et al. Hybrid SPECT/CT imaging in the management of differentiated thyroid carcinoma[J]. Asian Pacific Journal of Cancer Prevention: APJCP, 2018, 19(2): 303-308. |
| [5] | BARRACK F, SCUFFHAM J, MCQUAID S. Septal penetration correction in 131I imaging following thyroid cancer treatment[J]. Physics in Medicine & Biology, 2018, 63(7): 075012. |
| [6] | DEY N, BLANC-FERAUD L, ZIMMER C, et al. Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[J]. Microscopy Research and Technique, 2006, 69(4): 260-266. DOI:10.1002/jemt.20294 |
| [7] | LOENING A M, GAMBHIR S S. AMIDE: A free software tool for multimodality medical image analysis[J]. Molecular Imaging, 2003, 2(3): 131-137. DOI:10.1162/153535003322556877 |
| [8] | ZHAO W R, WANG X Q, LIANG Y K, et al. Initial evaluation of three novel quantiative imaging utilities from a newly developed SPECT camera in China[J]. The Journal of Nuclear Medicine, 2018, 59(S1): 219. |
| [9] | LANGE K, HUNTER D R, YANG I. Optimization transfer using surrogate objective functions[J]. Journal of Computational and Graphical Statistics, 2000, 9(1): 1-20. |
| [10] | WANG G B, QI J Y. Penalized likelihood PET image reconstruction using patch-based edge-preserving regularization[J]. IEEE Transactions on Medical Imaging, 2012, 31(12): 2194-2204. DOI:10.1109/TMI.2012.2211378 |
| [11] | HUTTON B F, BUVAT I, BEEKMAN F J. Review and current status of SPECT scatter correction[J]. Physics in Medicine & Biology, 2011, 56(14): R85-R112. |
| [12] | DEWARAJA Y K, LJUNGBERG M, KORAL K F. Characterization of scatter and penetration using Monte Carlo simulation in 131I imaging[J]. The Journal of Nuclear Medicine, 2000, 41(1): 123-130. |
| [13] | JAN S, SANTIN G, STRUL D, et al. GATE: A simulation toolkit for PET and SPECT[J]. Physics in Medicine & Biology, 2004, 49(19): 4543-4561. |
