【paper reading】Uncalibrated Photometric Stereo under Natural Illumination

经典的photometric stereo方法对于实验设置主要有两个假设：1.已标定方向和亮度的平行光；2.朗伯表面的物体反射。近些年来，很多研究都致力于放宽这两方面的假设使ps方法能够更加适用于实际应用。其中对于光照的假设，几年前很多研究使用未标定或者半标定（只标定光源位置，不标定光源亮度）的点光源，现在的模型则更多的使用二阶球谐函数近似的自然光照模型从而应对更加接近现实世界的场景。同样，本文的方法主要是为了对第一种假设进行放宽，提出了一种对朗伯物体在未标定自然光下的表面恢复方法。

1.简介

整个方法基于分治的思想，将图像上的整个物体分解成一个个的小邻域，在每个小邻域上对物体表面进行恢复得到有歧义的表面法向，然后将所有的块拼接起来得到最终完整的法向图同时消除歧义。之所以分解到小邻域来进行法向恢复是基于一个关键理由：在一个足够小的平整表面块上，为表面点提供照明的自然光照的半球面相差不大，可以近似认为它们的照明情况一样，从而通过将复杂的可见半球面上的自然光照转换成等价的平行光照，从而在每个小邻域上整个问题从未标定自然光的photometric stereo问题转换成了更加便于处理的未标定平行光的photometric stereo问题。不过这里的问题是在每个邻域上计算出的法线方向所选择的参考坐标系是不同的，即它们存在3自由度的旋转的歧义。直接将所有邻域通过它们重合的点拼接起来的话累积误差会越来越大，为了更加鲁棒的合并这些表面的法向块，作者通过类似计算最短路径的方法计算了整个表面上的法向角度距离矩阵，将所有的部分解合并成一个整体解。

2.等效方向光模型

整个方法是基于朗伯表面假设进行的，在忽略物体间自反射以及投射阴影这些因素的影响时，对于表面反射率为$\rho$，表面法向为$n=[n_x,n_y,n_z{]}^T\in {\mathbb{R}}^{3\times 1}$的表面点，它的光照模型可以表示成
$$b={\int }_{\Omega }\rho L(\omega )\max ((n^T\omega ),0)d\omega$$
$\omega$是在以$n$为正方向的自然照明上半球$\Omega$中的单位方向向量，$L(\omega )$则指的是这个方向的自然照明的强度。简单来说就是这个点最终的亮度是能够照到这个点的半球面上所有方向上的自然照明光对它照亮，积分起来的效果。
现在假设表面上有很多点它们的法向方向一样是$b_k$，那么它们的照明半球${\Omega }_k$也一样，因此都可以用上式来表示成
$$b_k=\rho n_k^T{\int }_{{\Omega }_k}L(\omega )\omega d\omega =\rho n_k^T{\bar{l}}_k$$
即整个自然光照可以等价成一个方向光${\bar{l}}_k$的形式，对于不同的法向方向$b_k$，有着不同的照明半球${\Omega }_k$，因此也有着不同的方向光${\bar{l}}_k$与它对应。
上面这些只是针对单个点来说的，现在假设在物体表面的一个小块${\mathcal{N}}_k$上，假设这个块足够小并且足够平整，上面的点的法向变化不大，即块上面任意两点的法向方向夹角小于一定值$⟨n_{k,i},n_{k,j}⟩<\delta$，那么这些点对应的照明上半球也几乎是一样的（如Fig.2中的黄色半球），也能转换为同一个方向光的形式。

为了验证这种想法，作者对一个标准球进行了实验，整个图像的大小是$256\times 256$，通过选取邻域大小为半径1，10，30的圆（在图上第一行用灰圈表示）。在三种环境照明下对球进行了渲染，然后计算不同大小的邻域中每个点对应的等效方向光，可以发现，在半径为10的情况下，三百多个像素中几乎所有的点的对应的方向光都很好的聚合在一起。

3. 法向估计

对于一个photometric stereo的图像序列来说，物体会在Q个不同的自然光照下拍摄Q张图像。我们将物体表面切割成彼此有重叠区域的很多个patch，$\mathcal{N}={\mathcal{N}}_1\cup {\mathcal{N}}_2\cup ...\cup {\mathcal{N}}_k$，假设一个patch中有K个点，将这些点的法向堆叠成一个大小为$K\times 3$的矩阵$N$，将这些点在Q张图像中的亮度堆叠成一个$K\times Q$的矩阵$B$，则有
$$B_{K\times Q}=N_{K\times 3}{L}_{3\times Q}$$
对于每一个patch，因为它很小，因此可以认为它中的像素点的albedo都是不变的，因此对于$\rho$，可以将它作为一个常数并入到后面的项中。

3.1 在每个patch中估计法向

经过之前的等效方向光转换，现在在每个patch中的自然光ps问题转换成了方向光下朗伯体的ps问题。将之前的矩阵B进行SVD分解得到$B=U\Sigma V^T$。在理想情况下，$\Sigma$是一个含有三个非0值的对角矩阵。在这里将$U\sqrt{\Sigma }$作为等效的法向矩阵$\tilde{N}$，将$\sqrt{\Sigma }{V}^T$作为等效的光照矩阵$\tilde{L}$。这里的法向和光照都不是真实的法向和光照，它们与真实的法向和光照存在一个$3\times 3$的歧义Q。之前也说过了在patch中可以假定$\rho$是不变的，因此这里作者将albedo归并到法向上，因此有$\parallel {\tilde{N}}_{i}Q\parallel ={\tilde{N}}_{i}Q{Q}^T{\tilde{N}}_i^T=\rho$，这里的${\tilde{N}}_i$是patch中第i个点的等效法向。经过证明添加了同一反照率这个条件之后，$3\times 3$的歧义可以转换为一个3秩的旋转歧义矩阵O，即
$$B=\hat{N}O{O}^T\hat{L}$$
因此我们就在每个patch中求出了所有点的法向方向以及反照率$\rho$，虽然是带有3自由度的歧义O的。而至于光照我们可以不去关心。

3.2 集成patch中的法向生成完整法向图

经过之前的计算我们得到了每个patch中每个点的带有旋转歧义的法向方向，虽然这个法向方向不是真实的法向方向，但是在同一个patch中由于它们的歧义矩阵都是一样的，因此两个点的法向夹角和在真实情况下的法向夹角是一样的。而patch与patch之间都有一些重合的点，因此下面就通过这个线索将许许多多的patch拼接起来。
为了达到这一目的，作者构建了一个角度距离矩阵$D\in {\mathbb{R}}^{P\times P}$，其中P是整个表面上所有像素点的数目，矩阵每一个值代表两个点的法向方向的距离。
首先对于每个patch中的所有点对，它们的法向距离可以直接计算它们的等效法向的反余弦值。而不同patch中的点由于它们的旋转歧义不同不能直接计算，而需要通过最短路径搜索的方法传播得到，如Fig.4所示。

在图4中，由于$n_1$和$n_2$在一个patch中，因此$⟨n_1,n_2⟩=arccos(n_1,n_2)$。但是由于$n_1$和$n_3$不在一个patch中，因此需要通过两个patch之间的公共部分，如$n_2$来传播，因此$⟨n_1,n_3⟩=⟨n_1,n_2⟩+⟨n_2,n_3⟩$。同样，经过这样的传播过程能够将角度距离越传越远，一直到$n_4$和$n_5$，但是可以看到所有的角度距离都是正值，越传播到远处数值会越大。为了解决这个问题并且尽量减小累积误差的影响。作者认为两个像素如果它们的intensity profile一样，就认为它们的法向方向一致。所谓的intensity profile就是将该像素处Q张图像中的RGB值堆叠成的一个矩阵。简单理解就是如果两个点在所有图像中的RGB值都一样，那么它们的法向也就一样。在图4中，假设$n_3$和$n_4$这两个点就是这种情况，于是$⟨n_3,n_4⟩=0$，因此$⟨n_1,n_5⟩=⟨n_1,n_3⟩+⟨n_3,n_4⟩+⟨n_4,n_5⟩=⟨n_1,n_3⟩+⟨n_4,n_5⟩$。除了同一patch中的点对的法向距离能够直接计算，其他点对都需要通过patch之间的传播来得到法向距离，这个传播路径很可能不唯一，并且得到很多不同的结果，最终的结果取所有可能传播路径中"最短"(最后的法向距离最小)的。

剩下的这一部分是我自己的理解。我觉得有必要贴一下原文，原文是这样的，

由于之前得到的距离矩阵D是稀疏的，需要将它进行完整。作者重新构建了一个同样大小的矩阵M，
$$M_{i,j}=\{\begin{array}{rl}\cos (⟨{\hat{n}}_i,{\hat{n}}_j⟩)& ,\mathrm{i}\mathrm{f}{D}_{i,j}\in \Omega \\ 0& ,\mathrm{i}\mathrm{f}{D}_{i,j}\in {\Omega }^c\end{array}$$
$\Omega$表示D中有值的部分，${\Omega }^c$表示D中无值的部分。这里注意两个式子的不同，我的式子和原论文中的式子在后面的条件中，我把$\Omega$和${\Omega }^c$对换了。我觉得，作者在文章中可能顺序写错了。要不然${\hat{n}}_i$和${\hat{n}}_j$指的是之前得到的每个patch中的带旋转歧义的法向量吗？如果是的话那这两个patch的旋转歧义都不同，其中的法向量进行余弦值计算有什么意义呢？
除此之外，作者定义了一个矩阵A作为理想情况下的没有空值的M矩阵，作者说理想情况下$A={\hat{N}}^T\hat{N}$且秩为3，这也证明了我的猜测。因此就通过下面的式子来求矩阵A，这里也贴一下原始文章：

通过求解式7来得到矩阵A，那么问题来了，式子中的矩阵E是什么，通篇没有看到任何矩阵E的描述，就在这个式子里冒出了一个矩阵E。这里我的理解是E是和M以及D同样大小的矩阵，里面的值可以是任何值，只要在$\Omega$区域中的值是0就行了。那么整个式子可以这样理解，这个A矩阵就是在M矩阵有非0值时尽可能接近M的取值，而在M矩阵中0值的部分不做任何要求，并且最后的秩要是3。
那么最后通过求解出$A^{\ast }$，将它进行SVD分解得到$A^{\ast }=U\Sigma V^T$。因此整个物体表面的等效法向值就是
$$\hat{N}=U\sqrt{\Sigma }=\sqrt{\Sigma }{V}^T$$
注意这里的$\hat{N}$仍然不是最后的法向结果，但是它已经把所有的patch上的法向拼合在一起了，即所有的法向量之间的夹角都是正确值了。$\hat{N}$和真实的法向量还差一个旋转参数，后面通过强制可积后可以将这个旋转歧义转换成一个二义的凹凸歧义，最后通过人为手工选择物体是凸的还是凹的来得到最终结果。

4.实验部分

这里的实验部分都比较简单。首先利用GT的法线和自然光照在球面模型上通过选取不同大小的patch来验证估计的等效方向光与真实的等价方向光的方向的夹角大小。可以看到patch越小，估计出的等效光源越准确。在后面的实验中，作者都是用$3\times 3$大小的patch进行实验。

在如下的合成数据上进行实验

结果如下。在物体表面法线变化不剧烈，即比较平滑时，如球和小熊，最后的平均误差和中值误差都在10度以内，而兔子由于表面较粗糙，对最后的结果影响很大。

而由于使用的patch是$3\times 3$大小的，因此对于大块的反照率变化来说，整个方法表现出了良好的鲁棒性。

另外作者在真实数据集上进行了实验，最后的结果如下：


上个月只更新了一篇，实在是惭愧惭愧，不过确实写这些文章能增加我自己对文献的理解，但是也需要耗费大量的时间，以后应该会坚持下去的，但是如果有事比较忙的话更新频率可能会慢些。