Skip to content

Commit

Permalink
[opencl] Hillshade renderer
Browse files Browse the repository at this point in the history
  • Loading branch information
elpaso committed Aug 8, 2018
1 parent 350829e commit d5dd50c
Show file tree
Hide file tree
Showing 4 changed files with 158 additions and 31 deletions.
53 changes: 53 additions & 0 deletions resources/opencl_programs/aspect_renderer.cl
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
#include "calcfirstder.cl"

// Aspect renderer for QGIS

__kernel void processNineCellWindow( __global float *scanLine1,
__global float *scanLine2,
__global float *scanLine3,
__global uchar4 *resultLine,
__global float *rasterParams
) {

// Get the index of the current element
const int i = get_global_id(0);

// Do the operation
//return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * mCellSizeX))
float derX = calcFirstDer( scanLine1[i], scanLine2[i], scanLine3[i],
scanLine1[i+1], scanLine2[i+1], scanLine3[i+1],
scanLine1[i+2], scanLine2[i+2], scanLine3[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3]
);
//return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * mCellSizeY));
float derY = calcFirstDer( scanLine1[i+2], scanLine1[i+1], scanLine1[i],
scanLine2[i+2], scanLine2[i+1], scanLine2[i],
scanLine3[i+2], scanLine3[i+1], scanLine3[i],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]
);


float res;
if ( derX == rasterParams[1] || derY == rasterParams[1] ||
( derX == 0.0f && derY == 0.0f) )
{
res = rasterParams[1];
}
else
{
// 180.0 / M_PI = 57.29577951308232
float aspect = atan2( derX, derY ) * 57.29577951308232;
if ( aspect < 0 )
res = 90.0f - aspect;
else if (aspect > 90.0f)
// 360 + 90 = 450
res = 450.0f - aspect;
else
res = 90.0 - aspect;
}

res = res / 360 * 255;

resultLine[i] = (uchar4)(res, res, res, 255);
}

88 changes: 58 additions & 30 deletions resources/opencl_programs/hillshade.cl
Original file line number Diff line number Diff line change
@@ -1,35 +1,63 @@
#include "calcfirstder.cl"

__constant sampler_t sampler =
CLK_NORMALIZED_COORDS_FALSE
| CLK_ADDRESS_CLAMP_TO_EDGE
| CLK_FILTER_NEAREST;

__kernel void processNineCellWindow(
__read_only image2d_t inputImage,
image2d_t outputImage,
__global float *rasterParams
) {


// 0=width, 1=height
int2 currentPosition = (int2)(get_global_id(0), get_global_id(1));
float4 currentPixel = (ufloat4)(0.0f);
float4 calculatedPixel = (float4)(1.0f);

currentPixel = read_imageuf(inputImage, sampler, currentPosition);
calculatedPixel = currentPixel;
write_imageuf(outputImage, currentPosition, calculatedPixel);
}


const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
__kernel void processNineCellWindow( __global float *scanLine1,
__global float *scanLine2,
__global float *scanLine3,
__global uchar4 *resultLine, // This is an image BGRA !
__global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY,
// azimuthRad, cosZenithRad, sinZenithRad ]

) {

// Get the index of the current element
const int i = get_global_id(0);

// Do the operation
//return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * mCellSizeX))
float derX = calcFirstDer( scanLine1[i], scanLine2[i], scanLine3[i],
scanLine1[i+1], scanLine2[i+1], scanLine3[i+1],
scanLine1[i+2], scanLine2[i+2], scanLine3[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3]
);
//return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * mCellSizeY));
float derY = calcFirstDer( scanLine1[i+2], scanLine1[i+1], scanLine1[i],
scanLine2[i+2], scanLine2[i+1], scanLine2[i],
scanLine3[i+2], scanLine3[i+1], scanLine3[i],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]
);


float res;
if ( derX == rasterParams[1] || derY == rasterParams[1] )
{
res = rasterParams[1];
}
else
{
// 180.0 / M_PI = 57.29577951308232
float aspect = atan2( derX, derY ) * 57.29577951308232;
if ( aspect < 0 )
aspect = 90.0f - aspect;
else if (aspect > 90.0f)
// 360 + 90 = 450
aspect = 450.0f - aspect;
else
aspect = 90.0 - aspect;
// aspect = aspect / 360 * 255; // for display
aspect = radians( aspect );

float slope = sqrt( derX * derX + derY * derY );
slope = atan( slope );
// res = slope * M_PI * 255; for display

// Final
res = 255.0f * ( rasterParams[6] * cos( slope )
+ rasterParams[7] * sin( slope )
* cos( rasterParams[5] - aspect ) );

}

resultLine[i] = (uchar4)(res, res, res, 255);

void kernel copy(__read_only image2d_t in, image2d_t out)
{
int x = get_global_id(0);
int y = get_global_id(1);
int2 pos = (int2)(x, y);
uint4 pixel = read_imageui(in, smp, pos);
write_imageui(out, pos, pixel);
}
3 changes: 2 additions & 1 deletion resources/opencl_programs/slope.cl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@ __kernel void processNineCellWindow( __global float *scanLine1,
__global float *scanLine2,
__global float *scanLine3,
__global float *resultLine,
__global float *rasterParams
__global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY ]

) {

// Get the index of the current element
Expand Down
45 changes: 45 additions & 0 deletions resources/opencl_programs/slope_renderer.cl
Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
#include "calcfirstder.cl"

// Slope renderer for QGIS

__kernel void processNineCellWindow( __global float *scanLine1,
__global float *scanLine2,
__global float *scanLine3,
__global uchar4 *resultLine, // This is an image BGRA !
__global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY ]

) {

// Get the index of the current element
const int i = get_global_id(0);

// Do the operation
//return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * mCellSizeX))
float derX = calcFirstDer( scanLine1[i], scanLine2[i], scanLine3[i],
scanLine1[i+1], scanLine2[i+1], scanLine3[i+1],
scanLine1[i+2], scanLine2[i+2], scanLine3[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3]
);
//return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * mCellSizeY));
float derY = calcFirstDer( scanLine1[i+2], scanLine1[i+1], scanLine1[i],
scanLine2[i+2], scanLine2[i+1], scanLine2[i],
scanLine3[i+2], scanLine3[i+1], scanLine3[i],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]
);


float res;
if ( derX == rasterParams[1] || derY == rasterParams[1] )
{
res = rasterParams[1];
}
else
{
float slope = sqrt( derX * derX + derY * derY );
slope = atan( slope );
res = slope * 255;
}

resultLine[i] = (uchar4)(res, res, res, 255);

}

0 comments on commit d5dd50c

Please sign in to comment.