openMVS-- RefineMesh （原理及代码解读）

0、参考论文

Large-scale and high-quality Multi-view stereo

论文阅读笔记见作者另一篇笔记：链接

1、参数理解

• nResolutionLevel — 与subsample有关
• nMaxViews ---- 采用neighboring view（可视性较高），避免无效计算 -
• max-face-area — subdivide 阈值
  default:9 —6
• fRatioRigidityElasticity --在计算smooth gradient时，用于平衡 Rigidity && Elasticity间的比例因子
• nEnsureEdgeSize --edge length threshold
• nAlternatePair – 计算pair score的参照图像选取方式
• scale— 对应 multi-scale images
• scale-step 与initialize image 时的gauss blur sigma有关
• gradient-step–步长
  default :45.05—45.1

2、step1 – Initialize images

1. Scale—图像下采样尺度，用于实现coarse -to -detail
2. Step – 与图像高斯模糊化处理的sigma(方差)参数有关
3. 核心实现：

void MeshRefine::ThInitImage(uint32_t idxImage, Real scale, Real sigma)

4. 重要：图像尺度变化，更新对应相机参数

imageData.UpdateCamera(scene.platforms);

5. initialize 步骤涉及图像梯度计算：

图像grad计算，采用DerivativeKernel(3x5)及其转置计算双梯度，再融合

6. 代码实现机制：

• step1：AddEVENT

• step2：多线程

FOREACHPTR(pThread, threads)
		pThread->start(ThreadWorkerTmp, this);

• step3：MeshRefine refine 作为接口，调用event Run函数

refine.ThreadWorker();

void MeshRefine::ThreadWorker(){
...
evt->Run(this);
...
}

3、step2 – subdivide mesh

3.1、decimate

1. 代码实现：

scene.mesh.Clean()

函数内部：
借助VCG mesh 实现mesh decimate ,存在数据传递 Scene Mesh – VGG mesh –Scene Mesh

2. 实质：

通过设置decimate 的相关参数，最终实现基于mesh face 数目减少的mesh
退化/简化，整个过程包括对mesh的各种优化与clean，并涉及mesh 的拓扑、non-manifold、edge size的控制等属性

3.2、simplify

1. 代码实现：

// make sure there are no edges too small or too long
CLN::EnsureEdgeSize(p, epsilonMin, epsilonMax, collapseRatio, degenerate_angle_deg, mode, max_iters, 0);

实质：借助CGAL 库实现 CGAL::Simple_cartesian<double>

3.3、subdivide

1. 代码实现：

scene.mesh.Subdivide(maxAreas, maxArea);
------
PS：注意
const uint32_t maxAreaTh(2*maxArea);// 阈值变为2倍

2. subdivide 方式：

• one -to -four：新增edge的中点作为vertex
• 注意！！！
  对于area 没有超过阈值的face,若其邻接face 采用4分割，也会导致共享的edge上有新增的vertex，因此，也会partially split

3.3.1、reproject

1. step1 – 确定相机的可视mesh

针对每个camera 可视的vertex（采用vertex octree 八叉树 + frustum ）–将vertex邻接的 face作为该camera 的可视face ,进一步进行重投影

2. step2 – mesh to pix 代码实现：

EVTProjectMesh(idxImage, arrCameraFaces[idxImage])  

void MeshRefine::ThProjectMesh(uint32_t idxImage, const CameraFaces& cameraFaces)
{}
// project mesh to the given camera plane
void MeshRefine::ProjectMesh()

//core-Class
rasterer.Project(facet);

3. step3 – reprojection area
   
4. step4 – re-map vertex and camera faces

Subdivide()函数内的操作均会使mesh 的 vertex发生改变，因此需要更新vertex – face 的邻接关系

// re-map vertex and camera faces
ListVertexFacesPre();

4、step3 – multi-scale iteration

4.1、pre

1. 不同scale下，surface迭代不同次数

实质：基于surface 离散化，本质为 vertex 沿梯度下降，实现能量函数（代码中对应为 mesh score）的最小化

2. 迭代次数的计算

• iter
• gstep ----同一个scale下，每次迭代，步长变为gstep *= 0.98;
• iterStart && iterStop
  以此作为 ratioRigidityElasticity && bAdaptMesh(具体见下)的操作指示

3. 代码核心函数：

// score mesh using photo-consistency
// and compute vertices gradient using analytical method
double MeshRefine::ScoreMesh(double* gradients)

void MeshRefine::ThProcessPair(uint32_t idxImageA, uint32_t idxImageB)
{} //！！！重要

PS： 函数ThProcessPair中，计算face 法向量，其中法向量的朝向，由输入的vertex顺序确定！！！代码如下：

// given a triangle defined by three 3D points,
// compute its normal (plane's normal oriented according to the given points order)
// given a triangle defined by three 3D points,
// compute its normal (plane's normal oriented according to the given points order)
template<typename TYPE>
inline TPoint3<TYPE> ComputeTriangleNormal(const TPoint3<TYPE>& v0, const TPoint3<TYPE>& v1, const TPoint3<TYPE>& v2) {
	return (v1-v0).cross(v2-v0);
} // ComputeTriangleNormal

4.2、Mesh score --重要

原理同参考论文

4.2.1、score-photo及score-photo梯度

1. score-photo与梯度同时计算，并基于image pair累加
2. 核心实现：

EVTProcessPair
void MeshRefine::ThProcessPair(uint32_t idxImageA, uint32_t idxImageB)//重投影，score-photo计算、梯度计算

3. 重要参数：nAlternatePair

根据nAlternatePair确定计算score-photo的 reprojection “方向”
【i—j // j –I //both】

4. 重投影

• 原理

ImageAB = A—X (3D point) –B (非整像素，线性内插)

• 核心函数

// project image from view B to view A through the mesh;
// the projected image is stored in imageA
// (imageAB is assumed to be initialize to the right size)
void MeshRefine::ImageMeshWarp()

• 应用：基于重投影结果，计算score-photo（ZNCC） && 梯度

PS：

• nReduceMemory减少内存:重新计算图像像素值local variance
  代码：
  ComputeLocalVariance(viewA.image, mask, _imageMeanA, _imageVarA);
• ZNCC
  window size: 7*7
  原理：见参考链接
  核心函数：
  ComputeLocalZNCC(imageA, *imageMeanA, *imageVarA, imageAB, imageMeanAB, imageVarAB, mask, imageZNCC, imageDZNCC)
  涉及DZNCC (ZNCC的梯度)–数学公式未看

5. 注意：计算photo-score 及梯度时，涉及vertex 的depth 的更新
   

4.2.1.1、score-photo

1.所有image pair 的score-photo以RegularizationScale为权值进行累加

scorePhoto += (float)RegularizationScale*score;

const Real RegularizationScale((Real)((REAL)(imageDataA.avgDepth*imageDataB.avgDepth)/(cameraA.GetFocalLength()*cameraB.GetFocalLength())));
//同参考论文：---引入scale factor，使得E_error & E_fair homogenized  
//Scale factor： 【avg(depth) /focal length】^2

2. 每个 image pair的score-photo计算：公共可视三维场景，重投影得到的imageI与imageI'的ZNCC && ZNCC window 的置信度 【与pix的local variance 有关】（见参考论文）

score += (float)(ReliabilityFactor*(Real(1)-ZNCC));

const Real ReliabilityFactor(minVAVB/(minVAVB+Real(0.0015)));

4.2.1.2、score-photo梯度

1. 核心函数及其细节：

ComputePhotometricGradient()

遍历Image pair 公共可视pix
Pix（A）所在的face 的vertex 梯度及梯度norm的计算
具体与face 法向、iamgeA 的视线方向、ZNCC的梯度以及比例因子有关

photoGrad[idxVert] += g;
++photoGradNorm[idxVert];

--------------------

const Grad g(N*(sg*(Real)b[v]));
const Real sg((gB*(xJac*(const TMatrix<Real,3,1>&)dA))(0)*dZNCC*RegularizationScale/Nd);
const Point3f& b(baryMapA(r,c));

4.2.2、score-smooth及score-smooth梯度

4.2.2.1、score-smooth – smooth1

1. 核心函数：

EVTSmoothVertices1
ComputeSmoothnessGradient1()

2. Smooth-score1：
   遍历多个vertex簇（VIndex idxStart, VIndex idxEnd）中vertex，累加score（vertex的smooth-score1）

scoreSmooth += score;//累加所有vertex簇的smooth-score 1
score += regularityScore;//某个vertex簇的smooth-score1
const float regularityScore((float)norm(grad));//单个vertex的smooth-score1
grad = grad/(Real)verts.GetSize() - Cast<Real>(vertices[idxV]);// grad/(Real)verts.GetSize()：邻接vertex的质心坐标
Grad& grad = smoothGrad1[idxV];

4.2.2.2、score-smooth梯度

1. smooth2 核心代码：

EVTSmoothVertices2

//实质：基于smooth1计算二阶smooth2，方法同smooth1 – 结合邻接vertex：
void MeshRefine::ComputeSmoothnessGradient2(）

2. Score-smooth -梯度

基于 smoothGrad2 || smoothGrad1 && smoothGrad2 （ratioRigidityElasticity）

4.2.3、final mesh score

(nAlternatePair ? 0.2f : 0.1f)*scorePhoto + 0.01f*scoreSmooth

4.3 、迭代更新-基于梯度下降

基于梯度（方向）+ 步长的surface 迭代更新（本质： vertex坐标更新 ）

vert -= Cast<Vertex::Type>(grad*gstep);// 沿梯度下降！！！

4.4、adaptMesh—vertex remove

1. bAdaptMesh 决定是否进行该步骤
2. Remove 条件：

• 非边界vertex边界vertex 的判定条件， 见以下函数:
• 梯度大小 && 阈值
• smoothGrad1 && 阈值

!refine.vertexBoundary[v] && (float)gn < th && norm(refine.smoothGrad1[v]) < th
remove planar vertices (small gradient and almost on the center of their surrounding patch) //!!!!