光度立体最小二乘计算：过滤阴影像素及奇异矩阵报错求助

问题背景 #

我正在实现Woodham 1980年提出的朗伯光度立体算法，目前遇到了阴影像素导致的奇异矩阵问题：

• 现有输入：self.L是8个光照方向向量堆叠的数组，形状为(8, 3)；self.M是8张灰度图堆叠重塑后的数组，形状为(8, 9082318)，每一列对应一个像素在8种光照下的灰度值。
• 无阴影场景下，用self.N = np.linalg.lstsq(self.L.T, self.M.T, rcond=None)[0].T可以正常得到形状为(9082318, 3)的法向量数组，归一化后即可使用。
• 当前痛点：M中存在位置不固定的阴影像素，需要在计算最小二乘时剔除这些无效值（L维度保持不变），但尝试的掩码处理方法触发了numpy.linalg.LinAlgError: Singular matrix错误。

原工作函数 #

def _solve_l2(self):
    """
    Lambertian Photometric stereo based on least-squares Woodham 1980
    :return: None
    Compute surface normal : numpy array of surface normal (p × 3)
    """
    self.N = np.linalg.lstsq(self.L.T, self.M.T, rcond=None)[0].T
    print(self.N.shape)
    self.N = normalize(self.N, axis=1) # normalize to account for diffuse reflectance

尝试的参考代码 #

Ma = self.M.copy()
thresh = 300
Ma[self.M <= thresh] = 0
Ma[self.M > thresh] = 1
Ma = Ma.T
self.M = self.M.T
self.L = self.L.T
print(self.L.shape)
print(self.M.shape)
print(Ma.shape)
A = self.L
B = self.M
M = Ma
# 尝试基于掩码的最小二乘计算
rhs = np.dot(A.T, M * B).T[:,:,None] # n x r x 1 tensor
T = np.matmul(A.T[None,:,:], M.T[:,:,None] * A[None,:,:]) # n x r x r tensor
self.N = np.squeeze(np.linalg.solve(T, rhs)).T # transpose to get r x n

解决思路 #

1. 先修正阴影掩码的合理性

你当前设置的阈值thresh=300明显有问题——8位灰度图的取值范围是0-255，这个阈值会让几乎所有像素都被标记为阴影，导致有效样本数严重不足，直接引发矩阵奇异。

• 调整阈值：改用符合图像实际灰度范围的阈值，比如基于图像统计的自适应阈值：

  # 取灰度值的5%分位数作为阴影阈值，可根据实际情况调整
  thresh = np.percentile(self.M, 5)
  shadow_mask = self.M > thresh  # True表示有效像素，False表示阴影
  

• 过滤无效像素：每个像素的有效光照样本数必须≥3（因为要解3个法向量分量的未知数），否则对应的求解矩阵必然奇异。提前过滤这类像素：

  valid_count = np.sum(shadow_mask, axis=0)
  # 只保留有效样本数≥3的像素
  valid_pixel_mask = valid_count >= 3
  

2. 重构带掩码的最小二乘计算逻辑

原参考代码的张量操作逻辑不够清晰，建议针对每个像素独立处理（可用向量化方式优化效率）：

def _solve_l2_with_shadow(self):
    # 1. 生成合理的阴影掩码
    thresh = np.percentile(self.M, 5)
    shadow_mask = self.M > thresh
    valid_count = np.sum(shadow_mask, axis=0)
    valid_pixel_mask = valid_count >= 3
    
    # 2. 初始化法向量数组，无效像素先设为0
    self.N = np.zeros((self.M.shape[1], 3))
    
    # 3. 处理有效像素
    for pixel_idx in np.where(valid_pixel_mask)[0]:
        # 提取该像素的有效光照向量和灰度值
        valid_lights = shadow_mask[:, pixel_idx]
        L_sub = self.L[valid_lights, :].T  # 形状(3, k)，k是有效光照数
        M_sub = self.M[valid_lights, pixel_idx]  # 形状(k,)
        
        # 用lstsq或伪逆求解，避免奇异矩阵问题
        try:
            normal = np.linalg.lstsq(L_sub.T, M_sub, rcond=None)[0]
        except np.linalg.LinAlgError:
            # 矩阵奇异时用伪逆求解
            normal = np.dot(np.linalg.pinv(L_sub.T), M_sub)
        
        self.N[pixel_idx] = normal
    
    # 4. 归一化有效像素的法向量
    self.N[valid_pixel_mask] = normalize(self.N[valid_pixel_mask], axis=1)

3. 奇异矩阵的通用处理方案

无论用哪种方式，遇到奇异矩阵时，不要直接用np.linalg.solve，改用以下两种方法：

• np.linalg.lstsq：本身就是最小二乘求解函数，能自动处理接近奇异的矩阵；
• np.linalg.pinv：计算矩阵的摩尔-彭罗斯伪逆，返回最小二乘意义上的最优解，完全兼容奇异矩阵场景。

4. 中间结果验证

在调试阶段，可以添加一些打印逻辑验证掩码和计算的正确性：

# 查看有效样本数的分布
print(f"有效样本数范围: {np.min(valid_count)} ~ {np.max(valid_count)}")
print(f"有效像素占比: {np.sum(valid_pixel_mask)/self.M.shape[1]:.2%}")
# 查看单个像素的有效数据
test_idx = np.where(valid_pixel_mask)[0][0]
print(f"测试像素有效光照数: {valid_count[test_idx]}")
print(f"测试像素有效光照向量:\n{self.L[shadow_mask[:, test_idx], :]}")

内容的提问来源于stack exchange，提问作者janb