Skip to content
Permalink
Browse files

[opencl] Fix nodata value in slope, aspect and hillshade

  • Loading branch information
elpaso committed Aug 8, 2018
1 parent 164bcc2 commit a6d5d479fa5e3b8f239ae2e26010bfebd0579586
@@ -4,42 +4,37 @@ __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
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]
);


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

float derY = calcFirstDer( scanLine3[i], scanLine2[i], scanLine1[i],
scanLine3[i+1], scanLine2[i+1], scanLine1[i+1],
scanLine3[i+2], scanLine2[i+2], scanLine1[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]);

if ( derX == rasterParams[1] || derY == rasterParams[1] ||
( derX == 0.0f && derY == 0.0f) )
{
resultLine[i] = rasterParams[1];
}
else
resultLine[i] = 90.0 - aspect;
{
resultLine[i] = 180.0f + atan2pi( derX, derY ) * 180.0f;
}
}
}
@@ -1,109 +1,56 @@
#include "calcfirstder.cl"

// Calculate the first derivative from a 3x3 cell matrix
float calcFirstDer( float x11, float x21, float x31, float x12, float x22, float x32, float x13, float x23, float x33,
float inputNodataValue, float zFactor, float mCellSize )
{
//the basic formula would be simple, but we need to test for nodata values...
//X: return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * cellSizeX));
//Y: return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * cellSizeY));


if ( x11 == inputNodataValue )
{
x11 = x22;
}
if ( x12 == inputNodataValue )
{
x12 = x22;
}
if ( x13 == inputNodataValue )
{
x13 = x22;
}
if ( x21 == inputNodataValue )
{
x21 = x22;
}
if ( x23 == inputNodataValue )
{
x23 = x22;
}
if ( x31 == inputNodataValue )
{
x31 = x22;
}
if ( x32 == inputNodataValue )
{
x32 = x22;
}
if ( x33 == inputNodataValue )
{
x33 = x22;
}
return ( ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * mCellSize )) * zFactor;
}


#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

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

) {
__global float *scanLine2,
__global float *scanLine3,
__global float *resultLine,
__global float *rasterParams // mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY, zenith_rad, azimuth_rad
)
{

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

float res;
if ( scanLine2[i+1] == rasterParams[1] )
if ( scanLine2[i+1] == rasterParams[0] )
{
res = rasterParams[2];
resultLine[i] = rasterParams[1];
}
else
{

// 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[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[2], rasterParams[4]
);
float derX = calcFirstDer( scanLine1[i], scanLine1[i+1], scanLine1[i+2],
scanLine2[i], scanLine2[i+1], scanLine2[i+2],
scanLine3[i], scanLine3[i+1], scanLine3[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3] );

float derY = calcFirstDer( scanLine3[i], scanLine2[i], scanLine1[i],
scanLine3[i+1], scanLine2[i+1], scanLine1[i+1],
scanLine3[i+2], scanLine2[i+2], scanLine1[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]);

if ( derX == rasterParams[1] || derY == rasterParams[1] )
{
res = rasterParams[2];
resultLine[i] = rasterParams[1];
}
else
{
float aspect = atan2( derX, - derY );
// aspect = aspect * M_PI * 255; // for display

float slope = sqrt( derX * derX + derY * derY );
slope = atan( slope );
// res = slope * M_PI * 255; for display
float slope_rad = sqrt( derX * derX + derY * derY );
slope_rad = atan( slope_rad );
float aspect_rad;
if ( derX == 0.0f && derY == 0.0f)
aspect_rad = rasterParams[7] / 2.0f;
else
aspect_rad = M_PI + atan2( derX, derY );

// Final
res = 255.0f * ( rasterParams[6] * cos( slope )
+ rasterParams[7] * sin( slope )
* cos( rasterParams[5] - aspect ) );
resultLine[i] = max(0.0f, 255.0f * ( ( rasterParams[5] * cos( slope_rad ) ) +
( rasterParams[6] * sin( slope_rad ) *
cos( rasterParams[7] - aspect_rad ) ) ) );

}
}

// Opacity
res = res * rasterParams[7];

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

}
@@ -9,13 +9,13 @@ __kernel void processNineCellWindow( __global float *scanLine1,
const int i = get_global_id(0);

float x11 = scanLine1[i];
float x12 = scanLine1[i+1];
float x13 = scanLine1[i+2];
float x21 = scanLine2[i];
float x21 = scanLine1[i+1];
float x31 = scanLine1[i+2];
float x12 = scanLine2[i];
float x22 = scanLine2[i+1];
float x23 = scanLine2[i+2];
float x31 = scanLine3[i];
float x32 = scanLine3[i+1];
float x32 = scanLine2[i+2];
float x13 = scanLine3[i];
float x23 = scanLine3[i+1];
float x33 = scanLine3[i+2];

if ( x22 == rasterParams[0] )
@@ -4,36 +4,39 @@ __kernel void processNineCellWindow( __global float *scanLine1,
__global float *scanLine2,
__global float *scanLine3,
__global float *resultLine,
__global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY ]

) {
__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]
);


if ( derX == rasterParams[1] || derY == rasterParams[1] )
if ( scanLine2[i+1] == rasterParams[0] )
{
resultLine[i] = rasterParams[1];
resultLine[i] = rasterParams[1];
}
else
{
float res = sqrt( derX * derX + derY * derY );
res = atanpi( res );
resultLine[i] = res * 180.0;
float derX = calcFirstDer( scanLine1[i], scanLine1[i+1], scanLine1[i+2],
scanLine2[i], scanLine2[i+1], scanLine2[i+2],
scanLine3[i], scanLine3[i+1], scanLine3[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3] );

float derY = calcFirstDer( scanLine3[i], scanLine2[i], scanLine1[i],
scanLine3[i+1], scanLine2[i+1], scanLine1[i+1],
scanLine3[i+2], scanLine2[i+2], scanLine1[i+2],
rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]);


if ( derX == rasterParams[1] || derY == rasterParams[1] )
{
resultLine[i] = rasterParams[1];
}
else
{
float res = sqrt( derX * derX + derY * derY );
res = atanpi( res );
resultLine[i] = res * 180.0f;
}
}
}

0 comments on commit a6d5d47

Please sign in to comment.
You can’t perform that action at this time.