基于OpenMP的并行图像相减算法实现与分析*

李继良，于 策，孙济洲，商朝晖，陈锦言，曹 玮，张旭明

(1.天津大学，天津 300072;2.天津师范大学，天津 300387)

图像相减是一种精确地获得密集星体区域中变星的相对光度变化的技术，它在凌日行星观测、微引力透镜观测和寻找超新星的数据处理中都起到了关键作用。如果星的亮度是变化的，而星所在的背景是基本不变的，那么亮度变化(而不是绝对亮度)就可以通过在纠正视宁度(seeing)和图像放缩的差异之后进行图像相减得到。在一段时间内对同一天区连续拍摄可以得到一系列图像，将这些图像与参考图像相减，得到的图像上就只剩变星了，它们相对于参考图像的光度变化可以通过更专业的孔径测光得到。2010年下半年在南极昆仑站安装的AST3［1］也是利用图像相减技术作为其中的一个核心计算过程。图像相减的最早尝试始于Tomaney和Crotts［2］，他们采用的方法是对每幅图像上的亮星做傅里叶变换，通过其比值确定卷积核，然而这并不能保证相减后图像的质量足够高。Alard和Luption［3］提出了一种新算法，通过把卷积核分解为一系列的基本操作，求卷积核的各个像素的值时可以直接在图像域上进行，而不用先转换到傅里叶域。Alard［4］又将其由空间固定卷积核推广到了空间可变卷积核。这篇论文将采用基于OpenMP多线程的技术对C Alard在2000年提出的Optimal Image Subtraction(OIS)算法进行优化，主要目的在于进一步提高图像相减的效率。

当前采用多核处理器的计算机已成为市场主流。但是，如果用户软件无法利用多个核，就会造成CPU资源的极大浪费。因为在主频、高速缓存等其他设备完全相同的条件下，没有针对多核优化的程序在多核CPU上运行的速度与它在单核CPU上运行的速度并不会有明显的提高。实际上，串行程序仅仅利用了多核处理器中的一个处理核心，而其他处理核心则处于空闲状态。因此，面对多核平台开展系统软件和应用软件的开发工作将是软件开发人员必须面对的问题。采用多线程技术是提高资源利用率的一个有效方法，在单核处理器时代，应用程序已经支持多线程技术。然而，单核内的多线程运行是串行的，在某一时刻，只能有一段代码被CPU执行，单核处理器只能将多个指令流交错执行，并不能真正将它们同时执行。在多核平台，各线程都是在相互独立的执行核上并行运行的，如果线程数目小于等于执行核数目，那么各行程在运行时相互之间不存在对CPU资源的竞争，在某一时刻，可以有多个线程同时运行，达到真正意义上的并行处理。天文观测往往会产生TB量级的海量数据，能否在短时间内对这些数据进行分析处理已成为影响天文研究进展速度的关键因素之一。在多核时代，充分利用多个CPU核心带来的额外计算资源来提高数据分析处理的效率，对天文学研究具有重要意义。

1 OpenMP多线程编程简介

OpenMP［5］是由OpenMP Architecture Review Board牵头提出的，并已被广泛接受，用于共享内存并行系统的多线程程序设计的一套指导性注释(Compiler Directive)。OpenMP支持的程序语言包括C语言、C++和Fortran;而支持OpenMP的编译器包括Sun Compiler、GNU Compiler和Intel Compiler等。OpenMP提供了对并行算法的高层的抽象描述，程序员通过在源程序中加入专用的pragma预处理指令(pragma pre-processor directive)来指明自己的意图，由此编译器可以自动将程序进行并行化，并在必要之处加入同步互斥以及通信。当选择忽略这些pragma指令，或者编译器不支持OpenMP时，程序又可退化为通常的程序(一般为串行)，代码仍然可以正常运作，只是不能利用多线程来加速程序执行。

OpenMP提供的这种对于并行描述的高层抽象降低了并行编程的难度和复杂度，这样程序员可以把更多的精力投入到并行算法本身，而非其具体实现细节。对基于数据分块的多线程程序设计，OpenMP是一个很好的选择。同时，使用OpenMP也提供了更强的灵活性，可以较容易地适应不同的并行系统配置。线程粒度和负载平衡等是传统多线程程序设计中的难题，但在OpenMP中，OpenMP库从程序员手中接管了部分这两方面的工作。

但是，作为高层抽象，OpenMP并不适合需要复杂的线程间同步和互斥的场合。OpenMP的另一个缺点是不能在非共享内存系统(如计算机集群)上使用。在这样的系统上，MPI使用较多。

OpenMP的主要优势有:相对简单，不需要显式设置互斥锁、条件变量、数据范围以及初始化;可扩展，主要是利用添加并行化指令到顺序程序中，由编译器完成自动并行化;移植性好，OpenMP规范中定义的指导指令、运行库和环境变量，能够使用户在保证程序的可移植性的前提下，按照标准将已有的串行程序逐步并行化，可以在不同的厂商提供的共享存储体系结构间比较容易地移植。OpenMP的主要缺点是程序的可维护性不够好，再者当程序比较复杂的时候，编程会显得比较困难。

2 使用OpenMP编译指导指令对图像相减程序的优化

2.1 问题描述

图像相减是将已知的标准图像与要检测的图像进行相减，从而找到其中的差异，并利用已知的领域知识判断差异代表的实际物理意义或者差异根源的方法。图像相减的方法在包括医学等在内的多个涉及稳定图像的领域都有应用，其基本原理为通过同一区域不同时间的图像间差异来获得该区域的时间变化信息。应用到天文学领域主要是通过同一区域不同时间的图像间差异发现流量变化的星。图像相减变源测光算法主要包括了图像配准(Image Registration)、图像相减和相减后的检测等几个步骤。

图像相减变源测光算法是通过图像差异来获得星的流量随时间变化信息的，所以要有一幅图像作为参考图像，其他图像与之相减。但是由于观测时间、观测仪器的运动会导致参考图像和要相减的图像指向的可能并不是一个完全重叠的天区，他们之间总会存在微小的误差，这些误差主要体现在平移和旋转两个方面。所谓图像配准，就是为了消除图像间的平移和旋转。

然而图像配准后的天文图像还是不能直接相减。受不同的大气状态和云层遮挡等因素的影响，不同时间拍摄图像的星像的PSF(Point Spread Function，点扩散函数)是不同的，通俗地说，就是不同时间拍摄的图像中，同一颗星在图像上的大小不同，而且不同的遮挡效果会导致星的亮度存在差异。不同的大小导致同一颗星相减后无法完全消除其轮廓，会有残留的边缘，从而导致不能准确判断是否有变星。不同的亮度，导致图像相减后每颗星都会留有与其大小相同的“圆斑”，于是会在检测变源时多出不存在的变源。这些都要求首先将图像进行处理，然后再进行图像相减。记参考图像为二维函数R(x，y)，处理图像为I(x，y)。那么可以通过卷积将参考图像的PSF转化为处理图像的PSF，也可以反向将处理图像卷积到参考图像。由于参考图像选择的是“最好”的图像，PSF最小，星象最锐利，噪音最小，因此将参考图像卷积到处理图像会引入较小的噪声，从而得到较好的差分图像。这里的难点是需要确定卷积核(Kernel)函数，记为K(x，y)，并且K(x，y)是空间可变的。为了得到K(x，y)，要求解方程组R(x，y)*K(x，y)=I(x，y)，其中*代表卷积运算。在对参考图像进行卷积运算之后，所得图像和处理图像就可以进行对应像素之间的相减操作了。相减后仍可能有部分轮廓没有完全消除，需要做一次平滑处理。需要在经过平滑处理的图像中检测是否存在亮星，如果仍然存在亮星，就说明两幅图像中某颗星的亮度发生了变化或者存在参考图像中不存在的星，即要寻找的变源。

2.2 开源图像相减程序ISIS分析

ISIS是使用图像相减算法处理FITS图像［6］的一个完整的软件包，它包含3个子程序包，能够从一系列CCD图像中生成星体的光变曲线。为了得到星体的光变曲线，必须进行下面的3个步骤:(1)图像配准。图像配准的目的是将每个图像映射到同一个坐标系上，通常取其中的一幅图像做参考系统。这一过程的输出是在参考坐标系上内插的一幅FITS图像，包含两个步骤，获取天文转换方程X=f(x_ref，y_ref)和图像内插(双三次曲线)过程;(2)图像相减。这是程序的主要部分，也是文［3］和文［4］提出的新算法的核心部分。在运行程序之前，需要通过将一些最好的图像叠加起来生成一个好的参考图像。然后使用图像相减程序调整参考图像与原始图像的视宁度(Seeing)，原始图像在前面已经做过配准和内插操作。图像相减程序可以将整个图像分块做处理，这在使用有限的内存资源处理非常大的天文图像时尤其有用。程序对变星有两个不同级别的抛弃:检查每颗星是否表现出光流量变化，检查每颗星的卡方(Chi-square)。图像相减的最终输出为原始图像与参考图像相减后光流量变化的被减后图像;(3)测光。软件包将使用被减后图像中的变星做测光。变星的光度将会采用在固定位置的PSF拟合测光(profile-fitting photometry)来测量。像图像相减过程一样，这个过程也可以分块进行。

图像相减变源测光以c shell脚本的方式组织，其中interp.csh和temp.csh完成图像配准过程，subtract.csh脚本负责图像相减过程，detect.csh和find.csh使用相减后图像发现变源并找到位置，phot.csh脚本负责进行变源测光。为了得到并行优化的重点所在，本文对每个脚本的运行时间做了测量(测试中使用了BATC［7］拍摄的2 K×2 K图像，测试时使用了20幅图像，对每幅图像测量5次，结果为平均值)，见表1。

表1 处理的各个流程运行时间及所占百分比Table 1 Processing time and percentage of each step

从表1可见，在未使用优化选项时，卷积相减这一过程的时间占据了整个程序运行时间的52.6%;即使在使用了－02优化选项后，这一过程占得时间百分比仍高达39.9%。因此卷积相减这一过程就是程序运行的瓶颈所在。卷积相减，就成为使用OpenMP并行处理的关键部分，具体到源代码，则是subtract程序包下面的源程序。深入分析后发现，时间密集部分主要集中在kernel_convolve函数调用上，主要包括求卷积核和卷积部分，占据的时间分别为subtract程序包的8%和80%。鉴于以上结果，本文将主要考虑ISIS软件包的第2个步骤，也就是图像相减过程。

图1 处理的各个流程百分比示意图Fig.1 Time percentages of processing steps

2.3 对原算法的优化

用于科学和工程应用的大量程序被表示为基于迭代构造的形式，即它们是基于循环的。通过严格集中于循环来优化这些程序，是对较老的向量巨型计算机的一种传统回溯。将这种方法扩展到现代的并行计算机中，得到一种并行算法策略，在该算法策略中，并行任务被识别为可并行循环的迭代。这种模式及解决问题的方式是构造基于循环的程序，以进行并行计算。当可以获得现成的代码时，目标是将一个串行程序“演化”为一个并行程序，所采用的方式是对循环进行一系列的转换。理想上，所有的改变局限于对循环的转换，这种转换需要消除循环的相关性，但并不改变整个程序的语义(这称为语义中立转换)。根据多核计算机的系统特点，对循环体的算法进行调整、改造和优化设计，这是提高程序并行化效率的关键之一［8］。

多核计算机的各个处理器都配备有自己的局部高速缓存(Cache)，处理器访问自身局部缓存的速度要比访问主存或访问其他处理器的缓存快得多。因此并行化设计必须考虑如何合理地使用缓存，以避免在读写数据时产生相关与冲突。为了实现快速访问缓存，对于多重嵌套循环体的并行化设计，必须遵循如下规则:第一，尽量并行化最外层循环，使得在并行区域内获得最大的计算工作量，增加其并行粒度;第二，最里层的循环变量应该是变化最快的可变下标变量，而且在最里层应该按照顺序访问数组元素，也就是应按数组列访问。这样可以提高缓存局部性。因此，数组的最左下标变量应当作为最里层的循环变量，可以保证从主存按数组列顺序把数组元素取出并放入缓存之中，从而达到快速访问缓存，提高并行效率的目的。

kernel_convolve.c的伪代码描述如下:

以上程序kernel_convolve中，由之前程序中已经确定的卷积核(kernel)大小为mesh_size×mesh_size的方阵，将参考图像分块，每块大小与kernel大小相同(conv_step=mesh_size)，即为conv_step×conv_step。从而参考图像被划分为nsteps_x×nsteps_y数量的大小为conv_step×conv_step的stamp块进行处理。

在对图像分块后，就进入了影响程序性能的多重循环中。前两重循环遍历图像的每一个stamp块。并根据不同的stamp块求解相应该块的卷积核(make_kernel)。完成对kernel的求解后，开始对stamp进行卷积操作，三四层循环为遍历stamp块，在对块中的点进行卷积之前，要判断点的位置和值是否符合卷积要求，如果符合，则进行五六层循环，即像素点的卷积值是由此像素点周围conv_step×conv_step的点与kernel进行卷积之后的结果。因此五六层循环即是conv_step×conv_step。

在求卷积核make_kernel中，根据之前程序由最小二乘拟合求得的kernel_sol来求解空间可变卷积核的系数kernel_coeffs，之后对卷积核(kernel)进行初始化，每个卷积核的系数kernel_coeffs乘以一个基向量kernel_vec，即求得相应stamp的卷积核kernel。

原ISIS程序忽略了高性能计算机缓存性能的发挥，在最内层使用了列主序而不是行主序来遍历数组，因此对并行或者非并行的执行都是极为不利的。特别是多处理时会影响使用缓存的性能，使运算效率受到严重影响。所以在对原串行程序进行并行化设计时，必须考虑充分地运用缓存性能的编程策略，尤其要将最内层的循环由列主序遍历转变为行主序遍历。

在对图像的各个stamp块分别求卷积核和做卷积的过程中，各个stamp块的处理彼此没有相关性，因此，可以运用“分而治之”(Divide and Conquer)策略划分并行子任务。具体的，在最外层循环添加OpenMP指导语句，即#pragma omp parallel default(shared)num_threads(THREAD_NUM)和#pragma omp for nowait，可以将stamp块的处理并行执行，提高程序的效率。

2.4 性能与效率测试

对并行算法的测试主要是测试加速比和并行效率，可以通过式(1)和式(2)完成［9］:

其中，S(p)表示加速比;ts表示使用单处理器系统执行算法的时间;tp表示使用具有p个处理器的多执行机执行所需的时间。

其中E表示并行效率;p表示处理器的个数，并行效率可定义为加速比与CPU个数的比值。

以BATC［9］拍摄的2 K×2 K的图像进行了测试，在Dell PowerEdge T300 System服务器上进行。Dell PowerEdge T300 System服务器上配置了Intel Xeon X3326四核中央处理器，每个核带有3072 KB的缓存，工作主频为2.50 GHZ，内存为DDR2@667 MHZ，内存大小为2 GB，更多关于该服务器的硬件信息请参考Dell网站;系统使用的操作系统为64位Fedora 9系统，其中Swap空间大小为2 GB。

在使用串行程序时程序的平均运行时间为8.855 s，在使用OpenMP加速后程序的平均运行时间为5.049 s，加速比S(p)=1.754，并行效率为E=0.439，具体的测试结果见表2和表3。根据Amdahl定律，并行程序加速的上限为串行部分所占比例的倒数，程序的加速效果还是比较理想的。

表2 使用不同的线程数时相减部分所用时间Table 2 Processing-time values of subtractions using different numbers of threads

表3 加速比和并行效率Table 3 Speedup ratios and parallel efficiencies

3 结论与展望

在多核计算时代，采用并行计算方法加速天文学数据的分析和处理对天文学研究具有重要的意义。OpenMP是一种优秀的基于共享内存的并行方法，通过采用OpenMP技术将OIS图像相减算法进行并行化，程序的效率得到了显著的提高，使用BATC的天文图像进行测试后发现加速比可达1.754。然而程序的结果跟加速的理想结果还有一定的差距。原因有两个，首先是因为程序中还存在许多串行执行部分，因此，最大限度的并行化源程序，是获得理想的OpenMP并行程序并行效率的前提。其次，如果加入不当的并行指导语句，由于缓存的一致性问题，会出现某些变量的假共享，致使程序运行速度变慢。

GPU计算是一种新兴的科学与工程计算技术，计算密集部分采用GPU(图形处理器)而不是CPU进行。后续工作将考虑把OIS移植成基于GPU的算法，程序的效率应该会得到更大幅度的提升。

［1］毛银盾，唐正宏，郑义劲，等.CCD漂移扫描的基本原理及在天文学上的应用 ［J］.天文学进展，2005，23(4):304－317.Mao Yindun，Tang Zhenghong，Zheng Yijin，et al.The Basic Principle and the Application in Astronomy of CCD Drift-Scan ［J］.Progress in Astronomy，2005，23(4):304 －317.

［2］Tomaney A B，Crotts A P S.Expanding the Realm of Microlensing Surveys with Difference Image Photometry ［J］.The Astrophysical Journal，1996(112):2872.

［3］C Alard，R H Lupton.A method for optimal image subtraction ［J］.The Astrophysical Journal，1998(503):325－331.

［4］C Alard.Image Subtraction Using a Space-varying Kernel［J］.Astronomy and Astrophysics Supplement，2000(144):363 －370.

［5］Open MP ARB.OpenMP ［EB/OL］.［2011－03－14］.http://openmp.org/wp/.

［6］胡新华，邓元勇，王先平.FITS图像处理技术荟萃及在太阳观测中的应用 ［J］.天文研究与技术——国家天文台台刊，2008，5(1):55－65.Hu Xinhua，Deng Yuanyong，Wang Xianping.Collecting of Techniques on Dealing with Image Files of FITS Format and Applications to the Solar Observation ［J］.Astronomical Research ＆Technology——Publications of National Astronomical Observatories of China，2008，5(1):55 －65.

［7］BATC课题组.BATC HomePage［EB/OL］.北京:国家天文台，(2010－12－20)［2011－03－14］.http://batc.bao.ac.cn/c-index.htm.

［8］Timothy G Mattson，Beverly A Sanders，Berna L Massingill.并行编程模式 ［M］.敖富江，译.北京:清华大学出版社，2005:113－125.

［9］Harry F Jordan，Gita Alaghband，迟利华.并行处理基本原理 ［M］.刘杰，译.北京:清华大学出版社，2004:33－37.