快速二维（共4篇）
快速二维 篇1
0 引言
自1974年N.Ahmed等人提出了离散余弦变换(Discrete Cosine Transform,DCT)算法后,很多学者对DCT算法的实现进行了研究。陈禾等人[1]总结了前人的研究,将DCT算法的结构归纳为5类,其中基于乘法器的DCT结构一直是很多学者研究的对象,具有代表性的有W.H.Chen、LLM[2]、AAN[3]等算法。LLM算法将1-D DCT所需的乘法器个数减小到理论最小值11个。
在编解码器应用方面,Xvid的MPEG-4编解码器使用LLM算法的12个乘法器形式,而很多JPEG编码器用AAN算法。李晗等人[4]利用加法和移位器代替乘法器,在提高频率上取得了比较好的效果,但是即使将其运用到LLM算法中,2-D DCT仍需要使用184个加法器。
本研究通过引入整数变换矩阵和增加缩放模块,设计面积小、速度快、频率高的且可复用多个标准的2-D DCT。
1 离散余弦变换DCT原理
离散余弦变换DCT,是指对一帧图像以块为单位(通常是8×8或4×4像素块),通过2维DCT变换将图像从空间域变换到频域的过程,目的在于去除像素间的相关性。
二维N×N图像块的DCT变换的定义[5]如下:
2 二维DCT的VLSI结构优化设计
2.1 算法优化
整数矩阵的引入,一方面有效地减少计算量,另一方面由于Ef的运算被放到量化端和量化同时考虑,其矩阵系数由量化步长QP的大小决定。所以H.264标准对于DCT的处理和别的标准不同,很难用于其他标准中,即不具有普遍性。笔者考虑将Cf矩阵引入,得到Ef矩阵:
其中,a,b,c,d,e,f都是浮点数,由于在VISI中运算不方便,将其转化为定点数,经过扩展放大216得到a=8 192,b=7 740,c=10 486,d=7 346,e=9 777,f=13 159。
如果直接乘以缩放矩阵式(3)的64个元素需要64个乘法,耗费时钟多,面积大。本研究考虑引入流水线技术,具体有以下两种方案:
(1)通过在乘法器前面增加选择器可以压缩到8个乘法器。例如,如果2-D DCT按列输出,则每一列相应的元素都乘以上面矩阵中对应列的相应元素,由于对称性,可以在让每个乘法器在3个数之间选择,送出乘数,进行相乘,如第1个乘法器在a,b,c之间切换,如果当前输入列为第1列或者第4列,则输出到第1个乘法器的一个乘数为a,另一个乘数为相应列的第
写成矩阵形式为F=AfAT,其中A矩阵为浮点数矩阵,在H.264 High Profile标准中,将A阵进行整数化,矩阵形式等效为:Y=(CfXCfT)Ef,其中Cf矩阵[6]为:
1个数据,如果是第2、4、6、8列,则选择相应的乘数为b,如果是第3、7列则选择乘数为c。同理,其他6个乘法器类似。
(2)如果对频率要求比较高,考虑到乘法器的延迟比较大,将乘法器用加法器和移位寄存器替代。根据3种乘数并行计算出乘法结果。例如,由于第2维DCT送出的某一列的第1个像素结果,有可能进行的乘法是“×a”、“×b”、“×c”。为了取代这个乘法器,先用加法器和移位器取代这些乘数,如:b=213-29+26-22。这样可以通过加法和移位来计算出当前数据经过3种乘法之后的结果。通过计数当前列数即可选出最后的结果。由于每个系数最多有3个加法器即可实现,因此,方案1需要的一级乘法器延迟可以变成2级加法的延迟,提高了工作频率。取代8个乘法器所需加法器数量如表1所示。
由表1可知,每组改成2级加法提高工作频率后,用66个加法器取代8个乘法器。针对不同的应用需求(频率、面积)可以考虑不同的方案。
2.2 算法验证
对以上算法进行C语言建模,并将其移植到Xvid的MPEG-4的开源代码[7]编码器的DCT中。在不改变其他模块(包括IDCT)的前提下,通过输入不同的YUV原始视频序列(格式为4:2:0),先用标准代码进行编码,再用改进后的代码进行编码,分别通过Xvid的MPEG-4解码器进行解码。最后分别将解码得到的YUV文件和原始YUV文件,在亮度Y、色度Cb和Cr分量上进行PSNR值对比。统计每个序列300帧中每种成分的均值和所有序列在Y、Cb、Cr上的均值。本研究所用YUV原始视频序列均为标准测试序列,图像大小为CIF,比特率设为512 Kbps时统计所得结果如表2、表3所示。
YUV序列1~9视频序列分别是akiyo,coastguard,container,mother_daughter,stefan,tempete,flower,vectra,footbal,帧数为300。表中Y0,Cb0,Cr0分别代表Xvid标准算法编码所得视频序列的亮度Y和色度Cb,Cr分量的PSNR值,而Y1,Cb1,Cr1,Y2,Cb2,Cr2分别为方案1和方案2在各分量的PSNR值。
两种方案相对于Xvid标准算法信噪比变化如表4所示,数值如果为正,表明PSNR值有所改善,如果为负,表明PSNR值有所损失。由表4可以看到本研究的两种方案图像质量损失不超过0.015 d B,而且方案2显示出本研究算法在色度分量上分别改善了0.004 d B和0.015 d B。经过大量统计,结果表明本研究提出的算法保证了图像质量,精度满足要求。
3 二维DCT的硬件实现
3.1 整体结构
本研究采用行列分解法进行二维DCT的设计,其硬件实现框图如图1所示。
整个2-D DCT采用流水线设计,每个时钟输入为8个数据,输出也是8个数据,中间除了流水深度需要的时钟外,不需要任何多余的等待,速度快。同时可以根据不同的应用需求以及对频率和速度要求更改流水节拍,实现速度和频率的最优化。
3.2 一维DCT的实现
这里一维DCT要实现的是式(2)中的Cf矩阵,其硬件实现如图2所示。硬件实现的时候可以设计成4级流水或2级流水。
8×8块的一行数据经过4级加法器之后输出。如果对频率要求比较高,考虑用4级流水,图2的1、2、3、4分别为1D-DCT的对应4级加法,在第3级的加法中上半部分的数据需要进行寄存以保证数据的同步。如果需要时钟周期少,则用2级流水,图中5和6分别对应2级加法。图中的“”为右移操作。
为了保证数据不溢出,对于第1级一维DCT输入的数据精度为9 bit,输出数据精度为13 bit。对于第2级一维DCT输入的数据精度为13 bit,输出为17 bit。
3.3 转置的实现
为了保证速度,流水作业不能被打断,如果用单口RAM乒乓操作则实际所需要的RAM数为16个深度为1、宽度为104 bit的RAM,所以无论从RAM的产生还是面积的角度来考虑都不现实。A.Aggoun等人[8]使用2组三角矩阵堆实现了矩阵转置的功能,通过控制8个8选1选择器来连接,这种结构需要72个寄存器,流水的深度为9 clk。李晗等人实现4×4块只用了16个寄存器,节省了面积和时钟。
本研究使用64个寄存器的寄存器组进行转置的硬件实现方案,如图3所示。图中每个方框(如r00)都是由一个选择器和一个寄存器组成。举例说明,如果输入的数据是按照8×8的列顺序输入,依次是x0,x1,…,x7。在第1个时钟将x0~x7依次存入寄存器r70,r60,…,r00中,第2个时钟,第2列的寄存器将第一列r70~r00的数据依次存入r71~r01中,而r70~r00的数据还是从输入端取数据存入。如此执行下去,直到第8个时钟,整个8×8块的数据已经存入寄存器组中,此时寄存器r77,r76…r70的值已经有效,依次从中取出数据,经过选择器按顺序排好输出,此时输出的是上个矩阵的第1列。在第9个时钟,r77~r70从r67~r60中取数据并且更新寄存,而r67~r60从r57~r50中取数据并且寄存,其他行的数据类似操作,唯独r07~r00分别从输入端取数据x0~x7并且寄存。第16个时钟时第2个8×8块的数据已经寄存完毕,可以依次从r77~r07的寄存器中取数据。如此流水下去,便完成了行列转换,流水深度为8 clk,中间不需要等待,每个时钟每个寄存器都处于工作状态,提高了寄存器的利用率,并且节省了面积,提高了速度。
3.4 缩放模块的实现
(1)方案1:使用乘法器。
采用方案1的缩放模块的硬件实现如图4所示,使用8个乘法器来实现。通过8个3选1选择器分别选择不同的乘法器系数,control通过计数确定当前输出的数据是8×8块中的第几列。如果是第1、5列,则8个选择器依次选择8 192,7 740,10 486,7 740,8 192,7 740,10 486,7 740。如果是第2、4、6、8列,即依次选7 740,7 346,9 777,7 346,7 740,73 46,9 777,7 346。如果是第2、4、6、8列,则依次选择10 486,9 777,13 159,9 777,10 486,9 777,13 159,9 777。本方案的延时主要是1级乘法器。
(2)方案2:使用加法器和移位器代替乘法器。
采用两级的加法器和移位器代替乘法器,提高电路的工作频率,如图5所示。第1个数据经过并行的3路输入,每一路都由移位器和加减法器组成,通过control模块对输入的数据所在8×8块的列数进行选择,以面积换速度,通过3路并行工作,最后控制选择输出。本方案的延时主要在2级加法器。
两种方案都需进行缩放,由于这些系数都经过放大216倍来保证中间计算值的精度,最后DCT结果需要经过右移16 bit,需要注意数据的截断。
3.5 性能对比
本研究用Verilog HDL语言描述电路,分别用VCS和DC进行了仿真和综合,下面将本研究与其他文献所设计的2-D DCT进行性能对比。2-D DCT算法使用乘法器和加法器个数对比如表5所示。
由表5可见,本研究提出的2种方案计算复杂度比其他文献都小。其中,将文献[4]应用于LLM的24个乘法器和64个加法器2D-DCT方案需要184个加法器。若应用于理论上最低乘法器(11个)的结构时,也需要168个加法器。AAN算法中1D-DCT蝶形中用了5个乘法器,但每次一维DCT后要乘以8个缩放系数需要额外的6个乘法器,或二维DCT后乘以一个缩放矩阵,对于除JPEG外的其他标准不适用。
ASIC设计实现2-D DCT性能对比如表6所示[9,10],本研究所设计的2-D DCT完成单个8×8像素块的周期被降到20个时钟频率能达到510 MHz。在90 nm CMOS工艺下,DC综合最高可以达到840 MHz。若数据保持流水送入,则只需要经过3.09μs可完成一帧高清1 080 pixels图像的数据处理,满足工程应用需求。
FPGA原型验证采用的器件为Xilinx Virtex-5XC5VLX220-1 FF1760,综合工具采用Xilinx ISE。对于2-D DCT,计算周期为20 clk时,方案1综合频率为180 MHz,方案2为200 MHz,周期为25 clk时方案2为290 MHz。
4 结束语
本研究旨在通过设计适用于MPEG-4,H.263,JPEG,H.264等编码器芯片的快速8×8 DCT的IP核,研究并设计了适合于VLSI实现的快速2-D DCT结构,利用行列分解法,通过引入整数变换矩阵、增加缩放模块实现,经算法验证满足精度要求。转置的实现方案提高了寄存器的利用率,并且减小了流水深度。提供了2种实现方案实现缩放模块,可以根据不同的应用需求在面积、速度、频率之间权衡。采用流水线设计,并减小了流水深度,提高了速度和频率。具有灵活性和实用性。
用Verilog HDL进行硬件设计,并用VCS进行了仿真,在90 nm CMOS工艺下电路综合频率达到840 MHz,FPGA模拟验证频率达到290 MHz,完全满足设计要求。与其他文献2-D DCT算法结构相比,本研究提出的结构具有面积小、速度快、频率高的特点。不仅可以很好地应用于MPEG-4和H.263和JPEG等标准,还将H.264 High Profile中的8×8整数DCT进行了复用,是一个基于多标准的2-D DCT的IP核,可以很好地应用于支持多标准的编码器中。
摘要：针对多标准图像视频编码器中二维离散余弦变换的复用问题,设计了适用于MPEG-4,H.263,JPEG,H.264 High Profile编码器芯片的快速8×8 2-D DCT IP核,并完成了RTL设计、仿真和FPGA原型验证。通过引入H.264 High Profile 8×8整数变换矩阵和增加缩放模块,完成了多个标准中2-D DCT的复用。经Xvid MPEG-4编解码器验证,满足精度要求。设计采用流水线技术,并优化了速度和频率,在90 nm CMOS工艺下频率达到840 MHz。研究结果表明,该技术能很好地应用于多标准编码器中。
关键词：离散余弦变换,多标准图像与视频编码,JPEG,MPEG-4,H.264
参考文献
[1]陈禾,毛志刚,叶以正.DCT快速算法及其VLSI实现[J].信号处理,1998,14(A12):62-70.
[2]LOEFFLER C,LIGHTENBERG A,MOSCHYTZ G S.PracticalFast 1-D DCT Algorithms with 11 Multiplications[C]//Int ConfAcoustics,Speech and Signal Processing.Glasgow:[s.n.],1989:988-991.
[3]ARAI Y,AGUI T,NAKAJIMA M.A fast DCT-SQ schemefor images[J].Transactions of Institute of Electronics,Information and Communication Engineers,1988,E71(11):1095-1097.
[4]李晗,孙义和,向采兰.二维离散余弦变换及其逆变换的VLSI实现[J].微电子学,2008,38(3):326-329.
[5]毕厚杰.新一代视频压缩编码标准—H.264/AVC[M].北京:人民邮电出版社,2005:112-114.
[6]GORDON S,MARPE D,WIEGAND T.Simplified use of8×8 transforms[S].JVT Doc,JVT-K028,Munich,2004.
[7]XVID MPEG-4 VIDEO CODEC[DB/OL].[日期不详].www.xvid.org.
[8]AGGOUN A,JALLOH I.Two-dimensional DCT/IDCT ar-chitecture[J].IEEE Proc.Computer Digital Technolo-gy,2003,150(1):2-10.
[9]KURODA T.A 0.9 V,150 MHz,10 mV,4 mm2,2-D dis-crete cosine transform core processor with variable threshold-voltage(VT)scheme[J].IEEE J.SSC,1996,31(11):1770-1779.
[10]KRISHNAN R,GANGWAL O P,EIJNDHOVEN J,et al.Design of a 2D DCT/IDCT Application Specific VLIW Pro-cessor Supporting Scaled and Sub-sampled Blocks[C].16th International Conference on VLSI,2003.
快速二维 篇2
关键词：三维仿真,碰撞检测,层次模型
0 引言
碰撞检测问题按运动物体所处的空间可分为二维平面碰撞检测和三维空间碰撞检测。由于平面物体的构造都可用多边形来表示, 故其检测算法相对要简单一些;而三维物体的构造比较复杂, 所以, 其碰撞检测算法也比较困难。而由于现实的应用, 需要对真实而且复杂的现实世界实现高质量的计算机模拟。因此, 开发高效的三维空间碰撞检测技术具有重要的现实意义。
而提高碰撞检测算法的效率, 不仅仅取决于基本干涉检验算法的效率, 也与基本检测算法使用的次数有很大关系。因此, 任何对于这两方面的改进, 都能使整个碰撞检测算法的效率提高。而本文就是运用平面碰撞检测算法结合动态八叉树的三维空间碰撞检测算法, 改进了基本检测算法, 并减少其使用次数, 以此获得一种新的三维空间的快速碰撞检测算法。
1 二维碰撞检测算法结合动态八叉树算法
1.1 平面碰撞问题
近10几年来, 许多专家学者对平面碰撞问题进行了深入的研究, 并取得一些很好的结果, 提出了许多算法。而由于要研究的碰撞对象处于同一平面内, 因此, 可得到有力和巧妙的技巧, 这使得平面碰撞问题已得到很深入的研究, 并提出了很多种最优算法, 这些成熟的理论很好地解决了平面碰撞问题, 本文在这方面就不多作论述。
1.2 传统八叉树算法及其缺点
以层次模型为基础的八叉树干涉检验算法, 是一个空间非均匀网格剖分算法。它是一种表示形体的分解方法, 在复杂三维体空间表达、求解体空间中的面模型和体模型之间的位置关系等领域得到了广泛的应用。八叉树是立方单元体的一种分层表示方法, 是一种空间分割的层次数据结构, 这种表示把三维空间 (通常为正立方体空间) 递归地分成8个单元或节点。目前在八叉树结构的研究上虽然有一些成果, 但所应用的八叉树方法大多是静态表示方法。静态八叉树在碰撞检测和干涉校验等方面的效果很显著, 但其缺点也很明显:在空间方面, 主要是占用的存储空间多, 只能近似地表示形体, 不易获取形体的边界信息, 不能实时修改等;在时间及计算复杂度方面, 由于场景中三角面片数量非常庞大, 因此, 对每个三角面片进行碰撞检测, 使得计算量巨大, 而且很大部分计算的并非将要发生碰撞的关键三角面片, 使得计算的效率低下。
1.3 二维碰撞检测算法结合动态八叉树算法
1.3.1 动态八叉树的立方体表示法
动态八叉树的基本组成是基元立方体, 采用立方体的中心坐标位置和立方体的边长来描述立方体的几何信息, 如图1所示。由自定义类Cube (x, y, z, a, k) 表示基元立方体, 其中x, y, z是立方体的中心点坐标, a是立方体的边长, k是基元裂变控制按钮, 则立方体的8个顶点坐标分别为:A (x-a/2, y+a/2, z+a/2) , B (x-a/2, y+a/2, z-a/2) , C (x+a/2, y+a/2, z-a/2) , D (x+a/2, y+a/2, z+a/2) , E (x-a/2, y-a/2, z+a/2) , F (x-a/2, y-a/2, z-a/2) , G (x+a/2, y-a/2, z-a/2) , H (x+a/2, y-a/2, z+a/2) 。
有了基元立方体之后, 以它为基体, 采用自相似分解规则, 对立方体进行多级层次分解, 而实现动态八叉树碰撞检测。基元立方体产生第一次分解后, 其8个子单元的中心点坐标为:O1 (x-a/4, y+a/4, z+a/4) , O2 (x-a/4, y+a/4, z-a/4) , O3 (x+a/4, y+a/4, z-a/4) , O4 (x+a/4, y+a/4, z+a/4) , O5 (x-a/4, y-a/4, z+a/4) , O6 (x-a/4, y-a/4, z-a/4) , O7 (x+a/4, y-a/4, z-a/4) , O8 (x+a/4, y-a/4, z+a/4) 。8个子单元的边长为a/2, k是每一单元的分解控制按钮, 其初始值为false, 当该单元检测到与碰撞检测对象发生碰撞且未达到精度要求、也不是完全位于检测对象内部或外部时设为true。图2是多层次几何结构实例, 其中被填充的单元表示达到碰撞检测要求的子立方体, 该单元的k值为false, 即不再对其进行分解。
1.3.2 动态八叉树的生成方法
当基元立方体或它的子立方体检测到与检测对象发生碰撞时, 该立方体进行自相似分解。为了满足虚拟环境基本实时性要求, 采用立方体模拟检测对象。初始状态下, 基元立方体表达为Cube (x, y, z, a, k) , k=false。根据检测对象的运动范围选择参加碰撞检测的基元立方体, 与检测对象进行碰撞检测, 构造动态八叉树。动态八叉树的构造方法如下:
(1) 进行碰撞检测。选取在检测对象运动范围内的基元立方体, 根据检测对象与该立方体的碰撞情况, 如果没有发生碰撞, k值不变;如果发生碰撞, 设置k=true。
(2) 当k=false时, 该基元立方体保持静态, 结束本次循环。
(3) 当k=true时, 将其均匀分成8个子立方体。在对每个子立方体进行碰撞检测之前, 要判断它与检测对象的位置关系:完全位于对象的外部或内部, 还是其他情况。若此子立方体完全位于检测对象的外部或内部, 停止对它的检测;否则检测它是否与对象发生碰撞。检测到碰撞后, 判断它的边长是否小于精度要求, 小于精度要求, 则发出碰撞报告, 否则将其继续分解。
(4) 对每一个子立方体重复以上步骤, 直至达到精度要求为止。
1.3.3 算法基本思想
二维空间的碰撞检测算法已经很完善, 而三维空间的八叉树碰撞检测算法却存在占有存储空间巨大和运算量大的缺点, 于是, 本文提出了一种将这两种算法结合在一起的新型算法。算法基本思想是, 将三维空间中的碰撞检测对象影射到二维平面上, 然后运用成熟完善的二维平面碰撞检测算法来进行碰撞检测, 在平面上得出发生碰撞的结果后, 再进行三维空间的动态八叉树碰撞检测计算。二维平面碰撞检测算法结合动态八叉树算法的基本思想如图3所示。
(1) 将待检测对象分别影射到非平行的几个平面上, 分别生成场景的平面影射图, 并在不同的平面内部, 使用平面碰撞检测算法对检测对象进行碰撞检测, 当几个平面内任何一个平面上未检测到碰撞, 则此次循环结束, 返回重新执行 (1) 。只有当所有平面上都检测到碰撞时系统才发出碰撞报告, 并执行 (2) 。
(2) 初始化待检测对象。为了用多个基元立方体的组合来代替检测对象进行碰撞检测, 并且尽量避免出现半个立方体, 首先将检测对象最大的长宽高数值找出来并整数化, 然后取它们的最大公约数作为立方体的边长。因为基元立方体很多, 为了减少运算量, 根据检测对象的运动轨迹、范围来选取参加碰撞检测的立方体, 只检测在检测对象运动范围内的立方体, 将那些不在检测对象运动范围内的立方体排除。
(3) 立方体由八叉树动态表示。在立方体未检测到与检测对象发生碰撞之前, 被静态地表示为一个立方体, 当检测到与检测对象发生碰撞时, 它将被均分为8个子立方体, 然后检测每个子立方体是否与检测对象发生碰撞, 再将发生碰撞的子立方体8均分, 未发生碰撞的子立方体停止分解。当出现下面3种情况之一时, 停止这一子单元的递归, 不再对其继续进行分解: (1) 当该子立方体已经裂变得足够小, 即它的边长达到所要求的精度时; (2) 当该子立方体完全位于检测对象的外部时; (3) 当该子立方体完全位于检测对象的内部时。当待检测对象的外形不太复杂, 或比较有规律时, 根据子立方体8个顶点的坐标很容易判断出它与待检测对象的具体位置关系。如此循环下去直到所有递归全部停止为止。在碰撞检测过程中, 根据检测对象的运动状态及系统的精度要求来设置时间段, 在每一时间段检测一次碰撞情况, 对八叉树结构刷新一次, 使得八叉树结构总处于变化中, 从而实现动态的碰撞检测。由上述算法, 得出算法流程图, 如图4所示。
2 实验及结果分析
碰撞检测问题是确定不同的物体在空间是否占有相同区域的问题。同时碰撞检测问题不仅仅归结为一般的求交问题, 针对不同的应用对象, 碰撞检测问题涉及到检测方法的复杂性、检测算法的可靠性和效率等, 不同的检测算法具有不同的特点和面向不同的应用对象。而本文提出的二维平面碰撞检测结合动态八叉树算法, 其本身应用的对象是外表面不太复杂的对象, 而由于运用影射到平面, 用平面上的碰撞检测算法先预处理三维碰撞检测, 这就极大地减小了计算的复杂性, 而且八叉树是在平面上得出碰撞的结论之后才生成, 且范围限制在对象的运动范围之内, 这样就大大减少了基元立方体的数量, 彻底地降低了运算规模, 大大降低了空间占用率。另外, 在判断两个基元立方体是否发生碰撞上, 本算法只通过各顶点坐标间的关系来判断, 而不需要向普通八叉树算法那样大量计算三角面片的距离, 这样更是降低了算法的计算复杂性, 大大提高了算法效率。
算法针对外表不太复杂的对象, 而算法在复杂性、空间占有率及效率方面, 都优于传统八叉树算法, 而在算法可靠性即判断碰撞的准确性上略差于它, 因此, 检验本算法的可靠性成了至关重要的一点。而下面的算例中, 本算法将与传统的碰撞检测算法在可靠性方面做一番比较。这种传统的算法是一种基于体模型的碰撞检测算法, 它在默认的情况下, 在某个视角中, 用规则的几何形体来逼近或替代检测对象。由于它的可靠性受视角的影响很大, 因此, 在实验中, 本文选取了3个不同视角分别进行测试。实验中一个检测对象为立方体, 边长固定为10mm, 另一检测对象的边长由200mm减小到5mm, 而动态八叉树的精度设定为5mm, 即子立方体分解到边长为5mm时会停止分解。实验结果如表1。表中a表示立方体边长, p表示误差, 即算法认为发生碰撞时两个立方体之间的实际距离, p值后面括号内的数字为相对误差, 即p/a。表中结果表明, 二维平面碰撞检测结合动态八叉树算法的可靠性要优于这种传统的碰撞检测算法。
3 结束语
八叉树干涉检验算法, 是一种表示形体的分解方法, 在复杂三维体空间表达、求解体空间中的面模型和体模型之间的位置关系等领域得到了广泛的应用。而目前在八叉树结构的研究成果中, 所应用的八叉树方法大多是静态表示方法。静态八叉树在碰撞检测和干涉校验等方面的效果很显著, 但也存在着占用存储空间巨大、计算复杂等缺点。本文提出的结合二维平面碰撞算法与动态八叉树的新型算法, 即在真正发生碰撞前, 用相对简单的二维平面检测算法代替复杂的三维空间检测算法, 并且在最终的碰撞检测时, 摒弃传统的需要大量计算三角面片距离的八叉树干涉检验算法, 而用简单的空间坐标间的相互关系来判断是否发生碰撞。新算法在空间占用, 计算规模、复杂度以及算法效率上, 比传统八叉树法有了很大提高。但是, 由新算法本身所决定的是, 算法应用的对象是外表面不复杂的对象, 而且在碰撞的判断准确性上要差于传统八叉树法。笔者对解决这两方面缺点的预期方法是, 在新算法执行到运用动态八叉树进行碰撞检测时, 方法回归传统, 即用计算三角面片距离的方法代替用坐标判断的方法, 当然, 这样改动后的算法效果还有待进一步研究。
注: (1) 传统算法部分数据引用参考文献[7]
参考文献
[1]王志强, 洪嘉振, 杨辉.碰撞检测问题研究综述[J].Journal Of Soft-ware, 1999 (10) .
[2]Tetsuya U, Toshiaki O, Mario T.Collision detection in motion simu-lation[J].Computer&Graphics, 1983 (2) .
[3]Chin F, Wang C A.Optimal algorithms for the intersection and the minimum distance problems between planar polygons[J].IEEE Transactions on Computers, 1983 (12) .
[4]David Baraff.Interactive simulation of solid rigid bodies[J].IEEE Computer Graphics&Applications, 1995 (5) .
[5]Ahuja N, Nash C.Octree representations of moving objects.Computer Vision, Graphics and Image Processing, 1984 (2) .
[6]吴明华, 余永翔, 周济.采用空间分割技术的八叉树干涉检验算法[J].计算机学报, 1997 (9) .
[7]齐晓松, 胡青泥, 刘晶.基于多视角的动态八叉树碰撞检测算法[J].东华大学学报, 2006 (5) .
快速二维 篇3
矩量法(MOM)被看作是数值"精确解"，被广泛应用于电磁散射问题的分析中，但每次计算只能得到一个频率点的电流分布和雷达散射截面(RCS)，而在日渐发展的遥感技术和雷达目标识别中，需要随机分布目标的宽带RCS以产生一维距离像。与其它频域方法一样，为了获得目标的宽带RCS，应用MOM就必须在频带内的每个频率间隔点上逐点计算，当目标的RCS随频率变化剧烈时，必须以很小的频率间隔计算才能得到精确的频率响应，这就意味着在整个频带内矩阵方程求解次数的增加。为了克服这个缺点，E.H.Newman和G.J.Burke分别通过内插阻抗矩阵[1]和使用基于模型参数估计[2]的方法获得了宽带数据。一种类似的技术--渐近波形估计(AWE)技术被提出,并首先用于超大规模集成电路的时域分析[3]。近年来，AWE技术被逐渐应用到电磁问题的分析中[4,5,6,7,8]。
本文将AWE技术应用到MOM中，计算了二维随机分布的理想导体柱的宽带RCS。数值计算表明，本方法与传统MOM的逐点计算结果吻合良好，而计算效率却显著提高。
1 渐近波形估计技术
对于任意形状二维理想导体柱的电磁散射问题，其电场积分方程(EFIE)为：
式中为导体柱截面的周线，为导体表面上的电流密度，为入射电场，为第二类零阶Hankel函数，和分别表示场点和源点的位置矢量，k、分别为自由空间波数和波阻抗。
将理想导体柱的周界划分成N段，选择脉冲函数作为基函数，函数作为检验函数，应用MOM可将EFIE化成矩阵方程Z(k)I(k)=V(k)(2)
式中Z为阻抗矩阵，V为激励向量，它们的元素表达式分别为：
为入射平面波与轴间的夹角。
应用矩量法求解式(2)，只能得到一个频率点的导体表面电流。为了得到给定频带内的导体表面电流分布，就必须以一定的频率间隔重复求解式(2)。AWE技术通过将展开成关于的泰勒级数得到频带内的导体表面电流分布，即
展开系数的表达式为:
式中Z(i)表示Z(k)的i阶导数，V(n)表示V(k)的n阶导数。
为了扩大泰勒级数的收敛半径，可通过Padé逼近将展开成有理函数,即
确定，将其代入式(8)，便得到给定频带内任意频率点的表面电流分布，进而计算出散射场和宽带RCS。显然，在AWE方法中，只需进行一次矩阵求逆运算，便可得到给定频带内电流密度的分布，这正是AWE方法提高计算效率的原因所在。
2 数值计算与比较
为了验证应用AWE方法分析二维随机分布导体目标宽带电磁散射特性的有效性，本文对如下两种情况的宽带RCS进行了计算。
(1)36个正方形理想导体柱随机分布在边长为0.0103米的正方形平面内，每个正方形理想导体柱的边长为λ0/20(λ0为中心频率对应的波长)，根据蒙特-卡罗方法[9]生成的二维理想导体柱分布状态如图1所示，选择入射波为TM波，E zi=e-jkx，中心频率0f=3 5GHz,Padé逼近中的L=4,M=3,20～50GHz范围内的RCS频率响应如图2所示。
(2)16个正方形理想导体柱随机分布在边长为0.0102米的正方形平面内，每个正方形理想导体柱的边长为λ0/10，分布状态如图3所示，其他参数同(1)，20～50GHz范围内的RCS频率响应如图4所示。
图3 N=16,a=0/1 0,l=0.0 10 2 m图4 RCS与频率间的关系(参见右栏)
在上述两种情况下，每个正方形理想导体柱离散单元尺寸选择为λ0/40，从以上两个算例的RCS频率响应曲线来看，用Padé逼近的AWE技术得到的结果与MOM逐点计算的结果有较好的一致性，而计算效率则显著提高，表1为上述两种情况所用的计算时间。
表1计算时间及与传统矩量法的比较*(参见右栏)
3 结论
本文给出了一种分析二维随机分布导体目标宽带电磁散射特性的方法，文中对两种情况下的随机分布理想导体柱的宽带RCS进行了计算，结果表明，本文的计算结果与矩量法逐点计算的结果吻合良好，而计算效率显著提高。对于任意给定的频带，要想获得精确的宽带RCS，需要使用多个频率展开点的AWE技术。
参考文献
[1]Newman E H.Generation of wide-band data from the method of moments by interpolating the impedance matrix[J].IEEE Trans.Antennas Propagation,1988,36(12):1820-1824.
[2]Burke G J,Miller E K,Chakrabarti S,et al.Using model-based parameter estimation to increase the efficiency of computing electromagnetic transfer functions[J].IEEE Trans.Magn.,1989,25(4):2807-2809.
[3]Pillage L T,Rohrer R A.Asymptotic waveform evaluation for timing analysis[J].IEEE Trans.Computer-Aided Design,1990,9(4):352-366.
[4]Polstyanko S V,Dyczij-Edlinger R,Lee J F.Fast frequency sweep technique for the efficient analysis of dielectric waveguides[J].IEEE Trans.Microwave Theory Tech.,1997,45(7):1118-1126.
[5]Li M,Zhang Q J,Nakhla M.Finite difference solution of EM fields by asymptotic waveform techniques[J].IEE Proc.Part H:Microwaves,antennas and propagation,1996,143(6):512-520.
[6]Erdemli Y E,Gong J,Reddy C J,et al.Fast RCS pattern fill using AWE technique[J].IEEE Trans.Antennas and Propagation,1998,46(11):1752-1753.
[7]孙玉发,徐善驾.渐近波形估计技术在三维电磁散射问题快速分析中的应用[J].电子学报,2002(60)794-796.
[8]施长海,孙玉发.二维电大导体目标宽带雷达散射截面的快速计算[J].电波科学学报,2004(6):325-328.
快速二维 篇4
目前, 超声成像技术无论是在医学成像还是工业探伤检测等诸多领域中均扮演了举足轻重的角色。特别在医学领域中, 当下二维超声成像具有图像显示直观, 声波穿透性强, 相对核磁共振成像和X-射线成像等更安全等优点。但是当前商业化的二维超声系统都是使用简化近似方法来描述超声场的, 这直接导致分辨率远低于理论可达水平。其主要原因是:
(I) 超声成像绝大多数算法是基于回波技术的, 即基于检测目标体的反射波实现图像重建;
(II) 算法的限制导致所成图像是定性的, 即重建的图像大多仅仅反映出检测目标体的空间位置 (或者边缘轮廓) 。
相对这些问题和不足, 超声CT的成像过程是基于检测目标体外围的散射场数据, 这些数据携含了更为全面的物体内部空间结构等方面的信息;因而超声CT成像的图像完全可以定量地真实地还原出包括在物体内的声速传播以及衰减系数等声学特性[1]-[3]。正是这种优势导致超声CT成像称为声学界一个研究的热点问题, 这也正是本文研究的原因;通过本文下述的理论推导和探究, 最终提出一个切实可行的快速高效的正演成像算法。
2、超声CT成像的积分方程
在构造算法的过程中, 我们基于声波在传播过程中的波动方程作为研究的出发点。基于波动方程的推导和演变, 我们可以根据研究的需要展开深层次的研究。基于参考文献[4]等的论述, 我们可知, 在二维空间下, 声波 (包括超声波, 电磁波等) 在穿越非均匀媒介时, 均会产生相应的“压力场”。该压力场又可以分解为两部分之和的形式。即为:
上式 (1) 中, P (r) 表示总场 (Total field) ;Pinc (r) 表示入射场 (Incident field) ;Psca (r) 表示散射场 (Scattered field) 。
在本次研究中, 我们关注讨论的是二维超声成像的问题。正如在图片1中描述的情形, 具有无反射边界的匹配材料通常被用于消除存在空气/组织的交界面的反射波。在这种条件下, 背景媒介是相匹配的液体, 该液体的参数特性是和正常组织是相似的。组织的声学特性主要由它的密度ρ, 声速c, 衰减常数a决定。因此, 当组织是各向异性时, 这些参数变量是关于组织的空间变量分布的函数。
图片1超声CT成像中针对正演问题的典型的几何示意图。发射和接受探头被安放在一个充满液体区域的外围, 该液体区域周围是吸收边界。
上图图片1描述了超声领域研究问题的普遍情形, 即一个非均匀体 (它的参数是由质量密度和衰减常数决定) 被浸泡在一个背景媒介中 (该媒介本身具有一致性的参数) 。从该图片中, 我们可以发现假如这个包含各向异性的组织的区域被定义为区域D, 那么由收发探头所围绕的整个完整区域被定义为区域S。因此, 在二维超声成像的情形下, 所有关于场和源的变量因子就是关于的r= (x, y) 函数。这个声波场振幅被假设是以频率w的条件下, 随着时间t变化出现震荡;相应的时间因子是exp (-iwt) 。声波的速度被定义为c。在区域D内部的声场的解可以被表示为一个积分方程的形式, 即:
其中, r表示声波源的具体空间位置;同时
是有固定密度ρb的背景媒介中的空间标量格林函数。此外, 各向异性区域中的复波数是k, 背景媒介中复波数是kb。是对比度函数。
在方程 (2) 中, 入射场强Pinc (r) 是来自于位于空间r处的发射源产生的场 (这是前提组织是各向同性的) 。因此, Pinc (r) 可以很容易被获得。在组织特性是多样性情形下, 这个各向异性的特性会产生一个感应的场强, 这个场就是散射场 (即在方程 (2) 第二部分) 。然而, 这个散射场是无法知道的;因为它基于分布在区域D中的未知的场强P (r) 。所以, 方程 (2) 确切上是一个积分方程定义了在区域D中的未知的场强P。这个正演问题的任务就是去求解区域D中未知场强P (r) 。一旦区域D中的场强P被求解, 区域D之外的其他区域的场强就可以通过格林函数很容易被得到了。诚然, 在正演问题中, 背景媒介和各向异性区域的参数特性均是前提已知的。考虑到散射的强度可以被假设是有一定范围有限的;所以, 这个散射的分布也就通常可以被定义包围在一个大型的圆区域D中。
3、正演成像的算法构造
在求解正演成像的方法中两个目的会被达到。首先, 该方法可以被用于模拟和校准相应的成像系统;其次, 该方法可以通过提供产生的模拟的“检测”数据用于检验反演算法的效果;此外, 很多的反演成像算法也许需要反复调用正演成像提供的数据结果。
3.1 扩展波恩近似算法
求解方程 (2) 的一个很有效的方法是采用波恩近似法 (BA) [5,6], 即假设方程 (2) 中声场可以被入射场Pinc近似。所以在波恩近似的假设下, 通过调用执行积分方程中格林函数可以很容易地计算出总场。因而方程 (2) 可以近似表述为:
然而, BA算法有着非常有限的适用范围, 仅仅当散射体的尺寸和 (或) 对比度比较小才成立。为了提高适用范围, 我可以采用EBA算法[7,8]。所以, 上述方程 (3) 表述的近似场带入积分方程 (2) 可以转化为:
其中:
3.2 采用FFT加速EBA:构造FFT-EBA算法
基于卷积定理, 我们可以把基于EBA方法求解的声场的方程转化为如下的形式:
其中, F和F-1分别表示二维的傅里叶正变换和逆变换。
3.3 EBA作为预处理器用于预处理CG-FFT (或BCG-FFT) 算法
我们这里的目的是求解方程 (2) 中在各向异性区域中的未知声场P。考虑到方程 (2) 是格林函数g (r) 和x (r) p (r) 之间的二维卷积的问题, 我们可以援引卷积定理把方程 (2) 从新写为如下算子形式:
上述方程然后可以采用CG-FFT算法[9,10]来进行迭代求解。我们此处采用EBA作为预处理器。因此, 不同于通常解法直接求解方程 (7) 那样, 我们变为求解如下一个等价的问题:
其中, 这个预处理器是EBA已经在方程 (6) 中定义过。既然这个预处理过程LEB-1是可以通过FFT算法 (方程 (7) 展示) 方式得到, 这个新的方程 (8) 同样可以采用CG-FFT有效求解出来。通过这种预处理方式, 这个设计方案倾向于比传统的CG-FFT更加快捷。方程 (9) 同样可以通过采用双共轭梯度-快速傅里叶算法 (BCG-FFT) [10,11]来求解。
值得指出的是:除了积分方程离散化这一过程外, EBA在预处理CG-FFT (或BCG-FFT) 时候并不存在近似计算的情况。换言之, EBA和CG-FFT (或BCG-FFT) 的组合是一种精确没有近似化的过程;该组合在对于高对比度和低对比度的散射体成像过程中均有效且保持良好的精度。此外, 我们可以定义一个可以接受的误差标准δmin去检查整个迭代过程。定义L2范数下的相对误差:
对于低对比度的问题, 假如EBA的初始值的结果误差已经小于δmin, 符合精度要求;那么CG-FFT (或BCG-FFT) 将不会被调用到。对于高对比度的问题, CG-FFT (或BCG-FFT) 会被调用到直到迭代过程中的相对误差小于δmin为止。相对比K空间下的传统的CG-FFT (或BCG-FFT) 方法, 我们这里的EBA预处理CG-FFT (或BCG-FFT) 的组合方法更加具高效性。
4、结论
本文构造出了一种快速迭代的预处理解析器和FFT组合的超声CT正演成像算法。相对经典的矩量法 (MOM) 等方法, 该方法在保持同等成像精度的条件下可以使得计算时间从O (N3) 或O (KN2) 减少到O (KNlog N) 以下;同时内存占用量从O (N2) 缩减到O (N) 以下, 此处的K表示迭代次数。该方法完全可被用于求解目标散射体是大尺度的复杂成像情况。
摘要：超声CT成像是采用散射场数据来重建物体内部结构的一种定量化的超声技术。本课题采用快速傅里叶变换 (FFT) 加速扩展波恩近似 (EBA) , 进而基于EBA来预处理共轭梯度-快速傅里叶变换 (CG-FFT) 的方法来实现正演成像的过程。该算法是在传统超声成像领域和商业超声软件中没有使用过的;同时, 在保持同等计算精度的条件下, 该算法对比传统超声成像正演算法具有速度更快等优点。