code: def _create_pairwise_gaussian_2d(sx, sy, shape):
      
      feat_size = 2
      feats = np.zeros((feat_size, shape[0], shape[1]), dtype=np.float32)
      for i in range(shape[0]):
          for j in range(shape[1]):
              feats[0, i, j] = i / sx
              feats[1, i, j] = j / sy
      return feats.reshape([feat_size, -1])
docs: A simple reference implementation for the 2D case. The ND implementation
      is faster.