基于FPGA的快速DCT医学图像的压缩算法

【摘要】 介绍了一种适用于医学图像压缩的二维DCT快速算法的FPGA实现结构，采用行列分解法来实现该算法，首先把8×8的二维DCT变换分解为两个一维DCT变换，通过对变换的系数矩阵进行化简，使加法器和乘法器数量减少到最少，并采用了FPGA特有的并行流水线技术，明显节省了计算时间，提高了图像处理速度。

【关键词】 FPGA;二维DCT;医学图像压缩;行列分解;并行流水线

　　Abstract：A FPGA design for fast Discrete Cosine transform(DCT) implementation architecture used in medical image compression is presented. The architecture is achieved using row-column decomposition. First， the 8×8 two-dimensional DCT is decomposed into two one-dimensional DCT， simplifying the coefficient matrix， we can reduce the number of adder and multiplier to minimum， and using the parallel pipeline technology of FPGA， the computation time is greatly saved and image processing speed is improved.

　　Key words：FPGA;2D DCT;Medical image compression;Row-column decomposition;Parallel pipeline

　　1 引 言

　　DCT变换是视频压缩编解码器中很重要的一部分，被广泛应用于各种视频格式的编码算法中，例如：JPEG，MPEG1，MPEG2，H.264等[1]。DCT变换虽然不能直接降低数据量，但是它可以利用图像的统计特性使得各种降低数据量的方法更为有效地工作。它能把图像的能量集中到少数的几个数据上，在很大程度上消除数据间的冗余性和相关性。

　　一般来讲，医学图像都比较大，例如腹部横断面的CT图像(512×512的bmp图像)，对于此类图像，如果要在整体上进行DCT变换，将耗费大量的时间，因此我们需要将整幅图像切割成若干的子块，对子块进行DCT变换，在本研究中采用的是8×8pixels的字块。

　　二维DCT变换需要进行大量的运算，大量的乘法及加法运算严重影响了变换速度，为减少运算次数，缩短运算时间，人们作了不懈的努力，并提出了多种快速算法，但处理过程中数学运算仍然相当复杂。我们在前人工作的基础上，对用行列分解法实现DCT变换的方法进行了优化，并充分利用了FPGA器件的嵌入式乘法器及存储器资源，明显缩短了变换时间。

　　2 DCT算法原理

　　2.1 一维DCT算法原理

　　设{x(n)}表示N个有限的一维实数信号序列集合，n=0，1，2，…，N-1，则一维DCT定义为：

　　y(k)=2Nc(k)∑N-1n=0x(n)cos((2n+1)kπ2N)，

　　其中c(k)=12 k=0

　　1 k=1，2，…，N-1

　　2.2 二维DCT算法

　　对于一个N×N的图像，x(i，j)表示图像样值，其二维DCT变换公式如下：

　　y(u，v)=2Ne(u)e(v)∑N-1i=0∑N-1j=0x(i，j)

　　cos(2i+1)uπ2Ncos(2j+1)vπ2N，

　　其中

　　e(u)，e(v)=12 u，v=0

　　1 u，v≠0(u，v=0，1，2，…，N1)

　　生 物 医 学 工 程 研 究 第28卷第3期 王 宁，等：基于FPGA的快速DCT医学图像的压缩算法 二维DCT变换具有很高的复杂度，从上式我们可以看出，二维DCT可以分解为两个一维的DCT变换：

　　y(u，v)=2Ne(v)∑N-1j=0

　　2Ne(u)∑N-1i=0x(i，j)cos(2i+1)uπ2Ncos(2j+1)vπ2N，

　　首先对二维DCT进行一维列变换：

　　x(u，j)=2Ne(u)∑N-1i=0x(i，j)cos(2i+1)uπ2N

　　转置之后再次进行一维列变换：

　　y(u，v)=2Ne(v)∑N-1j=0x(u，j)cos(2i+1)vπ2N

　　至此，我们就完成了一次完整的二维DCT变换。

　　3 DCT变换在FPGA上的实现

　　3.1 一维DCT的实现

　　在图像压缩处理中，通常将图像分为若干的子块，继而对每一个子块进行DCT变换，在本方案中我们采用8×8pixels的子块。根据一维DCT变换的原理，可以得到一列8点DCT的计算结果，如下：

　　y0

　　y1

　　y2

　　y3

　　y4

　　y5

　　y6

　　y7=d〖〗ddddddd

　　aceg-g-e-c-a

　　bf-f-b-b-ffb

　　c-g-a-eeag-c

　　d-d-ddd-d-dd

　　e-agc-c-ga-e

　　f-bb-f-fb-bf

　　g-ec-aa-ce-gx0

　　x1

　　x2

　　x3

　　x4

　　x5

　　x6

　　x7，

　　其中a

　　b

　　c

　　d

　　e

　　f

　　g=12cosπ/16

　　cosπ/8

　　cos3π/16

　　cosπ/4

　　cos5π/16

　　cos3π/8

　　cos7π/16

　　要完成一次列向量的DCT变换需要64次乘法和56次加法，根据三角函数性质和矩阵系数对称性，上式可以分解为两个矩阵表达式：

　　y0

　　y2

　　y4

　　y6=dddd

　　bf-f-b

　　d-d-dd

　　f-bb-fx0+x7

　　x1+x6

　　x2+x5

　　x3+x4，

　　y1

　　y3

　　y5

　　y7=aceg

　　c-g-a-e

　　e-a-gc

　　g-ec-ax0-x7

　　x1-x6

　　x2-x5

　　x3-x4，

　　现在只需32个乘法器和32个加法器[2]，上面矩阵还可以继续进行分解：

　　y0

　　y4=dd

　　d-dx0+x7+x3+x4

　　x1+x6+x2+x5

　　y2

　　y6=bf

　　f-bx0+x7-(x3+x4)

　　x1+x6-(x2+x5)

　　y1

　　y7=ga

　　-ag(x3-x4)-d(x2-x5)+d(x1-x6)

　　d(x2-x5)+d(x1-x6)+(x0-x7)

　　y5

　　y3=ce

　　-ec(x3-x4)+d(x2-x5)-d(x1-x6)

　　-d(x2-x5)-d(x1-x6)+(x0-x7)

　　经过上述分解运算，乘法器和加法器的数量分别减少到了16次和26次，极大地降低了计算量。

　　根据上面分解后的矩阵，我们可以得到图1的一维DCT流水线结构：

　　图1 一维DCT流水线结构

　　Fig 1 One-dimensional pipeline structure

　　如图1所示，一维DCT变换采用了3级流水线的并行处理结构，当有8点数据输入后，首先经过串并转换电路将其转换为并行数据，之后利用FPGA丰富的乘法器资源和并行处理特性，实现8点实序列的并行输出，只需1个时钟周期就可以完成一次列向量的DCT变换，一次完整的一维DCT变换需要8个时钟周期。

　　图2(左)是一维DCT流水线结构的基本单元[3]，由4个乘法器和2个加法器组成，若采用图2(右)的等价单元结构，则所需乘法器的数量将进一步减少，共需13个乘法器和29个加法器。

　　图2 基本结构单元

　　Fig 2 The basic structural unit

　　3.2 二维DCT的实现

　　实现了一维DCT变换后，要实现二维DCT就很容易，只需加上用于数据存储和转置的RAM存储器即可，二维DCT变换的结构框图见图3[4]：

　　图3 二维DCT变换的结构框图

　　Fig 3 The block diagram of 2D-DCT transform

　　行列变换RAM可以形象地看成8×8的阵列，RAM的结构有8块8×8的双端口同步RAM组成[5]。每一块RAM块中存放1/8的8×8DCT系数。RAM中的前8个字对应矩阵中的第一列，第二个8个字对应第二列，以此类推(一列中的8个数据分别放入8块RAM中)。

　　采用这样的RAM结构可以实现8个数据的同时读写。首先对8×8矩阵的每一列进行一维DCT变换，当一列数据完成第一次DCT变换后，列向量中的8个数据同时写入RAM存储器，等待8列数据全部写入RAM并进行转置，之后再次对列向量进行一维DCT变换，因此，完成一个8×8点的二维DCT变换总共需要16个时钟周期，并且在第9个时钟沿有数据输出。

　　若对前后两个一维DCT分别使用各自的DCT流水线结构，则可以在一个子块进行第二次DCT变换的同时，对下一个子块进行第一次DCT的计算，从而节约时钟周期。例如：第一个子块耗费了16个时钟周期，这时第二个子块也应经完成了一次DCT变换，只需再过8个时钟周期即可完成，因此两个子块总共耗费24个时钟，而不是16+16=32个时钟周期。对于N个子块来说，总共耗费的时钟周期数为(8+8N)个。

　　4 仿真结果与分析

　　本算法是在Quartus II8.1环境下进行的设计和仿真(verilog语言描述)，采用的是Cyclone III系列的EP3C120F780I7芯片，共消耗了3 308个logic element(3%)、388个pins(73%)、1 993个logic register(2%)。如图4所示，iclk为时钟信号输入;d0～d7为8×8矩阵的数据输入，分别代表矩阵的1～8行数据，列向量数据同时输入;F0S～F7S表示完成每一列一维DCT变换后的数据输出，数据输入输出在同一时钟内完成;m0～m7表示完成第二次一维DCT变换后的数据输出，即完成一次完整的二维DCT变换的数据输出。从图中可以看出，从第9个时钟沿开始有数据输出，并且在第16个时钟结束时完成一次完整的二维DCT变换(忽略误差不记)，与理论结果相符，验证了算法的正确性。

　　5 结论

　　本研究结合FPGA具有丰富的乘法器和RAM存储器资源的优势，对基于行列分解法的二维DCT变换进行了优化。优化后的DCT蝶形图只需13个乘法器和29个加法器，极大地降低了计算量，再加上FPGA高速的并行处理结构和流水线结构，实现了对列(行)向量中8个数据的同时读写，完成一次二维DCT变换只需耗费16个时钟周期，并且可以实现两个子块DCT的同时进行，进一步减少计算时间。本算法耗时极少，达到了优化的目的，适用于高速实时的图像数据传输。

参考文献

　[1]Conzalez R C，Woods R E著.阮秋琦，阮宇智译.数字图像处理[M].北京:电子工业出版社，2003.

　　[2]Patino A M， PerioM M，Bllester F，et al.2D-DCT on FPGA bvpolynomial transformation in two-dimensions[A].Circuits and System[C].Proceedings of the 2004 Internatinal symposium，2004.3:23-36.

　　[3]陈普跃，赵新璧，陈斌.二维DCT快速算法及FPGA实现[J].电子质量，2008，2:5-7.

　　[4]孙阳，余锋.基于脉动阵列的二维DCT算法及其VLSI设计[J].微电子技术，2003，31(5):21-26.

　　[5]山洪刚，郑南宁，杨晓衡.实时视频编码的二维DCT/IDCT的实现[J].数字电视与数字视频，2002，12:4-6.