openMVS ---dense point cloud（代码理解）

0 代码实现原理

1. 参考论文及阅读笔记：

• Bleyer M , Rhemann C , Rother C . PatchMatch Stereo - Stereo Matching with Slanted Support Windows[C]// British Machine Vision Conference 2011. 2011. 阅读笔记
• Goesele, Michael, Snaveley, Noah, Curless, Brian,等. Multi-View Stereo for Community Photo Collections[J]. Proc.int.conf.on Computer Vision, 2007:1-8.阅读笔记
• Shen S . Accurate Multiple View 3D Reconstruction Using Patch-Based Stereo for Large-Scale Scenes[J]. IEEE Transactions on Image Processing, 2013, 22(5):1901-1914.阅读笔记
• Merrell P , Akbarzadeh A , Wang L , et al. Real-Time Visibility-Based Fusion of Depth Maps[C]// IEEE 11th International Conference on Computer Vision, ICCV 2007, Rio de Janeiro, Brazil, October 14-20, 2007. IEEE, 2007.阅读笔记

2. 代码理解参考博客：链接

1 代码实现步骤

1.1 参数传递

1. option 设置：

Key -value
value 与外部变量绑定并设置初值

2. 命令行/config输入参数，解析形成 key-value 形式的option,并更新绑定的外部变量 OPT
3. OPT—更新至 OPTDENSE

1.2 loading .mvs数据

1. .mvs to Interface

if (!ARCHIVE::SerializeLoad(obj, fileName))

2. Interface to Scene

bool Scene::LoadInterface(const String & fileName)

1.3 稠密重建

接口函数：

scene.DenseReconstruction()

1.3.0 相关变量及概念

1. mvs数据传递
   scene-denseDepthMapData

DenseDepthMapData data(*this, nFusionMode);

2. 重要数据类型：

DenseDepthMapData

DepthMapsData // !!!structure used to compute all depth-maps

DepthMapsData成员函数

DepthDataArr arrDepthData;

DepthData  --- 对应一张深度图

ViewScoreArr neighbors; // array of all images seeing this depth-map (ordered by decreasing importance) 【IndexType idx;ScoreType score;】

BitMatrix mask; // mark pixels to be ignored

1.3.1 ComputeDepthMaps(data) --计算多张深度图

1.3.1.1. pre

1. image index\ image map ？？？
2. image solution

• 根据level（0表示原图） 进行降采样，减小分辨率，最终影像的大小（max(width,height)）需在参数 min_solution & max_solution 之间
• 重要！！！：Image solution 改变—camera参数K需改变

Normalized K 与图像分辨率无关，full K 与图像分辨率有关，K normalizing 的过程，除以scale=
max(width, height)

1.3.1.2 为每一张depth map Select view（s） ！！！

1. 接口函数：

data.depthMaps.SelectViews(depthData)

2. 为reference view选择 nNumViews 个对应视图，若nNumViews == 1，选择target view( 1:1 )，反之，选择N neighbors views（1：n ）

if (OPTDENSE::nNumViews == 1 &&
!data.depthMaps.SelectViews(data.images, imagesMap,
data.neighborsMap))

nNumViews == 1时，才会执行 &&后的函数，得到only one target view， 否则，一个reference
view 对应nNumViews个target view

1.3.1.2.1 N neighbors — global score + filter

1. 接口函数

scene.SelectNeighborViews

2. 相关指标：

OPTDENSE::nMinViews, 
OPTDENSE::nMinViewsTrustPoint
FD2R(OPTDENSE::fOptimAngle)

3. global score计算原理

Area_ratio * ( ∑_3d points( Wa*Ws))
4. global score代码细节

• P & X 计算depth！！！

REAL Camera::PointDepth(const Point3& X) const
{
	return P(2,0)*X.x + P(2,1)*X.y + P(2,2)*X.z + P(2,3);//参考P108
} // PointDepth

• Scale 计算：
  方法1 ：与原论文（Multi-View Stereo for Community Photo Collections）一致
  方法2 ：

camera.GetFocalLength()/camera.PointDepth(X)

5. 根据global score filter neighbors

score.points < 3 公共可视点个数

neighbors.GetSize() < MINF(nMinViews,nCalibratedImages-1)

Scene::FilterNeighborViews(depthData.neighbors, fMinArea, fMinScale, fMaxScale, fMinAngle, fMaxAngle, OPTDENSE::nMaxViews)
// remove invalid neighbor views
	const float fMinArea(OPTDENSE::fMinArea);
	const float fMinScale(0.2f), fMaxScale(3.2f);
	const float fMinAngle(FD2R(OPTDENSE::fMinAngle));
	const float fMaxAngle(FD2R(OPTDENSE::fMaxAngle));

6. 取score较高的前OPTDENSE::nMaxViews, 得到N

if (neighbors.GetSize() > nMaxViews)
		neighbors.Resize(nMaxViews);//取score较高的nMaxViews --- neighbors sorted

1.3.1.2.2 MRF ---- the only one target neighbor

1. 选择要求：

// trying in the same time the selected image pairs to cover the whole
scene;

2. 接口函数：

DepthMapsData::SelectViews(IIndexArr& images, IIndexArr& imagesMap, IIndexArr& neighborsMap)

3. 采用MRF优化，得到最终结果neighborsMap（MRF优化原理及代码实现未细看）！！！

向graph中添加nodes、edges ，并使energy最小

1.3.1.3 单张深度图计算–reference && target / N

1.3.1.3.1 pre

1. 选择data.idxImage代表的图像作为前处理的reference image
2. 单幅图像处理过程根据 event queue中的event 来进行

1.3.1.3.2 initial depth --EVT_PROCESSIMAGE

1.函数实现：

DepthMapsData::InitViews(DepthData& depthData, IIndex idxNeighbor, IIndex numNeighbors)

2. TODO：创建与reference image相对应的 depthData，主要包括depthData.image 的初始化（reference image+ other view(s)）

• 彩色转灰度
• 尺度：
  other view 的尺度统一到reference image ----图像尺度变化，相机内参变化

1.3.1.3.3 遍历像素值计算depth-- EVT_ESTIMATEDEPTHMAP

根据fusemode选择算法：

//本文主要分析算法1

//算法1
// extract depth-map using Patch-Match algorithm
data.depthMaps.EstimateDepthMap(data.images[evtImage.idxImage]);

//算法2 ---未看
// extract disparity-maps using SGM algorithm
data.sgm.Match()
data.sgm.Fuse()

论文

1. STEP1：initial depth

• 方法1：基于sparse points 插值
  代码基于CGAL实现

std::pair<float,float> TriangulatePointsDelaunay(const DepthData::ViewData& image, const PointCloud& pointcloud, const IndexArr& points, CGAL::Delaunay& delaunay)

initil 原理：得到sparse point深度值范围后，在x-y-depth空间，构建Delaunay三角网，并确定每个三角形对应的x-y范围，使得在x-y范围的像素根据该三角形在x-y-z空间所确定的平面，计算对应深度z。
本质：内插

重要：不直接在x-y-depth空间插值的原因：
像素间depth存在不连续性，而XYZ空间 (surface)depth连续性较强

关于corner point depth的计算

• 最开始设置为整个spare points的depth 均值，插入到x-y-depth 空间的vertex，并构成DT，在x-y-depth 空间的DT中找到与其相邻的face，并确定相邻face的vertex
• 在X-Y-Z真实三维空间中，根据vertex对应的x,y +depth 通过投影，得到真实XYZ（3组），并得到真实surface（plane）
• 将corner point 的可视光线ray 与该plane 相交，交点的depth作为corner point的待选depth,并将corner point 与空间三点_ XYZ（3组） 的重心的对应2D 的像素坐标距离的倒数作为权重
  最后，将3个待选depth 加权得到corner point的最终depth
• 其他像素depth initial 插值方法同上，DT xy范围内的像素，将可视光线与对应的surface-plane相交插值即可

方法2：较简单，略

2. ！！！STEP2：逐像素遍历
   1）三次空间遍历（左上-右下 // 右下-左上）
   2）参数空间传播、迭代更新（验证-- 采用NCC score）

• 优化迭代depth的核心函数：

void* STCALL DepthMapsData::EstimateDepthMapTmp(void* arg)

// run propagation and random refinement cycles;
// the solution belonging to the target image can be also propagated
estimator.ProcessPixel(idx);

• 优化迭代过程
  STEP1// check if any of the neighbor estimates are better then the current estimate

邻域像素Confidence 阈值约束
传播得到的conf 与原conf比较

STEP2:// try random values around the current estimate in order to refine it

第一次随机：
// try completely random values in order to find an initial estimate
第二次随机 ：
const Depth ndepth(rnd.randomMeanRange(depth,
depthRange*scaleRange));

STEP3:
迭代更新条件：

if (conf > nconf) {
			conf = nconf;
			depth = neighbor.depth;
			normal = neighbor.normal;
		}

• 迭代优化后处理：

// remove all estimates with too big score and invert confidence map
//核心函数：
void* STCALL DepthMapsData::EndDepthMapTmp(void* arg)
//remove条件：
depth <= 0 || conf >= OPTDENSE::fNCCThresholdKee

3） NCC计算

• 窗口对应—基于support plane homography matrix 引入local smoothness 进行改进
  
• 每个像素 对应 11* 11 的NCC 邻域窗口，取stepSize = 2, 实际参与计算的像素个数为6*6，对应的权重矩阵为 6 *6
  （权重定义：）

1.3.1.3.4 单张深度图优化 – EVT_OPTIMIZEDEPTHMAP

PS：论文参考

Accurate, Dense, and Robust Multi-View Stereopsis 拓展

1.3.1.4.1 remove small segment

1. 核心函数：

data.depthMaps.RemoveSmallSegments(depthData)

2. 遍历depthmap进行segment，在形成segments 的过程中判断是否需要remove
3. segmentation 过程：

• 从当前像素开始，遍历4邻域像素，计算二者depth的相似性

depth_neighbor>0 && IsDepthSimilar(depth_curr, depth_neighbor, fDepthDiffThreshold)

IsDepthSimilar()—核心实现：

inline T DepthSimilarity(T d0, T d1)

• 根据相似性，判断是否加入当前seg_list，并继续基于4邻域扩充

3. segment remove 条件：

// if segment NOT large enough => invalidate pixels
			if (seg_list_count < speckle_size) {

1.3.1.4.2 fill the gap

1. 核心函数：

data.depthMaps.GapInterpolation(depthData)

2. 实现过程：

分别按行遍历depthMap值，按列遍历depthMap值
当gap的大小，gap前后（行方向或者列方向）上像素的depth满足相似条件时，采用线性内插的方式得到gap的深度值，并内插normal(
基于direction),更新conf（较小值）

1.3.1.5 单张深度图融合 filter & adjust（参考neighbors）

1. 核心函数：
   DenseReconstructionFilter（）
2. 同上，采用event进行filter& adjust操作的管理
   

1.3.1.5.1 EVT_FILTERDEPTHMAP

1. 核心操作：

采用1 : numMaxNeighbors(8) ，实现基于target depthMaps的，对refer depthMap进行filter
PS :
target depthMaps降序排列，选择最重要的前8个neighbors【images seeing this depth-map (ordered by decreasing importance)】

2. 具体实现：

data.depthMaps.FilterDepthMap(depthData, idxNeighbors, OPTDENSE::bFilterAdjust)

• 将neighbors depthMap 投影到reference view，得到多个projection reference depthMap

// project all neighbor depth-maps to this image

- Neighbor view+ neighbor depthMap --- TransformPointI2W
- 3D point + P （reference image） -- TransformPointW2C 得到depth  TransformPointC2I-得到对应二维像平面亚像素位置：
- set depth on the 4 pixels around the image projection 投影到二维图像对应亚像素位置

• 根据bAdjust，设置不同的filter方式
  1） fuse更新reference depthMap

更新方式：Confidence-base fuse ： 和原论文基本一致？？？

更新的具体过程：

• depth_refer & depth_neighbor_project_2_refer 判断similar（是否相差阈值范围内）

±---------+
相似：
Depth_new: 原depth以conf为权重加权融合
Conf_new: 直接相加
±------------+
不相似：
// penalize confidence
根据disagree类型不同，设置negConf

• 判断最终的Depth_new、Conf_new是否是inlier(代码如下)
  若为oulier—diacard: Depth_new = 0、Conf_new = 0

if (nPosViews >= nMinViewsAdjust && posConf > negConf && ISINSIDE(avgDepth/=posConf, depthDataRef.dMin, depthDataRef.dMax))

oulier—diacard: Depth_new = 0、Conf_new = 0

2）直接remove
// remove depth if it does not agree with enough neighbors

衡量depth_refer & depth_neighbor_project_2_refer agree 是否enough：

• 首先
  // check if very similar with the neighbors projected to this pixel
• 进一步地
  // check if similar with the neighbors projected around this pixel ！！！
  【depth_refer （i,j）& depth_neighbor_project_2_refer （i,j）对应4邻域位置的depth】

IsDepthSimilar(depth, d, thDepthDiffStrict)
if (nGoodViews < nMinViews || nGoodViews < nViews*nMinGoodViewsProc/100)

1.3.1.5.2 EVT_ADJUSTDEPTHMAP

删除中间文件 filtered_damp 将最终的depth保存为 damp（将未filtered的damp直接覆盖掉）

1.3.2 data.depthMaps.FuseDepthMaps–多张深度图融合投影

1. 核心函数

// fuse all valid depth-maps in the same 3D point cloud;
// join points very likely to represent the same 3D point and
// filter out points blocking the view
void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, bool bEstimateNormal)

2. 具体过程：

• STEP1：每张image计算score，由高到低排序，依次处理所有images

connection.score = (float)scene.images[i].neighbors.GetSize();

• STEP2：当前image ，根据对应depthMap 得到3D point， 并设置权重（与conf & 深度值有关）

PS：// convert the ZNCC score to a weight used to average the fused points

inline float Conf2Weight(float conf, Depth depth) {
	return 1.f/(MAXF(1.f-conf,0.03f)*depth*depth);
}

• STEP3：// check the projection in the neighbor depth-maps

• 比较projection 得到的depth（测量值）pt.z 与 depthMap （计算值）depthB = depthMapB(xB)的差异
• check if normals agree（代码见下）

if (normal.dot(normalB) > normalError)

• STEP4：若差异较小：融合 X_refer & X_neighbor Normal color 通过权重融合！！！

• STEP5：最后，根据每个3D point 的可见视图的阈值，判断融合后的3D point 是否有效

• Remove： 与3D point相关view的depthMap\对应的2D 投影点
  PS：3D scene 中与点有关的量
  
• Store：
  // invalidate all neighbor depths that do not agree with it

• STEP6：Normal估计
  选择该3D点weight 最大的depthMap 对应的NormalMap 上的对应值

3. 后续步骤

• 若最终无normalMap，最终估计normalMap

// estimates the normals through PCA over the K nearest neighbors

采用CGAL实现基于PCA的法向估计
--------------------------------------+
PS:

KD tree 的相关概念 “

• K个维度，根据某维度中位值进行二分，并将中位值数据作为结点，直至无法二分（该区域只有一个数据点）

• KNN,基于K-Dtree 实现K邻近搜索：搜索路径 + 路径回溯（基于搜索点+d_min 确定的“球体”与父节点超平面的相交情况），更新搜索路径+ 距离_节点 -堆栈，回溯结束，得到K邻近

• K邻域法向估计： 该方向上，领域点+本身点的投影最集中，由此确定的平面，最能贴近该点群的分布
  原理：点群空间分布的连续性），即该点群的PCA方向（特征值最大）的垂直方向（特征值最小）！！！

// correct normal orientation
		ASSERT(!views.IsEmpty());
		const Image& imageData = images[views.First()];
		if (normal.dot(Cast<float>(imageData.camera.C)-point) < 0)
			normal = -normal;

• 若无colorMap，最终估计colorMap ： 采用3D 点与对应的视图中相机与3D点距离最近的视图的对应像素颜色值（投影到该视图可能对应亚像素—双线性插值）

1.3.3 其他

2 代码部分的疑惑

imageMap
MRF 优化选择target view
Initial depth –三角化插值
ZigzgMatric cords
weightMap –NCC score