基于定點DSP的ART算法實現研究

摘 要:數字信號處理器(DSP)在大量數據計算方面有明顯的速度優勢，為了提高CT ART算法的圖像重建速度，提出基于定點DSP TMS320C6416的算法計算方案。該方案結合ART算法的原理，以及USB 2.0協議在數據傳輸和定點DSP在數據計算上的速度優勢，通過CYPRESS的CY7C68001接口芯片，同時采用SDRAM作為數據緩存，進行DSP與PC機之間大量數據的正確高速傳輸;通過C語言編寫DSP片內算法計算程序并對其進行優化，以降低圖像重建時間。最終實驗以Shepp_Logan頭部剖析圖為原始圖像，進行重建算法計算，實驗結果證實了方案的可行性。關鍵詞:DSP; CT圖像重建; ART算法; SDRAM

中圖分類號:TN915-34文獻標識碼:A

文章編號:1004-373X(2010)18-0017-04

Implementation of ART Algorithm Based on Fixed Point DSP

SUN Yi-gang， WANG Qing-yong， ZHANG Hong-ying

(Aeronautical Automation College， Civil Aviation University of China， Tianjin 300300， China)

Abstract: DSP has high speed in the computation of mass data， a computation scheme based on fixed point DSP TMS320C6416 algorithm is proposed to increase the speed of ART image reconstruction in CT technology. This method based on the ART theory， and the merit in speed of the USB 2.0 in the data transmission and the fixed DSP in data processing. It takes the CYPRESS CY7C68001 as USB 2.0 interface chipset， the SDRAM as data buffer to make sure the high speed data communication between the DSP and the PC host; and runs and optimizes the ART algorithm with C language to reduce the time consumption in DSP computation. The experiment makes the Shepp_Logan model as the original image and performs the computation of reconstruction algorithm， and the result indicates that the method is feasible.Keywords: DSP; CT image reconstruction; ART algorithm; SDRAM

0 引 言

計算機層析成像技術(CT)已經廣泛應用在醫學及工業檢測領域。圖像重建是CT技術的關鍵內容，有兩類方法:解析算法和迭代算法。迭代算法中比較常用的是代數迭代算法(algebraic reconstruction technique，ART)。ART算法適合于不完全投影數據的圖像重建，抗噪聲干擾能力強，還可以結合一些先驗知識求解;但不足是計算量大，重建速度慢，已成為該算法應用發展的瓶頸[1]。若射線條數為I，ART收斂于最優解約需(3~8)I次的迭代，在ART重建方法中，對于一幅n×n的圖像，取m個投影，每個方向的投影有n條射線，如果直接用ART方法重建，則系數矩陣的元素個數約為o(n4);對于一幅256×256或者更大的圖像，計算量巨大。目前對ART算法的研究主要是針對算法本身，還不能把算法直接固化到硬件中實現，該算法在通用計算機上進行圖像重建時速度慢，耗費時間長[2-3]。

選用高性能的處理器能夠有效提高圖像重建的速度，目前DSP技術已經廣泛應用各種領域，如圖像處理[4-5]，這主要是因為DSP技術適合應用于有大量的數據，并且數據需要較為快速的數學計算場合[6]。DSP技術在ART算法計算中的應用能夠有效縮短該算法的圖像重建時間，從而有助于進一步擴大CT技術的應用范圍。本文選用TMS320C6416是一款高性能的定點DSP，主頻高達1 GHz，能力處理可達8 000 MIPS，并且有專門的硬件乘法器，可滿足CT圖像重建中大量數據需要高速計算的要求。

1 ART算法原理

ART算法首先是將要重建的圖像離散化，即假設在所要重建的未知圖像上疊加一個n×n方格網，方格網內的每個單元表示相應的像素，像素值為常量，斷層截面與投影的幾何布置見圖1所示。

圖1 斷層截面與投影的幾何布置

該圖像重建算法可用下列線性方程組表示:

r11f1+r12f2+…+r1NfN=p1

r21f1+r22f2+…+r2NfN=p2



rM1f1+rM2f2+…+rMNfN=pM(1)

式中:N=n×n為像素總數;M為投影射線總數;rij為投影系數，即加權因子，反映了第j個像素對i根射線積分的貢獻;fj表示第j(j=1，2，…，N)個方格(像素)內的常量值。式(1)寫成i，j=1，2，…，N矩陣形式為:

RF=P(2)

式中:R=(rij)M×N為投影系數矩陣;F=(f1，f2，…，fN)T為圖像矢量;P=(p1，p2，…，pM)T為投影矢量。對于式(1)的求解一般不直接通過矩陣求逆的方法求解，而采用迭代算法。迭代公式如式(3)所示:

Xi=Xi-1-αRiXi-1-piRi#8226;RiRi(3)

式中:i表示射線號;X表示需要重建的圖像向量;Ri表示投影系數矩陣的i行;α代表松弛系數，取值范圍為(0，2)。

2 CT圖像重建ART算法的DSP實現

2.1 ART算法數據的定點DSP表示

ART迭代算法中含有大量的實數數據，包含有整數部分和小數部分。結合定點DSP和ART算法的特點，本文采用了兩種數據表示方法:定標法和直接表示法。

定標法，通過小數點的不同位置來表示不同的實數，小數點前為整數部分，其后為小數部分。本文采用32位定點DSP，所能表示的整數數值范圍為[-232-1， 232-1-1]。定點DSP在表示算法數據時，需要考慮計算數據的數值范圍[7]。針對定標法，考慮了該算法的如下參量值[1]:

(1) 投影系數矩陣元素，見圖1。它表示所要重建圖像的像素對射線積分的貢獻，值等于射線穿過像素所占部分的面積，其數值取值范圍為[0，1];

(2) 投影矢量pi，pi=∑jrijxj。表示i根射線經過的所有像素總的投影之和。其中，rij為投影系數矩陣元素;xj為j號像素的像素值。為了說明情況，設定rij數值均為1，將xj設為灰度值255，pi為232-1-1，求得j的總數約為1.7×107(圖像數據矩陣約4 000×4 000)。換言之，當投影矢量元素為定標法所能表示32位定點數的最大值時，圖像數據矩陣約為4 000×4 000，每個像素的灰度值為最大值255;投影系數矩陣元素rij均為1。實際中，圖像數據矩陣一般小于4 000×4 000，并且投影系數矩陣是含有大量零元素的稀疏矩陣[1]，故投影矢量pi可用32位定點數表示。

(3) 松弛系數。該值范圍為(0，2)，目前已經有大量的試驗表明，選擇低的松弛因子，一般為[0.05，0.25]，往往能獲得較好的重建結果，在某些投影噪聲很大的情況下，松弛因子會選得更低，本文選擇0.05;

(4) 圖像灰度。每一個像素的灰度值由8位表示，值域為[0，255]。

由上述4點可知，采用定標法能夠滿足ART算法中參量數據的數值范圍要求。

直接表示法:編譯環境支持直接表示浮點數據，即將相應的數據類型定義為float型，其在內存中是按指數形式存放的。定點DSP對于直接法表示的浮點數據，沒有硬件單元來直接進行浮點運算，是以整型數據計算的形式，通過軟件算法進行浮點數據運算的。

作者在實現定點DSP算法的計算中，首先使用直接法表示DSP接受到的該算法參量數據，并編寫算法程序，即將接收到的數據類型定義為float型，以確保接受正確數據;其后的優化算法則將參量float型浮點數據轉化為定點形式，以定標法為手段[4，7]優化算法。

總之，在ART算法圖像重建計算中，定點DSP能夠滿足算法參量數值范圍的要求，并能實現算法計算。

2.2 算法的DSP計算實現

本文以DSP和PC機的USB2.0協議數據傳輸及DSP計算過程來實現ART算法的定點DSP計算方案。

2.2.1 數據傳輸過程[8]

(1) 數據傳輸硬件實現。DSP作為USB的設備端;PC作為USB的主機端。ART算法每次計算所需數據量大，故首先將探測器得到的數據，即投影系數矩陣數據從USB主機端通過USB 2.0接口傳給SDRAM緩存，DSP從SDRAM取數進行算法計算，計算所得的最終結果存放在SDRAM，然后再通過USB 2.0接口返回給主機端，顯示重建后的結果圖像。數據傳輸過程的硬件框圖如圖2所示。

圖2 數據傳輸過程的硬件框圖

在圖2中，SDRAM為大容量的64位緩存，配置在DSP的EMIFA(外部存儲器接口A，64位的數據寬度)的CE0空間，DSP映射地址范圍為0x80000000~0x8FFFFFFF。在USB 2.0接口芯片(CY7C68001)進行數據傳輸時[9]，數據寬度為16 b，配置在EMIFB(外部存儲器接口B，16位數據寬度)的CE3空間，DSP地址映射范圍為0x6C000000~0x6FFFFFFF。

(2) DSP及PC機的數據傳輸實現[10]:采用USB 2.0的bulk批傳送方式傳輸數據。對于數據傳輸，考慮兩種方法，并結合算法特點進行了取舍。

方法一:在發送端將實數小數點后移，即放大發送前的數據，再將放大后的數據轉化為shortint，使其滿足USB 2.0接口的16 b的數據寬度，然后在接收端將得到的相應數據小數點前移，以恢復發送端的數據。該方案只能處理整數類型的數據，對于含有小數部分的實數只能保證正確接收到數據的整數部分，傳輸數據誤差大。由于ART算法的參量數據大多為含有小數部分的實數，故本文沒有采用這種方法。

方法二:直接傳輸浮點型數據，以確保算法參量數據的正確傳輸。一方面，由于USB接口芯片CY7C68001寄存器均以8 b數據寬度操作，本文方法是一種利用標準C語言的共用體結構(結構內的成員共同占有同一段內存單元)，共用體內成員為一個浮點型數據(32 b)和4個字符型數據(4×8=32 b)，每次針對字符型數據(8 b)進行讀寫操作。另一方面，由于USB 2.0接口的數據寬度是16 b，決定了浮點型數據不能一次傳輸，本文采用2次傳輸，在DSP設備端，根據接口芯片(CY7C68001)對OUT型端口(端口號2，端口號4)的讀操作進行相應順序的整合后，完成一次浮點數據的接收，另外采用同樣的方法，將結果數據以浮點類型傳送給PC主機端，顯示重建圖像。

2.2.2 DSP計算過程

(1) TI CCS環境下對算法計算時間測定;

(2) PC端發送投影矢量數據，并緩存在SDRAM;

(3) PC端發送投影系數矩陣的第一行，同樣緩存在SDRAM;

(4) DSP針對投影系數矩陣的第一行進行算法計算，并將計算結果緩存在SDRAM中;

(5) DSP計算完成后，PC端發送投影系數矩陣的第二行數據，并將其緩存在SDRAM中;

(6) DSP針對投影系數矩陣的第二行進行算法計算，并再次將計算結果緩存在SDRAM中;主機端根據TI CCS環境下測得的算法計算時間暫停發送數據;

(7) DSP計算完成后，PC端發送投影系數矩陣的第三行數據，過程一直循環，直到PC端將投影系數矩陣數據傳送完畢，即DSP完成一次迭代計算;

(8) 下次迭代計算時，DSP直接從SDRAM中取數進行算法計算。

DSP將計算的最終結果緩存在SDRAM，并通過USB 2.0接口傳送給PC機進行圖像顯示。DSP計算過程流程圖如圖3所示。

圖3 DSP計算過程流程圖

3 實驗結果

實驗以Shepp_Logan的頭部剖析圖(80×80)作為原始圖像，設定射線投影次數為90次，每次的射線為70束，即投影矢量數據總數是6 300(90×70);投影系數矩陣為R (6 300×6 400)。實驗使用USB 2.0協議進行算法數據的傳輸，在TI CCS 2.2編譯環境進行ART算法編程及運行時間測試，首先采用標準C語言編程，直接將算法計算中的數據定義為浮點型，故ART迭代計算公式以投影矢量的6 300個浮點數據以及投影系數矩陣的1行6 400個浮點數完成公式的單次計算。

作者使用CCS profile功能測試ART算法迭代公式(見式(3))計算1次所需時間。測試結果見圖4。

圖4 迭代公式計算1次時間

圖4中，Incl.Total表示DSP在統計工程中剖析代碼段消耗的所有時鐘周期數，TMS320C6416的主頻為1 GHz即周期為1 ns?？傻?，迭代公式的單次計算時間為8.768 6 ms，進行1次迭代需要6 300次同樣的計算，花費的時間約為55 s。如果多次迭代，算法計算的耗費時間將不能夠滿足實時性要求，其代碼需要被進一步的優化。通過對算法代碼進行分段時間測定發現，浮點數除法運算耗費時間最多，其次為浮點數乘法運算。采用的優化方法主要有[11]:采用TI IQmath庫進行重新編寫乘除計算;將數組以指針的形式替代表示;將多重for循環展開，優化后將算法迭代公式單次計算時間縮短為0.735 6 ms。另外，本文在Microsoft Visual C++6.0環境中使用標準C語言同樣編寫算法程序，并測得其算法耗費時間0.793 7 ms，如表1所示。

表1 定點DSP與PC處理器計算時間比較

定點DSP(TMS320C6416)PC AMD處理器(Athlon 3000+，1.81 GHz)

標準C語言乘除運算代碼優化后

8.768 6 ms0.735 6 ms0.793 7 ms

通過對表1的分析，被優化后定點DSP算法的代碼計算時間比PC AMD處理器進行同樣的計算減少0.058 1 ms，即減少了7.32%。

4 結 語

CT ART重建算法所需數據量大，重建時間長。本文通過USB 2.0協議傳輸算法數據，以及對該算法進行定點DSP (TMS320C6416)的片上算法計算，使其成為該算法在當前PC環境下實現的另外一種可選方案。實驗結果表明，選用定點DSP能有效提高ART算法的CT圖像重建速度，從而驗證了本文方案的可行性。當然，本文算法本身還不是最優的，DSP片上的算法計算還有進一步的優化空間，采用多DSP并行計算也會提高算法的計算速度，這些都是作者今后進一步研究的方向。

參考文獻

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

[2]張順利，張定華，李山，等.ART算法快速圖像重建研究[J].計算機工程與應用，2006，42(24):1-3.

[3]高欣.新型迭代圖像重建算法的理論研究與實現[D].杭州:浙江大學，2004.

[4]NIKOLIU Zoran. Implementation considerations for single-camera steering assistance systems on a fixed point DSP [J].IEEE Intelligent Vehicles Symposium. Eindhover: Eindhoven University of Technology， 2008: 697-702.

[5]YAN Lei， ZHAO Gang， LEE Choon-Young， et al.The Platform of image acquisition and processing system based on DSP and FPGA[J]. International Conference on Smart Manufacturing Application，2008，9(11): 470-473.

[6]MAGAR S S， CAUDEL E R， LEIGH A W. A microcomputer with digital signal processing capability[C]//Proc. ISSCC. Piscataway: IEEE， 1982:32-36.

[7]GAN W S，KUO S M.Teaching DSP software development: from design to fixed-point implementations[J].IEEE Transactions on Education， 2006， 49(1): 122-130.

[8]Texas Instruments. TMS320C6000 EMIF to USB interfacing using cypress EZ-USB SX2[M]. Texas， USA: Texas Instsuments， 2004.

[9]北京合眾達公司.SEED-DEC6416用戶指南(Rev.B)[M].北京:北京合眾達公司，2005.

[10]北京合眾達公司.SEED-DEC系列模板USB驅動程序架構與使用說明[M].北京:北京合眾達公司，2005.

[11]馬斌.基于DSP的JPEG編碼器的實現與優化[D].西安:西安電子科技大學，2009.