Python 工程师对 3D 高斯溅射的介绍（第 2 部分）

理解并编码如何在 3D 高斯溅射中使用高斯

       欢迎来到雲闪世界现在开始讨论高斯！这是每个人最喜欢的分布。如果您刚刚加入我们，我们已经在第1 部分中介绍了如何根据相机的位置获取 3D 点并将其转换为 2D 。在本文中，我们将讨论高斯分布的高斯部分。我们将使用GitHub中的 part_2.ipynb 。

我们在这里要做出的一点小改动是，我们将使用透视投影，它利用的内部矩阵与上一篇文章中所示的不同。然而，当将点投影到 2D 时，这两种方法是等效的，我发现第 1 部分中介绍的第一种方法更容易理解，但是我们改变了我们的方法，以便在 Python 中尽可能多地复制作者的代码。具体来说，我们的“内部”矩阵现在将由此处显示的 OpenGL 投影矩阵给出，乘法顺序现在将是 points @ external.transpose() @ internal。

内部透视投影矩阵

       对于那些想了解这个新的内部矩阵的人（否则请随意跳过本段）r 和 l 是右侧和左侧的裁剪平面，本质上是相对于照片宽度可以看到哪些点，t 和 b 是顶部和底部裁剪平面。N 是近裁剪平面（点将被投影到该平面），f 是远裁剪平面。有关更多信息，及完整代码可联系博主。这还会返回标准化设备坐标（介于 -1 和 1 之间）中的点，然后我们将其投影到像素坐标。撇开题外话，任务保持不变，即取 3D 中的点并投影到 2D 图像平面上。但是，在本教程的这一部分中，我们现在使用高斯而不是点。

def getIntinsicMatrix(
    focal_x: torch.Tensor,
    focal_y: torch.Tensor,
    height: torch.Tensor,
    width: torch.Tensor,
    znear: torch.Tensor = torch.Tensor([100.0]),
    zfar: torch.Tensor = torch.Tensor([0.001]),,
) -> torch.Tensor:
    """
    Gets the internal perspective projection matrix
    
    znear: near plane set by user
    zfar: far plane set by user
    fovX: field of view in x, calculated from the focal length
    fovY: field of view in y, calculated from the focal length
    """
    fovX = torch.Tensor([2 * math.atan(width / (2 * focal_x))])
    fovY = torch.Tensor([2 * math.atan(height / (2 * focal_y))])
    
    tanHalfFovY = math.tan((fovY / 2))
    tanHalfFovX = math.tan((fovX / 2))

    top = tanHalfFovY * znear
    bottom = -top
    right = tanHalfFovX * znear
    left = -right

    P = torch.zeros(4, 4)

    z_sign = 1.0

    P[0, 0] = 2.0 * znear / (right - left)
    P[1, 1] = 2.0 * znear / (top - bottom)
    P[0, 2] = (right + left) / (right - left)
    P[1, 2] = (top + bottom) / (top - bottom)
    P[3, 2] = z_sign
    P[2, 2] = z_sign * zfar / (zfar - znear)
    P[2, 3] = -(zfar * znear) / (zfar - znear)
    return P

       3D 高斯图由 x、y 和 z 坐标以及相关的协方差矩阵组成。正如作者所指出的：“一种显而易见的方法是直接优化协方差矩阵 Σ 以获得表示辐射场的 3D 高斯。然而，协方差矩阵只有当它们是半正定的时候才具有物理意义。对于我们所有参数的优化，我们使用梯度下降法，这种梯度下降法不容易被限制来产生这样的有效矩阵，而更新步骤和梯度很容易产生无效的协方差矩阵。”¹

因此，作者使用协方差矩阵的分解，该分解将始终产生正半定协方差矩阵。具体来说，他们使用 3 个“尺度”参数和 4 个四元数，将其转换为 3x3 旋转矩阵 (R)。然后，协方差矩阵由下式给出

协方差矩阵方程，其中 R 表示由 4 个四元数导出的 3x3 旋转矩阵，S 是 3 个尺度参数

       请注意，在转换为旋转矩阵之前，必须对四元数向量进行归一化，才能获得有效的旋转矩阵。因此，在我们的实现中，高斯点由以下参数组成：坐标（3x1 向量）、四元数（4x1 向量）、比例（3x1 向量）以及与不透明度（splat 的透明度）相关的最终浮点值。现在我们需要做的就是优化这 11 个参数来获得我们的场景 — 很简单！

但实际上情况比这要复杂一些。如果你还记得高中数学，高斯在特定点的强度由以下公式给出：

点 x 处的高斯强度由平均值 (mu) 和协方差矩阵的逆给出

      然而，我们关心的是 2D 中 3D 高斯的强度，即在图像平面中。但你可能会说，我们知道如何将点投影到 2D！尽管如此，我们还没有讨论将协方差矩阵投影到 2D，所以如果我们还没有找到 2D 协方差矩阵，我们就不可能找到 2D 协方差矩阵的逆。

现在，这是最有趣的部分（取决于你如何看待它）。3D 高斯分布作者的一篇论文参考文献 EWA Splatting 精确地展示了如何将 3D 协方差矩阵投影到 2D。² 但是，这需要了解雅可比仿射变换矩阵，我们将在下面计算它。我发现代码在理解一个困难的概念时最有帮助，因此我在下面提供了一些代码，以举例说明如何从 3D 协方差矩阵转换为 2D。

def compute_2d_covariance(
    points: torch.Tensor,
    external_matrix: torch.Tensor,
    covariance_3d: torch.Tensor,
    tan_fovY: torch.Tensor,
    tan_fovX: torch.Tensor,
    focal_x: torch.Tensor,
    focal_y: torch.Tensor,
) -> torch.Tensor:
    """
    Compute the 2D covariance matrix for each gaussian
    """
    points = torch.cat(
        [points, torch.ones(points.shape[0], 1, device=points.device)], dim=1
    )
    points_transformed = (points @ external_matrix)[:, :3]
    limx = 1.3 * tan_fovX
    limy = 1.3 * tan_fovY
    x = points_transformed[:, 0] / points_transformed[:, 2]
    y = points_transformed[:, 1] / points_transformed[:, 2]
    z = points_transformed[:, 2]
    x = torch.clamp(x, -limx, limx) * z
    y = torch.clamp(y, -limy, limy) * z

    J = torch.zeros((points_transformed.shape[0], 3, 3), device=covariance_3d.device)
    J[:, 0, 0] = focal_x / z
    J[:, 0, 2] = -(focal_x * x) / (z**2)
    J[:, 1, 1] = focal_y / z
    J[:, 1, 2] = -(focal_y * y) / (z**2)

    # transpose as originally set up for perspective projection
    # so we now transform back
    W = external_matrix[:3, :3].T

    return (J @ W @ covariance_3d @ W.T @ J.transpose(1, 2))[:, :2, :2]

       首先，tan_fovY 和 tan_fovX 是半个视场角的切线。我们使用这些值来限制我们的投影，防止任何狂野的屏幕外投影影响我们的渲染。我们可以从 3D 到 2D 的变换中推导出雅可比矩阵，就像我们在第 1 部分中介绍的初始正向变换一样，但我为您省去了麻烦，并在上面展示了预期的推导。最后，如果您还记得，我们在上面转置了旋转矩阵以适应术语的重新排列，因此我们在返回最终协方差计算之前在倒数第二行转置回来。正如 EWA 溅射论文所述，我们可以忽略第三行和第三列，因为我们只关心 2D 图像平面。您可能会想，为什么我们不能从一开始就这样做？好吧，协方差矩阵参数将根据您从哪个角度观察它而变化，因为在大多数情况下它不会是一个完美的球体！现在我们已经转换到正确的视点，协方差 z 轴信息就没有用了，可以丢弃。

       假设我们有 2D 协方差矩阵，我们几乎能够计算出每个高斯对图像中任何随机像素的影响，我们只需要找到逆协方差矩阵。再次回想一下线性代数，要找到 2x2 矩阵的逆，你只需要找到行列式，然后对项进行一些重新排列。以下是一些代码，可帮助你完成该过程。

def compute_inverted_covariance(covariance_2d: torch.Tensor) -> torch.Tensor:
    """
    Compute the inverse covariance matrix

    For a 2x2 matrix
    given as
    [[a, b],
     [c, d]]
     the determinant is ad - bc

    To get the inverse matrix reshuffle the terms like so
    and multiply by 1/determinant
    [[d, -b],
     [-c, a]] * (1 / determinant)
    """
    determinant = (
        covariance_2d[:, 0, 0] * covariance_2d[:, 1, 1]
        - covariance_2d[:, 0, 1] * covariance_2d[:, 1, 0]
    )
    determinant = torch.clamp(determinant, min=1e-3)
    inverse_covariance = torch.zeros_like(covariance_2d)
    inverse_covariance[:, 0, 0] = covariance_2d[:, 1, 1] / determinant
    inverse_covariance[:, 1, 1] = covariance_2d[:, 0, 0] / determinant
    inverse_covariance[:, 0, 1] = -covariance_2d[:, 0, 1] / determinant
    inverse_covariance[:, 1, 0] = -covariance_2d[:, 1, 0] / determinant
    return inverse_covariance

        好了，现在我们可以计算图像中每个像素的像素强度了。但是，这样做非常慢而且没有必要。例如，除非协方差矩阵非常大，否则我们真的不需要浪费计算能力来弄清楚 (0,0) 处的 splat 如何影响 (1000, 1000) 处的像素。因此，作者选择计算他们所谓的每个 splat 的“半径”。如下面的代码所示，我们计算沿每个轴的特征值（请记住，特征值显示变化）。然后，我们取最大特征值的平方根以获得标准偏差测量值并将其乘以 3.0，这涵盖了 3 个标准偏差内的 99.7% 的分布。这个半径有助于我们找出 splat 接触的最小和最大 x 和 y 值。渲染时，我们只计算这些范围内像素的 splat 强度，从而节省了大量不必要的计算。很聪明，对吧？

def compute_extent_and_radius(covariance_2d: torch.Tensor):
    mid = 0.5 * (covariance_2d[:, 0, 0] + covariance_2d[:, 1, 1])
    det = covariance_2d[:, 0, 0] * covariance_2d[:, 1, 1] - covariance_2d[:, 0, 1] ** 2
    intermediate_matrix = (mid * mid - det).view(-1, 1)
    intermediate_matrix = torch.cat(
        [intermediate_matrix, torch.ones_like(intermediate_matrix) * 0.1], dim=1
    )

    max_values = torch.max(intermediate_matrix, dim=1).values
    lambda1 = mid + torch.sqrt(max_values)
    lambda2 = mid - torch.sqrt(max_values)
    # now we have the eigenvalues, we can calculate the max radius
    max_radius = torch.ceil(3.0 * torch.sqrt(torch.max(lambda1, lambda2)))

    return max_radius

       上述所有步骤都为我们提供了预处理场景，然后可以在渲染步骤中使用。总结一下，我们现在有了 2D 中的点、与这些点相关的颜色、2D 中的协方差、2D 中的逆协方差、排序的深度顺序、每个 splat 的最小 x、最小 y、最大 x、最大 y 值以及相关的不透明度。有了所有这些组件，我们终于可以开始渲染图像了！

感谢关注雲闪世界。（亚马逊aws和谷歌GCP服务协助解决云计算及产业相关解决方案）

 订阅频道(https://t.me/awsgoogvps_Host)
 TG交流群(t.me/awsgoogvpsHost)