   int blocksPerGrid = (size.z + threadsPerBlock -1) >> 8;
   dim3 DimBlock(threadsPerBlock);
   dim3 DimGrid(blocksPerGrid);
+  int nBlocks = blocksPerGrid;
 
   int threadsPerBlockDim = 16;
   int blocksPerGridX = ((size.x) + threadsPerBlockDim -1) >> 4;
   int blocksPerGridY = ((size.y) + threadsPerBlockDim -1) >> 4;
-  int nBlocks = blocksPerGridX*blocksPerGridY;
   dim3 TwoDimBlock(16,16,1);
   dim3 TwoDimGrid(blocksPerGridX,blocksPerGridY,1);
 

  cudaMalloc((void**) &d_modified, (pow2*sizeof(int)));
  CUT_CHECK_ERROR("Memory hysteresis image creation failed");

  kernel_hysteresis_glm3<<<TwoDimGrid,TwoDimBlock>>>(d_hys, size, d_modified, nBlocks, pow2);
  CUT_CHECK_ERROR("Hysteresis Kernel failed");

  /// bind a texture to the CUDA array
  cudaBindTexture (NULL ,hysTexRef, d_hys);


}

__device__ int reduceSum256(int tid, int sdata[]){
// data[] should be a 256 positios array

  if (tid < 128) { sdata[tid] += sdata[tid + 128]; }
  __syncthreads();
  if (tid < 64) { sdata[tid] += sdata[tid + 64]; }
  __syncthreads();
  if (tid < 32) { sdata[tid] += sdata[tid + 32]; }
  __syncthreads();
  if (tid < 16) { sdata[tid] += sdata[tid + 16]; }
  __syncthreads();
  if (tid < 8) { sdata[tid] += sdata[tid + 8]; }
  __syncthreads();
  if (tid < 4) { sdata[tid] += sdata[tid + 4]; }
  __syncthreads();
  if (tid < 2) { sdata[tid] += sdata[tid + 2]; }
  __syncthreads();
  if (tid == 0) { sdata[tid] += sdata[tid + 1]; }
  __syncthreads();
  
  return(sdata[0]);

}

__device__ int reduceSum(int tid, int data[], int N, int pow2){

  if (tid<N){

    for(unsigned int s=pow2/2; s>0; s>>=1) {
      if ((tid < s) && (tid + s < pow2)){
          data[tid] += data[tid+s];
      }
      __syncthreads();
    }

