基于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.