 // implements approximation to 8x8 bilateral filtering
 __global__ void adaptiveMeanGPU8 (float* D, float* D_copy, float* D_horiz, int32_t D_width, int32_t D_height) {
   
-  // Pixel id
+  // Global coordinates and Pixel id
   uint32_t u0 = blockDim.x*blockIdx.x + threadIdx.x + 4;
   uint32_t v0 = blockDim.y*blockIdx.y + threadIdx.y + 4;
   uint32_t idx = v0*D_width + u0;
+  //Local thread coordinates
+  uint32_t ut = threadIdx.x + 4;
+  uint32_t vt = threadIdx.y + 4;
+  
 
   if(u0 > (D_width - 4) || v0 > (D_height - 4))
     return;

  float weight_sum = 0;
  float factor_sum = 0;

  for(int32_t i=0; i < 8; i++){
    weight_sum0 = 4.0f - fabs(*(D_copy+idx+(i-4))-val_curr);
    weight_sum0 = max(0.0f, weight_sum0);
    weight_sum += weight_sum0;
    factor_sum += *(D_copy+idx+(i-4))*weight_sum0;
  }
  
  if (weight_sum>0) {
      float d = factor_sum/weight_sum;
      if (d>=0) *(D_horiz+idx) = d;
  }
  
  __syncthreads();

  // vertical filter

    *(D_horiz+idx)  = -10;
  }

  __shared__ float D_copy_S[32+8][32+8];
  if(threadIdx.x == blockDim.x-1){
      D_copy_S[ut+1][vt] = D[idx+1];
      D_copy_S[ut+2][vt] = D[idx+2];
      D_copy_S[ut+3][vt] = D[idx+3];
      D_copy_S[ut+4][vt] = D[idx+4];
  }
  if(threadIdx.x == 0){
      D_copy_S[ut-4][vt] = D[idx-4];
      D_copy_S[ut-3][vt] = D[idx-3];
      D_copy_S[ut-2][vt] = D[idx-2];
      D_copy_S[ut-1][vt] = D[idx-1];
  }
  if(threadIdx.y == 0){
      D_copy_S[ut][vt-4] = D[(v0-4)*D_width+u0];
      D_copy_S[ut][vt-3] = D[(v0-3)*D_width+u0];
      D_copy_S[ut][vt-2] = D[(v0-2)*D_width+u0];
      D_copy_S[ut][vt-1] = D[(v0-1)*D_width+u0];
  }
  if(threadIdx.y == blockDim.y-1){
      D_copy_S[ut][vt+1] = D[(v0+1)*D_width+u0];
      D_copy_S[ut][vt+2] = D[(v0+2)*D_width+u0];
      D_copy_S[ut][vt+3] = D[(v0+3)*D_width+u0];
      D_copy_S[ut][vt+4] = D[(v0+4)*D_width+u0];
  }

  if(D[idx] < 0){
      D_copy_S[ut][vt] = -10;
  }else{
      D_copy_S[ut][vt] = D[idx];
  }
  __syncthreads();
      
  // full resolution: 8 pixel bilateral filter width
