技术开发 频道

医学信号与图像处理算法中的并行化

  【IT168 技术】

  注:本文为IT168&NVIDIA联合举办的“如何并行化我的应用”方案征集活动参赛作品。本次方案征集活动详情见:http://cuda.itpub.net/thread-1299715-1-1.html。近期活动的大部分方案,将会逐步与大家分享,不可错过哦!

  CUDA ZONE专区:http://cuda.it168.com/

  CUDA技术论坛:http://cuda.itpub.net

  一、项目应用概述

  自从1895年伦琴(Ronggen)发现X射线以来,特别是它在医学上得到应用以来,现在医学影像学迅速形成和发展,这不仅是自然科学史上的一个重大里程碑,而且在相当程度上改变了医学科学尤其临床医学的进程,为人类的疾病防治做出了巨大贡献。近年来,我国医疗影像仪器制造行业发展迅速,一些中档设备技术的转移强化了本土研发和生产规模,相关海归人才进一步提升了国内企业的竞争力,部分医疗器械产品出口逐年增加,在国际市场上也占有相当份额。因此,医学信息处理日益成为学者,行业人士和普通大众关注的焦点。

  对于传统医学信号与图像处理大量依赖专用的硬件设备,不仅致使造价昂贵,而且很多性能瓶颈无法得到突破。众所周知,信号处理上一些重要的变换,如FFT,希尔伯托变换等,在硬件实现上的成本异常惊人,达到难以接受的程度,只有在高端仪器上才会出现。

  近年来,信号与图像处理无论在算法上、系统结构上,还是在应用上以及普及程度上取得了长足的进展。但同时也面临许多挑战,其中最主要的问题就是如何提高解决实际复杂问题的综合能力,就当前技术水平来说,这种综合能力包括图像处理的网络化、复杂问题的求解与处理速度的高速化。

  图像并行处理技术是图像处理中的一个重要方面,是提高图像处理速度最为有效的技术。图像并行处理技术的发展难度大,其原因不仅在于对系统的硬件及系统结构本身,以及对计算机技术和集成电路等技术的依赖关系,而且在于实际应用的复杂性和应用部门对系统价格的承受能力。但图像并行处理技术发展产生的效益也十分显著,它在处理速度上所获得的加速比是令人振奋的。

  多核并行处理的技术的快速发展,为并行信号与图像处理提供给了强大的助力。特别是随着显卡的发展,GPU越来越强大,而且GPU为显示图像做了优化,同时也克服了多核CPU的一些性能瓶颈,在计算上已经超越了通用的CPU。如此强大的芯片如果只是作为显卡就太浪费了,因此NVidia推出CUDA(Compute Unified Device Architecture),让显卡可以用于图像计算以外的目的。

  以上概述了医学信息处理的背景,应用现状和业务瓶颈。第二节将分析并行化医学信号和图像处理的一些思路,在第三节中会给出一些应用实例。

  二、 医学信号与图像处理的并行化设计思路

  随着带有并行处理功能计算机的发展,大型科学与工程计算成为可能,使得科学技术作为科学研究的一种有效手段,已上升为与科学理论和科学实验并重的三大科学方法之一。而且随着计算机硬件、软件及算法的进步,这种趋势正在加强。在并行计算领域,并行体系结构、并行软件和并行算法三者缺一不可,而其中并行算法则是核心和瓶颈技术。

  所谓并行算法是指在各种并行计算机上求解问题和处理数据的算法,其本质是把多任务映射到多处理器中执行,或将现实的多维问题映射到具有特定的拓扑结构的多处理器上求解。并行处理作为计算机科学长期研究的一个重大课题。它主要依据于计算机系统的体系结构的3个基本概念:时间重叠、资源重复和资源共享。

  并行算法的实现强烈的依赖于计算机硬件和软件环境,因此在设计并行算法时有必要使用一些并行计算模型(也称编程模型)将各种并行计算机的基本特征抽象出来,形成一个处于具体并行计算机之上的抽象并行计算机,用于设计并行处理算法。迄今为止已有多种并行计算模型,如针对共享存储并行计算机的PRAM模型、针对分布式并行计算机的LogP模型等等。那么设计并行处理算法就可以依据于这些并行计算模型进行,然后将其映射到具体的并行计算机中去。

  并行算法的一般设计方法大致包含以下3个途径:

  (1) 检测和开发现有串行算法中固有并行性而直接将其并行化,这种并行设计思路比较直观,就是挖掘现有算法中的计算无关性,将任务映射到多处理器上,但对于具体的硬件平台则需要有更多的考虑,例如CUDA平台下存储器的访问方式,下一节会以CUDA平台下的box filter运算为例,进行一些分析。

  (2) 从问题本身特征出发,设计一个新的并行算法,这类情况主要是指问题本身不具备计算无关性,或者迭代或者数据相关。如图像处理在的直方图处理,要产生一个直方图必定要遍历每个像素点,而且像素点之间的操作具有相关性。下一节会分析一个CUDA平台下的直方图运算。

  (3) 修改已有的并行算法使其可求解另一类相似问题,一些经典并行算法设计思路的借鉴,也是信号与图像处理并行化的重要手段。下一节会以Fan-in算法的分析与改进进行说明。

  三、 一些并行化技术应用实例

  (1)CUDA平台下的box filter

  Box filter是图像处理中非常常用的滤波处理方式,而且像素间是无关的。在并行算法设计中,主要有下面两种并行思想:细粒度和粗粒度。Box filter采用细粒度的并行方式,是和串行处理方式一样,将一帧图像细分为不同的处理过程。很明显,每一个点的处理需要n2-1步.这里可以有一个提高,分别沿着X和Y轴,完成两次一维的平均计算从而实现二维平均,这个并行过程可以在GPU上实现,这样就只需要n-1步了。

  而粗粒度的并行实现方式,让每一个线程处理一整行和一整列,这样的方式比起前面的细粒度方式,能获得更好的性能。这是因为,粗粒度的处理方式,在两个GPU遍里面,每个线程能处理一个单一的点。然而,沿着X方向遍历的时候,对于全局存储器的读取方式是不满足合并访问规则的的,这是因为每个线程读取的是不同行的数据。采取纹理存储器能很好的规避这个问题,而对于Y方向的遍历就是满足合并访问规则的了。

  (2)CUDA平台下的直方图生成

  传统的图像直方图算法是通过逐点统计实现的,基本程序算法如下:

for(int i = 0; i < BIN_COUNT; i++)

  histogram[i]
= 0;

  
for(int i = 0; i < dataNum; i++)

  histogram[image[i]]
++;

 

    由于要确定的是整个图像的灰度分布情况,所以即使在并行处理环境中,不进行迭代运算也是无法生成图像的直方图的。要利用并行处理环境进行算法优化可以考虑将输入数据分成多个块,让某一个线程去处理其中一块形成一个局部直方图,最后归并成一个全局直方图。但是对于GPU平台这个算法需要进行调整,这是因为读取图像数据是顺序的而访问结果直方图数组是数据依赖的,考虑到GPU平台的存储器结构最优的是采用共享存储器而非全局存储器。GPU硬件上共享存储器的大小为16KB,共享存储区没有原子操作的支持。还有就是一个线程块应该包含128~256个线程以获得高效的运行。那么如果采用一个线程直接对应于一个子直方图的话,一个明显的限制就是16KB平均到最大一个包含192条线程的块,每个线程只有85bytes,因此这种方法最大情况共享存储区能适合每线程子直方图可达64单字节针脚数。字节计数器也引入了255字节的限制于单个执行线程所处理的数据大小,在执行线程间分割数据时需要考虑。也就是单个线程不能生成针脚数大于64的子直方图。

  由于硬件提供了执行方式为SIMT的warp块(块内32 个并行线程为一组来创建、管理、调度和执行线程),利用这一重要性质就可以设计共享存储区的原子操作。使用的方法是线程间不再有独立的子直方图而是每个warp(32个线程)共享同一存储区域生成子直方图,用这种方式存在warp内部共享存储区冲突。但硬件平台并未提供相应原子操作,所以这里需要利用一点程序设计技巧来解决冲突。下面就是解决冲突的方法:

__device__ void addData256(volatile unsigned int *s_WarpHist,unsigned int data,unsigned int threadTag)

  {

  unsigned
int count;

  
do{

  count
= s_WarpHist[data] & 0x07FFFFFFU;

  count
= threadTag | (count + 1);

  s_WarpHist[data]
= count;

  }
while(s_WarpHist[data] != count);

  }

 

    这里图像数据值是从全局存储区读入的,取值范围是[0,255],每个warp线程根据输入的图像数据使数组s_WarpHist[]相应位置增加1,数组s_WarpHist[]是与块内当前warp对应的块子直方图s_Hist[]的一部分。为了解决warp内的线程冲突,通过最末一个执行写入的线程将直方图标识,标识为直方图频数变量的高5位,因为warp的大小为32,所以只需要5位。这样的话,对于每一个线程读取当前直方图针脚值高5位就被替代为当前线程的标识符,然后将增加的频数写回子直方图的共享存储区。如果线程间要写入的不是同一针脚位置,就没有额外的操作,但当两个或多个线程在写同一位置时硬件会执行共享存储区写结合,使线程中被标识的频数被接受而拒绝所有其它未定的线程修改。在写尝试后每个线程读取来自同一个共享内存的位置。而那个可以写回自己修改的频数的线程,则退出循环并空闲等待warp中其余的线程。Warp继续它的执行,直到所有线程退出循环。由于每个warp使用自己的子直方图并且warp线程总是同步的,所以不依赖于warp的调度顺序,这个循环执行的迭代次数不会超过32次,而且只有在warp内所有线程读到的灰度值相同时才会发生。

  子直方图(即warp直方图)计算算法中的第一步count = s_WarpHist[data] & 0x07FFFFFFU;是将当前该针脚的频数放入count变量中,并将高位清零。

  第二步增加相应针脚频数,并写入线程标识符,这里特别说明的是threadTag = threadIdx.x << (32 - WARP_LOG_SIZE);也就是将线程号经过移位操作放到高5位作为标识符。

  第三步则是将增加后的频数写回到数组s_WarpHist[]。然后判断当前写入的频数是不是最后修改频数的线程所写入的,是就接受,不是则通过循环重新执行上述操作。

  在完成了子直方图计算后就是归并每一个子直方图,这个过程应该分两步首先归并各个warp的子直方图到块直方图,然后将块直方图归并为全局直方图。但因为对于全局存储区存在原子操作,所以块内或块间的线程可同时更新全局直方图而不会写冲突,所以不需要每个块输出自己的子直方图。

  (3)CUDA平台下的Fan-in及优化分析

  所谓Fan-in算法求和就是利用多处理器实现树型求和,图3.1是一个树型求和的实例:

 一些并行化技术应用实例

图3.1 Fan-in算法处理实例

  可以看到该算法是利用多处理器的并行计算能力,将计算时间缩短,也可以很清楚的看到并行处理程序的特色“浅而宽”。

  在本算法中Fan-in算法求和方法主要应用在直方图的求和计算中,全局直方图的针脚数为256。使用Fan-in算法求和步骤如下:

  第一步将0和128,1和129,…,127和255求和,将结果保存在0,1,…,127中;

  第二步将0和64,1和65,…,63和127求和,将结果保存在0,1,…,63中;

  继续重复以上计算方式,直到只剩一个元素为止,最后结果保存在元素0中,这里要特别说明的是结合CUDA程序特性线程在取数据的执行方式是线程0→线程1→……→线程n,有局部性原理可以知道,最优的数据寻址方式是相邻寻址,所以这里采用的是0和128这样间隔求和而非0和1这样的相邻求和叠加模型,即图3.3

  这种Fan-in算法的另外一种实现形式是纵横加工算法[30],其核心思想与Fan-in算法一致,但在实现上略有不同。算法如下:

  假设,N个数a1,a2,…,an,其中N=P2

  第一步:a1+a2,…,ap+1+ap+2,…,aN-P+1+aN-P+2

  第二步:(a1+a2)+a3,…,(ap+1+ap+2)+aP+3,…,(aN-P+1+aN-P+2)+aN-P+3

  第P-1步:(a1+a2+…+aP-1)+ap,…,(aP+1+aP+2+…+a2P-1)+a2P,…,(aN-P+1+aN-P+2+…+aN-1)+aN

  第P步到第P-1+log2P步按Fan-in算法处理。

  由此我们对以上两种算法做简要的理论分析,Fan-in算法使用的处理线程个数为N/2,加速比为Sf和效率Ef分别为

  而使用纵横加工算法使用的线程个数为 ,加速比Sz和效率Ez分别为

  从上面的式子可以看到,传统Fan-in算法的加速比高于纵横加工算法,但后者为常数效率的并行算法,也就是效率不随问题规模和处理器个数的改变而改变,也就是说当问题规模较小,线程切换延时也比较小的时候采用Fan-in算法,程序运行时间更快,当超过CUDA的kernel函数承担线程极限,函数切换延时较大时则应采用纵横加工算法获得的执行效率更高。本文根据所要解决的实际问题采用的是传统Fan-in算法的基本改进方式。

 一些并行化技术应用实例

图3.2 CUDA程序中直方图求和的Fan-in算法应用

0
相关文章