離散化代數(shù)重建的全變差算法改進(jìn)

張 蕾，徐伯慶

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

?

離散化代數(shù)重建的全變差算法改進(jìn)

張 蕾，徐伯慶

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

文中提出一種離散化代數(shù)重建的全變差改進(jìn)算法，主要針對(duì)離散化的代數(shù)重建算法易受到噪聲等因素的影響而使圖像邊緣較為模糊的問題，利用全變差最小化的約束條件，提出一種改進(jìn)的DART重建算法。實(shí)驗(yàn)表明，該算法與傳統(tǒng)ART算法相比，能較快重建出圖像，與DART算法相比，改善圖像邊緣模糊的情況，具有較好的抗噪性。

圖像重建；DART；TV

計(jì)算機(jī)斷層成像技術(shù)(CT)無論在醫(yī)學(xué)放射診斷還是在工業(yè)領(lǐng)域均有著廣泛應(yīng)用[1]。根據(jù)Radon變換得到的投影數(shù)據(jù)來重建圖像的算法主要分為兩類：解析算法和迭代算法。解析算法中較有名的是濾波反投影(Filtered Back Projection, FBP)算法[2-3]，分辨率高，成像速度快，目前仍應(yīng)用廣泛，但要求投影數(shù)據(jù)采集密集。迭代算法主要是代數(shù)重建(Algebraic Reconstruction Algorithm，ART)算法[3-4]，其結(jié)構(gòu)簡(jiǎn)單；缺點(diǎn)是計(jì)算量較大，重建速度慢，對(duì)計(jì)算機(jī)的內(nèi)存和運(yùn)算速度要求較高。

ART算法改進(jìn)的關(guān)鍵是加快迭代收斂速度和減少迭代次數(shù)。影響ART算法的因素有以下幾點(diǎn)：投影系數(shù)的求取、投影次序的訪問順序、松弛因子以及先驗(yàn)知識(shí)的影響等[5-6]。實(shí)際醫(yī)學(xué)CT中的圖像通常僅由兩種或幾種已知的灰度值組成，利用這些先驗(yàn)知識(shí)便可降低重建條件，從更少的投影數(shù)據(jù),或相同投影數(shù)據(jù)下迭代次數(shù)少時(shí)，更精確地重建圖像。近年來，由K. J. Batenburg 等人[7-8]提出的離散迭代算法(Discrete Algebraic Reconstruction Technique, DART)就是利用該先驗(yàn)知識(shí)，在加快了迭代收斂速度和減少迭代次數(shù)方面取得了顯著效果，但DART算法受到噪聲影響會(huì)使得圖像邊緣比較模糊，本文提出的DART改進(jìn)算法就是基于有限灰度值的先驗(yàn)知識(shí)快速重建圖像，并利用全變差(Total Variation，TV)最小化約束[9-10]進(jìn)一步抑制噪聲，使得圖像邊緣更加清晰。

1 算法理論

1.1 ART算法原理

設(shè)f(x,y)為待重建圖像，將一個(gè)n×n的方形網(wǎng)格疊加到f(x,y)上，每個(gè)網(wǎng)格內(nèi)的f值均為常量，如圖1所示。

圖1 重建圖像坐標(biāo)示意圖

掃描射線可由式xcosθ+ysinθ=t表示，其中θ為射線與x軸所夾的角，t為原點(diǎn)到射線的距離。射線寬度為δ。投影值函數(shù)，即Radon變換可由式(1)表示

(1)

在ART重建過程中，設(shè)射線數(shù)目為M，所得到的投影數(shù)據(jù)p={p1,p2,…pM},pi為第i條射線的投影值。wij為第i條射線與第j個(gè)像素網(wǎng)格所交的面積，稱為投影系數(shù)

(2)

投影值向量和原圖像的關(guān)系可由式(3)的代數(shù)方程組表示

(3)

W(M×N)稱為投影系數(shù)矩陣。式(2)的矩陣形式如下

W·F=P

(4)

ART迭代算法最初于1937年由Kaczmarz[11]提出。其基本思想是松弛法，首先設(shè)定一列初值F(0)，然后代入式(5)進(jìn)行迭代

(5)

其中，i=0,1,2,…,M,k為當(dāng)前迭代次數(shù)，j∈(1,N)，λ為松弛因子。經(jīng)過反復(fù)迭代式(5)，不斷修正重建值來獲取較精確的重建圖像。

1.2 DART算法原理

由于實(shí)際醫(yī)學(xué)重建或工業(yè)無損檢測(cè)中，被測(cè)圖像通常僅由一個(gè)或有限個(gè)灰度值組成，因此圖像重建問題可簡(jiǎn)化為離散斷層成像(Discrete Tomography，DT)。近年來，由K. J. Batenburg 等人[7-8]研究的基于ART迭代的DART重建算法有顯著的效果。

設(shè)原圖像的灰度值取值范圍為V={v1,v2,…v1},定義臨界值υ={υ1,υ2,…υ1}，其中

(6)

通過上述定義設(shè)函數(shù)T對(duì)圖像像素值進(jìn)行離散化

(7)

以二值圖像為例，DART算法的具體步驟如下：

(1)首先對(duì)待測(cè)圖像進(jìn)行ART迭代獲得初始迭代值F(0)；

(2)通過式(7)對(duì)F(0)進(jìn)行離散化，得到離散像素值S(0)=TF(0)；

(3)根據(jù)8連通域方法將圖像分割為邊緣像素點(diǎn)集合B和平滑區(qū)域集合I(若像素點(diǎn)與其8鄰域中至少有一個(gè)像素點(diǎn)的像素值不同，則稱其為邊緣像素點(diǎn)，否則為平滑區(qū)域)；對(duì)F(0)做如下處理

(8)

(4)以處理后的F(0)作為初始值重新利用式(5)進(jìn)行ART迭代得到F(1)，更新F(1)，保持平滑區(qū)域I值不變，即

(9)

(5)重復(fù)步驟(2)到步驟(4)。

1.3 全變差(TV)原理

實(shí)際圖像重建過程中，會(huì)不可避免地受到噪聲影響，因此去噪問題非常重要。由Rudin等人提出的全變差(TV)圖像去噪模型[12-13]已在圖像復(fù)原領(lǐng)域得到廣泛的應(yīng)用和研究。目前，有諸多方法可用來最小化全變差圖像去噪模型的能量函數(shù)[14-15]。

去噪模型可表示為

min‖F(xiàn)‖TV,s.t.W·F=P

(10)

式(10)與式(4)是對(duì)應(yīng)的。由式(10)可知，求解min‖F(xiàn)‖TV是關(guān)鍵。而實(shí)際min‖F(xiàn)‖TV不能直接計(jì)算得到，一般均是通過用范數(shù)L1來求解

‖F(xiàn)‖TV=‖F(xiàn)‖l=∑xy|Fx,y|

(11)

(12)

L1范數(shù)是圖像所有像素點(diǎn)的像素值之和，可通過梯度下降法獲取，過程如下

(13)

其中，ε=10-8。通過上式可求得min‖F(xiàn)‖TV。

2 算法改進(jìn)

本文在上述算法的基礎(chǔ)上提出DART +TV算法，在DART算法迭代過程中加入TV約束，抑制噪聲，使得圖像邊緣更加清晰。

在DART+TV算法的實(shí)現(xiàn)過程中：首先利用DART對(duì)圖像像素迭代值進(jìn)行離散化，并通過比較八鄰域像素值區(qū)分邊緣區(qū)域B和平滑區(qū)域I,若像素點(diǎn)屬于平滑區(qū)域I，則其像素值賦值為離散化后的值，否則其像素值保持迭代值不變，將處理后的整幅圖像像素值重新進(jìn)行ART迭代，只更新邊緣區(qū)域B的值保持平滑區(qū)域I內(nèi)的值不變，然后用全變差(TV)最小化約束條件對(duì)圖像像素值進(jìn)行處理，然后對(duì)此圖像進(jìn)行下一步DART迭代處理，直至像素值滿足約束條件則循環(huán)結(jié)束。算法步驟如下：

圖2 本文算法流程圖

算法流程圖如圖2所示。(1)首先對(duì)待測(cè)圖像進(jìn)行ART迭代獲得初始迭代值F(0)；(2)通過式(7)對(duì)F(0)進(jìn)行離散化，得到離散像素值S(0)=TF(0)；(3)根據(jù)8連通域方法將圖像分割為邊緣像素點(diǎn)集合B和平滑區(qū)域集合I(若像素點(diǎn)與其8鄰域中至少有一個(gè)像素點(diǎn)的像素值不同，則稱其為邊緣像素點(diǎn)，否則為平滑區(qū)域)；對(duì)F(0)做如下處理;(4)以處理后的F(0)作為初始值重新利用公式進(jìn)行ART迭代得到F(1)，更新F(1))，保持平滑區(qū)域I值不變；(5)利用式(12)對(duì)F(1)通過梯度下降法求導(dǎo)，進(jìn)行最小化約束；；(6)重復(fù)步驟(2)～步驟(5)。

3 實(shí)驗(yàn)結(jié)果

接下來給出3組實(shí)驗(yàn)，用于對(duì)傳統(tǒng)ART算法和DART算法和利用全變差(TV)約束處理的DART改進(jìn)算法進(jìn)行比較，以說明改進(jìn)后的算法的優(yōu)勢(shì)，在每組實(shí)驗(yàn)中，投影角度均每隔 取值。

(1) 第1組實(shí)驗(yàn)。在這組實(shí)驗(yàn)中，使用的原圖像是分辨率為 的正方形二值模型，如圖3(a)所示。圖3中分別列出了傳統(tǒng)ART算法與DART算法和DART改進(jìn)算法迭代后的結(jié)果。

圖3 實(shí)驗(yàn)結(jié)果圖

如圖3 (a)是原始圖像，圖3(b)是直接進(jìn)行ART迭代的實(shí)驗(yàn)結(jié)果，圖3(c)是進(jìn)行DART迭代的實(shí)驗(yàn)結(jié)果，圖3(d)圖是進(jìn)行DART改進(jìn)算法迭代的實(shí)驗(yàn)結(jié)果。

實(shí)驗(yàn)中，本文用峰值信噪比(PSNR)值來定量評(píng)價(jià)重建圖像與原始圖像的誤差，客觀評(píng)價(jià)出重建圖像的質(zhì)量。表1顯示了3種算法的實(shí)驗(yàn)參數(shù)指標(biāo)對(duì)比，參數(shù)PSNR值越大說明圖像噪聲越小，圖像質(zhì)量越好，k為迭代次數(shù)。

表1 參數(shù)指標(biāo)對(duì)比圖

從表1的實(shí)驗(yàn)結(jié)果可看出，本文算法和傳統(tǒng)ART算法和DART算法相比均有所改進(jìn)。

(2) 第2組實(shí)驗(yàn)。在這組實(shí)驗(yàn)中，使用的原圖像是分辨率為 的正方形四灰度值模型。I=zeros (100,100);I(10:30,10:30)=1;I(10:30,60:80)=2;I(60:80,10:30)=3;I(60:80,60:80)=1。

圖4分別列出了傳統(tǒng)ART算法、DART算法和DART改進(jìn)算法迭代后的結(jié)果。

圖4 實(shí)驗(yàn)結(jié)果圖

圖5 PSNR曲線對(duì)比

圖5中畫出了3種算法的PSNR值隨迭代次數(shù)的變化曲線圖。

(3) 第3組實(shí)驗(yàn)。在這組實(shí)驗(yàn)中，使用的原圖像是分辨率為 的Shepp-Logan模型，如圖6(a)所示。

圖 6 中分別列出了傳統(tǒng)ART算法與DART算法和DART改進(jìn)算法迭代后的結(jié)果。

圖6 實(shí)驗(yàn)結(jié)果圖

圖7 PSNR曲線對(duì)比

圖7為3種算法的PSNR值隨迭代次數(shù)的變化曲線圖。

(4) 實(shí)驗(yàn)分析。在以上3組實(shí)驗(yàn)中，第1、2組采用的是灰度級(jí)較少的兩幅圖像，而第3組采用的是較為復(fù)雜的圖像。從這3組實(shí)驗(yàn)的PSNR值變化曲線可看出，使用改進(jìn)的DART算法能比傳統(tǒng)的ART算法和DART算法較快地達(dá)到較好的質(zhì)量，改進(jìn)后的算法在前幾次迭代時(shí)重建質(zhì)量提高得比較快，一般迭代5、6次后PSNR值趨于平穩(wěn)。

尤其從1、2組的比較中可明顯看出，當(dāng)原圖像的灰度值較少時(shí)，使用DART改進(jìn)后的算法在重建效果方面有大幅提高。

4 結(jié)束語

圖像重建領(lǐng)域追求的是重建質(zhì)量和重建速度這兩個(gè)指標(biāo)。本文利用一般待重建圖像灰度級(jí)較少的先驗(yàn)條件，對(duì)迭代初值進(jìn)行離散化分割，將重建問題簡(jiǎn)單化，并利用全變差約束條件，改善邊緣模糊的情況，且在仿真實(shí)驗(yàn)中證實(shí)了其相對(duì)于傳統(tǒng)ART算法和DART算法的優(yōu)勢(shì)，即重建質(zhì)量提高了、更少的迭代次數(shù)就能達(dá)到較好的質(zhì)量，這正契合了改進(jìn)的目標(biāo)。

[1] 莊天戈.CT原理與算法[M].上海:上海交通大學(xué)出版社,1992.

[2] 曾更生.醫(yī)學(xué)圖像重建[M].北京:高等教育出版社,2010.

[3] Chetih N,Messali Z.Tomographic image reconstruction using Filtered Back Projection (FBP) and algebraic reconstruction technique (ART)[C].Germany:3rd International Conference on Control, Engineering & Information Technology,IEEE,2015.

[4] Jiang M,Wang G.Development of iterative algorithms for image reconstruction[J].Journal of X-Ray Science and Technology,2002, 10(10):77-86.

[5] Andersen A H. Algebraic reconstruction in CT from limited views[J].IEEE Transactions on Medical Imaging,1989,8(1):50-55.

[6] 秦中元,牟軒沁,王平,等.一種內(nèi)存優(yōu)化的代數(shù)重建算法及其快速實(shí)現(xiàn)[J].電子學(xué)報(bào),2003,31(9):1327-1329.

[7] Batenburg K J,Sijbers J.Dart:a fast heuristic algebraic reconstruction algorithm for discrete tomography[C].Moskv:International Conference on Image Processing(ICIP),2007.

[8] Kees Joost B, Jan S. DART: a practical reconstruction algorithm for discrete tomography[J].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society, 2011, 20(9):2542-2553.

[9] Gopi V P,Fayiz T K,Palanisamy P. Regularization based CT image reconstruction using Algebraic techniques[C].Soul:International Conference on Electronics and Communication Systems, IEEE,2014.

[10] Sidky E Y,Kao C M,Pan X.Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J]. Journal of X-ray Science and Technology,2009,14(2):119-139.

[11] Kaczmarz B S.Angenherte auflsung von systemen linearer gleichungen[J].Bulletin International de l’Academie Polonaise des Sciences et des Lettres, Series A,1937(2):55-63.

[12] Chen G J,Leng S.Prior image constrained compressed sensing (PICCS):a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets [J].Medical Physics,2008,35(2):660-663.

[13] 郭靜鈺,齊宏亮,袁媛,等.先驗(yàn)圖像約束的有限角度CT圖像重建算法[J].核電子學(xué)與探測(cè)技術(shù),2014(12):1421-1424.

[14] Darbon J, Sigelle M. Exact optimization of discrete constrained total variation minimization problems[C].Auckland, New Zealand：Combinatorial Image Analysis, Internationalworkshop, 2004.

[15] Chambolle A.Total variation minimization and a class of binary MRF models[M].Springer Berlin Heidelberg: Energy Minimization Methods in Computer Vision and Pattern Recognition,2005.

An Improved Reconstruction Algorithm for Discrete Tomography

ZHANG Lei，XU Boqing

(School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

This paper presents an improved discretization Algebraic Reconstruction algorithm which bases on Total Variation, mainly for discretization Algebraic Reconstruction (DART) algorithm is easily affected by factors such as noise and make the problem of image edge blurred, using Total Variation (TV) to minimize the constraints, an improved DART Reconstruction algorithm. Experiments show that the algorithm is compared with the traditional ART algorithm, can quickly reconstruct image, compared with the DART algorithm, improved the condition of image edge blur and has good noise resistance.

image reconstruction; DART; TV

10.16180/j.cnki.issn1007-7820.2016.12.034

2016- 02- 01

張蕾(1991-)，女，碩士研究生。研究方向：信號(hào)與信息處理。徐伯慶(1958-)，男，博士，副教授。研究方向：通信及圖像處理。

TP391.41

A

1007-7820(2016)12-121-05