def METHOD_NAME(x1, x2, patch=(1, 1), shift=(0, 0), patch_step=(1, 1),
                      shift_step=(1, 1), padding=(0, 0, 0, 0)):
    # inputs and output are NHWC tensors
    oh = (x1.shape[1] + sum(padding[:2]) -
          patch[0] + patch_step[0]) // patch_step[0]
    ow = (x1.shape[2] + sum(padding[2:]) -
          patch[1] + patch_step[1]) // patch_step[1]
    och = int(np.ceil((1 + 2 * shift[0]) / shift_step[0]))
    ocw = int(np.ceil((1 + 2 * shift[1]) / shift_step[1]))
    out_array = np.empty((x1.shape[0], oh, ow, och, ocw), dtype=x1.dtype)
    out_index = np.ndindex(out_array.shape)
    x1 = np.pad(x1, ((0, 0), padding[:2], padding[2:], (0, 0)), 'constant')
    x2 = np.pad(x2, ((0, 0), padding[:2], padding[2:], (0, 0)), 'constant')
    x2 = np.pad(x2, ((0, 0), *zip(shift, shift), (0, 0)), 'constant')
    for n in range(x1.shape[0]):
        for y in range(0, x1.shape[1] - patch[0] + 1, patch_step[0]):
            for x in range(0, x1.shape[2] - patch[1] + 1, patch_step[1]):
                p1 = x1[n, y:y+patch[0], x:x+patch[1]]
                for yy in range(y, y + 1 + 2 * shift[0], shift_step[0]):
                    for xx in range(x, x + 1 + 2 * shift[1], shift_step[1]):
                        p2 = x2[n, yy:yy+patch[0], xx:xx+patch[1]]
                        out_array[next(out_index)] = np.sum(p1 * p2)
    return out_array