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山體陰影原理以及算法具體解釋