一種地理柵格數(shù)據(jù)任選多邊形區(qū)域的直方圖實時生成方法
【技術領域】
[0001] 本發(fā)明屬于地理信息處理技術領域，具體設及一種地理信息系統(tǒng)中地理柵格數(shù)據(jù) 任選多邊形區(qū)域的直方圖實時生成方法。
【背景技術】
[0002] 直方圖是一個統(tǒng)計學概念，它通過對選定地域某一方面特征（如海拔、顏色等)進 行統(tǒng)計來表征該地域的全局特征。直方圖廣泛應用于許多領域比如圖像處理、數(shù)據(jù)挖掘W 及數(shù)據(jù)庫應用。
[0003] 隨著地理空間數(shù)據(jù)采集技術的發(fā)展，人類獲取地理空間數(shù)據(jù)的效率不斷提高，用 于分析處理的地理空間數(shù)據(jù)的規(guī)模越來越大。與此同時，隨著地理數(shù)據(jù)分析需求的多樣化 發(fā)展，用戶往往需要通過交互的方式快速得到自己選定空間區(qū)域的直方圖數(shù)據(jù)。目前絕大 多數(shù)GIS軟件如ArcGIS、SuperMap、MapGIS等均具有生成直方圖的功能，但是都不支持任選 多邊形區(qū)域的直方圖生成，而且其直方圖生成耗時比較長，遠不能滿足與用戶實時交互的 需要。
[0004] 當前國內外研究直方圖快速生成算法主要有兩個方向。
[0005] 一個是利用CPlKCentral Processing Unit)進行直方圖生成。隨著技術的進步， CPU處理速度越來越快，計算的性能可W隨著CPU性能和數(shù)量的增加而不斷提高，相對而言 I/O性能的增長跟不上計算性能的增長。在地理柵格數(shù)據(jù)處理領域，當處理大規(guī)模數(shù)據(jù)時， 低性能的I/O將成為影響整體性能的瓶頸，運嚴重制約了生成直方圖運算效率的提高。因 此，利用CPU進行直方圖生成往往是通過減少數(shù)據(jù)的訪問量、提高訪問速度來提高生成速 度。
[0006] 另一個是利用GPU(Gra地ics Processing Unit)進行直方圖生成。GPU作為一種專 為計算密集型、高度并行化應用而設計的高性能計算平臺，其運算能力和存儲器帶寬上相 對于CPU有明顯的優(yōu)勢，通過計算統(tǒng)一設備構架CUDA(Compute Unified Device Architecture),GPU可W在單指令多數(shù)據(jù)SIMD(Single Instruction Multiple Data)編程 模型下發(fā)揮其強大的計算能力，因此，高效的直方圖生成可W依托運些并行處理器平臺?；?于GPU的直方圖生成算法的并行化難點在于如何減少數(shù)據(jù)分布的不規(guī)則所產生的寫沖突。 目前CUDA的軟件開發(fā)包已經支持32位色圖和256色圖兩種基礎的直方圖生成算法，兩者的 實現(xiàn)主要依靠于對直方圖進行多次復制，即為每個線程生成一個私有直方圖副本W(wǎng)減少寫 沖突。

【發(fā)明內容】

[0007] 本發(fā)明的目的是提供一種地理柵格數(shù)據(jù)任選多邊形區(qū)域的直方圖實時生成方法。 利用本發(fā)明為用戶提供服務的流程如圖1所示，用戶根據(jù)需求在客戶端屏幕自由框選任意 多邊形區(qū)域，客戶端將用戶輸入的請求發(fā)送到服務器，服務器經過計算后再將求得的直方 圖結果返回到客戶端，最后由客戶端將結果反饋給用戶，而且為了滿足實時交互的需求，從 用戶輸入完畢到得到直方圖結果，耗時不應超過0.5秒。
[0008] 為實現(xiàn)上述目的，本發(fā)明技術解決方案如下：
[0009] -種地理柵格數(shù)據(jù)任選多邊形區(qū)域的直方圖實時生成方法，包括W下步驟：
[0010] 1)進行預處理。首先將原始柵格數(shù)據(jù)分割成一定尺寸的瓦片（可將瓦片大小取為 256 X 256)，其分割示意圖如圖3所示，圖中陰影部分表示原始柵格數(shù)據(jù)，一個方格代表一個 瓦片。然后分別計算柵格數(shù)據(jù)中各瓦片的直方圖并將其W文件方式存儲到服務器。計算一 塊瓦片的直方圖的過程偽代碼如下：
[0012] 下面結合偽代碼對該過程進行簡要說明。其中tile_sizeX和tile_sizeY為瓦片的 長和寬，tile_data為瓦片中需要統(tǒng)計為直方圖的屬性數(shù)值，X和y表示當前讀取點在瓦片中 的位置坐標，interval對應直方圖的間隔，tile_hist為該瓦片的直方圖結果向量，i表示向 量tile_Mst的第i項。計算一塊瓦片的直方圖的過程為依此遍歷瓦片中的每一個點，對于 其中任意一個點(X，y)，判斷該點對應的瓦片屬性值ti 1 e_data(X，y)與直方圖間隔的關系， 如果存在i X interval < tile_data(x，y)<(i+l) X interval，則為結果直方圖向量tile_ hist的第i項加1。
[0013] 2)判定各瓦片與用戶選定多邊形區(qū)域的關系，分別確定處于選定區(qū)域內部的瓦片 和處于選定區(qū)域邊界的瓦片。如圖4所示，原始柵格數(shù)據(jù)被分成多個瓦片，虛線表示選定區(qū) 域的邊界。那么接下來需要判定哪些瓦片在選定區(qū)域內部哪些瓦片與選定區(qū)域邊界相交。 首先計算能夠將用戶選定多邊形區(qū)域完全包容進去的最小矩形即最小矩形外包框，將其表 示為(別，71，^2,72)，其中^1，71)和^2,72)分別為最小矩形外包框的左下角和右上角的 頂點坐標?；诳蛻舳孙@示的縮放級別可W確定分辨率R，那么客戶端顯示的一個像素點就 代表RXR個柵格單元。從而可利用下取整函數(shù)floor得到縱向填充的范圍跨度是從floor (71/1?巧化1〇〇1'(72/10+1，橫向填充的范圍跨度是從^〇〇1'(^/10至相1〇〇1'^2/10 + 1。然后利 用掃描線填充算法（參見Donald D.Hearn,M.Pauline Baker.Computer Graphics with 0penGL(3th Edition)[M].USA:Prentice 化11，2010.P187-P201.)在分辨率R下對選定區(qū) 域進行填充。假設一個瓦片有NXN個柵格單元，那么為每個瓦片建立一個(N/R) X (N/R)的 填充矩陣flag_tile便可W表示瓦片的所有點。在填充過程中，對于在選定區(qū)域內的點，將 flag_tile中與該點對應的位置置為1，對于選定區(qū)域外的點，flag_tile中與該點對應的位 置置為0。同時記錄各填充矩陣中被置為1的點的個數(shù)。對于一個填充矩陣，如果矩陣中的每 一個點都被置為1，那么可W判斷運個填充矩陣對應的瓦片完全在選定區(qū)域內。那些未全部 被置為1的填充矩陣對應的瓦片處在選定區(qū)域邊界。運些填充矩陣將被用于確定原始柵格 數(shù)據(jù)中的點的位置來計算直方圖數(shù)據(jù)。
[0014] 3)對于步驟2)確定的處于選定區(qū)域內部的瓦片，讀入步驟I)中生成的文件，讀取 對應瓦片的直方圖數(shù)據(jù)并將各瓦片的直方圖數(shù)據(jù)匯總求和得直方圖向量inner_hist，其計 算公式如下：
[0015] inner_hist(i) = 2tiie_xeinner_tiietile_hist_of_tile_x(i)
[0016] 其中inne;r_tile為處于選定區(qū)域內部的瓦片集合，i表示innerjiist的第i項， 1：；[16_11131：_0；1^_1：；[16_義為瓦片1：；[16_義對應的直方圖向量。
[0017] 4)對處于選定區(qū)域邊界的瓦片，采用并行方式讀取原始柵格數(shù)據(jù)并統(tǒng)計直方圖。
[0018] 首先根據(jù)步驟2)確定的處于選定區(qū)域邊界的瓦片，并根據(jù)瓦片對應的填充矩陣確 定瓦片處于選定區(qū)域內的部分。將每個處于選定區(qū)域邊界的瓦片作為一個子任務。如圖5所 示，程序的主進程PO將子任務分配給不同的工作進程，將處在選定區(qū)域邊界的瓦片對應的 填充矩陣flag_tile發(fā)送給不同的工作進程，各工作進程分別使用開源柵格空間數(shù)據(jù)轉換 庫GDAL(參見GDAUGeospatial Data Abstraction Library .http://www.gdal .org/)在原 始柵格數(shù)據(jù)中提取瓦片并根據(jù)瓦片對應的填充矩陣來確定瓦片處于選定區(qū)域內的部分，然 后統(tǒng)計運部分的直方圖。
[0019] 對每一個處在選定區(qū)域邊界的瓦片進行直方圖統(tǒng)計的過程偽代碼如下：
[0021]下面結合偽代碼對該過程進行簡要說明。其中tile_sizeX和tile_sizeY為瓦片的 長和寬，R為分辨率因子，flag_tile為該瓦片對應的填充矩陣，X和Y用于表示flag_tile中 的位置坐標，tile_da化為瓦片中需要統(tǒng)計為直方圖的屬性數(shù)值，X和y表示當前讀取點在瓦 片中的位置坐標，interva 1對應直方圖的間隔，bound_ti 1 e_hist為該瓦片的直方圖結果向 量，i表示向量6〇11扣1_^16_1113*的第1項。對每一個處在選定區(qū)域邊界的瓦片進行直方圖統(tǒng) 計的步驟為:遍歷該瓦片對應填充矩陣flag_tile中的每一個點，對于其中任意一個點(X， Y)如果flag_tile(X，Y) = l，則表示與該點對應的瓦片上的點處在選定區(qū)域內部，即瓦片上 橫向范圍從X X R到(X+1) X R-1縱向范圍從Y X R到(Y+1) X R-1的所有點均處于選定區(qū)域內 部。然后，遍歷瓦片上運部分區(qū)域內的點，對于任意一個點(x，y)，判斷該點對應的瓦片屬性 值tile_data(x，y)與直方圖間隔的關系，如果存在i X interval < tile_data(x，y)<(i+l) X interval，則為結果直方圖向量bound_tile_hist的第i項加1。
[0022] 然后將各工作進程求得的結果進行求和得直方圖向量boundjiist，計算公式如 下，其中bound_tile為處于選定區(qū)域邊界的瓦片集合，i表示第i項，bound_tile_hist_of_ tile_x為邊界處瓦片tile_x對應的直方圖向量。
[0023] bound_h i S t (i) = 2 ti ie_xeb〇imd_ti iebound_t i 1 e_h i S t_o f _t i 1 e_x (i)
[0024] 5)將步驟3)得到的直方圖向量inne;r_hist和步驟4)中得到的直方圖向量bound_ hist求和得到選定區(qū)域的直方圖histogram,計算公式如下，其中i表示向量histogram的第 i項：
[0025] histogram(i)=inner_hist(i)+bound_hist(i)
[0026] 本發(fā)明的有益效果是：
[0027] 通過對原始柵格數(shù)據(jù)進行一定的預處理，減少了計算直方圖時對原始柵格數(shù)據(jù)的 訪問量，并結合客戶端屏幕分辨率對生成方法進行了優(yōu)化，同時利用高性能集群進行并行 化，最終可W滿足實時獲取用戶選定任意多邊形區(qū)域直方圖的需求。
【附圖說明】
[0028] 圖1是本發(fā)明為用戶提供服務的流程示意圖
[0029] 圖2是本發(fā)明的流程示意圖
[0030] 圖3是柵格數(shù)據(jù)分割為瓦片示意圖
[0031] 圖4是選定區(qū)域內瓦片的判定示意圖
[0032] 圖5是使用GDAL在并行I/O模式下處理地理空間柵格數(shù)據(jù)的示意圖
[0033] 圖6是本發(fā)明方法與直接從原始柵格進行統(tǒng)計的直方圖生成方法的任務耗時對比 (兩種方法均采用256個工作進程進行并行化）
【具體實施方式】
[0034] 下面結合附圖和具體實例對本發(fā)明做進一步描述。
[0035] 圖1為本發(fā)明為用戶提供服務的流程示意圖。用戶根據(jù)需求在客戶端屏幕自由框 選任意多邊形區(qū)域，客戶端將用戶輸入的請求發(fā)送到服務器，服務器經過計算后再將直方 圖結果返回到客戶端，最后由客戶端將直方圖結果展示給用戶。
[