Pointclouds from rendered depth images do not align

Question

I am trying to render multiple RGB-D images of a given object in Blender 2.93 using the python interface. The scene contains a single object at position [0, 0, 0] and I am rotating a camera around it on the unit sphere. When I am generating the point clouds from the depth images, they do not align (see images below). Does Blender add an offset to the depth?

EDIT: I don't think the issue here lies with the inverse projection method. If this were the case, points would show bending artifacts, which is not the case. I have also overloaded the original mesh with the generated point cloud and this shows a clear downward/inward shift, which to me indicates a parameter that is internally computed in blender but not reflected in the intrinsic camera matrix that I am saving to file to compute the point cloud.

I have attached a minimal code example of loading the depth data and generating the point cloud from it below. Any help would be greatly appreciated!

import os
import glob
import numpy as np
import open3d as o3d
import cv2


def get_rays_np(H, W, K, c2w):
    i, j = np.meshgrid(
        np.arange(W, dtype=np.float32), np.arange(H, dtype=np.float32), indexing="xy"
    )
    dirs = np.stack(
        ((i - K[0, 2]) / K[0, 0], (j - K[1, 2]) / K[1, 1], np.ones_like(i)), axis=-1
    )
    rays_d = np.sum(
        dirs[..., np.newaxis, :] * c2w[:3, :3], -1
    )  # dot product, equals to: [c2w.dot(dir) for dir in dirs]
    # Translate camera frame's origin to the world frame. It is the origin of all rays.
    rays_o = np.broadcast_to(c2w[:3, -1], np.shape(rays_d))
    return np.concatenate((rays_o, rays_d), axis=-1)


base_path = "./images/30b9882d1d75be3537678474be485ca/"

img_files = sorted(glob.glob(os.path.join(base_path, "rgb", "rgb_*")))
mask_files = sorted(glob.glob(os.path.join(base_path, "rgb", "mask_*")))
depth_files = sorted(glob.glob(os.path.join(base_path, "rgb", "depth_*")))
pose_files = sorted(glob.glob(os.path.join(base_path, "pose", "*.npy")))

K = np.load(os.path.join(base_path, "intrinsics.npy"))
for ii in range(1, n_views):
    img = cv2.imread(img_files[ii]) / 255.0
    mask = (cv2.imread(mask_files[ii], cv2.IMREAD_GRAYSCALE) / 255.0).astype(int)
    depth = cv2.imread(depth_files[ii], cv2.IMREAD_ANYDEPTH)
    c2w = np.load(os.path.join(pose_files[ii]))
    H, W = depth.shape

    # Create the point cloud

    depth = depth.reshape(-1)
    rays = get_rays_np(H, W, K, c2w)

    u = np.arange(0, W)
    v = np.arange(0, H)
    u, v = np.meshgrid(u, v)
    u, v = u.reshape(-1), v.reshape(-1)
    fx, fy = K[0, 0], K[1, 1]
    cx, cy = K[0, 2], K[1, 2]

    uv = np.stack((u, v, np.ones_like(u)), axis=-1)
    xy = (np.linalg.inv(K) @ np.transpose(uv)) * depth[None, :]
    x = np.transpose(xy)[:, 0]
    y = np.transpose(xy)[:, 1]

    # Get distance along the ray from z-depth
    dist = np.sqrt(x ** 2 + y ** 2 + depth.reshape(-1) ** 2)

    rays = rays.reshape(-1, 6)
    rays0, raysd = rays[:, :3], rays[:, 3:]

    # Normalize ray direction
    raysd = raysd / np.linalg.norm(raysd, axis=-1, keepdims=True)
    pts = rays0 + raysd * (dist[:, None])

    indices = np.where(mask.reshape(-1))
    pts = pts[indices]
    colors = img.reshape(-1, 3)[indices]

    frame = o3d.geometry.TriangleMesh.create_coordinate_frame(0.1).transform(c2w)
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(pts)
    pcd.colors = o3d.utility.Vector3dVector(colors)
    pcds.append(pcd)
    pcds.append(frame)
    o3d.visualization.draw_geometries(pcds)