docs: def _convolve_spatial2(im, hs,
                        mode = "constant",
                        grid_dim = None,
                        pad_factor = 2,
                        plan = None,
                        return_plan = False):
      """
      spatial varying convolution of an 2d image with a 2d grid of psfs
  
      shape(im_ = (Ny,Nx)
      shape(hs) = (Gy,Gx, Hy,Hx)
  
      the input image im is subdivided into (Gy,Gx) blocks
      hs[j,i] is the psf at the center of each block (i,j)
  
      as of now each image dimension has to be divisible by the grid dim, i.e.
      Nx % Gx == 0
      Ny % Gy == 0
  
  
      mode can be:
      "constant" - assumed values to be zero
      "wrap" - periodic boundary condition
      """
  
      if grid_dim:
          Gs = tuple(grid_dim)
      else:
          Gs = hs.shape[:2]
  
  
      mode_str = {"constant":"CLK_ADDRESS_CLAMP",
                  "wrap":"CLK_ADDRESS_REPEAT"}
  
      Ny, Nx = im.shape
      Gy, Gx = Gs
  
  
      # the size of each block within the grid
      Nblock_y, Nblock_x = Ny // Gy, Nx // Gx
  
  
      # the size of the overlapping patches with safety padding
      Npatch_x, Npatch_y = _next_power_of_2(pad_factor*Nblock_x), _next_power_of_2(pad_factor*Nblock_y)
  
  
      prog = OCLProgram(abspath("kernels/conv_spatial2.cl"),
                        build_options=["-D","ADDRESSMODE=%s"%mode_str[mode]])
  
      if plan is None:
          plan = fft_plan((Gy, Gx, Npatch_y,Npatch_x), axes = (-2,-1))
  
      x0s = Nblock_x*np.arange(Gx)
      y0s = Nblock_y*np.arange(Gy)
  
  
      patches_g = OCLArray.empty((Gy,Gx,Npatch_y,Npatch_x),np.complex64)
  
      #prepare psfs
      if grid_dim:
          h_g = OCLArray.zeros((Gy,Gx,Npatch_y,Npatch_x),np.complex64)
          tmp_g = OCLArray.from_array(hs.astype(np.float32, copy = False))
          for i,_x0 in enumerate(x0s):
              for j,_y0 in enumerate(y0s):
                  prog.run_kernel("fill_psf_grid2",
                                  (Nblock_x,Nblock_y),None,
                          tmp_g.data,
                          np.int32(Nx),
                          np.int32(i*Nblock_x),
                          np.int32(j*Nblock_y),
                          h_g.data,
                          np.int32(Npatch_x),
                          np.int32(Npatch_y),
                          np.int32(-Nblock_x//2+Npatch_x//2),
                          np.int32(-Nblock_y//2+Npatch_y//2),
                          np.int32(i*Npatch_x*Npatch_y+j*Gx*Npatch_x*Npatch_y)
                              )
      else:
          hs = np.fft.fftshift(pad_to_shape(hs,(Gy,Gx,Npatch_y,Npatch_x)),axes=(2,3))
          h_g = OCLArray.from_array(hs.astype(np.complex64))
  
  
      #prepare image
      im_g = OCLImage.from_array(im.astype(np.float32,copy=False))
  
      for i,_x0 in enumerate(x0s):
          for j,_y0 in enumerate(y0s):
              prog.run_kernel("fill_patch2",(Npatch_x,Npatch_y),None,
                      im_g,
                      np.int32(_x0+Nblock_x//2-Npatch_x//2),
                      np.int32(_y0+Nblock_y//2-Npatch_y//2),
                      patches_g.data,
                      np.int32(i*Npatch_x*Npatch_y+j*Gx*Npatch_x*Npatch_y))
  
  
      #return np.abs(patches_g.get())
      # convolution
      fft(patches_g,inplace=True,  plan = plan)
      fft(h_g,inplace=True,  plan = plan)
      prog.run_kernel("mult_inplace",(Npatch_x*Npatch_y*Gx*Gy,),None,
                      patches_g.data, h_g.data)
      fft(patches_g,inplace=True, inverse = True, plan = plan)
  
  
      logger.debug("Nblock_x: {}, Npatch_x: {}".format(Nblock_x, Npatch_x))
      #return np.abs(patches_g.get())
      #accumulate
      res_g = OCLArray.empty(im.shape,np.float32)
  
      for j in range(Gy+1):
          for i in range(Gx+1):
              prog.run_kernel("interpolate2",(Nblock_x,Nblock_y),None,
                              patches_g.data,res_g.data,
                              np.int32(i),np.int32(j),
                              np.int32(Gx),np.int32(Gy),
                              np.int32(Npatch_x),np.int32(Npatch_y))
  
      res = res_g.get()
  
      if return_plan:
          return res, plan
      else:
          return res