利用高频光源快速分离直接光照和全局光照

直接光照与全局光照的定义

如上图所示，在场景中光源发出光线经不同传播途径最终在 \(P\) 点反射进入到相机，这里 \(P\) 处进入相机的光照成分有如下几种：

光线 \(A\) ，光源直接照射到 \(P\) 处，然后反射到相机中。这种光照被称作直接光照（Direct）。
光线 \(B\) ，光源先直接照射到表面的其他地方后，反射到 \(P\) 处，然后反射到相机中。这种光照被称作互反射光照（Interrelection）。如桌上的几何体表面相互反射光照。
光线 \(C\) ，光源先直接照射到表面的其他地方后，因为光线穿透到表面下方，经散射后，从 \(P\) 处射出，最后到相机中。这种光照被称作亚表面散射光照（Subsurface）。如人类的皮肤看到被光照打亮而显得红润。
光线 \(D\) ，光源先直接照射到参与介质中，经介质散射到 \(P\) 处，然后反射到相机中。这种光照被称作体散射光照（Volumetric）。如水中的鱼因为水的散射光而看得清。
光线 \(E\)，光源先直接穿过半透明表面，然后投射到\(P\) 处，最后反射到相机中。这种光照被称作漫射光照（Diffusion）。如磨砂玻璃后看到的人影。

下面是一些案例：

在上述光照中，只有光线 \(A\) 这种被称作直接光照，而其余的都称作全局光照。

一般来说只有被光源直接照射才会产生直接光照，没有被照射则可以认为不产生直接光照。而全局光照始终都有，全局光照来自场景中其他非当前点所贡献的光照。而对于场景中其他非当前点，其自身照向当前点的光照成分同样可以继续分为直接光照和全局光照。

全亮光照模型

现在我们对上述场景进行建模，假设光源照亮整个表面，对于表面上的一小块面片 \(i\) ，其被相机 \(c\) 拍摄到，那么进入相机光照可记作 \(L[c,i]\)，令该光照拆分为直接光照分量 \(L_{d}[c,i]\) 和全局光照分量 \(L_{g}[c,i]\)，则有：

\[L[c,i] = L_{d}[c,i] + L_{g}[c,i] \]

根据前面提到的光照类型，我们知道全局光照分量包含多种，为了简化研究，我们假设表面不透明，不存在亚表面散射，仅有互反射。那么面片 \(i\) 的全局光照分量来自于其他所有面片对他互反射。

假设整个表面被划分为 \(N\) 个小面片，被光源直接照亮的面片个数为 \(M\) 。那么可以定义除面片 \(i\) 之外的其他面片的集合为 \(P = \{j | 1 \le j \le N, j\ne i\}\) 。

取集合 \(P\) 中的任意一个面片 \(j\) 研究其对面 \(i\) 所产生的互反射光照成分，令 \(L[i,j]\) 为面片 \(j\) 在指向面片 \(i\) 方向上的光照，\(A[i,j]\) 为包含面片 \(i\) 和面片 \(j\) 的集合位置关系与面片 \(i\) 的BRDF属性的复杂函数，描述面片 \(j\) 对面片 \(i\) 的光照贡献权重。于是全局光照分量 \(L_{g}[c,i]\) 可表示为：

\[L_g[c,i] = \sum_{j \in P} A[i,j] \cdot L[i,j] \]

而对于面片 \(j\) 所产生的光照 \(L[i,j]\) 又可以根据直接光照和全局光照的概念，继续拆分为直接光照分量 \(L_{d}[i,j]\) 和全局光照分量 \(L_{g}[i,j]\)：

\[L[i,j] = L_{d}[i,j] + L_{g}[i,j] \]

将上式代入前面的式子，则有：

\[\begin{align*} L_g[c,i] &= \sum_{j \in P} A[i,j] \cdot L[i,j]\\ &= \sum_{j \in P} A[i,j] \cdot(L_{d}[i,j] + L_{g}[i,j])\\ & = \sum_{j \in P} A[i,j] \cdot L_{d}[i,j] + \sum_{j \in P} A[i,j] \cdot L_{g}[i,j]\\ \end{align*} \]

为了方便表达，令：

\[\begin{align*} L_{gd}[c,i] &= \sum_{j \in P} A[i,j] \cdot L_{d}[i,j]\\ L_{gg}[c,i] &= \sum_{j \in P} A[i,j] \cdot L_{g}[i,j] \end{align*} \]

则：

\[L_g[c,i] = L_{gd}[c,i] + L_{gg}[c,i] \]

于是完整的光照模型可以表示为：

\[\begin{align*} L[c,i] &= L_{d}[c,i] + L_{g}[c,i] \\ &= L_{d}[c,i] + L_{gd}[c,i] + L_{gg}[c,i] \end{align*} \]

其中：

\(L_{d}[c,i]\) ，光源直接照射面片 \(i\) ，反射至相机 \(c\) 的光照。
\(L_{gd}[c,i]\)，场景中其他面片 \(j(j\ne i)\)，被光源直接照射产生 \(L_d[i,j]\) 后，再间接照射到面片 \(i\)，反射到相机 \(c\) 的光照。
\(L_{gg}[c,i]\)，场景中其他面片 \(j(j\ne i)\)，自身接收全局光（高阶互反射 / 散射） \(L_g[i,j]\) 后，再间接照射到面片 \(i\)，反射到相机 \(c\) 的光照。

高频光源下的光照模型

在前面的光照模型中，我们的光源可以认为是一个低频光源，场景中的表面被全部光源直接照射点亮。若换成可以控制光源某些位置点亮或者不点亮，那么场景中的表面被光源照射的位置才产生直接光照，没被照射的位置则不产生直接光照。以条纹投影为例（即投射高频图案），一般来说亮条纹投射的表面被认为是光源直接照射的，即点亮，而暗条纹则没有被光源直接照射，即未点亮。

在高频光源下，假设整个表面被划分为 \(N\) 个小面片，当光源全亮时（即低频），被光源直接照亮的面片个数为 \(M\) 。而当控制光源投射高频图案后，被点亮的面片比例为 \(\alpha\) ，那么此时被光源直接照亮的面片个数为 \(\alpha M\) 个。令被点亮的面片集合为 \(Q=\{k|k\in N, lit(k) = 1\}\) ，其中 \(lit(k) = 1\) 表示该面片被照亮。

此外因为投影仪通常有漏光的情况发生，即使光源未直接照射面片，也会漏出部分光照射到面片上，假设这部分的亮度是正常照射的亮度的 \(b(0\le b\le1)\) 倍。

高频图案

对于任意面片 \(i\) ，假设其被高频光源直接点亮，那么其直接光照分量 \(L_{d}[c,i] \ne 0\) ，与全量时并无区别：

\[L_{d}^{+}[c,i] = L_{d}[c,i] \]

对于分量 \(L_{gd}[c,i]\) ，其表示的是场景中其他面片 \(j(j\ne i,j\in P)\)，被光源直接照射产生 \(L_d[i,j]\) 后，再间接照射到面片 \(i\)，反射到相机 \(c\) 的光照。因此只有被高频光源直接照射时（\(j\in Q\)），才有 \(L_d^{+}[i,j]=L_d[i,j]\) ，即被点亮，若未被点亮（\(j \notin Q\)），则因为漏光有 \(L_d^{+}[i,j] = bL_d[i,j]\) 。

因此若将高频光照下的该分量记作 \(L_{gd}^{+}[c,i]\)，则有：

\[\begin{align*} L_{gd}^{+}[c,i] &= \sum_{j \in Q} A[i,j] \cdot L_{d}[i,j] + \sum_{j \notin Q} A[i,j] \cdot bL_d[i,j] \\ &= \sum_{j \in Q} A[i,j] \cdot L_{d}[i,j] + b\sum_{j \notin Q} A[i,j] \cdot L_d[i,j] \\ \end{align*} \]

若 \(A[i,j]\) 和 \(L_d[i,j]\) 在空间平滑，高频均匀采样，被点亮的面片比例为 \(\alpha\) ，且面片全亮时的分量 \(L_{gd}[c,i]\) 表达式为：

\[L_{gd}[c,i] = \sum_{j \in P} A[i,j] \cdot L_{d}[i,j] \]

于是：

\[\begin{align*} \sum_{j \in Q} A[i,j] \cdot L_{d}[i,j] = \alpha \sum_{j \in P} A[i,j] \cdot L_{d}[i,j]\\ \sum_{j \notin Q} A[i,j] \cdot L_{d}[i,j] =(1-\alpha ) \sum_{j \in P} A[i,j] \cdot L_{d}[i,j]\\ \end{align*} \]

所以上式可进一步改写为：

\[\begin{align*} L_{gd}^{+}[c,i] &= \alpha L_{gd}[c,i] + b(1-\alpha) L_{gd}[c,i]\\ &= (\alpha + b(1-\alpha)) L_{gd}[c,i] \end{align*} \]

对于分量 \(L_{gg}[c,i]\) ，其表示的是场景中其他面片 \(j(j\ne i)\)，自身接收全局光（高阶互反射 / 散射） \(L_g[i,j]\) 后，再间接照射到面片 \(i\)，反射到相机 \(c\) 的光照。而 \(L_g[i,j]\) 是全局光照分量，全局分量与场景中的总光照功率成正比，在考虑漏光的条件下，总光照功率来自两部分，即点亮区域和未点亮区域，点亮区域的面片比例为 \(\alpha\) 亮度为 1，而未点亮区域的面片比例为 \(1-\alpha\)，亮度为点亮亮度的 \(b\) 倍，于是有：

\[\begin{align*} L_g^{+}[i,j] &= \alpha L_g[i,j] + b(1-\alpha)\alpha L_g[i,j]\\ & = (\alpha + b(1-\alpha)) L_g[i,j] \end{align*} \]

对应的高频下的分量为：

\[\begin{align*} L_{gg}^{+}[c,i] &= \sum_{j \in P} A[i,j] \cdot L_{g}^{+}[i,j]\\ &= \sum_{j \in P} A[i,j] \cdot (\alpha + b(1-\alpha)) L_g[i,j]\\ &= (\alpha + b(1-\alpha))\sum_{j \in P} A[i,j] \cdot L_g[i,j]\\ &= (\alpha + b(1-\alpha)) L_{gg}[c,i] \end{align*} \]

因此对于面片 \(i\) 来说，在高频光源下，来自其他面片产生的全局光照为

\[\begin{align*} L_g^{+}[c,i] &= L_{gd}^{+}[c,i] + L_{gg}^{+}[c,i]\\ &= (\alpha + b(1-\alpha)) L_{gd}[c,i] + (\alpha + b(1-\alpha)) L_{gg}[c,i] \\ &= (\alpha + b(1-\alpha))( L_{gd}[c,i] + L_{gg}[c,i]) \\ &= (\alpha + b(1-\alpha)) L_g[c,i] \end{align*} \]

因此相机在高频光源下获得的面片 \(i\) 处的全部光照为：

\[\begin{align*} L^{+}[c,i] &= L_{d}^{+}[c,i] + L_{g}^{+}[c,i]\\ &= L_{d}[c,i] + (\alpha + b(1-\alpha)) L_{g}[c,i] \\ \end{align*} \]

至此，面片 \(i\) 被高频光源的图案点亮时的光照方程推导完毕。

高频互补图案

若继续投射一副互补的图案，场景中的倍点亮的面片个数为 \((1-\alpha)M\) 个，那么此时面片 \(i\) 将不被点亮，但因为投影仪漏光，此时面片 \(i\) 的直接光照分量 :

\[L_{d}^{-}[c,i] = bL_{d}[c,i] \]

类似地，全局分量只需将前述公式中的 \(\alpha\) 替换为 \(1-\alpha\) ：

\[\begin{align*} L_{gd}^{+}[c,i] &= (\alpha + b(1-\alpha)) L_{gd}[c,i]\\ L_{gg}^{+}[c,i] &= (\alpha + b(1-\alpha)) L_{gg}[c,i] \end{align*} \]

则有：

\[\begin{align*} L_{gd}^{-}[c,i] &= (1-\alpha + b\alpha) L_{gd}[c,i]\\ L_{gg}^{-}[c,i] &= (1-\alpha + b\alpha) L_{gg}[c,i] \end{align*} \]

于是：

\[\begin{align*} L_g^{-}[c,i] &= L_{gd}^{-}[c,i] + L_{gg}^{-}[c,i]\\ &= (1-\alpha + b\alpha) L_{gd}[c,i] + (1-\alpha + b\alpha) L_{gg}[c,i] \\ &= (1-\alpha + b\alpha)( L_{gd}[c,i] + L_{gg}[c,i]) \\ &= (1-\alpha + b\alpha) L_g[c,i] \end{align*} \]

因此相机在高频光源下获得的面片 \(i\) 处的全部光照为：

\[\begin{align*} L^{-}[c,i] &= L_{d}^{-}[c,i] + L_{g}^{-}[c,i]\\ &= bL_{d}[c,i] + (1-\alpha + b\alpha) L_g[c,i] \end{align*} \]

其他类型全局光照

上述推导过程同样适用于其他类型的全局光照：

直接光照与全局光照分离

联立上述两幅互补高频光源投射所得到的方程：

\[\left\{ {\matrix{ {L^{+}[c,i] = L_{d}[c,i] + (\alpha + b(1-\alpha)) L_{g}[c,i]} \cr {L^{-}[c,i] = bL_{d}[c,i] + (1-\alpha + b\alpha) L_g[c,i]} \cr } } \right. \]

于是可推导出：

\[\left\{ \begin{align*} L_{d}[c,i] &= \frac{\bigl((1-\alpha)+b\alpha\bigr)L^{+}[c,i] - \bigl(\alpha + b(1-\alpha)\bigr)L^{-}[c,i]}{(1-b)(1+b)(1-\alpha)}, \\[8pt] L_{g}[c,i] &= \frac{L^{-}[c,i] - b\,L^{+}[c,i]}{(1-b)(1+b)(1-\alpha)}. \end{align*} \right. \]

为分离直接光照 \(L_d\) 与全局光照 \(L_g\) ，通常取 \(\alpha = \frac{1}{2}\) ，可使得两幅图亮度均衡，采样最稳定。而 \(b\) 通常根据光源点亮和未点亮的功率比来确定。此时：

\[\left\{ \begin{align*} L_{d}[c,i] &= \frac{L^{+}[c,i] - L^{-}[c,i]}{1-b}, \\[8pt] L_{g}[c,i] &= \frac{2}{1+b}\cdot\frac{L^{-}[c,i] - b\,L^{+}[c,i]}{1-b}. \end{align*} \right. \]

通常来说 \(L^+ \ge L^-\)，所以有：

\[\left\{ \begin{align*} L_{max} &=L^+\\[8pt] L_{min} &=L^- \end{align*} \right. \]

若不考虑漏光，即 \(b=0\) ，且 \(\alpha = \frac{1}{2}\) ，可得：

\[\left\{ \begin{align*} L_{d} &= L_{max} - L_{min} \\[8pt] L_{g} &= 2L_{min} \end{align*} \right. \]

工程处理

在实际的面结构光中，通常会投射多张不同频率的条纹投影，理论上应该从最高频率的图像中去分离直接光照和全局光照，但实际处理时，这样做会导致解码图像出现明显的伪影，所以通常会取第二高频率的图像序列来计算，若采用不同的投影图案精度会更高。

假设采集到的结构光投影图案图片集合为 \(S=\{I_1,\cdots,I_k\}\)，那么对于图像中的的任意位置的像素 \(p\) ，其直接光照为分量为 \(L_d(p)\)，全局光照分量为 \(L_g(p)\)，则根据前述理论，当 \(\alpha = \frac{1}{2}\) 时，有：

\[\begin{align*} L_p^{+} = \max_{0\lt i \le k} I_i(p),\\ L_p^{-} = \min_{0\lt i \le k} I_i(p). \end{align*} \]

于是：

\[\begin{align*} L_d(p) = \dfrac{L_p^{+} - L_p^{-}}{1-b},\\[6pt] L_g(p) = 2\dfrac{L_p^{-} - bL_p^{+}}{1-b^2}. \end{align*} \]

其中 \(b \in [0,1)\) ，通常取 \(b = 0.3\) 。

可用的光源类型

上述分离方法可以有多种光源类型，包括：

效果展示：

算法验证

import cv2
import numpy as np

def estimate_direct_light(images, b):
    """
    @brief 直接光照分量估计
    @param[in]   images  投影仪投射的一组图像中第二高频的编码图像（列表，元素为单通道灰度图）
    @param[in]   b       用于模拟关闭状态下的投影仪像素发出的光量(b∈[0,1])
    @return      2通道的图像,分别为直接光照分量图和全局光照分量图
    """
    # 最多使用10张图计算

    count = len(images)
    COUNT = count

    # 无图像输入，返回空矩阵
    if count < 1:
        return np.array([])
    
    print(" --- estimate_direct_light START ---")
    
    # 限制最大使用图像数量
    if count > COUNT:
        count = COUNT
        print(f"WARNING: Using only {COUNT} of {len(images)} images")
    
    # 检查所有图像是否为8位单通道灰度图
    for img in images:
        if len(img.shape) != 2 or img.dtype != np.uint8:
            print("Gray images required")
            return np.array([])
    
    # 获取图像尺寸
    height, width = images[0].shape
    # 初始化2通道结果图：通道0=直接光照，通道1=全局光照
    direct_light = np.zeros((height, width, 2), dtype=np.uint8)
    
    # 计算系数
    b1 = 1.0 / (1.0 - b)
    b2 = 2.0 / (1.0 - b * b)
    
    # 逐像素计算（等效C++逻辑）
    for h in range(height):
        # 取出每一行的像素指针
        rows = [img[h, :] for img in images[:count]]
        
        for w in range(width):
            # 获取当前像素在所有图像中的亮度值
            pixels = [row[w] for row in rows]
            Lmax = max(pixels)
            Lmin = min(pixels)
            
            # 计算直接光照和全局光照（四舍五入取整）
            Ld = int(b1 * (Lmax - Lmin) + 0.5)
            Lg = int(b2 * (Lmin - b * Lmax) + 0.5)
            
            # 完全复刻C++赋值逻辑
            if Lg > 0:
                direct_light[h, w, 0] = Ld
                direct_light[h, w, 1] = Lg
            else:
                direct_light[h, w, 0] = Lmax
                direct_light[h, w, 1] = 0
    
    print(" --- estimate_direct_light END ---")
    return direct_light

def read_video_all_frames_rgb(video_path):
    """
    读取视频所有帧，不转灰度，分离出 R、G、B 三个通道的图像列表
    :param video_path: 视频路径
    :return: r_frames, g_frames, b_frames  三个分别存放 R、G、B 通道的灰度图列表
    """
    cap = cv2.VideoCapture(video_path)
    r_frames = []
    g_frames = []
    b_frames = []

    if not cap.isOpened():
        print("错误：无法打开视频！")
        return r_frames, g_frames, b_frames

    print("正在读取视频所有帧并分离 RGB 通道...")
    while True:
        ret, frame = cap.read()
        if not ret:
            break  # 视频读取完毕

        # OpenCV 读取的图像顺序是 BGR
        b, g, r = cv2.split(frame)

        r_frames.append(r)
        g_frames.append(g)
        b_frames.append(b)

    cap.release()
    print(f"视频读取完成！共载入 {len(r_frames)} 帧，已分离 R、G、B 三通道\n")
    return r_frames, g_frames, b_frames

def calc_single(gray_list,b = 0.1):
    if len(gray_list) == 0:
        print("未读取到任何帧")
        exit()
    
    # 2. 调用直接光照估计函数
    result = estimate_direct_light(gray_list, b)
    
    # 3. 拆分输出
    direct_img = result[:, :, 0]
    global_img = result[:, :, 1]
    return direct_img, global_img

def calc_rgb(r_frames, g_frames, b_frames,b=0.1):
    r_direct_img, r_global_img = calc_single(r_frames,b)
    g_direct_img, g_global_img = calc_single(g_frames,b)
    b_direct_img, b_global_img = calc_single(b_frames,b)
    color_direct_img = cv2.merge([b_direct_img, g_direct_img, r_direct_img])
    color_global_img = cv2.merge([b_global_img, g_global_img, r_global_img])
    return color_direct_img,color_global_img

def run_rgb():
    VIDEO_PATH = r".\validation_400x400.avi"
    b = 0.3
    r_frames, g_frames, b_frames = read_video_all_frames_rgb(VIDEO_PATH)
    color_direct_img,color_global_img = calc_rgb(r_frames, g_frames, b_frames,b)
    cv2.imwrite("direct_light.png", color_direct_img)
    cv2.imwrite("global_light.png", color_global_img)
    
if __name__ == "__main__":
    run_rgb()