梁亚星,张 萌,陈 燕,桂志国,, 张 权,
(1.中北大学 电子测试技术国家重点实验室,山西 太原 030051;2.中北大学 生物医学成像与影像大数据山西省重点实验室,山西 太原 030051;3.中北大学 信息与通信工程学院,山西 太原 030051)
计算机断层成像技术(Computed Tomography,CT)是一种重要的疾病诊断手段.在医学临床诊断中,由于病变部位往往位于一个局部区域,同时诊断也仅需病变位置的结构信息[1].而在目前常规全局扫描的检查方式中,因扫描区域覆盖了病变部位之外的正常组织,X射线会对其造成不必要的辐射伤害.
相比之下,基于感兴趣区域的局部断层扫描成像[2]技术则显得更具临床应用意义.该技术通过局部扫描方式,即位于感兴趣区域外的发射源和探测器不工作,在一定程度上降低了辐射剂量和设备的硬件损耗,但同时也造成了投影数据的截断,在重建的局部图像中形成截断伪影,呈现为重建图像边界高亮的圆环,且含有灰度移位伪影[2-3].为消除截断伪影,目前多采用投影数据外插延拓的方法,依据扩充缺损数据方法和重建过程的不同,衍生了不同的算法[3-7].其中,陈云斌[3]、王克军[4]和罗戎蕾[7]利用每个角度下处于边缘的两个投影值分别向两端扩充,直至填充整个数据区域;薛少峰[5]利用截断数据边界已知的若干连续投影数据进行曲线拟合,并根据该曲线插值采样,获得缺损投影数据.此外,王浩[6]则利用级比生成法获取截断投影数据相邻部分的缺损投影数据.
根据文献[8]基于投影域校正方法去除CT图像环形伪影,文献[9]利用霍夫变换检测井壁图像的平面地质特征,以及文献[10]引入权重函数模型融合不同滤波手段处理后的CT图像的启发,区别于依据某一投影角度的已知投影数据进行截断投影数据扩充的方法,本文在分析截断投影数据特点的基础上,根据投影数据的正弦结构边缘信息,提出一种新的局部重建算法.
图1(a)是对Shepp-Logan体模进行360°全局扫描得到的投影数据,对其截断得到图1(b)所示的局部扫描数据.
图1 仿真数据
投影数据为正弦曲线的叠加,在局部扫描条件下,截断投影数据仍保留正弦结构信息.因此,本文利用从截断投影数据提取的正弦边缘信息作为扩充边界,将投影数据按照图1(b)箭头指示方向向外扩充至白色正弦曲线边界处.显然正弦边缘检测是其中的关键,本文算法则采用霍夫变换进行检测.
算法先对截断数据进行边缘提取,接着基于霍夫变换对其进行正弦曲线检测,利用确定的特征参数修复缺失的正弦曲线部分.
正弦曲线可表示为
y=Asin(wx-θ)+y0,
(1)
式中:A为振幅;ω为角速度;θ为初始相位;y0为正弦曲线的基线;T为周期.
假定P1(x1,y1)为图2(a)中正弦曲线上任意一点,图中P2(x2,y2),P3(x3,y3),其中x2=x1+T/2,x3=x1+T/4,则有
y0=(y1+y2)/2,
(2)
(y1-y0)/(y3-y0)=tan(wx1-θ),
(y1-y0)/(y2-y0)=tan(wx2-θ),
(3)
(4)
对任意的x1,依据式(2)~式(4)依次求解基线y0,相位θ,振幅A.从而获取图1(b)中的正弦边界.在此基础上,将截断投影数据上、下边界处的投影数据原值外推,得到一次扩充的投影数据,如图2(b)所示.但由于霍夫变换检测的局限,发现图2(b)中仍存在箭头所示的零值区域.为此,算法采用线性插值法进行二次投影数据修复.
图2 正弦曲线和修复后投影数据
针对图2(b)中修复不完全问题,基于两侧已知的投影数据进行线性插值,进行二次修复.
若n1、n2分别是第i个探测器上缺失零值区域左右边界的非零值位置,则第i行缺失区域[n1+1,n2-1]内任一位置n处的投影值为
P(i,n+n1)=P(i,n1)+((P(i,n2)-
P(i,n1))/(n2-n1-1))*n.
(5)
扩充后的投影数据如图2(c)所示.从中可见,截断投影数据得到了较好的扩充.
算法可总结为投影数据基于正弦边界的一次修复、基于线性插值的二次修复,以及依据修复后的投影数据经滤波反投影重建得到CT局部重建图像.算法具体步骤为:
1)一次修复:首先采用Canny算子对截断投影数据P1进行边缘检测,获得二值图像;采用霍夫变换对二值图像进行正弦曲线检测、计算获得正弦曲线特征参数:基线、相位和振幅;将截断投影数据边界处的投影数据原值外推至检测出的正弦边界,得到第一次修复结果P2;
2)二次修复:依据式(5),利用未修复完全的零值区域的左右邻近的非零值数据,进行线性插值,得到第二次修复结果P3;
3)对经1)和2)两次修复后的投影数据P3,采用滤波反投影算法重建,得到CT局部重建图像.
实验中,采用Shepp-Logan头部体模和实际骨盆数据,验证算法进行局部重建的可行性.实验编程环境:64位Window 10操作系统,处理器Intel(R)Core(TM)i7-5500U CPU @2.4 GHz,内存4 G.
为了作对比验证,将截断数据分别进行常数延拓、局部均值延拓、对称延拓,并重建图像.其中,常数延拓将截断正弦图上、下边缘处的投影数据作为常数向两侧扩充;局部均值延拓首先求取截断投影数据上、下边界处的正下方、正上方3*3邻域内已知的投影数据的均值,再将求取的均值进行外推扩充;对称延拓则是以截断数据边界处的数据作为对称点,以对称的方式,分别对上、下边界之外的缺失数据进行扩充.
此外,为了定量评价算法的有效性,采用归一化均方距离判据NMSD、归一化平均绝对距离判据NAAD、均方误差MSE、峰值信噪比PSNR和结构相似性SSIM进行进一步验证.
(6)
(7)
(8)
(9)
(10)
式中:ui,j和vi,j表示原始图像和重建图像第i行、j列的像素灰度值;μu和μv表示原始图像和重建后图像的平均灰度值;σu和σv表示原始图像和重建后图像的标准差,图像大小N×N.n是图像灰度级数.
实验采用大小为256×256的Shepp-Logan仿真模型.选取以图像中心为圆心、半径为50像素的圆形区域为待重建的感兴趣区域.采用等角度扇束扫描,投影角度360°,探元数目246个,选取中间68个探元数据为截断后投影数据.
图3(a)、(b)分别是原始体模及其感兴趣区域放大显示图.图3(c)为未截断的完整投影数据的FBP重建结果的感兴趣区域,图3(d)则为截断数据的FBP重建结果图.显然,图3(d)边界处产生了环状的高亮截断伪影,且图像整体模糊,对比度下降.图3(e)、图3(f)分别给出了图3(c)与图3(b)、图3(d)与图3(b)的差值图,从中也观察到截断投影重建结果边界处的高亮环状伪影及损失的边缘细节.
图3 仿真模型、FBP重建图和ROI差值图
图4 为常数延拓、局部均值延拓、对称延拓和本文算法四种算法的CT局部重建结果对比图.左列为FBP重建结果,右列为四种方法重建结果分别与原图相应的感兴趣区域图3(b)的差值图.可以观察到图4(b)、(d)、(f)中矩形框标注区域内存在比较明显的边界特征,相比而言,图4(h)边界信息较少.这说明,本文算法重建图像与原图最为相似,抑制了局部重建图像中边界处的环状截断伪影,且边缘细节信息损失较少、对比度较高.
图4 仿真数据FBP局部重建图及对应差值图
表1 仿真数据各局部重建算法定量评价参数
表1 列出了各局部重建算法定量评价参数.分析比较可见,本文算法的NMSD、NAAD和MSE值最低,PSNR值和SSIM值最高.这进一步说明了本文重建图像与原始图像差异程度小,失真程度小.同时也验证了利用截断数据中提取的正弦边缘信息进行数据修复的有效性.
实验采用实际骨盆数据,选取以图像中心为圆心、半径为50像素的圆形区域为待重建的感兴趣区域,重建图像大小为256×256.投影角度360°,探元数目512个,选取中间180个探元数据为截断后投影数据.
图5(a)、(b)分别为完整投影数据和截断投影数据.图5(c)、(d)则分别为正弦边缘修复和线性插值两次扩充的投影数据.
图5 实际数据
图6 实际数据FBP重建图和差值图
图6(a)、(b)分别是完整投影数据FBP全重建及其感兴趣区域放大显示图.图6(c)为截断投影数据FBP重建结果的感兴趣区域,图6(d)则为图6(c)与图6(b)的差值图.图6(c)、(d)的边界处都存在环状亮伪影,这与仿真数据局部重建结果边界处存在有高亮环状伪影一致.
图7 左列依次为四种算法FBP重建结果图,右列依次为四种算法重建结果分别与原图相应的感兴趣区域图6(b)的差值图.本文算法重建结果图7(g)与原图更加相似,视觉效果较好.对比图7(b)、(d)和(f),图7(h)矩形框标注区域内存在更少的边界信息.
图7 实际数据FBP局部重建图及对应差值图
图8 给出了原图相应的感兴趣区域图6(b)分别与图6(b)、图6(c)和四种算法局部重建结果图的结构相似索引图.索引图中的像素越亮,则其SSIM值越接近于1,图像质量越好.比较图8(b)~(f),本文算法处理的图8(f)中大部分像素亮的程度比较高,这说明本文算法的结构相似性数值也最高.
图8 结构相似索引图
表2 列出了各局部重建算法定量评价参数,分析比较可见,本文算法NMSD、NAAD和MSE值最低,PSNR值和SSIM值最高.
表2 实际数据各局部重建算法定量评价参数
可见,实际数据实验结果的主观及客观评价均得到了与仿真数据实验一致的结论.算法的有效性进一步得到了验证.
本文提出了一种新的投影数据恢复扇束CT局部重建方法.利用截断投影数据的正弦结构边缘信息及线性插值对其进行修复,改善了直接使用截断数据重建图像造成的截断伪影问题.对比其他截断投影数据扩充方法,从局部重建图像的视觉效果及客观评价参数等方面,评估了算法的可行性.