DEM山體陰影原理以及算法具體解釋
山體陰影原理以及算法具體解釋
山體陰影基本原理:
山體陰影是假想一個光源在某個方向和某個太陽高度的模擬下。用過臨近像元的計算來生成一副0-255的灰度圖。
一、山體陰影的主要參數:
1、 太陽光線的入射角度:這個角度的量算起點是正北方向,依照順時針的方向,角度的範圍是0到360度。例如以下圖所看到的,默認的角度是315度,西北方向,例如以下圖所看到的:
2、 太陽高度角:太陽高度角也簡稱太陽高度。是太陽光線和當地地平面之間的夾角,範圍是0-90度,默認的太陽高度是45度,例如以下圖所看到的:
二、山體陰影計算方法
山體陰影的計算公式例如以下
(1) Hillshade = 255.0 * ((cos(Zenith_rad
(sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad- Aspect_rad)))
假設Hillshade < 0, 則設Hillshade=0.
當中,Zenith_rad是太陽天頂角的的弧度數。Slope_rad是某一點的坡度弧度數,Azimuth_rad是指太陽光線方向角的弧度數。Aspect_rad是某一點的坡向弧度數。
計算山體陰影的照明光源的角度默認是太陽高度角,可是真正計算時。須要用到太陽天頂角,太陽天頂角的計算方法是90°-太陽高度角。所以有例如以下計算公式:
計算太陽天頂角弧度數
(2) Zenith_deg = 90 - Altitude
轉換為弧度數:
(3) Zenith_rad = Zenith * pi / 180.0
計算照明的方向:
照明的方向角是指定的角度數,山體陰影的計算公式須要弧度數。
首先,須要將地理上的指南針方向轉換為數學上的向右的方向。即向右為起算的方向;其次。須要將角度轉換為弧度。
轉為數學上的方向角:
(4) Azimuth_math = 360.0 - Azimuth + 90
註意假設 Azimuth_math >=360.0, 那麽:
(5) Azimuth_math = Azimuth_math - 360.0
轉換為弧度:
(6) Azimuth_rad= Azimuth_math * pi / 180.0
計算坡度和坡向
坡度和坡向是利用一個3*3的窗體在輸入影像中訪問每一個像素,9個像素從左到右、從上到下分別用a-i表示,如圖所看到的:
a | b | c |
d | e | f |
g | h | i |
E像元X方向上的變化率採用例如以下的算法:
(7) [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / (8 * cellsize)
E像元Y方向上的變化率採用例如以下的算法:
(8) [dz/dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * cellsize)
坡度的弧度計算公式,考慮了Z因子(協調Z方向的單位和Xy平面上單位的一個系數):
(9) Slope_rad = atan (z_factor * sqrt ([dz/dx]2 + [dz/dy]2))
坡向通過以下的方法進行計算:
If [dz/dx] is non-zero:
Aspect_rad= atan2 ([dz/dy], -[dz/dx])
if Aspect_rad< 0 then
Aspect_rad= 2 * pi + Aspect_rad
If [dz/dx] iszero:
if [dz/dy] >0 then
Aspect_rad= pi / 2
else if [dz/dy]< 0 then
Aspect_rad= 2 * pi - pi / 2
else
Aspect_rad = Aspect_rad
山體陰影計算演示樣例:出自arcgis10.0幫助文檔。
最後,奉上OpenCL實現的代碼:
__kernel void hillshade_kernel( __global const float *pSrcData, __global float *pDestData,const int nWidth,const int nHeight , struct HillshadeOption hillOption) { int j = get_global_id(0); int i = get_global_id(1); if (j >= nWidth || i >= nHeight) return; int nTopTmp = i-1; int nBottomTmp = i+1; int nLeftTep = j-1; int nRightTmp = j+1; //處理邊界情況 if (0 == i) { nTopTmp = i; } if (0 == j) { nLeftTep = j; } if (i == nHeight-1) { nBottomTmp = i; } if (j == nWidth-1) { nRightTmp = j; } __local float afRectData[9]; afRectData[0] = pSrcData[nTopTmp*nWidth+nLeftTep]; afRectData[1] = pSrcData[nTopTmp*nWidth+j]; afRectData[2] = pSrcData[nTopTmp*nWidth+nRightTmp]; afRectData[3] = pSrcData[i*nWidth+nLeftTep]; afRectData[4] = pSrcData[i*nWidth+j]; afRectData[5] = pSrcData[i*nWidth+nRightTmp]; afRectData[6] = pSrcData[nBottomTmp*nWidth+nLeftTep]; afRectData[7] = pSrcData[nBottomTmp*nWidth+j]; afRectData[8] = pSrcData[nBottomTmp*nWidth+nRightTmp]; const float degreesToRadians = (M_PI_F / 180); float dx = ((afRectData[2]+ afRectData[5]*2 + afRectData[8]) - (afRectData[0] + afRectData[3]*2 + afRectData[6])) / (8 * hillOption.dbEwres); float dy = ((afRectData[6] + afRectData[7]*2 + afRectData[8]) - (afRectData[0]+ afRectData[1]*2 + afRectData[2])) / (8 * hillOption.dbNsres); //計算坡度(弧度) float key = sqrt(dx *dx + dy * dy); float dfSlope = atan( hillOption.dfzScaleFactor * key); //計算坡向(弧度) float dfAspect = 0; if (dx != 0) { dfAspect = atan2(dy,-dx); if (dfAspect < 0) { dfAspect += 2* M_PI_F; } } if (dx == 0) { if (dy > 0.0f) { dfAspect = M_PI_F / 2; } else if (dy < 0.0f) dfAspect = M_PI_F + M_PI_F / 2; } //將太陽高度和太陽光線角度轉換為要求的格式 float dfZenithDeg = hillOption.dfAltitude; float dfAzimuthRad = hillOption.dfAzimuth; //最後計算山體陰影值 float dfHillshade = 255 * (cos(dfZenithDeg)*cos(dfSlope) + sin(dfZenithDeg)*sin(dfSlope)*cos(dfAzimuthRad-dfAspect)); if (dfHillshade < 0) { dfHillshade = 0; } if (dfHillshade >= 255) { dfHillshade = 255; } pDestData[i*nWidth+j] = (int)(dfHillshade+1/2); }
DEM山體陰影原理以及算法具體解釋