[opencl] Hillshade renderer

qgis · Aug 8, 2018 · d5dd50c · d5dd50c
1 parent 350829e
commit d5dd50c
Show file tree

Hide file tree

Showing 4 changed files with 158 additions and 31 deletions.
diff --git a/resources/opencl_programs/aspect_renderer.cl b/resources/opencl_programs/aspect_renderer.cl
@@ -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);
+}
+
diff --git a/resources/opencl_programs/hillshade.cl b/resources/opencl_programs/hillshade.cl
@@ -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);
 }
diff --git a/resources/opencl_programs/slope.cl b/resources/opencl_programs/slope.cl
@@ -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

diff --git a/resources/opencl_programs/slope_renderer.cl b/resources/opencl_programs/slope_renderer.cl
@@ -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);
+
+}