def METHOD_NAME(keypoints: np.ndarray, heatmaps: np.ndarray,
                          blur_kernel_size: int) -> np.ndarray:
    """Refine keypoint predictions using distribution aware coordinate
    decoding. See `Dark Pose`_ for details. The operation is in-place.
    Note:
        - instance number: N
        - keypoint number: K
        - keypoint dimension: D
        - heatmap size: [W, H]
    Args:
        keypoints (np.ndarray): The keypoint coordinates in shape (N, K, D)
        heatmaps (np.ndarray): The heatmaps in shape (K, H, W)
        blur_kernel_size (int): The Gaussian blur kernel size of the heatmap
            modulation
    Returns:
        np.ndarray: Refine keypoint coordinates in shape (N, K, D)
    .. _`Dark Pose`: https://arxiv.org/abs/1910.06278
    """
    N, K = keypoints.shape[:2]
    H, W = heatmaps.shape[1:]
    # modulate heatmaps
    heatmaps = gaussian_blur(heatmaps, blur_kernel_size)
    np.maximum(heatmaps, 1e-10, heatmaps)
    np.log(heatmaps, heatmaps)
    for n, k in product(range(N), range(K)):
        x, y = keypoints[n, k, :2].astype(int)
        if 1 < x < W - 2 and 1 < y < H - 2:
            dx = 0.5 * (heatmaps[k, y, x + 1] - heatmaps[k, y, x - 1])
            dy = 0.5 * (heatmaps[k, y + 1, x] - heatmaps[k, y - 1, x])
            dxx = 0.25 * (
                heatmaps[k, y, x + 2] - 2 * heatmaps[k, y, x] +
                heatmaps[k, y, x - 2])
            dxy = 0.25 * (
                heatmaps[k, y + 1, x + 1] - heatmaps[k, y - 1, x + 1] -
                heatmaps[k, y + 1, x - 1] + heatmaps[k, y - 1, x - 1])
            dyy = 0.25 * (
                heatmaps[k, y + 2, x] - 2 * heatmaps[k, y, x] +
                heatmaps[k, y - 2, x])
            derivative = np.array([[dx], [dy]])
            hessian = np.array([[dxx, dxy], [dxy, dyy]])
            if dxx * dyy - dxy**2 != 0:
                hessianinv = np.linalg.inv(hessian)
                offset = -hessianinv @ derivative
                offset = np.squeeze(np.array(offset.T), axis=0)
                keypoints[n, k, :2] += offset
    return keypoints