一、简介

之前无意中看到这样一种点云孔洞修补方法，感觉很有意思，就下载了作者的源码看看效果，具体流程如下所述：
1、首先，我们需要先定位网格上的孔洞。
2、然后计算孔洞边界顶点的中心位置(孔洞边界顶点位置的平均值)。边界上的顶点与中心位置上新创建的顶点进行连接，从而创建一个填充孔洞的表面补丁。在这个补丁中的三角形会被上采样(细分)，这样当我们进行补丁平滑时，我们可以有更多的顶点来操作。

注意，这里对原始网格的三角形进行了上采样，因为这使它非常容易将原始补丁和原始几何图形的顶点融合在一起。所谓“融合”，就是指确保原始网格的边界，和补丁的边界共享相同的顶点。

3、接下来，通过在补丁上求解方程 $Δ^{2}f$ =0来创建一个平滑的补丁。这需要离散拉普拉斯算子来实现这一点，而拉普拉斯算子在网格表面上的离散被称为拉普拉斯-贝尔特拉米算子。这里使用libigl库创建了离散拉普拉斯- beltrami运算子，然后使用该运算子求解方程 $Δ^{2}f$ =0。

4、最后，我们使用网格抽取算法对网格进行下采样，得到最终修补之后的结果。

二、实现代码

CMakeLists.txt

cmake_minimum_required(VERSION 3.1)
project(hole_fixer)

set(CMAKE_BUILD_TYPE "Release")
set(CMAKE_CXX_FLAGS_RELEASE "-Zi -Od")
set(CMAKE_EXE_LINKER_FLAGS "-Debug")
message("Build Type:" ${CMAKE_BUILD_TYPE})
message("Release mode:" ${CMAKE_CXX_FLAGS_RELEASE})
message("Linker mode:" ${CMAKE_EXE_LINKER_FLAGS})

include_directories(libigl/ libigl/eigen/)
add_executable(hole_fixer main.cpp)

这里对原始代码做了一下改动，让它更符合我的习惯。

main.cpp

#include <igl/boundary_loop.h>
#include <igl/cat.h>
#include <igl/colon.h>
#include <igl/slice.h>
#include <igl/upsample.h>
#include <igl/decimate.h>
#include <igl/shortest_edge_and_midpoint.h>
#include <igl/harmonic.h>
#include <igl/readOFF.h>
#include <igl/writeOFF.h>

using Eigen::VectorXi;
using Eigen::MatrixXd;
using Eigen::MatrixXi;
using Eigen::VectorXd;
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> Tri;

void printHelpExit() {
    printf("Invalid command line arguments specified!\n\n");

    printf("USAGE: hole_fixer [options]\n\n");

    printf("OPTIONS: \n");

    printf("  -in\t\t\ttarget mesh file in .off-format, with a hole\n");
    printf("  -out\t\t\toutput mesh file in .off-format\n");
    printf("  -outfaces\t\tHow many faces to decimate the mesh to\n");
    printf("  -upsample\t\tHow much upsampling to use when creating the patch\n");

    exit(1);
}

const char* findToken(const char* param, int argc, char* argv[]) {
    const char* token = nullptr;
    for (int i = 0; i < argc; ++i) {
        if (strcmp(argv[i], param) == 0) {
            if (i + 1 < argc) {
                token = argv[i + 1];
                break;
            }
        }
    }

    if (token == nullptr) {
        printf("Could not find command-line parameter %s\n", param);
        return nullptr;
    }

    return token;
}

const char* parseStringParam(const char* param, int argc, char* argv[]) {
    const char* token = findToken(param, argc, argv);
    return token;
}

bool parseIntParam(const char* param, int argc, char* argv[], unsigned int& out) {
    const char* token = findToken(param, argc, argv);
    if (token == nullptr)
        return false;

    int r = sscanf(token, "%u,", &out);
    if (r != 1 || r == EOF) {
        return false;
    }
    else {
        return true;
    }
}

int main(int argc, char *argv[])
{
    
    //
    // command line parsing.
    //
    const char* inFile = parseStringParam("-in", argc, argv);
    if (inFile == nullptr) printHelpExit();

    const char* outFile = parseStringParam("-out", argc, argv);
    if (outFile == nullptr) printHelpExit();

    unsigned int outFacesN;
    if (!parseIntParam("-outfaces", argc, argv, outFacesN)) printHelpExit();

    unsigned int upsampleN;
    if (!parseIntParam("-upsample", argc, argv, upsampleN)) printHelpExit();
    
    // original mesh vertices and indices. This is the original mesh, which has a hole.
    MatrixXd originalV;
    MatrixXi originalF;

    if (!igl::readOFF(inFile, originalV, originalF)) {
        printHelpExit();
    }

    VectorXi originalLoop; // indices of the boundary of the hole. 
    igl::boundary_loop(originalF, originalLoop);

    if (originalLoop.size() == 0) {
        printf("Mesh has no hole!");
        printHelpExit();
    }

    // upsample the original mesh. this makes fusing the original mesh with the patch much easier.
    igl::upsample(Eigen::MatrixXd(originalV), Eigen::MatrixXi(originalF), originalV, originalF, upsampleN);

    // compute boundary center.
    VectorXd bcenter(3);
    {
        VectorXi R = originalLoop;
        VectorXi C(3); C << 0, 1, 2;
        MatrixXd B;
        MatrixXd A = originalV;
        igl::slice(A, R, C, B);
        bcenter = (1.0f / originalLoop.size()) * B.colwise().sum();
    }

    // a flat patch that fills the hole.
    MatrixXd patchV = MatrixXd(originalLoop.size() + 1, 3); // patch will have an extra vertex for the center vertex.
    MatrixXi patchF = MatrixXi(originalLoop.size(), 3);

    {
        VectorXi R = originalLoop;
        VectorXi C(3); C << 0, 1, 2;
        MatrixXd A = originalV;
        MatrixXd temp1;
        igl::slice(A, R, C, temp1);

        MatrixXd temp2(1, 3);
        temp2 << bcenter(0), bcenter(1), bcenter(2);

        // patch vertices will be the boundary vertices, plus a center vertex. concat these together.
        igl::cat(1, temp1, temp2, patchV);

        // create triangles that connect boundary vertices and the center vertex.
        for (int i = 0; i < originalLoop.size(); ++i) {
            patchF(i, 2) = (int)originalLoop.size();
            patchF(i, 1) = i;
            patchF(i, 0) = (1 + i) % originalLoop.size();
        }

        // also upsample patch. patch and original mesh will have the same upsampling level now
        // making it trivial to fuse them together.
        igl::upsample(Eigen::MatrixXd(patchV), Eigen::MatrixXi(patchF), patchV, patchF, upsampleN);
    }

    // the final mesh, where the patch and original mesh has been fused together.
    std::vector<std::vector<double>> fusedV;
    std::vector<std::vector<int>> fusedF;

    int index = 0; // vertex index counter for the fused mesh.

                   // first, we add the upsampled patch to the fused mesh.
    {
        for (; index < patchV.rows(); ++index) {
            fusedV.push_back({ patchV(index, 0), patchV(index, 1), patchV(index, 2) });
        }

        int findex = 0;
        for (; findex < patchF.rows(); ++findex) {
            fusedF.push_back({ patchF(findex, 0), patchF(findex, 1), patchF(findex, 2) });
        }
    }

    // now, we fuse the patch and the original mesh together.
    {
        // remaps indices from the original mesh to the fused mesh.
        std::map<int, int> originalToFusedMap;

        for (int itri = 0; itri < originalF.rows(); ++itri) {

            int triIndices[3];
            for (int iv = 0; iv < 3; ++iv) {

                int triIndex = originalF(itri, iv);

                int ret;

                if (originalToFusedMap.count(triIndex) == 0) {
                    int foundMatch = -1;

                    // the vertices at the boundary are the same, for both the two meshes(patch and original mesh).
                    // we ensure that these vertices are not duplicated.
                    // this is also how we ensure that the two meshes are fused together.
                    for (int jj = 0; jj < patchV.rows(); ++jj) {
                        VectorXd u(3); u << fusedV[jj][0], fusedV[jj][1], fusedV[jj][2];
                        VectorXd v(3); v << originalV(triIndex, 0), originalV(triIndex, 1), originalV(triIndex, 2);

                        if ((u - v).norm() < 0.00001) {
                            foundMatch = jj;
                            break;
                        }
                    }

                    if (foundMatch != -1) {
                        originalToFusedMap[triIndex] = foundMatch;
                        ret = foundMatch;
                    }
                    else {
                        fusedV.push_back({ originalV(triIndex, 0), originalV(triIndex, 1), originalV(triIndex, 2) });
                        originalToFusedMap[triIndex] = index;
                        ret = index;
                        index++;
                    }
                }
                else {
                    ret = originalToFusedMap[triIndex];
                }

                triIndices[iv] = ret;
            }

            fusedF.push_back({
                triIndices[0],
                triIndices[1],
                triIndices[2] });


        }

    }

    MatrixXd fairedV(fusedV.size(), 3);
    MatrixXi fairedF(fusedF.size(), 3);
    // now we shall do surface fairing on the mesh, to ensure
    // that the patch conforms to the surrounding curvature.
    {

        for (int vindex = 0; vindex < fusedV.size(); ++vindex) {
            auto r = fusedV[vindex];

            fairedV(vindex, 0) = r[0];
            fairedV(vindex, 1) = r[1];
            fairedV(vindex, 2) = r[2];
        }

        for (int findex = 0; findex < fusedF.size(); ++findex) {
            auto r = fusedF[findex];

            fairedF(findex, 0) = r[0];
            fairedF(findex, 1) = r[1];
            fairedF(findex, 2) = r[2];
        }

        VectorXi b(fairedV.rows() - patchV.rows());
        MatrixXd bc(fairedV.rows() - patchV.rows(), 3);
        // setup the boundary conditions. This is simply the vertex positions of the vertices not part of the patch.
        for (int i = (int)patchV.rows(); i < (int)fairedV.rows(); ++i) {
            int jj = i - (int)patchV.rows();

            b(jj) = i;

            bc(jj, 0) = fairedV(i, 0);
            bc(jj, 1) = fairedV(i, 1);
            bc(jj, 2) = fairedV(i, 2);
        }

        MatrixXd Z;
        int k = 2;
        // surface fairing simply means that we solve the equation
        // Delta^2 f = 0
        // with appropriate boundary conditions.
        // this function igl::harmonic from libigl takes care of that.

        // note that this is pretty inefficient thought.
        // the only boundary conditions necessary are the 2-ring of vertices around the patch.
        // the rest of the mesh vertices need not be specified.
        // we specify the rest of the mesh for simplicity of code, but it is not strictly necessary,
        // and pretty inefficient, since we have to solve a LARGE matrix equation as a result of this.
        igl::harmonic(fairedV, fairedF, b, bc, k, Z);
        fairedV = Z;
    }

    // finally, we do a decimation step.
    MatrixXd finalV(fusedV.size(), 3);
    MatrixXi finalF(fusedF.size(), 3);
    VectorXi temp0; VectorXi temp1;
    igl::decimate(fairedV, fairedF, outFacesN, finalV, finalF, temp0, temp1);
    
    igl::writeOFF(outFile, finalV, finalF);
    std::cout << "执行成功!" << std::endl;
    system("pause");

    return 0;
}

全部的代码可以点击网址：https://github.com/Erkaman/hole_fixer。