基于正弦圖恢復的CT局部重建算法

梁亞星，張 萌，陳 燕，桂志國，, 張 權，

(1.中北大學 電子測試技術國家重點實驗室，山西 太原 030051；2.中北大學 生物醫(yī)學成像與影像大數(shù)據(jù)山西省重點實驗室，山西 太原 030051；3.中北大學 信息與通信工程學院，山西 太原 030051)

0 引 言

計算機斷層成像技術(Computed Tomography，CT)是一種重要的疾病診斷手段.在醫(yī)學臨床診斷中，由于病變部位往往位于一個局部區(qū)域，同時診斷也僅需病變位置的結(jié)構(gòu)信息[1].而在目前常規(guī)全局掃描的檢查方式中，因掃描區(qū)域覆蓋了病變部位之外的正常組織，X射線會對其造成不必要的輻射傷害.

相比之下，基于感興趣區(qū)域的局部斷層掃描成像[2]技術則顯得更具臨床應用意義.該技術通過局部掃描方式，即位于感興趣區(qū)域外的發(fā)射源和探測器不工作，在一定程度上降低了輻射劑量和設備的硬件損耗，但同時也造成了投影數(shù)據(jù)的截斷，在重建的局部圖像中形成截斷偽影，呈現(xiàn)為重建圖像邊界高亮的圓環(huán)，且含有灰度移位偽影[2-3].為消除截斷偽影，目前多采用投影數(shù)據(jù)外插延拓的方法，依據(jù)擴充缺損數(shù)據(jù)方法和重建過程的不同，衍生了不同的算法[3-7].其中，陳云斌[3]、王克軍[4]和羅戎蕾[7]利用每個角度下處于邊緣的兩個投影值分別向兩端擴充，直至填充整個數(shù)據(jù)區(qū)域；薛少峰[5]利用截斷數(shù)據(jù)邊界已知的若干連續(xù)投影數(shù)據(jù)進行曲線擬合，并根據(jù)該曲線插值采樣，獲得缺損投影數(shù)據(jù).此外，王浩[6]則利用級比生成法獲取截斷投影數(shù)據(jù)相鄰部分的缺損投影數(shù)據(jù).

根據(jù)文獻[8]基于投影域校正方法去除CT圖像環(huán)形偽影，文獻[9]利用霍夫變換檢測井壁圖像的平面地質(zhì)特征，以及文獻[10]引入權重函數(shù)模型融合不同濾波手段處理后的CT圖像的啟發(fā)，區(qū)別于依據(jù)某一投影角度的已知投影數(shù)據(jù)進行截斷投影數(shù)據(jù)擴充的方法，本文在分析截斷投影數(shù)據(jù)特點的基礎上，根據(jù)投影數(shù)據(jù)的正弦結(jié)構(gòu)邊緣信息，提出一種新的局部重建算法.

1 本文算法

1.1 算法基本思想

圖1(a)是對Shepp-Logan體模進行360°全局掃描得到的投影數(shù)據(jù)，對其截斷得到圖1(b)所示的局部掃描數(shù)據(jù).

圖1 仿真數(shù)據(jù)

投影數(shù)據(jù)為正弦曲線的疊加，在局部掃描條件下，截斷投影數(shù)據(jù)仍保留正弦結(jié)構(gòu)信息.因此，本文利用從截斷投影數(shù)據(jù)提取的正弦邊緣信息作為擴充邊界，將投影數(shù)據(jù)按照圖1(b)箭頭指示方向向外擴充至白色正弦曲線邊界處.顯然正弦邊緣檢測是其中的關鍵，本文算法則采用霍夫變換進行檢測.

1.2 霍夫變換

算法先對截斷數(shù)據(jù)進行邊緣提取，接著基于霍夫變換對其進行正弦曲線檢測，利用確定的特征參數(shù)修復缺失的正弦曲線部分.

正弦曲線可表示為

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，依據(jù)式(2)～式(4)依次求解基線y0，相位θ，振幅A.從而獲取圖1(b)中的正弦邊界.在此基礎上，將截斷投影數(shù)據(jù)上、下邊界處的投影數(shù)據(jù)原值外推，得到一次擴充的投影數(shù)據(jù)，如圖2(b)所示.但由于霍夫變換檢測的局限，發(fā)現(xiàn)圖2(b)中仍存在箭頭所示的零值區(qū)域.為此，算法采用線性插值法進行二次投影數(shù)據(jù)修復.

圖2 正弦曲線和修復后投影數(shù)據(jù)

1.3 線性插值

針對圖2(b)中修復不完全問題，基于兩側(cè)已知的投影數(shù)據(jù)進行線性插值，進行二次修復.

若n1、n2分別是第i個探測器上缺失零值區(qū)域左右邊界的非零值位置，則第i行缺失區(qū)域[n1+1,n2-1]內(nèi)任一位置n處的投影值為

P(i,n+n1)=P(i,n1)+((P(i,n2)-

P(i,n1))/(n2-n1-1))*n.

(5)

擴充后的投影數(shù)據(jù)如圖2(c)所示.從中可見，截斷投影數(shù)據(jù)得到了較好的擴充.

1.4 算法描述

算法可總結(jié)為投影數(shù)據(jù)基于正弦邊界的一次修復、基于線性插值的二次修復，以及依據(jù)修復后的投影數(shù)據(jù)經(jīng)濾波反投影重建得到CT局部重建圖像.算法具體步驟為：

1)一次修復：首先采用Canny算子對截斷投影數(shù)據(jù)P1進行邊緣檢測，獲得二值圖像；采用霍夫變換對二值圖像進行正弦曲線檢測、計算獲得正弦曲線特征參數(shù)：基線、相位和振幅；將截斷投影數(shù)據(jù)邊界處的投影數(shù)據(jù)原值外推至檢測出的正弦邊界，得到第一次修復結(jié)果P2；

2)二次修復：依據(jù)式(5)，利用未修復完全的零值區(qū)域的左右鄰近的非零值數(shù)據(jù)，進行線性插值，得到第二次修復結(jié)果P3；

3)對經(jīng)1)和2)兩次修復后的投影數(shù)據(jù)P3，采用濾波反投影算法重建，得到CT局部重建圖像.

2 實驗結(jié)果與分析

實驗中，采用Shepp-Logan頭部體模和實際骨盆數(shù)據(jù)，驗證算法進行局部重建的可行性.實驗編程環(huán)境：64位Window 10操作系統(tǒng)，處理器Intel(R)Core(TM)i7-5500U CPU @2.4 GHz，內(nèi)存4 G.

為了作對比驗證，將截斷數(shù)據(jù)分別進行常數(shù)延拓、局部均值延拓、對稱延拓，并重建圖像.其中，常數(shù)延拓將截斷正弦圖上、下邊緣處的投影數(shù)據(jù)作為常數(shù)向兩側(cè)擴充；局部均值延拓首先求取截斷投影數(shù)據(jù)上、下邊界處的正下方、正上方3*3鄰域內(nèi)已知的投影數(shù)據(jù)的均值，再將求取的均值進行外推擴充；對稱延拓則是以截斷數(shù)據(jù)邊界處的數(shù)據(jù)作為對稱點，以對稱的方式，分別對上、下邊界之外的缺失數(shù)據(jù)進行擴充.

此外，為了定量評價算法的有效性，采用歸一化均方距離判據(jù)NMSD、歸一化平均絕對距離判據(jù)NAAD、均方誤差MSE、峰值信噪比PSNR和結(jié)構(gòu)相似性SSIM進行進一步驗證.

(6)

(7)

(8)

(9)

(10)

式中：ui,j和vi,j表示原始圖像和重建圖像第i行、j列的像素灰度值；μu和μv表示原始圖像和重建后圖像的平均灰度值；σu和σv表示原始圖像和重建后圖像的標準差，圖像大小N×N.n是圖像灰度級數(shù).

2.1 仿真數(shù)據(jù)實驗

實驗采用大小為256×256的Shepp-Logan仿真模型.選取以圖像中心為圓心、半徑為50像素的圓形區(qū)域為待重建的感興趣區(qū)域.采用等角度扇束掃描，投影角度360°，探元數(shù)目246個，選取中間68個探元數(shù)據(jù)為截斷后投影數(shù)據(jù).

圖3(a)、(b)分別是原始體模及其感興趣區(qū)域放大顯示圖.圖3(c)為未截斷的完整投影數(shù)據(jù)的FBP重建結(jié)果的感興趣區(qū)域，圖3(d)則為截斷數(shù)據(jù)的FBP重建結(jié)果圖.顯然，圖3(d)邊界處產(chǎn)生了環(huán)狀的高亮截斷偽影，且圖像整體模糊，對比度下降.圖3(e)、圖3(f)分別給出了圖3(c)與圖3(b)、圖3(d)與圖3(b)的差值圖，從中也觀察到截斷投影重建結(jié)果邊界處的高亮環(huán)狀偽影及損失的邊緣細節(jié).

圖3 仿真模型、FBP重建圖和ROI差值圖

圖4 為常數(shù)延拓、局部均值延拓、對稱延拓和本文算法四種算法的CT局部重建結(jié)果對比圖.左列為FBP重建結(jié)果，右列為四種方法重建結(jié)果分別與原圖相應的感興趣區(qū)域圖3(b)的差值圖.可以觀察到圖4(b)、(d)、(f)中矩形框標注區(qū)域內(nèi)存在比較明顯的邊界特征，相比而言，圖4(h)邊界信息較少.這說明，本文算法重建圖像與原圖最為相似，抑制了局部重建圖像中邊界處的環(huán)狀截斷偽影，且邊緣細節(jié)信息損失較少、對比度較高.

圖4 仿真數(shù)據(jù)FBP局部重建圖及對應差值圖

表1 仿真數(shù)據(jù)各局部重建算法定量評價參數(shù)

表1 列出了各局部重建算法定量評價參數(shù).分析比較可見，本文算法的NMSD、NAAD和MSE值最低，PSNR值和SSIM值最高.這進一步說明了本文重建圖像與原始圖像差異程度小，失真程度小.同時也驗證了利用截斷數(shù)據(jù)中提取的正弦邊緣信息進行數(shù)據(jù)修復的有效性.

2.2 實際數(shù)據(jù)實驗

實驗采用實際骨盆數(shù)據(jù)，選取以圖像中心為圓心、半徑為50像素的圓形區(qū)域為待重建的感興趣區(qū)域，重建圖像大小為256×256.投影角度360°，探元數(shù)目512個，選取中間180個探元數(shù)據(jù)為截斷后投影數(shù)據(jù).

圖5(a)、(b)分別為完整投影數(shù)據(jù)和截斷投影數(shù)據(jù).圖5(c)、(d)則分別為正弦邊緣修復和線性插值兩次擴充的投影數(shù)據(jù).

圖5 實際數(shù)據(jù)

圖6 實際數(shù)據(jù)FBP重建圖和差值圖

圖6(a)、(b)分別是完整投影數(shù)據(jù)FBP全重建及其感興趣區(qū)域放大顯示圖.圖6(c)為截斷投影數(shù)據(jù)FBP重建結(jié)果的感興趣區(qū)域，圖6(d)則為圖6(c)與圖6(b)的差值圖.圖6(c)、(d)的邊界處都存在環(huán)狀亮偽影，這與仿真數(shù)據(jù)局部重建結(jié)果邊界處存在有高亮環(huán)狀偽影一致.

圖7 左列依次為四種算法FBP重建結(jié)果圖，右列依次為四種算法重建結(jié)果分別與原圖相應的感興趣區(qū)域圖6(b)的差值圖.本文算法重建結(jié)果圖7(g)與原圖更加相似，視覺效果較好.對比圖7(b)、(d)和(f)，圖7(h)矩形框標注區(qū)域內(nèi)存在更少的邊界信息.

圖7 實際數(shù)據(jù)FBP局部重建圖及對應差值圖

圖8 給出了原圖相應的感興趣區(qū)域圖6(b)分別與圖6(b)、圖6(c)和四種算法局部重建結(jié)果圖的結(jié)構(gòu)相似索引圖.索引圖中的像素越亮，則其SSIM值越接近于1，圖像質(zhì)量越好.比較圖8(b)～(f)，本文算法處理的圖8(f)中大部分像素亮的程度比較高，這說明本文算法的結(jié)構(gòu)相似性數(shù)值也最高.

圖8 結(jié)構(gòu)相似索引圖

表2 列出了各局部重建算法定量評價參數(shù)，分析比較可見，本文算法NMSD、NAAD和MSE值最低，PSNR值和SSIM值最高.

表2 實際數(shù)據(jù)各局部重建算法定量評價參數(shù)

可見，實際數(shù)據(jù)實驗結(jié)果的主觀及客觀評價均得到了與仿真數(shù)據(jù)實驗一致的結(jié)論.算法的有效性進一步得到了驗證.

3 結(jié) 論

本文提出了一種新的投影數(shù)據(jù)恢復扇束CT局部重建方法.利用截斷投影數(shù)據(jù)的正弦結(jié)構(gòu)邊緣信息及線性插值對其進行修復，改善了直接使用截斷數(shù)據(jù)重建圖像造成的截斷偽影問題.對比其他截斷投影數(shù)據(jù)擴充方法，從局部重建圖像的視覺效果及客觀評價參數(shù)等方面，評估了算法的可行性.