 #define BLOCKX 16
 #define BLOCKY 16
 #define THREADS_PER_BLK (BLOCKX * BLOCKY)
+#define SCAN_BLOCK_DIM  (BLOCKX * BLOCKY) // needed by sharedMemExclusiveScan implementation
+
 __constant__ float  cuConstColorRamp[COLOR_MAP_SIZE][3];
 
 

//
// Each thread renders a pixel
__global__ void kernelRenderBlocks() {
    __shared__ int support[THREADS_PER_BLK]; 

    int imageX = blockIdx.x * blockDim.x + threadIdx.x;
    int imageY = blockIdx.y * blockDim.y + threadIdx.y;

    float boxB = (float)(blockIdx.y * blockDim.y)/height;
    float boxT = (float)((blockIdx.y+1) * blockDim.y)/height;

    int linearThreadIndex = threadIdx.y * blockDim.x+threadIdx.x;

    //The loop guard ensures all threads will enter the for loop, no matter whether it corresponds to a circle
    //Example of possible thread_iter value: 0, 256, 512,...
    for (int thread_iter= 0; thread_iter<cuConstRendererParams.numCircles; thread_iter+=THREADS_PER_BLK){

        //Each thread is assigned to a circleIndex, whether exists or not
        int circleIndex = thread_iter + linearThreadIndex;

        if (circleIndex<cuConstRendererParams.numCircles){ //check circle index in-bound
            int index3_circ = 3 * circleIndex;
            float3 p_circ = *(float3*)(&cuConstRendererParams.position[index3_circ]);
            float rad = cuConstRendererParams.radius[circleIndex];
            scanInput[linearThreadIndex] = (circleInBoxConservative(p_circ.x, p_circ.y, rad, boxL,boxR,boxT,boxB)) ? 1:0;
        } else {
            scanInput[linearThreadIndex] = 0;
        }

        __syncthreads();

        // Exclusive scan
        sharedMemExclusiveScan(linearThreadIndex, scanInput, scanOutput, scanScratch, THREADS_PER_BLK);
        __syncthreads();
        // See if last circle in thread block is in block region
        int num_circle_in_block = scanInput[THREADS_PER_BLK-1] ? scanOutput[THREADS_PER_BLK-1]+1 : scanOutput[THREADS_PER_BLK-1];
        
        
        if (scanInput[linearThreadIndex]){
            // Scan output index is off by 1
            int circleInBlockIdx = (linearThreadIndex == THREADS_PER_BLK-1)?  scanOutput[THREADS_PER_BLK-1]: scanOutput[linearThreadIndex+1]-1; 
            // Push all circle in block region forward
            circleInBlock[circleInBlockIdx] = linearThreadIndex;
        }
        __syncthreads();

        //printf("%d find %d\n", threadIdx.y*blockDim.x+threadIdx.x, support[0]);*/
        for (int i=0; i < num_circle_in_block; i++){               
                int index3 = 3 * (circleInBlock[i] + thread_iter);

                // read position and radius
                float3 p = *(float3*)(&cuConstRendererParams.position[index3]);
                float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(imageX) + 0.5f),
                                                     invHeight * (static_cast<float>(imageY) + 0.5f));
                shadePixel(circleInBlock[i]+thread_iter,pixelCenterNorm, p, imgPtr);
        }
        __syncthreads();
    }
