（无线电物理专业论文）基于cuda的实时mip重建及其立体显示.pdf

、 2 0 1 0s c h o o ly e a rm a s t e r sd e g r e et h e s i s i n s t i t u t i o nc o d e ：1 0 2 6 9 n u m b e r ：51 0 7 0 6 0 2 0 5 4 屁s 兹砀面忍一砌 c u d a b a s e dr e a l t i m er e c o n s t r u c t i o n a n d s t e r e o s c o p i cd i s p l a yo f m i p d e p a r t m e n t ：d e p t p h y s i c s m a j o r ： r a d i op h y s i c s r e s e a r c hd i r e c t i o n ：n u c l e a r m a g n e t i c r e s o n a n c es o f t w a r e t e c h n i q u e s u p e r v i s o r ：y a n gg u a n g g r a d u a t e ：x i 髓i m a y , 2 0 1 0 i i 8川9 舢7m 3舢禾 舢7iiii-哪y 鹰嚯麓 华东师范大学学位论文原创性声明 郑重声明：本人呈交的学位论文幽洲鼬么如惴及寡刍钵璺互 ，是 在华东师范大学攻读砸匆博士( 请勾选) 学位期间，在导师的指导下进行的研究工作及 取得的研究成果。除文中已经注明引用的内容外，本论文不包含其他个人已经发表或撰 写过的研究成果。对本文的研究做出重要贡献的个人和集体，均已在文中作了明确说明 并表示谢意。 作者签名： 日期：f ，年, - - - n 柙 华东师范大学学位论文著作权使用声明 碗刁洲泌亥町化产弛否美易锛多而，系本人在华东师范大学攻读学位期间 在导师指导下完成的砚左博士( 请勾选) 学位论文，本论文的研究成果归华东师范大学 所有。本人同意华东师范大学根据相关规定保留和使用此学位论文，并向主管部门和相 关机构如国家图书馆、中信所和“知网”送交学位论文的印刷版和电子版；允许学位论 文进入华东师范大学图书馆及数据库被查阅、借阅；同意学校将学位论文加入全国博士、 硕士学位论文共建单位数据库进行检索，将学位论文的标题和摘要汇编出版，采用影印、 缩印或者其它方式合理复制学位论文。 本学位论文属于( 请勾选) ( ) 1 经华东师范大学相关部门审查核定的“内部 或“涉密学位论文宰， 于年月同解密，解密后适用上述授权。 ( i ) 2 不保密，适用上述授权。 导师签名拿岳喜一 幸“涉密”学位论文应是已经华东师范大学学位评定委员会办公室或保密委员 会审定过的学位论文( 需附获批的华东师范大学研究生申请学位论文“涉密审 批表方为有效) ，未经上述部门审定的学位论文均为公开学位论文。此声明栏不 填写的，默认为公开学位论文，均适用上述授权) 。 n i l 墨吖 名 鹳 人本 墨笠硕士学位论文答辩委员会成员名单 姓名职称单位备注 套组欺 截稻 安乐，胁欠压 磊丐余历j 彼蓰牟每，切大 毛蛆霍翻程疟拜呼- 切欠 i v 答辩时间：芦阜生善j d 摘要 医学成像中的3 d 物体的可视化是观察和分析人体器官和组织的重要工具。然而可 视化的计算速度限制了可视化的实时交互性，从而影响了其应用。本文对最大强度投影 ( m i p ) 进行加速，达到了实时计算的要求。在此基础上，利用最新的立体显示技术， 实现了m i p 图像的立体显示，使用户达到了很真实的显示效果。 我们首先用o p e n m p ，利用现代c p u 多核的特点，对m i p 进行加速，加速的倍数 接近c p u 的内核数。然后使用n v i d i a 的c u d a 平台，充分利用g p u 的计算性能， 实现了共享内存和利用纹理两个版本的m i p 。在实验条件下，纹理版本的实现达到了实 时显示的效果。本文还实现了子体积重建功能，即先估计静止组织信号强度，然后用一 个灰度映射方法改善对比度。 为了更加真实地显示物体，我们使用n v i d i a3 d v i s i o n 技术对物体进行3 d 显示。 对血管造影成像的数据利用d i r e c t 3 d 的立体纹理，用基于层的光线投射法进行立体显 示。通过1 2 0 h z 的显示器和配套的眼镜，可以观察到具有真实立体感的3 dm i p 图像并 可实时交互式地操作图像。 关键词：c u d am 口实时计算立体显示 a b s t r a c t m e d i c a lv i s u a l i z a t i o ni sap o w e r f u lt o o lt h a to f f e r sa ni n s i g h ti n t ob o d yo r g a n sa n d t i s s u e s h o w e v e r , i t sa p p l i c a t i o ni so f t e nr e s t r i c t e db yi t sl a c k i n go fr e a l - t i m ei n t e r a c t i o n b e c a u s eo fl o n gc o m p u t i n gt i m e i nt h i sp a p e r , w ea c h i e v e dr e a l - t i m em a x i m u mi n t e n s i t y p r o j e c t i o n ( m i p ) c o m p u t i n gb yo p t i m i z i n go r i g i n a la l g o r i t h m f u r t h e r m o r e ，w er e a l i z e d i n t e r a c t i v es t e r e o s c o p i cd i s p l a yo ft h em i pi m a g e sw i t ht h el a t e s t3 dv i s i o nt e c h n i q u e f i r s t l y , o p e n m pi su s e dt oe x p l o i t i n gt h ec o m p u t a t i o nc a p a c i t yo fm o d e mm u l t i - c o r e c p u st oa c c e l e r a t em i p a na c c e l e r a t i o nf a c t o rn e a rt h en u m b e ro fc o r e so ft h ec p ui s a c h i e v e dt h o u g h t h e nw eu s e dc u d at e c h n o l o g yd e v e l o p e db yn v i d i a t oe x p l o i tg p u s c o m p u t i n gp e r f o r m a n c e t w oa l g o r i t h m s w e r ep r e s e n t e dh e r e ，b a s e do nt e x t u r ea n d s h a r e d m e m o r yr e s p e c t i v e l y a m o n gt h e m ，t h et e x t u r e - b a s e da l g o r i t h ma c h i e v e sr e a l - t i m e c o m p u t i n gt a r g e t t h e nw ei m p l e m e n t e das u b - v o l u m e r e c o n s t r u c t i o na l g o r i t h mw h i c h e s t i m a t e st h es t a t i ct i s s u es i g n a li n t e n s i t y , a n da d j u s t sm i pi m a g ec o n t r a s ta c c o r d i n g l y f i n a l l y , w e u s e dn v i d i a3 d v i s i o nt e c h n i q u et od i s p l a yi m a g e m r ad a t ai s r e c o n s t r u c t e dw i t hs l i c e - b a s e dt e x t u r er a yc a s t i n ga l g o r i t h m w i t ha1 2 0 h zd i s p l a ya n dap a i r o fs p e c i f i cg l a s s e s ，t r u e3 dm i pi m a g e sc a nb ev i e w e da n dm a n i p u l a t e di n t e r a c t i v e l y k e yw o r d s ：c u d a ； m i p ； r e a l t i m ec o m p u t i n g ；s t e r e od i s p l a y 2 目录 摘要1 a b s t r a c t 2 第一章绪论4 1 磁共振成像4 2 血管造影成像4 3 数据的可视化5 4 最大强度投影5 5 本文工作7 第二章o p e n m p 加速8 1 2 3 第三章 1 2 3 4 第四章 1 2 3 4 第五章 光线投射算法介绍。8 s h e a r w a r p 算法介绍8 o p e n m p 的实现10 c u d a 与m i p 加速算法实现1 3 c u d a 平台简介1 3 使用c u d a 共享内存的实现2 0 使用c u d a 纹理缓存的实现2 4 子体积重建2 6 立体显示。2 8 透视投影2 8 视差与立体显示2 8 o p e n g l 四缓冲立体显示的实现2 9 n v i d i a3 d v i s i o n 立体显示3 0 工作展望3 4 3 第一章绪论 1 磁共振成像 磁共振成像( m a g n e t i cr e s o n a n c ei m a g i n g ，m r i ) ，是一种安全的成像技术【1 】，已 经广泛运用于医学影像领域。1 9 4 6 年，美国哈佛大学的p u r c e l l 及斯坦福大学的b l o c h 分别独立地发现了磁共振现象，于1 9 5 2 年获得诺贝尔物理奖。1 9 6 6 年e r n s t 发展了脉 冲傅里叶变换n m r 测谱方法。1 9 7 3 年，l a u t e r b u r 用线性梯度场对空间进行编码，首次 从实验上得到核磁共振图像。1 9 8 2 年，多层面技术的引入，使m r i 得以应用在临床。 随后m r i 技术在此基础上飞速发展。 m r i 成像的成像参数和对比度机制很多，即可以对人体组织成像，也可以对组织的 某些物理特性成像。与c t 相比，m r i 具有对软组织更高的分辨力，对人体无放射性伤 害等特点。并且，m r i 可以任意方位成像，也可以三维旋转成像，为临床医学提供了更 好的应用。m r i 的发展极大促进了对各种疾病的研究及其临床诊断。2 0 0 3 年，l a u t e r b u r 和m a n s f i e l d 因此获得了诺贝尔生理学或医学奖。 2 血管造影成像 血管造影成像( m a g n e t i cr e s o n a n c e a n g i o g r a p h y ，m r a ) 【1 】就是血管成像，是磁共 振成像技术中的一种，可以使人的心血管系统可视化。在m r i 中，生理运动是伪影的 重要来源。但是，生理运动产生的信号变化又可被反过来分析运动的情况。血液的流动， 就使得血管部分的信号强度和其他组织有很大的差别。 血管造影成像技术可以分为两大类：流入敏感法或飞行时间( t o f ) 法和相位对比 度血管造影( p cm r a ) 法。飞行时问法的原理比较简单，一个层面被射频脉冲重复激 发，静止的组织将被饱和，而流入这个层面的新鲜血液，将产生较大的信号。使用最大 强度投影算法可以在二维图像上显示血管。 与其他成像方法相比，磁共振血管造影成像有几个优点：一，是非侵入性成像方法， 对身体无危害，二，对比度高，分辨率高，三，可以提供血管的生理信息。 m r a 成像采集的都是三维的数据，不能被直接显示出来，目前将m r a 三维数据 转换成二维图片的方法有以下几种：最大强度投影，子体积重建，时域滤波【1 】。 4 3 数据的可视化 对于一些三维数据( 如m r a 和d t i 数据) ，需要将三维体数据适当地显示出来， 这样专业人员才能从数据中得到有用的信息。将一系列二维断面数据重构，并且显示物 体的三维结构，是医学信息可视化的研究内容【2 】。现在，医学信息可视化已经广泛应用 在诊断医疗，手术规划等方面。通过可视化技术，对数据进行处理，构造三维物体，定 位各个组织，可以帮助医生对病情做出正确的判断。 对体数据进行显示可以分为两种方法，面绘制和体绘带u 3 1 。面绘制是对三维物体的 表面进行显示，通过生成等值面，重构物体的表面，加上光照来绘制图像。而体绘制是 直接利用光和三维体数据的相互作用，表现物体的结构。从图片质量和信息上来说，体 绘制优于面绘制，但是体绘制一般计算速度较慢，不能满足实时交互的要求。体绘制常 用的算法有光线投射法( r a yc a s t i n g ) 【4 】，抛雪球算法( s p l a t t i n g ) 5 】，错切变形法 ( s h e a r - w a r p 6 1 ) ，三维纹理贴图法等等。 4 最大强度投影 最大强度投影法( m a x i m u mi n t e n s i t yp r o j e c t i o n ，m i p ) 7 ，是一种3 d 数据场的可视化 方法，把一个三维数据投影显示到二维图像上。对于二维图像上的每一个点，都有一根 视线穿过三维数据。最大强度投影算法，就是将二维图像中的每个像素的值设置成为对 应于该像素的视线经过的所有体素的强度的最大值。而根据投影算法，又有正投影和透 视投影。在正投影中，每根投影的视线平行，而在透视投影中，所有投影视线相交于同 一点( 我们的观察位置) 。本文中，若没特别指出，最大强度投影算法采用的都是正投 影。在磁共振血管造影成像中，血管呈高信号( 见图1 ) ，可视化时一般采用最大强度投 影算法。 5 图1m r a 成像原始数据图，其中的亮点为高信号的血管。 图2 m r a 成像数据m i p 图像 最大强度投影算法由于在投影时得到的是每根视线经过的体素的最大值，得到的投 影图像就缺少了深度信息。因此需要从多个视角进行观察，才能够获得立体感。对观察 角度实时交互旋转，这需要每秒显示2 0 幅左右的图像，才能流畅地显示。但m i p 重建 的原始数据量较大，重建速度较慢，对于典型的2 5 6 2 5 6 2 5 6 普通电脑上仅能达 到每秒1 帧的渲染速度，因此目前的商用软件中都采用事先设定视角变动的范围和步进， 渲染完成后再进行电影播放。为了实现m i p 的实时交互显示，必须对m i p 算法进行加 速。 对计算速度的要求使我们要对最大强度投影算法进行加速，一般可以分为两类。 第一种，软件方法【6 】【8 儿9 】【1 0 】【1 1 】。其中一些方法在得到加速效果的同时影响了图 像质量。由于我们需要的信息都在那些信号比较高的体数据上，k e u nh ok i m 8 对数据 用一个阈值进行了处理，但是使用阈值时对图像质量有一定影响。p h i l i p p el a c r o u t 6 使 用了s h e a r - w a r p 算法，并结合了数据分类方法，由于s h e a r - w a r p 算法中有两次采样，同 时采样的时候使用的双线性采样，对图像质量有影响。s 6 r e ng r i m m 9 1 使用了查询表。 y u nt i a n 1 0 q a x 寸体数据进行了分类，并且比较各个插值方法对m i p 结果的影响。虽然 在采样时使用最近点采样或双线性采样会比三线性采样快很多，同时图像质量会下降。 l u k a sm r o z 1 1 预先去掉了那些对m i p 没影响的体数据，采用了一个新的存储方法。 第二种，硬件加速方法【1 2 】 1 3 】【1 4 】【1 5 】【1 6 】【1 7 】。k e u nh og i m 8 对数据阈值处理的 同时，使用了计算服务器来计算。h a n s p e t e rp f i s t e r 1 2 】中使用特定的硬件加速卡 v o l u m e p r o ，价格昂贵。s i m o ns t e g m a i e r 1 5 q h 提出了一系列可视化方法的实现，包括 m i p ，但是它使用了n v i d i a 公司的特定o p e n g l 扩展，并且使用了当时的低级的着色 6 语言，在当时o p e n g l 着色语言( o p e n g ls h a d i n gl a n g u a g e ，g l s l ) n v i d i a 公司并 没有完全支持，移植性很差。d a n i e lr u i j t e r s 1 6 分析了在g p u 上实现的体绘制方法的瓶 颈，然后提出了在特定硬件上的解决这些瓶颈的方法。这些方法大都需要硬件的价格和 普及性的限制。另外它们采用的可编程的着色语言当时正在发展过程中，这些都影响的 了它们的通用性。 需要指出的是，以上这些方法虽然都有一定的加速效果，但是都没有达到实时计算 的要求。 5 本文工作 本文的目标是追求更快的m i p 算法实现。良好的用户体验需要实时交互的功能，而 实时交互需要实时计算的支持；快速的算法实现使我们可以考虑得到更优解的算法；更 快的计算速度使立体显示等新的交互模式具备实用性。这些都需要我们更好的优化算法 实现，更好的利用电脑的计算性能。 本文采用多种技术对m i p 算法进行加速。对于一般大小的三维体数据，利用n v i d i a 公司推出的通用计算平台c u d a ( c o m p u t eu n i f i e dd e v i c ea r c h i t e c t u r e ) 实现了实时计 算显示的目标，再次基础上进一步实现了m i p 图像的交互式三维立体显示。 第二章先介绍了s h e a r - w a r p 算法和光线投射法，介绍了o p e n m p ( o p e n m u l t i p r o c e s s i n g ) 并行计算，并且实现了m i p 的o p e n m p 加速，加速倍率和c p u 的核 的数量接近。 第三章详细介绍了c u d a 平台。利用c u d a 的共享存储器实现了s h e a r - w a r p 版本 的m i p ，然后又在c u d a 平台上利用纹理缓冲功能实现了采用光线投射法的m i p ，实 现了实时计算的目标。本章还比较了两个m i p 实现版本的优缺点。 第四章将m i p 进行立体显示。人在观察物体时，两眼看到的图像有视差，大脑会根 据两眼的得到的图像合成一幅立体的图像。而一般的设备上，左右眼观察到的是显示器 上的同一幅图像，没有真实的立体感觉。利用o p e n g l 四缓冲功能或者是n v i d i a 提供 的3 dv i s i o n 功能，可以将血管真实立体的显示。 第五章对未来的工作进行了展望。 7 第二章 o p e n m p 加速 光线投射算法介绍 光线投射是项运用广泛的3 d 数据场可视化的技术。对于观察点和最后图像中的 每个像素点，可以确定一条射线( v i e w i n gr a y s ) 。找到与这条射线相交的体素点，根据 这些相交的体素点计算图像上该点的像素值 1 7 1 1 7 ，如图3 。光线投射法有两种计算 方法，物体空间序和图像空间序。物体空间序，就是顺序遍历物体来计算。图像空间序， 就是一根根光线的顺序遍历计算。光线投射法，可以应用在最大强度投影，等值面重建 等重建方法上。 图像序的光线投射法，用遍历一条条投射光线的顺序，计算图像，有两个明显的缺 点，一，计算过程中有大量的三维坐标的变换，计算量大，二，对体数据的数据读取不 是有序的。这两点对算法的计算速度有很大影响。 e 图3 基本的光线投射，观察点和图像像素点确定一条视线 2 s h e a r w a r p 算法介绍 l a c r o u t e 和l e v o y 提出了s h e a r - w a r p 算法【6 】【3 】，做到了对数据的有序读取，并减 少了计算量。s h e a r - w a r p 算法并不是沿着投射光线计算每个像素点的值，而是将数据一 层一层的投影，这就做到了有序的访问数据。具体做法就是先将三维体数据沿片层方向 错切，投影到一个扭曲的中间图像，然后通过w a r p 变形得到最后的图像。如图4 。 8 v , e w m gr a y s e s h e a r 一1 - 图4s h e a r - w a r p 算法示意图a ，物体经过错切后，视线方向与层面垂直， 图5 s h e a r - w a r p 算法示意图b 错切之后，光线的方向与z 方向平行，将三维体数据按片层的顺序，逐层投影到中 间图像上。正投影时，一层数据投影到中间图像时插值系数不变。 三维体数据的自然坐标系为( x o ， y o ，z o ) ，最终图像的坐标系为( x ，y ， z ) ， 三维坐标变换用数学表达式表示为： x y z w | m 西 y o z 0 w o ( 2 1 ) 从公式。每一层的数据都投影到中间图像的不同位置，就像是移动每一层数据，然 后垂直投影到中间图像。这就是对层数据进行错切，错切变换的数学表示为 m s 蛔一 1o 01 oo 0 0 墨0 s i 0 10 01 9 ( 2 2 ) m 订= m 。f 。 1 0 - s i 0 01一s ：0 0010 0 oo1 剩下的部分就是对中间图像w a r p 变形 m w 呻= m 咖 1 0 邑0 o1s ：0 ， 0010 oo 0l 1 0 - s i 0 01 一s ：0 0o1o 000 1 ( 2 3 ) ( 2 4 ) s h e a r - w a r p 算法中，若视线与z 轴夹角很大，那么中间图像就很大甚至无穷大，为 了避免这个情况，在错切时要先选定一个主轴，x ， y ，z 中的一个轴为错切后的视线 方向( 选与视线夹角最小的轴) ，根据主轴方向，预先旋转好数据，然后再进行错切变 换。 综上，不考虑主轴方向，s h e a r - w a r p 可以表示为 m m 。一m 。m 曲。， ( 2 5 ) 图6s h e a r 。w a r p 算法流程图 3 o p e n m p 的实现 现在多核处理器已经成为了主流，但是对于单线程程序，多核电脑并不能发挥它的 优势；而对多线程程序，几个线程可以同时在多个核上执行，从而速度得到提升。o p e n m p 是一个支持共享内存并行编程的应用程序接口f 1 8 ，包括一套编译指令语句和运行库例 程和环境变量。它能很方便地将我们的程序改编成多线程的。o p e n m p 使用f o r k j o i n 并行执行模型。程序的有一个主线程( 单线程) ，运行到某些编译指导语句部分，程序 就根据编译指导语句指定的方式，分成几个线程执行，所有线程执行完之后再执行主线 程中随后的部分。程序以主线程、多线程并行、主线程、多线程并行的形式执行。 1 0 具体o p e n m p 编程细节参考o p e n m p 编程手册。 f j f 3 o000 主线程一 riri kn k n 并行域并行域 图7o p e n m p 执行模型，程序以主线程、多线程并行、主线程、多线程并行的形式执行 看下面的编译制导语句，指示要将f o r 循环并行执行。设置电脑的环境变量，双核 电脑中，将o m p n u m t h r e a d s 设置为2 ，这段f o r 循环就会分成两个线程执行。 在第- d , 节提到的s h e a r - w a r p 算法中，计算几乎集中在投影到中间图像和中间图像 向最终图像的变换中。我们可以将这两部分计算并行执行。对于投影计算部分，我们以 层为最小单位。我们要将k 层数据投影到中间图像，将这k 层数据投影分成几个部分并 行执行。如果电脑有n 个核，就将k 层数据分成n 组，一个核执行一组投影计算。需要 注意的是，在将k 成数据分成n 组之后，应该产生n 幅中问图像，这n 幅中间图像再以 最大强度投影类似做法合成一幅图像。然后再w a r p 变换。对于w a r p 变换部分，以行为 最小单位，分成n 个部分并行执行。 图8 0 p e n m p 加速流程图 在m i p 时，计算时间和投影图像的大小成正比。试验时，原始数据的相邻体素 距离对应于最终图像大小中为一个像素的大小。之后测试中图像大小均以此计算。 测试电脑为i n t e lc o r e 26 3 0 01 8 6 g h z ，内存为2 g ，操作系统为w i n d o w sx p 。 望堡堡堕囹堕塑塑墅堕囹婴兰虫 未使用o p e n m p 9 2 9 8 1 m s9 1 9 4 9 m s1 0 3 0 m s 堡旦iq 旦! 璺丛!三三q ! 鱼堡垒至! 墨至! 里兰兰笪堡璺 表格1 使用了o p e n m p 之后的实验数据，数据大小2 8 8 x 3 8 4 x 1 5 6 ，重复计算1 0 0 次 2 5 6 2 5 6 2 5 62 8 8 3 8 4 1 5 6 未使用o p e n m p9 1 5 2 4 m s ( 1 0 9 f p s ) 9 6 8 0 0m s ( 1 0 3 f p s ) 堡塑iq 望呈璺! 堡主q ! q ! 里兰! ! ：2 1 垒! ! ! 兰塑兰望堕! ! ：竺兰尘! ! 表格2 未使用和使用了o p e n m p 加速的计算时间，数据大小2 8 8 3 8 4 1 5 6 ，重复计算1 0 0 次 从以上数据可以看出，s h e a r - w a r p 算法的计算量集中在错切投影部分。采用o p e n m p 并行计算后，m i p 的计算速度有了提升。实验电脑为双核，理论加速就是2 倍。从实验 数据上可以看出，使用o p e n m p 的版本比未使用o p e n m p 的版本快了1 7 8 倍，加速效 果比较明显，但是距离实时计算的要求还有明显的差距。 1 2 第三章c u d a 与m i p 加速算法实现 1 c u d a 平台简介 c u d a 平台是n v i d i a 公司推出的通用计算平台，利用它，可以用c 语言进行编 程，使用显卡上的图形处理器( g p u ) 进行计算【1 9 1 【2 0 】。g p u 与c p u 相比，计算性能 和性价比上都有很大的优势。g p u 是为了满足显示需要设计的，而3 d 建模和显示需要 大量的计算，这就要求g p u 有很强的计算性能。g p u 经过发展，已经变成高度并行化、 多线程、多核的处理器，计算能力和存储器带宽已经远远地超越了c p u ( 图9 ) 。在c u d a 推出之前，已经有了g p g p u ( g e n e r a l p u r p o s eg p u ) 2 1 ，可以利用图形处理器进行通用 计算。但是由于g p g p u 编程的接口不友好，编程的难度大，因此n v i d i a 公司推出了 c u d a 平台，以充分发挥g p u 的性能。 图9 g p u 和c p u 计算- 性能的发展 不过并不是所有计算都适合在g p u 上运算。一般来说，适合在g p u 上运算的，计 算强度要大，还有，数据要有很高的并行度。有一些任务，如图像去噪、三维数据处理， 矩阵运算等，天生具有数据并行的特点，非常适合用c u d a 编程。c u d a 还包括几个 应用函数库，c u f f t 和c u b l a s ，实现了傅里叶变换和矩阵运算，使c u d a 可以方便 应用在一般科学计算上。图1 0 给出了些图像处理算法在c p u 和g p u 计算时间对比， 可以看到g p u 有很大的优势。 c u d a 在推出初期，只支持单精度浮点数运算。2 0 0 8 后，计算能力为1 3 或更高的 显卡支持双精度运算，但是执行双精度运算时性能会将会降低到原来的1 1 0 左右。 1 3 图1 0g p u 和c p u 在图像处理的计算性能对e t l 2 2 ，在各个算法中g p u 的计算时间均比c p u 的实 现短，最快的只有c p u 的1 7 5 1 2 2 】 舢c u d a 的硬件结构 为了编写出高效率的c u d a 程序，我们有必要了解一下它的硬件结构【1 9 】。 如图1 1 所示，一块n v i d i ag e f o r c e9 8 0 0g t 显卡上有全局存储器和1 6 个多处 理器( m u l t i p r o c e s s o r ) ，而每个多处理器上又有8 个处理器( p r o c e s s o r ) 。那么一块这样 的显卡上就是有1 2 8 个核，可以有1 2 8 个活动的线程在同时执行。多处理器上除了有处 理器外，还有共享存储器，常量缓存，纹理缓存和指令单元。每个处理器上还有一组寄 存器。其中最重要的就是共享存储器，它的数据是被运行在该处理器上的线程所共享的。 全局存储器的数据是所有线程共享的。 c u d a 的模型为s i m t ( 单指令，多线程) ，类似于s i m d ( 单指令，多数据) 架构。 但是c u d a 有所不同。c u d a 3 2 个线程( 一组3 2 个线程组织成一个w a r p ) 执行同一条 指令，但是它允许指定单一线程执行和允许分支。在程序运行时，线程按照顺序被组织 成一个个w a r p 。在w a r p 中若执行路径有分支，连续执行各个分支，同时禁用未在此分 支上的线程，总的执行时间为各个分支执行时间的总和。分支只是在w a r p 内出现，不 同的w a r p 独立执行。 1 4 ：o 掣眦 m ”憾印饿硼喇r 舛 图1 1c u d a 的硬件结构，g p u 由n 个多处理器组成，每个多处理器上有m 个处理器1 1 9 】 b 1c u d a 的计算模型 g p u 是高度并行化，多线程的多核处理器，c u d a 是可编程g p u 发展而成的可伸 缩的编程模型，它是c 语言的扩展。c u d a 的核心有三个重要抽象概念：线程组层次 结构、共享存储器、屏蔽同步。 进行c u d a 编程，首先要熟悉一个概念，什么是线程。线程是一串连续的指令。 一个线程内，执行的代码是顺序的。而线程之间，程序是并行执行的。 在c u d a 编程中，程序员定义一种称为内核( k e r n e l ) 的函数，在调用时，还要指 定运行该内核的线程数量。一个n 线程的c u d a 内核将在g p u 上并行执行n 次。 c u d a 的线程组织有两个层次。第一层次为b l o c k ，第二层次为t h r e a d 。一个内核 函数分成多个b l o c k 执行，总的线程数量为b l o c k 的数量乘以每个b l o c k 的线程数量。 线程调度时以b l o c k 为单位。联系前面提到的g p u 上的多处理器和处理器的结构，c u d a 把一个b l o c k 分配给一个多处理器去执行，只要多处理器上的资源足够，几个b l o c k 可 以分配到同一个多处理器上，一个多处理器上最多并发执行8 个b l o c k 。 每一个线程，它们都有着线程编号索引和一个唯一的i d ( 从o 到n 1 的整数) 。线 程编号索引用内置的变量类型表示，b l o c k i d x ，t h r e a d i d x ，在线程执行时，由c u d a 对 这两个变量赋值，程序运行时每个线程可以读取这两个变量获得自己的索引。b l o c k i d x 和t h r e a d i d x 变量有三个数据成员：x ，y ，z ，对应到现实世界的三维，可以方便的构建 1 5 一维，二维，三维的线程b l o c k ( 如要计算一幅二维图像，可以将每个像素视为一个线 程进行编程) 。根据这个线程索引的不同，我们就可以使它们执行不同的操作，操作在 不同的数据上。 c u d a 加入了一个新的操作符，三个尖括号的操作符。上面的例子中出现了一个变 量b l o c k d i m ，为d i m 3 类型，有x ，y ，z 三个数据成员，是用来表示b l o c k 的大小的。上 例中第一个数b l o c k n u m 指定了总共有多少个b l o c k ，第二个数t h r e a d _ n u m 指定了每个 b l o c k 有多少个线程，这两个变量都是d i m 3 类型。 以图1 2 为例，对一个二维的( d i m x = 1 2 ，d i m y = 6 ) 的线程组织，我们把线程如 下分块 那么b l o c k ( 1 ，o ) 中的t h r e a d ( 3 ，2 ) 在二维空间的索引就为 对应的线程的i d 是y 木d i m x + x = 3 1 。 1 6 f 霉茹”一一 ； ；黧麟黧鬻溅溺 ! 粼粼黼黧溺 。 、0 。 毒： 、。一 、t 一 瞄船矗f 戗l 氛托灌连赶o l l 秘州啦d f 工a i l l k 鬯a 巷f 置村饲 图1 2c u d a 中的线程层次结构，图中总共有3 x 2 个b l o c k ，每个b l o c k 有4 x 3 个线程【1 9 】 c ) c u d a 是c 语言的扩展 c u d a 对c 编程语言的扩展共有四种： a )函数类型限定符，指定函数是在主机上还是设备上执行，以及函数是可通过 主机还是可通过设备调用。有d e v i c e ， g l o b a l 。h o s t d e v i c e 限定的函数被设备调用，在设备上运行。 菩o b a l限定的函数被主机调用，在设备上运行。 h o s t 限定的函数被主机调用，在主机上运行。 b )变量类型限定符，指定一个变量在设备上的存储器位置，有d e v i c e 。 c o n s t a n t _ _ _ , s h a r e d ； d e v i c e 限定符用于声明在设备上的变量。 c o n s t a n t 漫定符用于声明在固定存储器上的变量。 s h a r e d _ _ 限定符用于声明在共享内存中的变量。 c )一个新的操作符 ，指定内核调用是的执行配置； d ) 四个内置变量，用于指定网格和块维度以及块和线程索引，有g r i d d i m ， b l o c k i d x 。b l o c k d i m 。t h r e a d i d x 。 多线程运行时的数据共享和同步 c u d a 提供了一个轻量级的同步函数u s y n c t h r e a d s 0 。在一个线程块中，通过这个 1 7 函数指定同步点，同一个块中的线程执行到这句时必须等待其他线程执行到这句，才可 以继续往下执行。 在c u d a 中，每个线程都有各自的本地存储器。每一个b l o c k 有一个共享存储空间， 共享存储的生命周期与b l o c k 相同，可以被b l o c k 上的所有线程访问。 同一个b l o c k 的线程可以通过在共享存储中交换数据来彼此协作。在交换数据时必 须先同步数据，以保证每个线程都完成了数据写入。在下面的代码t d a t a 在分配在共享内存上，每个线程读入一个数据放在d a t a 数组。在进行线程同步 之后，块上的线程就可以见到其他线程读取的数据，协同完成计算。 d 1c u d a 的数据读取和性能优化 c u d a 编程的性能瓶颈是对全局存储器的读取，频繁的读取会极大的降低程序的性 能。在编写c p u 程序时，有时也会遇到数据读取影响程序性能的情况。c p u 采用了多 级缓存来减少数据读取的延迟。g p u 则不同，它并没有采用缓存的设计，它采取的是一 种读取延迟隐藏的设计。在一个多处理器上有很多个w a r p ，线程调度程序会选择一个活 动的w a r p 执行，当一个w a r p 需要读取全局存储器上的数据时，它就让出处理器，让其 他已准备好的w a r p 执行。不同w a r p 的执行和读取数据可以同时进行，这样就掩盖了数 据读取的延迟。事实上，一次读取全局存储器需要几百个时钟周期，因此，频繁的读取 全局存储器对程序性能影响很大。为了读取的效率和减少对全局存储器的访问，c u d a 对半个w a r p ( 也就是1 6 个线程) 的读取操作接合，减少对全局的读取。若满足读取条 件，半个w a r p 的读取只需一个或两个存储器事务就可以完成。具体的接合条件参见 c u d a 编程手册【1 9 】，各个计算能力的c u d a 显卡有所区别。 另一个有效的方法是先将数据读取到共享存储中。由于共享存储器的读取速度比起 全局存储器快很多，并且可以数据共享，可以先将数据读进共享存储再计算。块内线程 之间就可以协作读取数据，重复利用，减少读取。 1 8 若很难利用共享内存和做到接合读取，可以使用缓存，尽管速度比共享存储器慢， 对已经读取过的数据可以在缓存中直接得到，减少对全局存储器的访问。同时缓存提供 了一些很实用的功能，比如线性插值等。 共享存储器的访问速度远远大于全局存储器，但还是应避免访问共享存储时的存储 体冲突。共享内存被设计成1 6 个大小相等的存储块，称为存储体( b a n k ) 。访问1 6 个 不同的存储体可以同步实现。如果存储器请求同时访问一个存储体，则需要序列化，降 低了有效带宽。 寄存器是最快的存储器，编译器和程序调度程序负责寄存器的访问冲突问题，但是 在线程数为6 4 的倍数时最优。 c u d a 编程时，基本遵循以下几个规则： a ) 在程序设计时，尽量使算法并行化 b ) 优化全局存储器的读取 曲减少共享存储器访问冲突 d ) 在一个w a r p 中尽量避免分支 e ) c u d a 程序的结构： a ) 将内存中的数据复制到g p u 存储器上； 执行c u d a k e m e l 函数，将结果写入g p u 存储器上； 将结果复制回内存中； b ) 一般使用共享内存的c u d a k e r n e l 函数结构： 分配共享内存，将全局存储器上的数据读入共享内存； 调用同步函数，保证读入数据操作已经完成； 使用共享内存中的数据进行计算； 将结果写入到g p u 存储器上； nc u d a 的应用 c u d a 开始使用在一些商用软件、科研计算【2 1 】【2 2 】【2 3 】【2 4 】【2 5 】【2 6 儿2 7 】【2 8 】和可视 化【2 9 】上，得到了很好的加速效果。y u s u k eo k i t s u 3 0 在图像重建中使用了c u d a ，得到 了2 4 倍的加速比。基于g p u 视线的f f r 与基于多核c p u 的f f r 快了1 倍 3 1 1 1 3 1 1 ，很 多日常的应用程序在g p u 上都有很好的加速效果【3 2 】。 c u d a 非常适合在m r i 数据处理上的应用【2 1 】【2 1 】，某些m r i 成像数据量( 如m r a ) 1 9 巨大，原始数据需经过处理才成为一幅幅图像，数据处理的时间一般都较长。同时一些 优秀的算法的计算强度很大，计