  * Code adapted from https://github.com/andmax/gpufilter
  */
 
+#include <thrust/device_vector.h>
+
+texture<float, 2, cudaReadModeElementType> textureReadDepth;
+
+inline unsigned int divUp(unsigned int a, unsigned int b) {
+  return (a % b != 0) ? (a / b + 1) : (a / b);
+}
+
+#define WS 32 // Warp size (defines b x b block size where b = WS)
+#define HWS 16 // Half Warp Size
+#define DW 8 // Default number of warps (computational block height)
+#define CHW 7 // Carry-heavy number of warps (computational block height for some kernels)
+#define OW 6 // Optimized number of warps (computational block height for some kernels)
+#define DNB 6 // Default number of blocks per SM (minimum blocks per SM launch bounds)
+#define ONB 5 // Optimized number of blocks per SM (minimum blocks per SM for some kernels)
+#define MTS 192 // Maximum number of threads per block with 8 blocks per SM
+#define MBO 8 // Maximum number of blocks per SM using optimize or maximum warps
+#define CHB 7 // Carry-heavy number of blocks per SM using default number of warps
+#define MW 6 // Maximum number of warps per block with 8 blocks per SM (with all warps computing)
+#define SOW 5 // Dual-scheduler optimized number of warps per block (with 8 blocks per SM and to use the dual scheduler with 1 computing warp)
+#define MBH 3 // Maximum number of blocks per SM using half-warp size
+
 /**
  *  @ingroup gpu
  *  @brief Algorithm SAT stage 1

      g_ysum += c_m_size;
    }

        // fix ybar -> y -------------------------

    g_ybar += c_width;
    y = *g_ybar += y;

 *  @param[out] g_ysum All \f$s(P_{m,n}(Y))\f$
 */
__global__
void algSAT_stage2(const int c_n_size, const int c_m_size, const int c_width, float *g_ybar,
    float *g_ysum )
{
  const int tx = threadIdx.x, ty = threadIdx.y,
