基于定点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 断层截面与投影的几何布置

该图像重建算法可用下列线性方程组表示:

r11f1+r12f2+…+r1NfN=p1

r21f1+r22f2+…+r2NfN=p2

rM1f1+rM2f2+…+rMNfN=pM(1)

式中:N=n×n为像素总数;M为投影射线总数;rij为投影系数,即加权因子,反映了第j个像素对i根射线积分的贡献;fj表示第j(j=1,2,…,N)个方格(像素)内的常量值。┦(1)写成i,j=1,2,…,NЬ卣笮问轿:

RF=P(2)

式中:R=(rij)M×N为投影系数矩阵;F=(f1,f2,…,fN)T为图像矢量;P=(p1,p2,…,pM)T为投影矢量。对于式(1)的求解一般不直接通过矩阵求逆的方法求解,而采用迭代算法。迭代公式如式(3)所示:

Xi=Xi-1-αRiXi-1-piRi•RiRi(3)

式中:i表示射线号;X表示需要重建的图像向量;Ri表示投影系数矩阵的i行;α代表松弛系数,取值范围为(0,2)。

2 CT图像重建ART算法的DSP实现

2.1 ART算法数据的定点DSP表示

ART迭代算法中含有大量的实数数据,包含有整数部分和小数部分。结合定点DSP和ART算法的特点,本文采用了两种数据表示方法:定标法和直接表示法。

定标法,通过小数点的不同位置来表示不同的实数,小数点前为整数部分,其后为小数部分。本文采用32位定点DSP,所能表示的整数数值范围为[-232-1, 232-1-1]。定点DSP在表示算法数据时,需要考虑计算数据的数值范围[7]。针对定标法,考虑了该算法的如下参量值[1]:

(1) 投影系数矩阵元素,见图1。它表示所要重建图像的像素对射线积分的贡献,值等于射线穿过像素所占部分的面积,其数值取值范围为[0,1];

(2) 投影矢量pi,pi=∑jrijxj。表示i根射线经过的所有像素总的投影之和。其中,rij为投影系数矩阵元素;xj为j号像素的像素值。为了说明情况,设定rij数值均为1,将xj设为灰度值255,pi为232-1-1,求得j的总数约为1.7×107(图像数据矩阵约4 000×4 000)。换言之,当投影矢量元素为定标法所能表示32位定点数的最大值时,图像数据矩阵约为4 000×4 000,每个像素的灰度值为最大值255;投影系数矩阵元素rij均为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.