光度立体法估计法线与反射率重建场景

1 从明暗恢复形状

从明暗恢复形状（Shape from Shading，SfS）是指从图像的明暗信息推断出物体表面几何形状的过程。这个问题假设光照条件已知，目标表面是光滑且均匀的，并且照明是单向的。其基本思想是根据目标表面对光照的反应，推断出表面法线，从而得到表面的三维形状。

1.1 渲染方程

渲染方程为：

我们将其简化： \(L_{\text{out}} = \rho \cdot L_{\text{in}} \cdot n^\top s = R(n)\)

• $ L_{\text{out}} $：表示出射光线
• $ \rho $：表示该点的漫反射率
• $ L_{\text{in}} $：表示入射光线
• $ n $：表示表面法线

1.2 梯度空间表示法

梯度空间表示法 Gradient Space Representation，使用梯度信息来表示物体表面的几何形状。

在梯度空间表示法中，给定光线 $ s $ 和观察到的反射率 $ R $，我们可以通过以下方式来计算法线 $ n $：

1. 计算 $ n^\top s = \cos(\theta) $，得到 $ \theta $ —— 法线 $ n $ 与光线 $ s $ 之间的夹角
2. 按照给定的角度 $ \theta $，从光线 $ s $ 开始，投影这个集合到 $ z = 1 $ 平面上，我们可以得到一个锥面曲线即反射率曲线 Iso-Reflectance Contour，它表示出了所有与 $ s $ 的夹角为 $ \theta $​ 的可能法线

通过不同的反射率 $ R $ ，我们可以得到反射率图（Reflectance Map）：

因此，在梯度空间表示法中，我们可以通过结合反射率图和其他几何信息来有效地参数化法线，并利用这些信息来推断物体表面的形状和材质。

2 光度立体法

光度立体法（Photometric Stereo）是通过在相同视点但采用不同（已知）点光源的多张图像来实现三维重建的。其需要对每个像素的法线和反射率进行估计。

2.1 采集 K 张图像

在相同的相机视点下，采集 K 张图像，每张图像使用不同的已知点光源。这些光源的位置和方向应该是已知的，并且在不同的图像中有所变化。

通过如下的反射率图可以看到，对于每个像素，我们必须要有3个反射率图才能确定一个真正的法线，所以，我们至少要有从3个不同的方向来的光源

但是要避免共线光源：当光度立体设置中使用的所有光源都是共线时（位于同一直线或平面上），所得的线性系统将变得秩亏。这使得不可能唯一地确定每个像素的表面法线。因此，光度立体无法提供准确的重建。

2.2 光度法线与反射率估计

对于每个像素，利用 K 张图像中的光照信息，通过解光度立体方程组来估计法线和反射率：

使用兰伯反射，并且入射光强度为 $ L_{\text{in}} = 1 $，那么图像的亮度 $ I $ 可以表示为：

对于给定的三个观测（相同的 $ v $，不同的 $ s $），我们可以将其表示为矩阵形式如下：

其中，$ I_1, I_2, I_3 $ 分别是三个观测得到的图像亮度，$ s_1, s_2, s_3 $ 分别是对应的三个光源方向向量，$ \rho $ 是漫反射率，$ n $​ 是法线向量。

通过使用更多的光源可以获得更好的结果（通过对测量进行平均）。通过最小二乘解法我们得到：

简单实现代码如下：

def compute_normals_albedo_map(imgs, mask, light_positions):
    """
    imgs np.array [k,h,w] np.float32 [0.0, 1.0]
    mask np.array [h,w] np.bool
    light_positions np.array [k,3] np.float32
    ---
    dims:
    k: number of images
    h: image height (num rows)
    w: image width (num cols)
    """

    S = light_positions
    I = imgs.reshape(imgs.shape[0], -1)

    # rho n = (S^T S)^-1 S^T I
    rho_n = np.linalg.inv(S.T @ S) @ S.T @ I
    # rho = ||rho_n||
    rho = np.linalg.norm(rho_n, axis=0)
    # n = rho_n / rho
    n = np.divide(rho_n, rho, out=np.zeros_like(rho_n), where=rho != 0)

    # mask out
    mask_flat = mask.flatten()
    n[:, ~mask_flat] = 0

    normals_unit = n.T.reshape(imgs.shape[1], imgs.shape[2], 3)
    rho = rho.reshape(imgs.shape[1], imgs.shape[2])

    assert normals_unit.shape == (imgs.shape[1], imgs.shape[2], 3)
    assert rho.shape == (imgs.shape[1], imgs.shape[2])

    rho = np.clip(rho, 0.0, 1.0)
    normals_unit = np.clip(normals_unit, 0.0, 1.0)

    return normals_unit, rho, mask

输出是：

• normals_unit：三维数组，表示每个像素点的单位法线向量。它的形状是 (imgs.shape[1], imgs.shape[2], 3)，其中 imgs.shape[1] 和 imgs.shape[2] 分别是图像的高度和宽度，3 表示每个像素点的法线有三个分量（x、y、z）
• rho：二维数组，表示每个像素点的漫反射率。它的形状是 (imgs.shape[1], imgs.shape[2])，与图像的大小相同，对应每一个像素的反射率

2.3 重新照亮场景

我们现在知道整个图像的像素法线和反照率，这使我们能够重新照亮场景（即人为地改变灯光位置）

def relight_scene(light_pos, normals_unit, albedo, mask):
    """
    light_pos np.array [k,3] np.float32
    mask np.array [h,w] np.bool
    ----
    dims:
    h: image height (num rows)
    w: image width (num cols)
    ----
    returns:
        imgs np.array [h,w] np.float32 [0.0, 1.0]
    """
    assert light_pos.shape == (3,)
    assert np.allclose(1.0, np.linalg.norm(light_pos))
    assert normals_unit.shape[-1] == 3
    assert len(normals_unit.shape) == 3

    img = albedo * (normals_unit @ light_pos)
    # mask out
    img[~mask] = 0
    img_norm = np.clip(img, 0.0, 1.0)

    assert np.all(
        np.logical_and(0.0 <= img_norm, img_norm <= 1.0)
    ), "please normalize your image to interval [0.0,1.0]"
    return img_norm

其接受灯光位置、像素法线、反射率和掩码作为输入，并返回重新照亮后的图像。

其实其中核心就是：根据之前的简化的渲染公式 $L_{\text{out}} = \rho \cdot L_{\text{in}} \cdot n^\top s$ ，其中设置入射光线为1，即可得：

img = albedo * (normals_unit @ light_pos)