基础知识

稠密点云重建的目的是在相机位姿已知的前提下，逐像素的计算图像中每一个像素点对应的三维点，得到场景物体表面密集的三维点云。

稠密点云重建的前提相机的姿态已知,且待重建的三维点具有图像一致性,即三维点位于物体的表面上而非物体内部，那么这个三维点根据相机的内外参投影到不同视角的图像上，以投影点为中心取出两个方形的小patch，那么它们所包含的场景应该是比较类似的.

极线搜索

两个视角的图像中，A视角中的一点在不同的深度的三维点投影到B视角中形成一条线，称为极线。当三维点的坐标未知，A图像中的一点，在B图像中可能匹配上对应的极线上的某点，在极线上找到B点的具体位置与相机的姿态之后，可以利用三角量量测的得到三位点的深度信息。

如何判断极线上的某点是否为A点的同名点，主要利用以下约束：光度一致性约束、可视性约束。

光度一致性约束

同一个空间的点在不同的视角的投影应当具有相同的光度，重建的核心在于恢复空间中具有光度一致性的点。（要求物体的表面为朗伯面，玻璃、镜子、水面不属于朗伯面）

光度一致性假设的度量方法

Sum of squared Differences(SSD):

ρSSD(f,g)=∣∣f−g∣∣2\rho_{SSD}(f, g)=||f-g||^{2}ρSSD​(f,g)=∣∣f−g∣∣2

Sum of Absolute Differences(SAD):

ρSAD(f,g)=∣∣f−g∣∣1\rho_{SAD}(f,g)=||f-g||_{1}ρSAD​(f,g)=∣∣f−g∣∣1​

Normalized Cross Correlation(NCC):

ρNCC(f,g)=(f−f‾)(g−g‾)δfδg\rho_{NCC}(f,g)=\frac{(f-\overline f)(g-\overline g)}{\delta_{f}\delta_{g}}ρNCC​(f,g)=δf​δg​(f−f​)(g−g​)​

其中f代表以A视图中某点的领域的像素组成的向量，g代表B视图中对应同名点邻域像素组成的向量

可视性约束

1、图像上的点不能被遮挡。 2、重建的点前面不能出现点。 3、不能出现在物体的内部。

多视角立体技术

基于体素的方法

规则划分

基于体素的方法等价于3D空间中的像素标定问题，在物体内部的点标记为1，在物体外部的点，标记为0，介于0-1之间的点就是物体的表面。

标记完之后，使用光度一致性约束和可视性约束来对物体表面的点进行计算，评估标定质量。基于体素的方法的缺点是计算量大，在分辨率很大的时候，计算量会增长很多。

不规则划分

在点云较为稀疏的区域采用分辨率较小的体素，在点云较为稠密与细节比较精确的的区域采用分辨率较高的体素，可以保证物体重建的精度，另一方面可以自适应调整四面体的大小，减少计算量。

基于体素的方法的优化

典型的MRF离散化优化问题

参考图像上的每个像素都分配一个标签（内部或者外部）（无向图的优化）

E(f)=Edata(f)+Esmoothness(f)+Evisibility(f)E(f)=E_{data}(f) + E_{smoothness}(f) + E_{visibility}(f)E(f)=Edata​(f)+Esmoothness​(f)+Evisibility​(f)

其中E(f)E(f)E(f)代表能量模型。

f表示无向图的所有顶点的向量，需要对每个顶点赋值

EdataE_{data}Edata​是光度一致性假设，在不同视图上观察到该点的纹理或者光度尽量相似。

EsmoothnessE_{smoothness}Esmoothness​平滑项约束，相邻的体素（四面体）的标签趋于一致。

EvisibilityE_{visibility}Evisibility​不能被别的东西遮挡，不能遮挡别的体素。

基于体素的方法的优小结

优点：1、生成规则的点云

 2、便于提取物体的平面

缺点：1、精度受到空间划分分辨率的影响

 2、难以处理精度高、规模大的场景

基于空间patch扩散的方法(pmvs)

1、假设空间中的3d矩形patch

2、通过一定规则的扩张方法，使得pathch覆盖物体表面

过程：

初始种子点生成采用sift、HOG等稀疏特征点，对这些点进行深度和法向量的重建扩张过程对已重建三维点的邻域进行匹配滤波过程采用光度一致性约束和可视性约束来去除噪点

3D patch定义：

位置、法向量（每个patch具有25个点， pathch的中心为点的位置，法向量为三维点邻域的法向量）

对三维空间上的一个pathch投影到不同视角中，计算NCC

1.在图像上均匀计算HOG/Harris特征

2.沿极线进行搜索找到匹配特征点

3.对匹配对进行三角化建立patch,

中心: 三角化确定

法向量: 指向参考图像

可视图像V(p)V(p)V(p): 法向量和视角的夹角足够小,通过光度一致性约束对可视图像进行筛选V∗(p)={I∣I∈V(P),h(p,I,R(P))≤a}V^{*}(p)=\{I|I\in V(P), h(p,I,R(P))\le a\}V∗(p)={I∣I∈V(P),h(p,I,R(P))≤a}

4.对patch位置和法向量进行优化maxc(p),n(p)g∗(p)=1∣V∗(p)\R(p)∣∑I∈V∗(p)\R(p)h(p,I,R(p))\mathop{max}\limits_{c(p), n(p)} g^{*}(p) = \frac{1}{|V^{*}(p) \backslash R(p)|}\sum_{I\in V^{*}(p) \backslash R(p)}h(p,I,R(p))c(p),n(p)max​g∗(p)=∣V∗(p)\R(p)∣1​I∈V∗(p)\R(p)∑​h(p,I,R(p))

深度图融合(重点)

基本步骤：

为图像选择邻域图像构成立体图像对(通过图像pair进行深度图重建)计算每一幅参考图像的深度图，得到每幅参考图像的点云对所有深度图进行融合

全局视角的选择(建立候选视角池)

初步筛选出一些适合用于重建的图像,局部视角从全局视角中选择.

1、图像具有相同的内容、外观和尺度

2、图像具有较大的视差（宽基线）

其中，FVF_{V}FV​表示V视角， FRF_{R}FR​表示R视角，wN(f)w_{N}(f)wN​(f)代表相机的三角测量夹角大小，ws(f)w_{s}(f)ws​(f)代表两幅图像的尺度相似性。

图像分辨率的估计

用于衡量图像的分辨率，图像的一个像素的宽度对应的三维空间中的物体尺寸。

相关视角的选择(针对每个像素)

全局视角的选择时固定的，局部视角的选择是针对每个像素的，每个像素在重建过程中不需要用到全部的全局视角，有些视角存在遮挡。局部视角选择是从全局视角中再一次进行筛选得到。一般从全局视角中选择四个当作局部视角进行优化。

局部视角的选择：

1、NCC值确定候选视角

2、实现要足够分散（不共面）

lR(V)=gR∗∏V′∈Awe(V,V′)l_{R}(V)=g_{R}*\prod_{V^{'}\in A} w_{e}(V, V^{'})lR​(V)=gR​∗V′∈A∏​we​(V,V′)

深度图融合算法整体思路

区域生长法扩张（宽度优先搜索）：

1、根据重建的置信度（NCC）建立优先级队列（排序）

2、从初始的稀疏特征点开始深度估计

3、对每个种子点进行非线性深度优化（核心）

4、每次优化完成判断两种情况，满足的话将邻域像素添加到队列中当作种子点，并将当前像素patch的信息作为邻域像素的初始值：（1）邻域没有深度值（2）当前像素的置信度高于邻域像素一定范围（当前点因为重建优化完，置信度得到提升，和邻域拉开很大差距，说明这个区域很需要优化，有很大的优化空间）

新建patch的可视图像和法向量的初始值等同于邻域patch

n(p′)=n(p)V(p′)=V(P)n(p^{'})=n(p)\\V(p^{'})=V(P)n(p′)=n(p)V(p′)=V(P)

深度值非线性优化

假设参考视角为R,k∈NR, k\in NR,k∈N为与参考视角相关的局部邻域视角,PkP_{k}Pk​为三维点投影到第k个视角的投影矩阵. 每个像素提供了一个初始的深度和法向量和局部视角后，对深度值进行非线性优化。建立几何模型如下：

其中，oRo_{R}oR​表示参考视角的相机中心，(s, t)表示参考视角像素点,其初始深度为h(s,t)h(s,t)h(s,t)，已知相机姿态, 则对应的三维空间中的射线的单位向量为r→R(s,t)\overrightarrow{r}_{R}(s, t)rR​(s,t), 则该像素对应的三维空间的点为:

XR(s,t)=OR+h(s,t)r→R(s,t)X_{R}(s, t)=O_{R}+h(s,t)\overrightarrow{r}_{R}(s, t)XR​(s,t)=OR​+h(s,t)rR​(s,t)

为了对像素点(s,t)(s,t)(s,t)的深度和三维中的法向量进行优化,以该像素为中心，建立nxn的patch, pathch对应三维中间中的一个很小的平面.同时引入两个变量来帮助表示patch中每个像素的三维坐标. 则像素(s+i,t+j)(s+i, t+j)(s+i,t+j)对应的深度为:

h(s+i,t+j)=h(s,t)+i∗hs(s,t)+j∗ht(s,t)h(s+i, t+j)=h(s,t)+i*h_{s}(s,t)+j*h_{t}(s,t)h(s+i,t+j)=h(s,t)+i∗hs​(s,t)+j∗ht​(s,t)

假设(s+i,t+j)(s+i, t+j)(s+i,t+j)处的射线方向近似为r→R(s,t)\overrightarrow{r}_{R}(s, t)rR​(s,t),则像素(s+i,t+j)(s+i, t+j)(s+i,t+j)对应的三维坐标为:

XR(i,j)=OR+r→R(s,t)(h(s,t)+i∗hs+j∗ht)X_{R}(i,j)=O_{R}+\overrightarrow{r}_{R}(s, t)(h(s,t)+i*h_{s}+j*h_{t})XR​(i,j)=OR​+rR​(s,t)(h(s,t)+i∗hs​+j∗ht​)

光度模型：

对朗伯反射进行简单假设,以提升颜色匹配的准确度. 为邻域内的每个视角K分配一个颜色尺度ck=[ckr,ckg,ckb]T∈R3×1c_{k}=[c_{k}^{r}, c_{k}^{g}, c_{k}^{b}]^{T} \in{R^{3\times 1}}ck​=[ckr​,ckg​,ckb​]T∈R3×1,如果深度h(s,t)h(s,t)h(s,t)估计准确，则有：

IR(s+i,t+j)=ck⋅Ik(Pk(XR(s+i,t+j)))I_{R}(s+i,t+j)=c_{k}\centerdot I_{k}(P_{k}(X_{R}(s+i, t+j)))IR​(s+i,t+j)=ck​⋅Ik​(Pk​(XR​(s+i,t+j)))

对所有的邻域视角和所有的patch中的patch中的三维点成立,其中IR(.)I_{R}(.)IR​(.)表示取参考视角对应像素处的颜色(3通道,为向量形式), Ik(.)I_{k}(.)Ik​(.)表示取第k个视角中图像上对应像素处的颜色,表示参考视角的一点和其他视角在对应同一点上的颜色一致.

结合几何模型和光度模型,可以得到深度值和法向量的优化数学模型如下：

E=∑k∈N∑i=−n−12n−12∑j=−n−12n−12[IR(s+i,t+j)−ck.Ik(Pk(XR(s+i,t+j)))]2E=\sum_{k\in N} \sum_{i=-\frac{n-1}{2}}^{\frac{n-1}{2}} \sum_{j=-\frac{n-1}{2}}^{\frac{n-1}{2}}[I_{R}(s+i,t+j)-c_{k}.I_{k}(P_{k}(X_{R}(s+i, t+j)))]^{2}E=k∈N∑​i=−2n−1​∑2n−1​​j=−2n−1​∑2n−1​​[IR​(s+i,t+j)−ck​.Ik​(Pk​(XR​(s+i,t+j)))]2

其中XR(s+i,t+j)X_{R}(s+i, t+j)XR​(s+i,t+j)表示将参考图像中的(s+i,t+k)(s+i, t+k)(s+i,t+k)映射到三维空间中的坐标,乘以PkP_{k}Pk​表示将该点的三维坐标投影到第k个视角图像上的像素坐标,对这个像素坐标取IkI_{k}Ik​表示将第k个视角图像上这个点的像素颜色值取出来,再乘以前面的ckc_{k}ck​,表示一个调节因子,增加鲁棒性

为了书写方便IR(s+i,t+j)I_{R}(s+i,t+j)IR​(s+i,t+j)简记为IR(i,j)I_{R}(i,j)IR​(i,j), 将Ik(Pk(X(s+i,t+j)))I_{k}(P_{k}(X(s+i,t+j)))Ik​(Pk​(X(s+i,t+j)))简记为Ik(i,j)I_{k}(i,j)Ik​(i,j),将hs(s,t),ht(s,t)h_{s}(s,t), h_{t}(s,t)hs​(s,t),ht​(s,t)分别简记为hs,hth_{s}, h_{t}hs​,ht​于是,上式可以简记为:

E=∑ijk[IR(i,j)−ck.Ik(i,j)]2E=\sum_{ijk}[I_{R}(i,j)-c_{k}.I_{k}(i,j)]^{2}E=ijk∑​[IR​(i,j)−ck​.Ik​(i,j)]2

i,j∈[−n−12,n−12]i,j\in[-\frac{n-1}{2}, \frac{n-1}{2}]i,j∈[−2n−1​,2n−1​]是patch中的采样点，k表示视角的个数，需要优化的变量为三维坐标的深度h(s,t),hs,hth(s,t), h_{s},h_{t}h(s,t),hs​,ht​和颜色尺度因子{ck}\{c_{k}\}{ck​}

颜色尺度的优化：

∂E∂ckr=∑ijrkr(i,j)IRr(i,j)+ckr∑ij(Ikr(i,j))2∂E∂ckg=∑ijrkg(i,j)IRg(i,j)+ckg∑ij(Ikg(i,j))2∂E∂ckb=∑ijrkb(i,j)IRr(i,j)+ckb∑ij(Ikb(i,j))2ck=[∑ijIKr(i,j)IRr(i,j)∑ij(IKr(i,j))2,∑ijIKg(i,j)IRg(i,j)∑ij(IKg(i,j))2,∑ijIKb(i,j)IRb(i,j)∑ij(IKb(i,j))2]T\frac{\partial E}{\partial c_{k}^{r}} = \sum_{ij}r_{k}^{r}(i,j)I_{R}^{r}(i,j) + c_{k}^{r}\sum_{ij}(I_{k}^{r}(i,j))^{2} \\\frac{\partial E}{\partial c_{k}^{g}} = \sum_{ij}r_{k}^{g}(i,j)I_{R}^{g}(i,j) + c_{k}^{g}\sum_{ij}(I_{k}^{g}(i,j))^{2} \\\frac{\partial E}{\partial c_{k}^{b}} = \sum_{ij}r_{k}^{b}(i,j)I_{R}^{r}(i,j) + c_{k}^{b}\sum_{ij}(I_{k}^{b}(i,j))^{2} \\c_{k}=[\frac{\sum_{ij}I_{K}^{r}(i,j)I_{R}^{r}(i,j)}{\sum_{ij}(I_{K}^{r}(i,j))^2}, \frac{\sum_{ij}I_{K}^{g}(i,j)I_{R}^{g}(i,j)}{\sum_{ij}(I_{K}^{g}(i,j))^2}, \frac{\sum_{ij}I_{K}^{b}(i,j)I_{R}^{b}(i,j)}{\sum_{ij}(I_{K}^{b}(i,j))^2}]^{T}∂ckr​∂E​=ij∑​rkr​(i,j)IRr​(i,j)+ckr​ij∑​(Ikr​(i,j))2∂ckg​∂E​=ij∑​rkg​(i,j)IRg​(i,j)+ckg​ij∑​(Ikg​(i,j))2∂ckb​∂E​=ij∑​rkb​(i,j)IRr​(i,j)+ckb​ij∑​(Ikb​(i,j))2ck​=[∑ij​(IKr​(i,j))2∑ij​IKr​(i,j)IRr​(i,j)​,∑ij​(IKg​(i,j))2∑ij​IKg​(i,j)IRg​(i,j)​,∑ij​(IKb​(i,j))2∑ij​IKb​(i,j)IRb​(i,j)​]T

颜色尺度的优化：

分别对上式求三个偏导，并令结果为0，得：

ckr=∑ijIKr(i,j)IRr(i,j)∑ij(IKr(i,j))2ckg=∑ijIKg(i,j)IRg(i,j)∑ij(IKg(i,j))2ckb=∑ijIKb(i,j)IRb(i,j)∑ij(IKb(i,j))2]Tc_{k}^{r} = \frac{\sum_{ij}I_{K}^{r}(i,j)I_{R}^{r}(i,j)}{\sum_{ij}(I_{K}^{r}(i,j))^2} \\c_{k}^{g}= \frac{\sum_{ij}I_{K}^{g}(i,j)I_{R}^{g}(i,j)}{\sum_{ij}(I_{K}^{g}(i,j))^2} \\c_{k}^{b}=\frac{\sum_{ij}I_{K}^{b}(i,j)I_{R}^{b}(i,j)}{\sum_{ij}(I_{K}^{b}(i,j))^2}]^{T}ckr​=∑ij​(IKr​(i,j))2∑ij​IKr​(i,j)IRr​(i,j)​ckg​=∑ij​(IKg​(i,j))2∑ij​IKg​(i,j)IRg​(i,j)​ckb​=∑ij​(IKb​(i,j))2∑ij​IKb​(i,j)IRb​(i,j)​]T

h(s,t),hs,hth(s,t), h_{s},h_{t}h(s,t),hs​,ht​的优化：

待优化的模型为最小二乘的非线性优化问题，无法得到该函数的解析解，采用迭代的方式求取其近似解。首先引入IR(i,j)I_{R}(i,j)IR​(i,j)关于h(s,t),hs,hth(s, t), h_{s}, h_{t}h(s,t),hs​,ht​的线性表达式:

IR(i,j)≈ck⋅Ik(Pk(OR+r→R(s,t)(h(s,t)+i∗hs+j∗ht)))+ck⋅∂Ik(i,j)∂h(s,t)(dh(s,t)+i∗dhs+j∗dht)=ck⋅Ik(i,j)+ck⋅∂Ik(i,j)∂h(s,t)(dh(s,t)+i∗dhs+j∗dht)\begin{aligned} \boldsymbol{I}_{R}(i, j) \approx \boldsymbol{c}_{k} \cdot & \boldsymbol{I}_{k}\left(\boldsymbol{P}_{k}\left(\boldsymbol{O}_{R}+\overrightarrow{\boldsymbol{r}}_{R}(s, t)\left(h(s, t)+i * h_{s}+j * h_{t}\right)\right)\right)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \\ =& \boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \end{aligned}IR​(i,j)≈ck​⋅=​Ik​(Pk​(OR​+rR​(s,t)(h(s,t)+i∗hs​+j∗ht​)))+ck​⋅∂h(s,t)∂Ik​(i,j)​(dh(s,t)+i∗dhs​+j∗dht​)ck​⋅Ik​(i,j)+ck​⋅∂h(s,t)∂Ik​(i,j)​(dh(s,t)+i∗dhs​+j∗dht​)​

由此,我们可以得到能量函数E的近似表达式为:

E=∑ijk(IR(i,j)−ck⋅Ik(i,j)+ck⋅∂Ik(i,j)∂h(s,t)(dh(s,t)+i∗dhs+j∗dht))2E=\sum_{i j k}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right)\right)^{2}E=ijk∑​(IR​(i,j)−ck​⋅Ik​(i,j)+ck​⋅∂h(s,t)∂Ik​(i,j)​(dh(s,t)+i∗dhs​+j∗dht​))2

令:

bijk=IR(i,j)−ck⋅Ik(i,j)∇x=[dh(s,t),dhs,dht]TAijk=(ck⋅∂Ik(i,j)∂h(s,t))[1ij]T\boldsymbol{b}_{i j k}=\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\\\;\\ \nabla \boldsymbol{x}=\left[d h(s, t), d h_{s}, d h_{t}\right]^{T}\\ \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T}bijk​=IR​(i,j)−ck​⋅Ik​(i,j)∇x=[dh(s,t),dhs​,dht​]TAijk​=(ck​⋅∂h(s,t)∂Ik​(i,j)​)⎣⎡​1ij​⎦⎤​T

则上式可以表达为:

E=∑ijk(Aijk∇x+bijk)2E=\sum_{i j k}\left(A_{i j k} \nabla x+b_{i j k}\right)^{2}E=ijk∑​(Aijk​∇x+bijk​)2

求EEE关于∇x\nabla x∇x的梯度，可以得到：

∂E∂∇x=∑ijkAijkT(Aijkx+bijk)\frac{\partial E}{\partial \nabla \boldsymbol{x}}=\sum_{i j k} \boldsymbol{A}_{i j k}^{T}\left(\boldsymbol{A}_{i j k} \boldsymbol{x}+\boldsymbol{b}_{i j k}\right)∂∇x∂E​=ijk∑​AijkT​(Aijk​x+bijk​)

令∂E∂∇x=0\frac{\partial E}{\partial \nabla x}=0∂∇x∂E​=0得：

(∑ijkAijkTAijk)∇x+∑ijkAijkTbijk)=0\left.\left(\sum_{i j k} A_{i j k}^{T} A_{i j k}\right) \nabla x+\sum_{i j k} A_{i j k}^{T} b_{i j k}\right)=0⎝⎛​ijk∑​AijkT​Aijk​⎠⎞​∇x+ijk∑​AijkT​bijk​⎠⎞​=0

令A=∑ijkAijkTAijk,b=∑ijkAijkbijkA=\sum_{ijk}A^{T}_{ijk}A_{ijk}, b=\sum_{ijk}A_{ijk}b_{ijk}A=∑ijk​AijkT​Aijk​,b=∑ijk​Aijk​bijk​，则有：

∇x=A−1b\nabla x=A^{-1}b∇x=A−1b

补充:

Aijk=(ck⋅∂Ik(i,j)∂h(s,t))[1ij]TAijkTAijk=[1ij](ck⋅∂Ik(i,j)∂h(s,t))T(ck⋅∂Ik(i,j)∂h(s,t))[1ij]T=∥ck⋅∂Ik(i,j)∂h(s,t)∥2[1,i,ji,i2,ijj,ij,j2]AijkTbijk=[1ij](ck⋅∂Ik(i,j)∂h(s,t))T(IR(i,j)−ck⋅Ik(i,j))=sijk[1ij]\begin{array}{c} \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial I_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T}\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T}=\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}\right]\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)=s_{i j k}\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right] \end{array}Aijk​=(ck​⋅∂h(s,t)∂Ik​(i,j)​)⎣⎡​1ij​⎦⎤​TAijkT​Aijk​=⎣⎡​1ij​⎦⎤​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(ck​⋅∂h(s,t)∂Ik​(i,j)​)⎣⎡​1ij​⎦⎤​T=∥∥∥​ck​⋅∂h(s,t)∂Ik​(i,j)​∥∥∥​2⎣⎡​1,i,j,​i,i2,ij,​jijj2​⎦⎤​AijkT​bijk​=⎣⎡​1ij​⎦⎤​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(IR​(i,j)−ck​⋅Ik​(i,j))=sijk​⎣⎡​1ij​⎦⎤​​

其中sijk=(ck⋅∂Ik(i,j)∂h(s,t))T(IR(i,j)−ck⋅Ik(i,j))s_{i j k}=\left(c_{k} \cdot \frac{\partial I_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(I_{R}(i, j)-c_{k} \cdot I_{k}(i, j)\right)sijk​=(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(IR​(i,j)−ck​⋅Ik​(i,j))是标量，因此可以得到：

A=∑ijkAijkTAijk=∑ijk∥ck⋅∂Ik(i,j)∂h(s,t)∥2[1,i,ji,i2,ijj,ij,j2]b=∑ijkAijkTbijk=∑ijk(ck⋅∂Ik(i,j)∂h(s,t))T(IR(i,j)−ck⋅Ik(i,j))[1ij]\begin{array}{l} \boldsymbol{A}=\sum\limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\sum\limits_{i j k}\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}\right]\\ \boldsymbol{b}=\sum \limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\sum\limits_{i j k}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right] \end{array}A=ijk∑​AijkT​Aijk​=ijk∑​∥∥∥​ck​⋅∂h(s,t)∂Ik​(i,j)​∥∥∥​2⎣⎡​1,i,j,​i,i2,ij,​jijj2​⎦⎤​b=ijk∑​AijkT​bijk​=ijk∑​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(IR​(i,j)−ck​⋅Ik​(i,j))⎣⎡​1ij​⎦⎤​​

算法流程：

只进行h(s,t)h(s,t)h(s,t)的优化(迭代4次)While(iteration <20)每间隔5次迭代，进行h(s,t)，hs,hth(s,t)， h_{s}, h_{t}h(s,t)，hs​,ht​的优化，否则仅仅优化深度优化完成后跟踪判断每个视角的置信度如果视角的置信度太小，或者超过迭代14次尚未收敛，则从局部视角中移除该视角，并重新进行视角选择