    // threads : used for features
 
    const int xOut = blockIdx.x*blockDim.y+threadIdx.y;
-   const bool withinImageBounds = xOut < gradOutput_width;
-   const bool withinGridBounds = blockIdx.x*blockDim.y + threadIdx.x / 2 < gradOutput_width;
-
    const int yOut = blockIdx.y;
    const int width = inputImages_width;
    const int height = inputImages_height;
-
    const int b = blockIdx.z;
 
-   float yf,xf;
-
-   __shared__ float gridData[32];
-   if (threadIdx.y==0 && withinGridBounds)
+   const bool withinImageBounds = xOut < gradOutput_width;
+   const bool withinGridBounds = blockIdx.x*blockDim.y + threadIdx.x / 2 < gradOutput_width;
+
+   if(withinImageBounds && withinGridBounds)
    {
-      gridData[threadIdx.x] = grids_data[b*grids_strideBatch + yOut*grids_strideHeight + xOut*grids_strideWidth + threadIdx.x];
-   }
-   __syncthreads();
-
-   if(withinImageBounds)
-   {
-      yf = gridData[threadIdx.y*2];
-      xf = gridData[threadIdx.y*2+1];
-
-      int yInTopLeft, xInTopLeft;
-      float yWeightTopLeft, xWeightTopLeft;
-      getTopLeft(xf, inputImages_width, xInTopLeft, xWeightTopLeft);
-      getTopLeft(yf, inputImages_height, yInTopLeft, yWeightTopLeft);
-
-      xWeightTopLeft = 0.5;
-      yWeightTopLeft = 0.5;
-
-      const int inTopLeftAddress = inputImages_strideBatch * b + inputImages_strideHeight * yInTopLeft + inputImages_strideWidth * xInTopLeft;
-      const int inTopRightAddress = inTopLeftAddress + inputImages_strideWidth;
-      const int inBottomLeftAddress = inTopLeftAddress + inputImages_strideHeight;
-      const int inBottomRightAddress = inBottomLeftAddress + inputImages_strideWidth;
-
-      const int gradInputImagesTopLeftAddress = gradInputImages_strideBatch * b + gradInputImages_strideHeight * yInTopLeft + gradInputImages_strideWidth * xInTopLeft;
-      const int gradInputImagesTopRightAddress = gradInputImagesTopLeftAddress + gradInputImages_strideWidth;
-      const int gradInputImagesBottomLeftAddress = gradInputImagesTopLeftAddress + gradInputImages_strideHeight;
-      const int gradInputImagesBottomRightAddress = gradInputImagesBottomLeftAddress + gradInputImages_strideWidth;
-
-      const int inTopLeftMaskAddress = masks_strideBatch * b + masks_strideHeight * yInTopLeft + masks_strideWidth * xInTopLeft;
-      const int inTopRightMaskAddress = inTopLeftMaskAddress + masks_strideWidth;
-      const int inBottomLeftMaskAddress = inTopLeftMaskAddress + masks_strideHeight;
-      const int inBottomRightMaskAddress = inBottomLeftMaskAddress + masks_strideWidth;
-
-      const int gradMasksTopLeftAddress = gradMasks_strideBatch * b + gradMasks_strideHeight * yInTopLeft + gradMasks_strideWidth * xInTopLeft;
-      const int gradMasksTopRightAddress = gradMasksTopLeftAddress + gradMasks_strideWidth;
-      const int gradMasksBottomLeftAddress = gradMasksTopLeftAddress + gradMasks_strideHeight;
-      const int gradMasksBottomRightAddress = gradMasksBottomLeftAddress + gradMasks_strideWidth;
-
-      const int gradOutputAddress = gradOutput_strideBatch * b + gradOutput_strideHeight * yOut + gradOutput_strideWidth * xOut;
-
-      float topLeftDotProduct = 0;
-      float topRightDotProduct = 0;
-      float bottomLeftDotProduct = 0;
-      float bottomRightDotProduct = 0;
-
-      bool topLeftIsIn = between(xInTopLeft, 0, width-1) && between(yInTopLeft, 0, height-1);
-      bool topRightIsIn = between(xInTopLeft+1, 0, width-1) && between(yInTopLeft, 0, height-1);
-      bool bottomLeftIsIn = between(xInTopLeft, 0, width-1) && between(yInTopLeft+1, 0, height-1);
-      bool bottomRightIsIn = between(xInTopLeft+1, 0, width-1) && between(yInTopLeft+1, 0, height-1);
-
-      float v = 0;
-      float inTopLeft=0;
-      float inTopRight=0;
-      float inBottomLeft=0;
-      float inBottomRight=0;
-
-      float c = 0;
-
-      float m = 0;
-      float inTopLeftMask=0;
-      float inTopRightMask=0;
-      float inBottomLeftMask=0;
-      float inBottomRightMask=0;
-
-      if(topLeftIsIn) inTopLeftMask = masks_data[inTopLeftMaskAddress];
-      if(topRightIsIn) inTopRightMask = masks_data[inTopRightMaskAddress];
-      if(bottomLeftIsIn) inBottomLeftMask = masks_data[inBottomLeftMaskAddress];
-      if(bottomRightIsIn) inBottomRightMask = masks_data[inBottomRightMaskAddress];
-
-
-      m = xWeightTopLeft * yWeightTopLeft * inTopLeftMask
-        + (1 - xWeightTopLeft) * yWeightTopLeft * inTopRightMask
-        + xWeightTopLeft * (1 - yWeightTopLeft) * inBottomLeftMask
-        + (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * inBottomRightMask;
-
-      /*
-         In that loop we accumulate
-         - gradients into the gradInputImages array with atomic adds
-         - we compute the dot product that we need for the grid gradient
-      */
-
-      for(int t=threadIdx.x; t<inputImages_channels; t+= blockDim.x)
-      {
-
+     int xOut_l = xOut > 0 ? (xOut - 1) : xOut;
+     int yOut_t = yOut > 0 ? (yOut - 1) : yOut;
+     int xOut_r = xOut < (output_width - 1) ? (xOut + 1) : xOut;
+     int yOut_b = yOut < (output_height - 1) ? (yOut + 1) : yOut;
+
+     float yf = grids_data[b*grids_strideBatch + yOut*grids_strideHeight + xOut*grids_strideWidth];
+     float xf = grids_data[b*grids_strideBatch + yOut*grids_strideHeight + xOut*grids_strideWidth + 1];
+
+     // get coordinates for four corners of target response map in the source response map
+     float yf_tl = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_l*grids_strideWidth];
+     float xf_tl = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_l*grids_strideWidth + 1];
+
+     float yf_tr = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_r*grids_strideWidth];
+     float xf_tr = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_r*grids_strideWidth + 1];
+
+     float yf_bl = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_l*grids_strideWidth];
+     float xf_bl = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_l*grids_strideWidth + 1];
+
+     float yf_br = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_r*grids_strideWidth];
+     float xf_br = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_r*grids_strideWidth + 1];
+
+     // compute the bottom left coordinates of source response map
+     float xs_tl = min(min(min(xf_tl, xf_tr), xf_bl), xf_br);
+     float ys_tl = min(min(min(yf_tl, yf_tr), yf_bl), yf_br);
+
+     float xs_br = max(max(max(xf_tl, xf_tr), xf_bl), xf_br);
+     float ys_br = max(max(max(yf_tl, yf_tr), yf_bl), yf_br);
+
+     int yi_t, xi_l;
+     float yWeight_t, xWeight_l;
+
+     int yi_b, xi_r;
+     float yWeight_b, xWeight_r;
+
+     // compute the top left, top right, bottom left, bottom right corner for current grid
+     // compute the nearest bottom right coordiate in source map of top left grid coordiate
+     getBottomRight(xs_tl, inputImages_width, xi_l, xWeight_l);
+     getBottomRight(ys_tl, inputImages_height, yi_t, yWeight_t);
+     // compute the nearest top left coordiate in source map of bottom right grid coordiate
+     getTopLeft(xs_br, inputImages_width, xi_r, xWeight_r);
+     getTopLeft(ys_br, inputImages_height, yi_b, yWeight_b);
+
+     bool topLeftIsIn = between(xi_l, 0, width-1) && between(yi_t, 0, height-1);
+     bool topRightIsIn = between(xi_r, 0, width-1) && between(yi_t, 0, height-1);
+     bool bottomLeftIsIn = between(xi_l, 0, width-1) && between(yi_b, 0, height-1);
+     bool bottomRightIsIn = between(xi_r, 0, width-1) && between(yi_b, 0, height-1);
+
+     const int gradOutputAddress = gradOutput_strideBatch * b + gradOutput_strideHeight * yOut + gradOutput_strideWidth * xOut;
+     if (!topLeftIsIn && !topRightIsIn && !bottomLeftIsIn && !bottomRightIsIn) {
+       for(int t=threadIdx.x; t<inputImages_channels; t+= blockDim.x)
+       {
+         float gradOutValue = gradOutput_data[gradOutputAddress + t];
+         gradCanvas_data[gradOutputAddress + t] = gradOutValue;
+       }
+     }
+     else { // if there are some points inside the freground image and mask
+       float weights[32];
+       float bias_x[32];
+       float bias_y[32];
+       float bias_x_weighted = 0;
+       float bias_y_weighted = 0;
+       float m = 0;
+       float weight_sum = 0;
+       int id_point = 0;
+       for (int y = yi_t; y <= yi_b; ++y) {
+         if (!between(y, 0, height-1)) continue;
+         for (int x = xi_l; x <= xi_r; +=x) {
+           if (!between(x, 0, width-1)) continue;
+           bias_x[id_point] = 2 * (xf - x); // multiply 2 since it is derivative of square
+           bias_y[id_point] = 2 * (yf - y); // multiply 2 since it is derivative of square
+           float weight = __expf(bias_x[id_point] * bias_x[id_point] + bias_y[id_point] * bias_y[id_point]);
+           bias_x_weighted += weight * bias_x[id_point];
+           bias_y_weighted += weight * bias_y[id_point];
+           weight_sum += weight;
+           weights[id_point] = weight;
+           int address = masks_strideBatch * b + masks_strideHeight * y + masks_strideWidth * x;
+           m += weight * masks_data[address];
+           ++id_point;
+         }
+       }
+       m /= weight_sum;
+       float c = 0;
+       float grad_xf = 0;
+       float grad_yf = 0;
+       for(int t=threadIdx.x; t<inputImages_channels; t+= blockDim.x)
+       {
          float gradOutValue = gradOutput_data[gradOutputAddress + t];
          float gradOutValue_fg = gradOutValue * m;
          float gradOutValue_bg = gradOutValue * (1 - m);
          // bool between(int value, int lowerBound, int upperBound)
-         if(topLeftIsIn)
-         {
-            float inTopLeft = inputImages_data[inTopLeftAddress + t];
-            topLeftDotProduct += inTopLeft * gradOutValue_fg;
-            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesTopLeftAddress + t], xWeightTopLeft * yWeightTopLeft * gradOutValue_fg);
-         }
-
-         if(topRightIsIn)
-         {
-            float inTopRight = inputImages_data[inTopRightAddress + t];
-            topRightDotProduct += inTopRight * gradOutValue_fg;
-            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesTopRightAddress + t], (1 - xWeightTopLeft) * yWeightTopLeft * gradOutValue_fg);
-         }
-
-         if(bottomLeftIsIn)
-         {
-            float inBottomLeft = inputImages_data[inBottomLeftAddress + t];
-            bottomLeftDotProduct += inBottomLeft * gradOutValue_fg;
-            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesBottomLeftAddress + t], xWeightTopLeft * (1 - yWeightTopLeft) * gradOutValue_fg);
-         }
-
-         if(bottomRightIsIn)
-         {
-            float inBottomRight = inputImages_data[inBottomRightAddress + t];
-            bottomRightDotProduct += inBottomRight * gradOutValue_fg;
-            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesBottomRightAddress + t], (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * gradOutValue_fg);
-         }
-
-         // jw2yang: copy the gradients outside the object region to canvas, and inside region
-         if (!topLeftIsIn && !topRightIsIn && !bottomLeftIsIn && !bottomRightIsIn) {
-            gradCanvas_data[gradOutputAddress + t] = gradOutValue;
-         }
-         else {
-            gradCanvas_data[gradOutputAddress + t] = gradOutValue_bg;
-         }
-
-         // jw2yang: compute the gradient mask value
-         if(topLeftIsIn) inTopLeft = inputImages_data[inTopLeftAddress + t];
-         if(topRightIsIn) inTopRight = inputImages_data[inTopRightAddress + t];
-         if(bottomLeftIsIn) inBottomLeft = inputImages_data[inBottomLeftAddress + t];
-         if(bottomRightIsIn) inBottomRight = inputImages_data[inBottomRightAddress + t];
-         v = xWeightTopLeft * yWeightTopLeft * inTopLeft
-           + (1 - xWeightTopLeft) * yWeightTopLeft * inTopRight
-           + xWeightTopLeft * (1 - yWeightTopLeft) * inBottomLeft
-           + (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * inBottomRight;
-
+         v = 0;
+         id_point = 0;
+         for (int y = yi_t; y <= yi_b; ++y) {
+           if (!between(y, 0, height-1)) continue;
+           for (int x = xi_l; x <= xi_r; +=x) {
+             if (!between(x, 0, width-1)) continue;
+             int address = inputImages_strideBatch * b + inputImages_strideHeight * y + inputImages_strideWidth * x + t;
+             v += weight[id_point] * inputImages_data[address];
+           }
+         }
+         v /= weight_sum;
          c = canvas_data[gradOutputAddress + t];
-
          float gradMaskValue = gradOutValue * (v - c);
 
-
-         // update gradient on mask map
-         if(topLeftIsIn)
-         {
-            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksTopLeftAddress], xWeightTopLeft * yWeightTopLeft * gradMaskValue);
-         }
-
-         if(topRightIsIn)
-         {
-            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksTopRightAddress], (1 - xWeightTopLeft) * yWeightTopLeft * gradMaskValue);
-         }
-
-         if(bottomLeftIsIn)
-         {
-            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksBottomLeftAddress], xWeightTopLeft * (1 - yWeightTopLeft) * gradMaskValue);
-         }
-
-         if(bottomRightIsIn)
-         {
-            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksBottomRightAddress], (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * gradMaskValue);
-         }
-      }
-      /*
-         Here we reduce the dot product and compute the grid gradient before writing it.
-      */
-
-      /* could do shuffles and use no shmem at all but cuda arch is 2.0 */
-      __shared__ volatile float __shmem[16][32];
-      __shmem[threadIdx.y][threadIdx.x] = topLeftDotProduct;
-      sumReduceShMem(__shmem[threadIdx.y]);
-      topLeftDotProduct = __shmem[threadIdx.y][0];
-
-      __shmem[threadIdx.y][threadIdx.x] = topRightDotProduct;
-      sumReduceShMem(__shmem[threadIdx.y]);
-      topRightDotProduct = __shmem[threadIdx.y][0];
-
-      __shmem[threadIdx.y][threadIdx.x] = bottomLeftDotProduct;
-      sumReduceShMem(__shmem[threadIdx.y]);
-      bottomLeftDotProduct = __shmem[threadIdx.y][0];
-
-      __shmem[threadIdx.y][threadIdx.x] = bottomRightDotProduct;
-      sumReduceShMem(__shmem[threadIdx.y]);
-      bottomRightDotProduct = __shmem[threadIdx.y][0];
-
-      yf = - xWeightTopLeft * topLeftDotProduct + xWeightTopLeft * bottomLeftDotProduct - (1-xWeightTopLeft) * topRightDotProduct + (1-xWeightTopLeft) * bottomRightDotProduct;
-      xf = - yWeightTopLeft * topLeftDotProduct + yWeightTopLeft * topRightDotProduct - (1-yWeightTopLeft) * bottomLeftDotProduct + (1-yWeightTopLeft) * bottomRightDotProduct;
-
-      if(threadIdx.x==0)
-      {
-         gridData[threadIdx.y*2] = yf * (inputImages_height-1) / 2;
-         gridData[threadIdx.y*2+1] = xf * (inputImages_width-1) / 2;
-      }
-   }// must put a big if condition in order not to hang at __syncthreads()...
-   __syncthreads();
-
-   if(threadIdx.y==0 && withinGridBounds)
-       gradGrids_data[b*gradGrids_strideBatch + yOut*gradGrids_strideHeight + xOut*gradGrids_strideWidth + threadIdx.x] = gridData[threadIdx.x];
-}
-
-
-
-
+         id_point = 0;
+         for (int y = yi_t; y <= yi_b; ++y) {
+           if (!between(y, 0, height-1)) continue;
+           for (int x = xi_l; x <= xi_r; +=x) {
+             if (!between(x, 0, width-1)) continue;
+             int address = inputImages_strideBatch * b + inputImages_strideHeight * y + inputImages_strideWidth * x;
+             atomicAdd(&gradInputImages_data[address + t], weights[id_point] * gradOutValue_fg);
+             atomicAdd(&gradMasks_data[address], weights[id_point] * gradMaskValue);
+             ++id_point;
+           }
+         }
+         gradCanvas_data[gradOutputAddress + t] = gradOutValue_bg;
+
+         // compute gradient on grid coordinates
+         id_point = 0;
+         for (int y = yi_t; y <= yi_b; ++y) {
+           if (!between(y, 0, height-1)) continue;
+           for (int x = xi_l; x <= xi_r; +=x) {
+             if (!between(x, 0, width-1)) continue;
+             int address = inputImages_strideBatch * b + inputImages_strideHeight * y + inputImages_strideWidth * x;
+             grad_yf += gradOutValue_fg * inputImages_data[address + t]
+                        * weights[id_point]
+                        * (bias_y[id_point] - bias_y_weighted / weight_sum)
+                        / weight_sum;
+             grad_xf += gradOutValue_fg * inputImages_data[address + t]
+                        * weights[id_point]
+                        * (bias_x[id_point] - bias_x_weighted / weight_sum)
+                        / weight_sum;
+           }
+         }
+       }
+       gradGrids_data[b*gradGrids_strideBatch + yOut*gradGrids_strideHeight + xOut*gradGrids_strideWidth] = grad_yf * (inputImages_height-1) / 2;
+       gradGrids_data[b*gradGrids_strideBatch + yOut*gradGrids_strideHeight + xOut*gradGrids_strideWidth + 1] = grad_xf * (inputImages_width-1) / 2;
+     }
+   }
+}
 
 static int cunn_SubSamplerBHWD_updateGradInput(lua_State *L)
 {

      if (!topLeftIsIn && !topRightIsIn && !bottomLeftIsIn && !bottomRightIsIn)
        output_data[outAddress + t] = canvas_data[outAddress + t];

      if(topLeftIsIn) inTopLeft = inputImages_data[inTopLeftAddress + t];
      if(topRightIsIn) inTopRight = inputImages_data[inTopRightAddress + t];
      if(bottomLeftIsIn) inBottomLeft = inputImages_data[inBottomLeftAddress + t];
      if(bottomRightIsIn) inBottomRight = inputImages_data[inBottomRightAddress + t];

      v = xWeightTopLeft * yWeightTopLeft * inTopLeft
        + (1 - xWeightTopLeft) * yWeightTopLeft * inTopRight
        + xWeightTopLeft * (1 - yWeightTopLeft) * inBottomLeft
        + (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * inBottomRight;

      // we do not replace the canvas region with foreground, instead, we add value together.
      output_data[outAddress + t] = (1 - m) * canvas_data[outAddress + t] + m * v;

   const bool withinImageBounds = xOut < output_width; // check whether x exceed the boundary
   const bool withinGridBounds = blockIdx.x*blockDim.y + threadIdx.x / 2 < output_width; // check whether

   if(!withinImageBounds || !withinGridBounds) return;

   const int width = inputImages_width;
   const int height = inputImages_height;

   int xOut_l = xOut > 0 ? (xOut - 1) : xOut;
   int yOut_t = yOut > 0 ? (yOut - 1) : yOut;
   int xOut_r = xOut < (output_width - 1) ? (xOut + 1) : xOut;
   int yOut_b = yOut < (output_height - 1) ? (yOut + 1) : yOut;

   float yf = grids_data[b*grids_strideBatch + yOuy*grids_strideHeight + xOut*grids_strideWidth];
   float xf = grids_data[b*grids_strideBatch + yOut*grids_strideHeight + xOut*grids_strideWidth + 1];

   // get coordinates for four corners of target response map in the source response map
   float yf_tl = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_l*grids_strideWidth];
   float xf_tl = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_l*grids_strideWidth + 1];

   float yf_tr = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_r*grids_strideWidth];
   float xf_tr = grids_data[b*grids_strideBatch + yOut_t*grids_strideHeight + xOut_r*grids_strideWidth + 1];

   float yf_bl = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_l*grids_strideWidth];
   float xf_bl = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_l*grids_strideWidth + 1];

   float yf_br = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_r*grids_strideWidth];
   float xf_br = grids_data[b*grids_strideBatch + yOut_b*grids_strideHeight + xOut_r*grids_strideWidth + 1];

   // compute the bottom left coordinates of source response map
   float xs_tl = min(min(min(xf_tl, xf_tr), xf_bl), xf_br);
   float ys_tl = min(min(min(yf_tl, yf_tr), yf_bl), yf_br);

   float xs_br = max(max(max(xf_tl, xf_tr), xf_bl), xf_br);
   float ys_br = max(max(max(yf_tl, yf_tr), yf_bl), yf_br);

   int yi_t, xi_l;
   float yWeight_t, xWeight_l;

   int yi_b, xi_r;
   float yWeight_b, xWeight_r;

   // compute the top left, top right, bottom left, bottom right corner for current grid
   // compute the nearest bottom right coordiate in source map of top left grid coordiate
   getBottomRight(xs_tl, inputImages_width, xi_l, xWeight_l);
   getBottomRight(ys_tl, inputImages_height, yi_t, yWeight_t);
   // compute the nearest top left coordiate in source map of bottom right grid coordiate
   getTopLeft(xs_br, inputImages_width, xi_r, xWeight_r);
   getTopLeft(ys_br, inputImages_height, yi_b, yWeight_b);

   // compute average mask value in region [topleft, bottom right]
   float weights[32];
   float m = 0;
   float weight_sum = 0;
   int id_point = 0;
   for (int y = yi_t; y <= yi_b; ++y) {
     if (!between(y, 0, height-1)) continue;
     for (int x = xi_l; x <= xi_r; +=x) {
       if (!between(x, 0, width-1)) continue;
       float weight = __expf((x - xf) * (x - xf) + (y - yf) * (y - yf));
       weight_sum += weight;
       weights[id_point] = weight;
       int address = masks_strideBatch * b + masks_strideHeight * y + masks_strideWidth * x;
       m += weight * masks_data[address];
       ++id_point;
     }
   }
   m /= weight_sum;

   float v=0;
   float inTopLeft=0;
