Dem Hillshade principle and algorithm detailed

Hillshade principle and algorithmic explanation

Hillshade Fundamentals:

A hillshade is a hypothetical light source that simulates a certain direction and the height of a certain sun, using the calculations of neighboring cells to generate a 0-255 grayscale image.

First, the main parameters of the mountain Shadow:

1, the angle of incidence of the sun rays: this angle of the acreage starting point is North direction, in accordance with the clockwise direction, the angle of the range is 0-360 degrees, as shown, the default angle is 315 degrees, northwest direction, as shown:

2, the Sun height angle: The Sun height angle also referred to as the sun height. is the angle between the sun light and the local ground plane, the range is 0-90 degrees, the default sun height is 45 degrees, as shown in:

Second, the method of calculating the mountain shadow

The Hillshade formula is calculated as follows

(1) Hillshade = 255.0 * ((cos (zenith_rad) * cos (slope_rad)) +

(Sin (zenith_rad) * sin (slope_rad) * cos (azimuth_rad- Aspect_rad)))

If Hillshade < 0, then set hillshade=0.

Where Zenith_rad is the number of radians of the Sun's zenith, Slope_rad is the number of slope radians of a point, Azimuth_rad is the number of radians in the direction of the sun's rays, and Aspect_rad is the number of slope radians at a point.

The angle of the illumination light source that calculates the hillshade is the sun height angle by default, but the sun zenith angle is used when the real calculation is made, and the solar Zenith angle is calculated by 90°-the sun height angle. Therefore, the following formula is calculated:

calculate the number of radians of the sun Zenith Angle :

zenith_deg Altitude

Convert to radians number:

Zenith_rad Zenith * pi/180.0

To calculate the direction of illumination:

The direction angle of the illumination is the number of angles specified, and the calculation formula for the Hillshade requires radians. First, you need to convert the direction of the compass from the geography to the mathematically right direction, which is the direction to the right, and second, you need to convert the angle to radians.

To a mathematically oriented angle:

Azimuth_math Azimuth + 90

Note If azimuth_math >=360.0, then:

Azimuth_math Azimuth_math -360.0

Convert to radians:

Azimuth_rad Azimuth_math * pi/180.0

Calculate slope and aspect

Slope and aspect is the use of a 3*3 window to access each pixel in the input image, 9 pixels from left to right, from top to bottom, respectively, with the A-i representation:

 
   
  

  
    
   
 
    
 
     
A
 
     
B
 
     
C
 
    
 
    
 
     
D
 
     
E
 
     
F
 
    
 
    
 
     
G
 
     
H
 
     
I
 
    
 
   
 
   
  

The rate of change in the x direction of the E cell is based on the following algorithm:

(7) [DZ/DX] = ((c + 2fi)-(a + 2dgcellsize)

The rate of change in the y direction of e cells is based on the following algorithm:

(8) [Dz/dy] = ((g + 2hi)-(a + 2bccellsize)

The arc of the slope is calculated by considering the Z-factor (a factor that coordinates the units in the Z-direction and the units on the XY plane):

Slope_rad = atan (z_factor

The aspect is calculated by the following method:

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

Hillshade Calculation Example: From the arcgis10.0 help document.

Finally, the code implemented in OpenCL:

__kernel void Hillshade_kernel (__global const float *psrcdata, __global float *pdestdata,const int nwidth,const int NHeig HT, struct hillshadeoption hilloption) {int j = get_global_id (0); int i = get_global_id (1); if (j >= nwidth | | I >= nH Eight) Return;int ntoptmp = I-1;int nbottomtmp = i+1;int Nlefttep = J-1;int nrighttmp = j+1;//processing boundary condition 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 = ((afrectd ATA[6] + afrectdata[7]*2 + afrectdata[8])-(afrectdata[0]+ afrectdata[1]*2 + afrectdata[2])/(8 * hilloption.dbnsres);// Calculate slope (radians) float key = sqrt (dx *dx + dy * dy), float dfslope = atan (Hilloption.dfzscalefactor * key);//calculate aspect (radians) float Dfaspec t = 0;if (dx! = 0) {dfaspect = atan2 (DY,-DX), if (Dfaspect < 0) {dfaspect + = 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;} Convert the sun height and the sun light angle to the required format float dfzenithdeg = hilloption.dfaltitude;float Dfazimuthrad = hilloption.dfazimuth;// The last computed Hillshade value 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 Hillshade principle and algorithm detailed