【笔记】《Subpixel Photometric Stereo》的思路

这段时间真的好忙，周更啊什么的都停滞了。前几天又看了一圈谭平的关于如何提高光度立体成像法线分辨率的这个论文，看完也写了长长的笔记。

因此这篇只是笔记，而且由于我能力有限里面还有很多没有彻底搞懂的地方，所以里面自然也可能有很多理解错了地方。由于文章本身难度就很大，再加上笔记的语无伦次，可能会很难懂，就此发上来吧。

论文的封面就是下面的这张图，最后的效果是很棒的，图贴在最后面了。

PartA 像素子法线分布和凹凸性的计算依据

1. 利用阴影明暗关系求出每个面片(对应像素)的法线GMM分布(2维，xy，z固定是1，GMM及后面会用到的EM算法的相关理论可以在https://zhuanlan.zhihu.com/p/30483076了解)，这个GMM代表了表面上各个方向的法线分布概率。之所以这么写是由于分辨率低的像素内可能含有复杂的不符合高斯分布的法线结构，这些结构导致了最终光线对这个像素区域的反射不是镜面的，这个大像素的法线分布用GMM来代替，这被称为微表面理论。

1. 假设像素是由内部更多的子像素构成，每个子像素就是一个微表面，其法线分布符合高斯分布且使得每个子像素的反射都比较接近于镜面反射从而合成为大的GMM，我们的目标就是找到这个大的像素法线的GMM分布，从中分离出子像素可能的法线分布，每个子像素的法线分布则是普通的高斯分布
2. 首先基于上面对像素集合结构的设计，每个像素的反射方程就变成下面这样(式子即来自BRDF)：

1. l和v是光照方向和视角，p是法线分布，F是Fresnel菲涅尔反射，Ks是微结构里的阴影因子，Km是遮蔽因子，h是视角与光照方向的平分向量，ρ是表面反射微分辐射率和入射微分辐照度的比率，与图像的对于像素强度有关，这条反射方程是来自于基于物理着色BRDF(有关BRDF模型的简单介绍可以在https://zhuanlan.zhihu.com/p/21376124找到)的。目标就是恢复法线分布p(h)
2. 常用的菲涅尔函数

1. 在这里由于视角v不变只有光照方向l在变，所以首先Km就是不变的（遮蔽与视角有关）。而为了简化对GMM的估计且平时(光线角度不刁钻时)Ks和Km对反射的影响比法线本身小很多，所以在这把F和Ks都先当作是常量
2. 于是将把像素的法线分布用GMM代替，代替后就是下面的样子，显然G就是GMM的概率函数，为N个子法线分布组成：

1. 按照前面所说的常量假设，在这里A就先认为是一个常量，此模型来解法线分布GMM
2. 但是仅有法线分布还不够，虽然在平时Ks和Km对反射的影响比法线本身小很多，但是当光线聚焦方向大幅改变时Ks的影响就很大了。例如两个像素如下图的凹凸性相反时，它们仍可能是由相同的法线构成所以法线本身对反射的影响相同，但是由于凹凸性的关系它们对光照的反应不同，通过改变光的方向可以来计算其凹凸性。例如下图中大角度光照会使得凹表面的很多部分没有被照射。这个特性用来计算出表面的凹凸性以在下面优化排列。

PartB 恢复像素的子法线GMM分布

1. 前面说到可以暂时把公式里的A当作常量，由于我们将会有很多个不同方向光照下得到的图片，所以原则上是可以利用这么多图，视角方向，光照方向为参数通过最小二乘解得到GMM的参数的

1. rk是辐射强度，n是使用传统PS算出的主法线方向，O是各个图的反射信息，A·G是新预测，尽可能使得新预测误差低
2. 然而，如同[22]中所说，由于这个反射方程的非线性，最小二乘解是很不稳定不可靠的
3. 由于之前的推导，我们知道A与G是独立的，认为A是已知的，自然就可以算下面的式子

1. 上面的G就是一个定义在半球上的概率分布函数pdf，问题是如何通过不确定的ρ函数来估算G
2. 如果我们将像素当作镜面反射，通过不同的光照方向，可以算得当前光线与视角的确切平分线hk的值，由于镜面反射的亮度最大情况就是平分线h正好与法线n方向平行的情况，所以通过采样半球上的多个平分线，就可以进行一个确切的计算，缓缓逼近几个法线方向，这个就叫做hk的一个采样，得到的是GMM的一个个采样值。在很多次的采样下我们就可以还原出这个GMM

1. 但是由于采样时的不均匀性，稀疏区域的采样可能产生的误差会很大，我们需要对这个采样O^进行1/t(h)加权，稀疏区域误差加权比较小，稠密区域加权大。经过加权得到的才是想要的GMM。对样本加权后可以通过EM算法迭代计算这个GMM。这时分为两种情况：
2. 简单的情况即光线(平分线)在半球上均匀采样时，t(h)会变成恒定的t0(每块的面积都相同)，所以采样的权值就是O^k/t0，那加不加权也就无所谓了。然后使用下面的EM算法迭代求解

1. 由于接下来我们会认为这GMM里的子高斯分布概率都是相同的(是平分出来的法线)，所以αi为1/N
2. 解出来后A可以使用最小二乘解算得

1. 而复杂的情况即不均匀采样时得到的结果会有偏差，这是很普遍的问题，采样越密的地方对误差校准贡献越大。在非均匀时权值自然是O^k/t(h)，为了优化这个问题我们在半球上计算这个值时会让权值1/t(h)由此处h在Voronoi图(冯洛诺伊图/泰森多边形https://ww2.mathworks.cn/help/matlab/ref/voronoi.html)上的对应面积来设置。

PartC 恢复像素的子结构凹凸度

1. 在这一部分里，前面说到Km是不变的，但是在这里的计算下Ks就需要改变了
2. 前面说到了凸结构会比凹结构得到更多的照射，也就是说凹结构当光照角度提高时Ks会随之增加(斜率越大像素越暗)，而凸结构则没有那么大的影响
3. 由于具体数值不好解出，只想要了解结构凹凸性的我们需要利用光照变化下亮度变化的比值。由于前面已经得到了像素的法线GMM，又可以按照光照角度例如设置为80和60用O^k来代替ρ，相除得到A的比值。
4. 由于A中Km的变化忽略不计，F的变化也不大，A的比值就可以约等于Ks的比值，即：Agrazing/Anongrazing ≈ Ks(lnongrazing)/Ks(lgrazing)，再使用一个阈值来处理这个比值就可以知到像素微结构的凹凸度。在此基础上计算出凹凸度图r

PartD 构造每个像素所对应的子法线排列结构

1. 由于上一部分求得了每个像素的子法线GMM，我们可以在这个法线分布中切割出想要的子法线成分。那么总的目标就是如何将这些切割出来的子法线安排到恰当的位置上
2. 直接安排是不现实的，所以需要简化问题。主要想法是解出一些固定的2*2法线排列，这些排列称为基元(texton)，然后直接安放在适合的像素上，相似的像素会对应相同的基元，划分像素并找到的可对应某一类像素的基元就称为解基元IL(solution texton)
3. 具体来说，首先前面已经给每个像素分割出了2*2个法线，这样的法线可以有(2*2)！=24种排列，由于法线都是2维的，将每一种排列拼成一个1*8(2*2*2)的基元向量
4. 将所有这样的向量放置在分布空间中，同样使用EM算法拟合出它们的基元GMM
5. 在以一定的阈值过滤GMM删去一些过于接近的高斯成分后，ＧＭＭ的每个单高斯分布成分就看做一类基元簇，我们使用这些基元的平均值来代表这个簇
6. 可想而知这样每个基元簇就是由很多排列相近而分属于不同像素的基元组成
7. 其中有一些基元是原先同个像素中被划分出来的，也就是有相同的子法线但排列不同，如今属于不同的基元簇
8. 为了最小化这个几何描述，我们使用投票机制，两个不同基元簇间的基元进行投票，在一个阈值的过滤下来自相同像素达到一定数量的基元簇被划分为一个等价类
9. 这样处理后来自相同像素且排列相近的基元就被划入了同个等价类中，这些基元对于它们自身的像素比较有代表性且有很多共通相似之处，所以我们可以用能代表这个等价类的基元(解基元)来直接代替图像的像素，问题就转换为如何在等价类中找到可以代表这个等价类的解基元

PartE 寻找各个像素对应的解基元用到的能量函数

1. 首先法线的排列(基元)选择有三个相关约束：各个像素的凹凸度，可积性限制和基元描述距离
2. 凹凸度：对于一个闭合曲线，其法线积分大于0则曲面是凸的，积分小于0则曲面是凹的。所以这里取每个像素的解基元(abcd是这四个子法线，z是1)，将它们认为是这个闭合曲线表面的法线来积分，公式如下：

1. 然后利用partC算得的凹凸度图r来选择下面的映射函数计算凹凸度的能量函数：

1. 最大化E0就使得法线排列更偏向于符合凹凸度状态，λ是凹凸的阈值
2. 可积性：再由于闭合曲面的积分为0，可积性限制也就是像素区域看作一个闭合曲面所以积分也应该是0，即：

1. 最大化下面的能量函数E1使得法线排列更偏向于符合曲面闭合：

1. 基元描述距离：对于每一对相邻的像素和它们各自的基元s和基元t，它们的交界处也是闭合曲面，所以可积性限制一样要在两个像素的交界处发生作用，即：

1. 而且交界处的排列也需要尽量符合求得的解基元，所以还要得到交界处最符合的基元的概率，先将相邻的四个子像素进行串联得到新排列，然后求得这个排列在之前的基元GMM中对应的概率最高的基元成分的概率maxP
2. 这样得到下面的能量函数E2，最大化E2就可以尽可能使得交界处的子像素排列一样符合可积性和解基元限制

1. 原则上这个交界判断还需要发生在例如四个角的重叠处，但是为了防止更大的复杂度而不考虑了（这里似乎只判断了水平方向）

1. 将这些能量函数串起来得到我们需要优化的总能量如下，目标就是找到一组能使这个能量函数平均值最大化的解基元IL，接下来用MCMC来计算这些IL

PartF 使用置信传播来初始化MCMC的状态

1. 有了上面的能量函数后，接下来要用MCMC(马尔科夫链蒙特卡洛方法)来进行优化。在优化之前，我们知道MCMC的效果和速度与初始值有关，所以首先我们可以用置信传播来为MCMC初始化一个合适的初值以加快收敛提高效果
2. 对于这个置信传播，一开始的时候要构造一张无向图，图中的每个节点是一个等价类，包含了很多基元，两个节点只有当其包含了某像素相互是四邻居状态才会相连
3. 而每个节点需要遵守前面的可积性限制和闭合限制E0和E1，得到节点能量函数，其中IL是等价类目前的代表基元，v是等价类，s就是等价类中包含的像素

1. 相应的，有边的节点中的像素需要符合前面提到的交界处能量函数E2，即下面的边能量函数。这里的E2'的意思是先抛去maxP也就是暂不考虑交界处是否符合现有解基元来计算的E2

1. 然后利用这个无向图中的能量函数置信传播就可以初始化MCMC的参数

PartG 使用MCMC计算和优化

1. 到了这一步我们可以有一个总的能量函数了(MCMC算法我并不明白，这段可能有很多问题)

1. 这个是前面的三个能量的积，目标是使得其IL的平均值最大，也就是PartE结尾的那个函数，那么在这里把Etotal简写为π，参数为基元，一种基元的选择为马尔可夫链的状态x，于是状态转移便写为这样：

1. 由于马尔可夫链的状态转移依赖于概率转移提议函数q(x*|xk)，在这里构造了两个提议函数，本地提议ql，全局提议qw，混合得到的函数如下：

1. αl和αw是权值，总和为1。以上两个提议函数都是基于首要独立采样器MIS的
2. 由于实际上之前初始化中让像素属于一个等价类的算法由于诸如阈值选择和图像噪声的原因是不可靠的，它们很可能都属于一个拥有最多基元的等价类，所以使用本地提议函数一轮随机选择一个像素改变其等价类关联并在新的等价类中找到适合的基元，以求达到更好的结果
3. 而这个改变关联的概率算法中，首先需要按照下面的概率，通过一个像素的能量概率选择一个像素

1. 接着随机选择一个等价类，由于目前的基元都属于一个基元簇GMM中，也就是这个基元本身也是有出现概率的，相乘后再计算采纳新基元的概率，采纳后再计算总能量看是否有改进

1. 又为了使马尔可夫链更容易跳出局部最优，全局提议原理类似局部提议，先是按照下面的能量概率选择一个等价类：

1. 再按照概率改变这个等价类的基元为一个新的基元

1. 解了两种提议后，它们的权值理应有条件发生作用。全局的改变需要在局部改变了多次却仍然无法进行进一步优化时采用，也就是αl可以为(其中iter是局部提议采用却没有提高结果的次数计数器，σ就是前面说到的局部提议被采用的固定概率(?)，设为100)：

1. αw与αl的和为1
2. 到此为止所有步骤都结束了，剩下的就是多次迭代求解。这种分辨率提升的方法计算量很大，所以无法提升到很高的程度，2*2或3*3是比较合适的

PartH 效果