四元数

在这篇文章中,我们处理上篇文章遗留的四元数解析。


从“四元数”到“3x3 矩阵”

在图形学渲染中,我们最习惯使用的是 3 \times 3 的旋转矩阵(Rotation Matrix)。然而,很多数据集(或 Colmap 这样的三维重建软件)为了避免“万向节死锁”并压缩存储,通常会使用四元数 (Quaternion) 来表示旋转。

四元数由四个值组成:x, y, z, w。其中 x, y, z 表示旋转所围绕的方向向量,而 w 则与旋转的角度相关。为了能让我们的相机工具矩阵正常工作,我们此处手写一个转换函数。

将四元数转换为矩阵时,我们应用以下数学公式:

R = \begin{bmatrix} 1 - 2y^2 - 2z^2 & 2xy - 2zw & 2xz + 2yw \\ 2xy + 2zw & 1 - 2x^2 - 2z^2 & 2yz - 2xw \\ 2xz - 2yw & 2yz + 2xw & 1 - 2x^2 - 2y^2 \end{bmatrix}

  1. def quat_to_rotmat(quat):
  2.     """
  3.     Convert quaternion [x, y, z, w] to a 3x3 rotation matrix.
  4.     Supports batch dimensions (..., 4) -> (..., 3, 3).
  5.     """
  6.     # Normalize quaternion to ensure valid rotation matrix
  7.     quat = quat / torch.norm(quat, dim=-1, keepdim=True)
  8.  
  9.     x, y, z, w = torch.unbind(quat, dim=-1)
  10.  
  11.     # Precompute products
  12.     x2 = x * x
  13.     y2 = y * y
  14.     z2 = z * z
  15.     xy = x * y
  16.     xz = x * z
  17.     yz = y * z
  18.     xw = x * w
  19.     yw = y * w
  20.     zw = z * w
  21.  
  22.     # Calculate R components
  23.     r00 = 1.0 - 2.0 * (y2 + z2)
  24.     r01 = 2.0 * (xy - zw)
  25.     r02 = 2.0 * (xz + yw)
  26.  
  27.     r10 = 2.0 * (xy + zw)
  28.     r11 = 1.0 - 2.0 * (x2 + z2)
  29.     r12 = 2.0 * (yz - xw)
  30.  
  31.     r20 = 2.0 * (xz - yw)
  32.     r21 = 2.0 * (yz + xw)
  33.     r22 = 1.0 - 2.0 * (x2 + y2)
  34.  
  35.     # Stack along last dimension and reshape
  36.     rot_matrix = torch.stack([r00, r01, r02, r10, r11, r12, r20, r21, r22], dim=-1)
  37.     rot_matrix = rot_matrix.reshape(*quat.shape[:-1], 3, 3)
  38.     return rot_matrix

输入 quat 的形状可能是 (N, 4)。所以我们不能简单地用索引提取标量,因为每一列会包含 N 个相机的数据。unbind(-1) 会沿着最后一个维度将张量拆分成 4 个形状为 (N,) 的张量,这保证了后续运算都是高度并行的矩阵计算。


Colmap:从“世界到相机”逆转为“相机到世界”

Colmap 软件输出的矩阵是“世界到相机 (World-to-Camera, W2C)”的矩阵,而我们的渲染管线需要的是“相机到世界 (Camera-to-World, C2W)”的矩阵。

我们需要进行一次逆变换。变换的原理我们在之前的文章中已经了解过了。

假设原始的 W2C 旋转矩阵为 R_{w2c},平移向量为 T_{w2c}。转换为 C2W 时:

1. 旋转矩阵的逆等于它的转置:

R_{c2w} = R_{w2c}^T

2. 平移向量结合新的旋转矩阵进行反向推导:

T_{c2w} = R_{c2w} \times (-T_{w2c})

  1. def w2c_to_c2w(q_w2c, t_w2c):
  2.     """
  3.     Convert world-to-camera parameters to camera-to-world matrix.
  4.     q_w2c: quaternion of shape (..., 4)
  5.     t_w2c: translation of shape (..., 3)
  6.     Returns: c2w matrix of shape (..., 4, 4)
  7.     """
  8.     R_w2c = quat_to_rotmat(q_w2c)
  9.  
  10.     # R_c2w = R_w2c.T (transpose the rotation matrix)
  11.     R_c2w = R_w2c.transpose(-1, -2)
  12.  
  13.     # t_c2w = -R_c2w @ t_w2c
  14.     t_c2w = -torch.matmul(R_c2w, t_w2c.unsqueeze(-1)).squeeze(-1)
  15.  
  16.     # Construct 4x4 matrix
  17.     shape = q_w2c.shape[:-1]
  18.     if len(shape) == 0:
  19.         c2w = torch.eye(4, dtype=q_w2c.dtype, device=q_w2c.device)
  20.         c2w[:3, :3] = R_c2w
  21.         c2w[:3, 3] = t_c2w
  22.     else:
  23.         c2w = torch.eye(4, dtype=q_w2c.dtype, device=q_w2c.device).repeat(*shape, 1, 1)
  24.         c2w[..., :3, :3] = R_c2w
  25.         c2w[..., :3, 3] = t_c2w
  26.  
  27.     return c2w

内参缩放

搞定了四元数,我们解析得到相机外参。

  1. # 3. Load camera poses and create camera pose
  2. cameras_path = "datasets/out_colmap/bonsai/cameras.npy"
  3. print(f"Loading camera parameters from {cameras_path}...")
  4. cameras = np.load(cameras_path, allow_pickle=True)
  5.  
  6. # Select the first camera pose (cam_id = 0)
  7. cam_id = 0
  8. cam = cameras[cam_id]
  9. print(f"Using camera {cam_id} ({cam['name']}): q={cam['q']}, t={cam['t']}")
  10.  
  11. q_w2c = torch.tensor(cam['q']).float()
  12. t_w2c = torch.tensor(cam['t']).float()
  13. c2w = w2c_to_c2w(q_w2c, t_w2c)

但是最终渲染的结果不是很清楚,只能隐约看到一盆“盆景”的点云轮廓。

原因是此处使用的本就是稀疏点云,加上原图又非常大,每个点落在单一像素上,广袤的屏幕中只有零星的像素被点亮,肉眼根本看不清。

为了可视化更清晰,我们可以在较小的分辨率上(例如缩小 4 倍或 8 倍)进行渲染。

这就引出了一个非常关键的操作:如果改变了图像的渲染分辨率,必须等比例地修改相机的内参(Intrinsics)。

相机的内参包括焦距(f_x, f_y)和光心偏移(c_x, c_y)。它们与原图分辨率是绑定的。当我们缩小渲染画布时,这些参数必须同比缩放。

  1. def scale_intrinsics(W_target, H_target, W_source, H_source, fx, fy, cx, cy):
  2.     # 1. 计算宽高的缩放比例 (Scale Factors)
  3.     scale_x = W_target / W_source
  4.     scale_y = H_target / H_source
  5.  
  6.     # 2. 对焦距进行等比例缩放
  7.     scaled_fx = fx * scale_x
  8.     scaled_fy = fy * scale_y
  9.  
  10.     # 3. 对光心(主点)位置进行等比例缩放
  11.     scaled_cx = cx * scale_x
  12.     scaled_cy = cy * scale_y
  13.  
  14.     return scaled_fx, scaled_fy, scaled_cx, scaled_cy