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

【發(fā)明內(nèi)容】

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