【三维重建】OpenMVS代码解析——点云稠密化：详细过程
（2）
上⼀节内容过了⼀遍DensifyPointCloud.cpp中的main()函数，其中的每⼀步的具体操作在本篇⽂章中进⾏介绍。本节介绍的函数都是独⽴的辅助功能，不是点云稠密化的核⼼功能。
1. 功能⼀：在曲⾯上进⾏点云采样的功能
1 //option1: enable uniformly samples points on a mesh
2 Scene scene(OPT::nMaxThreads);
3 if (OPT::fSampleMesh != 0) {
4 // sample input mesh and export the obtained point-cloud
伏羲与女娲5 //load .ply of .obj pointcloud
6 if (!sh.Load(MAKE_PATH_SAFE(OPT::strInputFileName)))
7 return EXIT_FAILURE;
8 TD_TIMER_START();
9
10 PointCloud pointcloud;
11 if (OPT::fSampleMesh > 0) // t points density
12 sh.SamplePoints(OPT::fSampleMesh, 0, pointcloud);
13 el // t points number
14 sh.SamplePoints((unsigned)ROUND2INT(-OPT::fSampleMesh), pointcloud);
15
16 VERBOSE("Sample mesh completed: %u points (%s)", pointcloud.GetSize(), TD_TIMER_GET_FMT().c_str());
17 pointcloud.Save(MAKE_PATH_SAFE(Util::getFileFullName(OPT::strOutputFileName))+_T(".ply"));
18 Finalize();
19 return EXIT_SUCCESS;
20 }
OPT::fSampleMesh是进⾏点云采样的操作指令，OPT::fSampleMesh != 0意味着要进⾏这部分功能，return EXIT_SUCCESS意味着进⾏完这部分内容后不再执⾏后⾯的代码了。这⼀⼤段代码⾥，核⼼内容就两句
话，sh.Load(MAKE_PATH_SAFE(OPT::strInputFileName))和sh.SamplePoints((unsigned)ROUND2INT(-OPT::fSampleMesh), pointcloud)。
1.1 导⼊数据
Load()函数在libs/MVS/Mesh.cpp下，这部分代码很直观，从这⾥可以看出来，OpenMVS的作者他是⾃⼰去实现输⼊输出，点
云，octotree，甚⾄数据类型都是⾃定义的，⽽不是调库。
站在代码阅读的⾓度来看，读到这⾥已经可以了，如果感兴趣的话可以进⼀步读到LoadOBJ()和LoadPLY()⾥⾯去，可以帮助理解.obj 和.ply⽂件是怎么描述mesh的，这⾥就不展开介绍了。这两个函数的作⽤是读取对应的mesh模型⽂件，并把对应的数据保存到内存中去。
1// import the mesh from the given file
2bool Mesh::Load(const String& fileName)
3{
4 TD_TIMER_STARTD();
5 const String ext(Util::getFileExt(fileName).ToLower());
6 bool ret;
7 if (ext == _T(".obj"))
8 ret = LoadOBJ(fileName);
9 el
10 ret = LoadPLY(fileName);
11 if (!ret)
12 return fal;
13 DEBUG_EXTRA("Mesh loaded: %u vertices, %u faces (%s)", vertices.GetSize(), faces.GetSize(), TD_TIMER_GET_FMT().c_str());
14 return true;
15}
1.2 在mesh上进⾏点云采样
SamplePoints()函数也是在libs/MVS/Mesh.cpp下，
1void Mesh::SamplePoints(REAL samplingDensity, unsigned mumPointsTheoretic, PointCloud& pointcloud) const 2{
3 //allocate memory
4 ASSERT(!IsEmpty());
至强e3
5 pointcloud.Relea();
6 if (mumPointsTheoretic > 0) {
7 ve(mumPointsTheoretic);
8 if (HasTexture())
9 ve(mumPointsTheoretic);
10 }
11
12 // for each triangle
13 std::mt19937 rnd((std::random_device())()); //generate random number
14 std::uniform_real_distribution<REAL> dist(0,1);
15 FOREACH(idxFace, faces) {
16 const Face& face = faces[idxFace];
17
18 // vertices (OAB)
19 const Vertex& O = vertices[face[0]];
20 const Vertex& A = vertices[face[1]];
21 const Vertex& B = vertices[face[2]];
22
23 // edges (OA and OB)
24 const Vertex u(A - O);
25 const Vertex v(B - O);
26
27 // compute triangle area
28 const REAL area(ss(v)) * REAL(0.5));
29
30 // deduce the number of points to generate on this face
31 const REAL fPointsToAdd(area*samplingDensity);
32 unsigned pointsToAdd(static_cast<unsigned>(fPointsToAdd));
33
34 // take care of the remaining fractional part;
35 // add a point with the same probability as its (relative) area
36 const REAL fracPart(fPointsToAdd - static_cast<REAL>(pointsToAdd));
37 if (dist(rnd) <= fracPart)
38 pointsToAdd++;
39
40 for (unsigned i = 0; i < pointsToAdd; ++i) {
41 // generate random points as in:
42 // "Generating random points in triangles", Greg Turk;
43 // in A. S. Glassner, editor, Graphics Gems, pages 24-28. Academic Press, 1990
44 REAL x(dist(rnd));
45 REAL y(dist(rnd));
46
47 // test if the generated point lies on the right side of (AB)
48 if (x + y > REAL(1)) {
49 x = REAL(1) - x;
50 y = REAL(1) - y;
51 }
52
53 // compute position
54 place_back(O + static_cast<Vertex::Type>(x)*u + static_cast<Vertex::Type>(y)*v);
55
56 if (HasTexture()) {
57 // compute color
58 const FIndex idxTexCoord(idxFace*3);
59 const TexCoord& TO = faceTexcoords[idxTexCoord+0];
60 const TexCoord& TA = faceTexcoords[idxTexCoord+1];
61 const TexCoord& TB = faceTexcoords[idxTexCoord+2];
61 const TexCoord& TB = faceTexcoords[idxTexCoord+2];
62 const TexCoord xt(TO + static_cast<TexCoord::Type>(x)*(TA - TO) + static_cast<TexCoord::Type>(y)*(TB - TO));
63 place_back(textureDiffu.sampleSafe(Point2f(xt.x*textureDiffu.width(), (1.f-xt.y)*textureDiffu.height())));
64 }
65 }
66 }
67}
⾸先给点云开辟内存空间，再遍历每⼀个mesh，根据采样分辨率和计算在当前的mesh上应该采样多少个点。这⾥计算的点数可能是⼩数，⽐如说100.6个点，此时随机数函数就排上了⽤场，在0-1之间创建⼀个随机数，如果这个随机数⼩于那个0.6，那么就创建101个点，否则就创建100个点。然后假设当前mesh有三个顶点OAB，那么采样获得的点就在这三个点围成的三⾓形上，颜⾊取该位置到三个点的距离作为权值，融合颜⾊给到当前的采样点上。
磁铁相吸相斥原理
2. 导⼊数据
这部分内容是导⼊输⼊数据，
1 // load and estimate a den point-cloud
2 // input images and pos
3 if (!scene.Load(MAKE_PATH_SAFE(OPT::strInputFileName)))
4 return EXIT_FAILURE;
5 // input mesh
6 if (!OPT::pty())
7 sh.Load(MAKE_PATH_SAFE(OPT::strMeshFileName));
8 if (scene.pointcloud.IsEmpty() && OPT::pty() && sh.IsEmpty()) {
爱是一道光9 VERBOSE("error: empty initial point-cloud");
10 return EXIT_FAILURE;
11 }
scene.Load(MAKE_PATH_SAFE(OPT::strInputFileName))导⼊的是相机⾥程计和图
⽚，sh.Load(MAKE_PATH_SAFE(OPT::strMeshFileName))导⼊的是经过SfM获得的mesh模型，对于后⾯的内容，以上数据缺⼀不可。
2.1 导⼊相机数据
Load()函数在libs/MVS/Scene.cpp下，
1bool Scene::Load(const String& fileName, bool bImport)
2{
3 TD_TIMER_STARTD();
4 Relea();
5
6 #ifdef _USE_BOOST
7 // open the input stream
唐山海怎么死的8 std::ifstream fs(fileName, std::ios::in | std::ios::binary);
9 if (!fs.is_open())
10 return fal;
11 // load project header ID
12 char szHeader[4];
ad(szHeader, 4);
14 if (!fs || _tcsncmp(szHeader, PROJECT_ID, 4) != 0) {
15 fs.clo();
16 if (bImport && Import(fileName))
17 return true;
18 if (LoadInterface(fileName))
19 return true;
20 VERBOSE("error: invalid project");
21 return fal;
22 }
23 // load project version
24 uint32_t nVer;
ad((char*)&nVer, sizeof(uint32_t));
ad((char*)&nVer, sizeof(uint32_t));
26 if (!fs || nVer != PROJECT_VER) {
27 VERBOSE("error: different project version");
28 return fal;
29 }
30
31 // load stream type
32 uint32_t nType;
ad((char*)&nType, sizeof(uint32_t));
34
35 // skip rerved bytes
36 uint64_t nRerved;
ad((char*)&nRerved, sizeof(uint64_t));
38
39 // rialize in the current state
40 if (!SerializeLoad(*this, fs, (ARCHIVE_TYPE)nType))
41 return fal;
42
43 // init images
44 nCalibratedImages = 0;
45 size_t nTotalPixels(0);
46 FOREACH(ID, images) {
47 Image& imageData = images[ID];
48 if (imageData.poID == NO_ID)
49 continue;
512汶川大地震
50 imageData.UpdateCamera(platforms);
51 ++nCalibratedImages;
52 nTotalPixels += imageData.width * imageData.height;
53 }
54 DEBUG_EXTRA("Scene loaded (%s):\n"
55 "\t%u images (%u calibrated) with a total of %.2f MPixels (%.2f MPixels/image)\n"
56 "\t%u points, %u vertices, %u faces",
57 TD_TIMER_GET_FMT().c_str(),
58 images.GetSize(), nCalibratedImages, (double)nTotalPixels/(1024.0*1024.0), (double)nTotalPixels/(1024.0*1024.0*nCalibratedImages),
59 pointcloud.points.GetSize(), mesh.vertices.GetSize(), mesh.faces.GetSize());
60 return true;
61 #el
62 return fal;
63 #endif
64} // Load
3.功能⼆：降采样
1 //option2: split the scene in sub-scenes such that each sub-scene surface does not exceed the given maximum sampling area (0 - disabled)
2 if (OPT::fMaxSubsceneArea > 0) {
3 // split the scene in sub-scenes by maximum sampling area
4 Scene::ImagesChunkArr chunks;
5 // downsampling depth pixels and images and get pintcloud
6 scene.Split(chunks, OPT::fMaxSubsceneArea);
7 // save each chunk as a .mvs
8 scene.ExportChunks(chunks, Util::getFilePath(MAKE_PATH_SAFE(OPT::strOutputFileName)), (ARCHIVE_TYPE)OPT::nArchiveType);
9 Finalize();
10 return EXIT_SUCCESS;
相机成像原理
11 }
OPT::fMaxSubsceneArea是执⾏降采样的操作指令，>0意味着执⾏这部分代码。scene.Split()是降采样的操作过
程，scene.ExportChunks()是输出结果。
3.1 降采样
Split()函数在libs/MVS/Scene.cpp下，这⼀长串原⽣注释介绍了这个函数是⼲什么⽤的：限制读取到内存⾥的照⽚数量和照⽚分辨率，⽬的就是让内存⼀次就能把当前的数据都加载了。
1// split the scene in sub-scenes such that each sub-scene surface does not exceed the given
2// maximum sampling area; the area is compod of overlapping samples from different cameras
3// taking into account the footprint of each sample (pixels/unit-length, GSD inver),
4// including overlapping samples;
5// the indirect goals this method tries to achieve are:
6// - limit the maximum number of images in each sub-scene such that the depth-map fusion
7// can load all sub-scene's depth-maps into memory at once
8// - limit in the same time maximum accumulated images resolution (total number of pixels)
9// per sub-scene in order to allow all images to be loaded and procesd during mesh refinement
10unsigned Scene::Split(ImagesChunkArr& chunks, float maxArea, int depthMapStep) const {}
3.1.1 定义⼀些变量
1 const float areaScale(0.01f);//resolution
2 typedef cList<Point3f::EVec,const Point3f::EVec&,0,4096,uint32_t> Samples;
3 typedef TOctree<Samples,float,3> Octree;
4 Octree octree;//octotree in world
5 FloatArr areas(0, images.size()*4192);//area of each surface piece
6 IIndexArr visibility(0, (IIndex)areas.capacity());//for each points in Samples, the idx of image it belongs to
7 Unsigned32Arr imageAreas(images.size());//for each image, how many points can be obrved on it
8 Samples samples(0, (uint32_t)areas.capacity());//pointcloud in world
3.1.2 对每⼀张image进⾏降采样
1 FOREACH(idxImage, images) {
2 const Image& imageData = images[idxImage];
3 if (!imageData.IsValid())
4 continue;
困难的反义词是什么
5 DepthData depthData;
6 depthData.Load(CompoDepthFilePath(imageData.ID, "dmap"));
7 if (depthData.IsEmpty())
8 continue;
9
10 //1. limit maximum accumulated images resolution per sub-scene
11 const IIndex numPointsBegin(visibility.size());
12 const Camera camera(imageData.GetCamera(platforms, depthData.depthMap.size()));
13 for (int r=(ws%depthMapStep)/2; r<ws; r+=depthMapStep) {
14 for (int c=(ls%depthMapStep)/2; c<ls; c+=depthMapStep) {
15 const Depth depth = depthData.depthMap(r,c);
16 if (depth <= 0)
17 continue;
18 const Point3f& X = place_back(Cast<float>(camera.TransformPointI2W(Point3(c,r,depth))));
19 place_back(Footprint(camera, X)*areaScale);
20 place_back(idxImage);
21 }
22 }
23 imageAreas[idxImage] = visibility.size()-numPointsBegin;
24 }
3.1.3 构造⼋叉树
