Skip to content

Commit 41e7e71

Browse files
elpasoelpaso
committed
[opencl] Hillshade fix alpha and nodata
1 parent eaa7982 commit 41e7e71

File tree

2 files changed

+156
-44
lines changed

2 files changed

+156
-44
lines changed

‎resources/opencl_programs/hillshade_renderer.cl‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
__kernel void processNineCellWindow( __global float *scanLine1,
2+
__global float *scanLine2,
3+
__global float *scanLine3,
4+
__global uchar4 *resultLine, // This is an image!
5+
__global float *rasterParams
6+
7+
)
8+
{
9+
10+
// Get the index of the current element
11+
const int i = get_global_id(0);
12+
13+
float x11 = scanLine1[i];
14+
float x12 = scanLine1[i+1];
15+
float x13 = scanLine1[i+2];
16+
float x21 = scanLine2[i];
17+
float x22 = scanLine2[i+1];
18+
float x23 = scanLine2[i+2];
19+
float x31 = scanLine3[i];
20+
float x32 = scanLine3[i+1];
21+
float x33 = scanLine3[i+2];
22+
23+
if ( x22 == rasterParams[0] )
24+
{
25+
float alpha = rasterParams[8] * rasterParams[16];
26+
resultLine[i] = (uchar4)(rasterParams[13] * alpha,
27+
rasterParams[14] * alpha,
28+
rasterParams[15] * alpha, 255 * alpha);
29+
}
30+
else
31+
{
32+
if ( x11 == rasterParams[0] ) x11 = x22;
33+
if ( x12 == rasterParams[0] ) x12 = x22;
34+
if ( x13 == rasterParams[0] ) x13 = x22;
35+
if ( x21 == rasterParams[0] ) x21 = x22;
36+
// Note: skip central cell x22
37+
if ( x23 == rasterParams[0] ) x23 = x22;
38+
if ( x31 == rasterParams[0] ) x31 = x22;
39+
if ( x32 == rasterParams[0] ) x32 = x22;
40+
if ( x33 == rasterParams[0] ) x33 = x22;
41+
42+
float derX = ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * rasterParams[3] );
43+
float derY = ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -rasterParams[4]);
44+
45+
if ( derX == rasterParams[0] ||
46+
derX == rasterParams[0] )
47+
{
48+
float alpha = rasterParams[5] * rasterParams[16];
49+
resultLine[i] = (uchar4)(rasterParams[13] * alpha,
50+
rasterParams[14] * alpha,
51+
rasterParams[15] * alpha, 255 * alpha);
52+
}
53+
else
54+
{
55+
float res;
56+
if ( rasterParams[17] ) // Multi directional?
57+
{
58+
// Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
59+
// Fast formula from GDAL DEM
60+
const float xx = derX * derX;
61+
const float yy = derY * derY;
62+
const float xx_plus_yy = xx + yy;
63+
// Flat?
64+
if ( xx_plus_yy == 0.0f )
65+
{
66+
res = clamp( 1.0f + rasterParams[9], 0.0f, 255.0f );
67+
}
68+
else
69+
{
70+
// ... then the shade value from different azimuth
71+
float val225_mul_127 = rasterParams[10] +
72+
( derX - derY ) * rasterParams[11];
73+
val225_mul_127 = ( val225_mul_127 <= 0.0f ) ? 0.0f : val225_mul_127;
74+
float val270_mul_127 = rasterParams[10] -
75+
derX * rasterParams[12];
76+
val270_mul_127 = ( val270_mul_127 <= 0.0f ) ? 0.0f : val270_mul_127;
77+
float val315_mul_127 = rasterParams[10] +
78+
( derX + derY ) * rasterParams[11];
79+
val315_mul_127 = ( val315_mul_127 <= 0.0f ) ? 0.0f : val315_mul_127;
80+
float val360_mul_127 = rasterParams[10] -
81+
derY * rasterParams[12];
82+
val360_mul_127 = ( val360_mul_127 <= 0.0f ) ? 0.0f : val360_mul_127;
83+
// ... then the weighted shading
84+
const float weight_225 = 0.5f * xx_plus_yy - derX * derY;
85+
const float weight_270 = xx;
86+
const float weight_315 = xx_plus_yy - weight_225;
87+
const float weight_360 = yy;
88+
const float cang_mul_127 = (
89+
( weight_225 * val225_mul_127 +
90+
weight_270 * val270_mul_127 +
91+
weight_315 * val315_mul_127 +
92+
weight_360 * val360_mul_127 ) / xx_plus_yy ) /
93+
( 1 + rasterParams[8] * xx_plus_yy );
94+
res = clamp( 1.0f + cang_mul_127, 0.0f, 255.0f );
95+
}
96+
}
97+
else
98+
{
99+
res = ( rasterParams[9] -
100+
( derY * rasterParams[6] -
101+
derX * rasterParams[7] )) /
102+
sqrt( 1 + rasterParams[8] *
103+
( derX * derX + derY * derY ) );
104+
res = res <= 0.0f ? 1.0f : 1.0f + res;
105+
}
106+
res = res * rasterParams[5];
107+
resultLine[i] = (uchar4)(res, res, res, 255 * rasterParams[5]);
108+
}
109+
}
110+
}

‎resources/opencl_programs/hillshade_renderer_multidirectional.cl‎

Lines changed: 46 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,34 @@
11
__kernel void processNineCellWindow( __global float *scanLine1,
2-
__global float *scanLine2,
3-
__global float *scanLine3,
4-
__global uchar4 *resultLine, // This is an image!
5-
__global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY,
6-
// cos_az_mul_cos_alt_mul_z_mul_254, sin_az_mul_cos_alt_mul_z_mul_254, square_z, sin_altRadians_mul_254]
2+
__global float *scanLine2,
3+
__global float *scanLine3,
4+
__global uchar4 *resultLine, // This is an image!
5+
__global float *rasterParams
76

8-
) {
7+
)
8+
{
99

10-
// Get the index of the current element
11-
const int i = get_global_id(0);
10+
// Get the index of the current element
11+
const int i = get_global_id(0);
1212

13-
float x11 = scanLine1[i];
14-
float x12 = scanLine1[i+1];
15-
float x13 = scanLine1[i+2];
16-
float x21 = scanLine2[i];
17-
float x22 = scanLine2[i+1];
18-
float x23 = scanLine2[i+2];
19-
float x31 = scanLine3[i];
20-
float x32 = scanLine3[i+1];
21-
float x33 = scanLine3[i+2];
13+
float x11 = scanLine1[i];
14+
float x12 = scanLine1[i+1];
15+
float x13 = scanLine1[i+2];
16+
float x21 = scanLine2[i];
17+
float x22 = scanLine2[i+1];
18+
float x23 = scanLine2[i+2];
19+
float x31 = scanLine3[i];
20+
float x32 = scanLine3[i+1];
21+
float x33 = scanLine3[i+2];
2222

23-
float res;
24-
if ( x22 == rasterParams[1] )
25-
{
26-
res = rasterParams[2];
27-
}
28-
else
29-
{
23+
if ( x22 == rasterParams[0] )
24+
{
25+
float alpha = rasterParams[5] * rasterParams[16];
26+
resultLine[i] = (uchar4)(rasterParams[13] * alpha,
27+
rasterParams[14] * alpha,
28+
rasterParams[15] * alpha, 255 * alpha);
29+
}
30+
else
31+
{
3032
if ( x11 == rasterParams[0] ) x11 = x22;
3133
if ( x12 == rasterParams[0] ) x12 = x22;
3234
if ( x13 == rasterParams[0] ) x13 = x22;
@@ -40,36 +42,41 @@ __kernel void processNineCellWindow( __global float *scanLine1,
4042
float derX = ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * rasterParams[3] );
4143
float derY = ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -rasterParams[4]);
4244

43-
if ( derX == rasterParams[1] || derY == rasterParams[1] )
45+
if ( derX == rasterParams[0] ||
46+
derX == rasterParams[0] )
4447
{
45-
res = rasterParams[2];
48+
float alpha = rasterParams[5] * rasterParams[16];
49+
resultLine[i] = (uchar4)(rasterParams[13] * alpha,
50+
rasterParams[14] * alpha,
51+
rasterParams[15] * alpha, 255 * alpha);
4652
}
4753
else
4854
{
55+
float res;
4956
// Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
5057
// Fast formula from GDAL DEM
5158
const float xx = derX * derX;
5259
const float yy = derY * derY;
5360
const float xx_plus_yy = xx + yy;
54-
// Flat? -> white
61+
// Flat?
5562
if ( xx_plus_yy == 0.0f )
5663
{
57-
res = clamp( 1.0f + rasterParams[12], 0.0f, 255.0f );
64+
res = clamp( 1.0f + rasterParams[9], 0.0f, 255.0f );
5865
}
5966
else
6067
{
6168
// ... then the shade value from different azimuth
62-
float val225_mul_127 = rasterParams[13] +
63-
( derX - derY ) * rasterParams[14];
69+
float val225_mul_127 = rasterParams[10] +
70+
( derX - derY ) * rasterParams[11];
6471
val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
65-
float val270_mul_127 = rasterParams[13] -
66-
derX * rasterParams[15];
72+
float val270_mul_127 = rasterParams[10] -
73+
derX * rasterParams[12];
6774
val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
68-
float val315_mul_127 = rasterParams[13] +
69-
( derX + derY ) * rasterParams[14];
75+
float val315_mul_127 = rasterParams[10] +
76+
( derX + derY ) * rasterParams[11];
7077
val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
71-
float val360_mul_127 = rasterParams[13] -
72-
derY * rasterParams[15];
78+
float val360_mul_127 = rasterParams[10] -
79+
derY * rasterParams[12];
7380
val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
7481
// ... then the weighted shading
7582
const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
@@ -81,16 +88,11 @@ __kernel void processNineCellWindow( __global float *scanLine1,
8188
weight_270 * val270_mul_127 +
8289
weight_315 * val315_mul_127 +
8390
weight_360 * val360_mul_127 ) / xx_plus_yy ) /
84-
( 1 + rasterParams[11] * xx_plus_yy );
85-
91+
( 1 + rasterParams[8] * xx_plus_yy );
8692
res = clamp( 1.0f + cang_mul_127, 0.0f, 255.0f );
93+
res = res * rasterParams[5];
94+
resultLine[i] = (uchar4)(res, res, res, 255 * rasterParams[5]);
8795
}
8896
}
8997
}
90-
91-
// Opacity
92-
res = res * rasterParams[7];
93-
94-
resultLine[i] = (uchar4)(res, res, res, 255 * rasterParams[7]);
95-
9698
}

0 commit comments

Comments
 (0)