最小二乘网格优化

基于meshlab的最小二乘网格优化

参考资料：http://blog.csdn.net/hjimce/article/details/46505863

基本思路

简单的说就是每一个点的坐标是它邻接的点的坐标的平均值：
v = （v1 + v2 + v3 +……+vn） / n
v1到vn就是v点的n个邻接点，所有的关系用矩阵表示就是：

每一行乘以每一列就是上述的等式，可以得到aij为：

wij就是1/n，n为改点的邻接点个数

这个方法需要一些固定点以及所有点的邻接关系，那么固定点一般就是需要优化的一块网格的边界点

实现代码

基于meshlab的插件实现，参考了参考资料中的代码

//优化选中面
bool Coptimize::OptimizeSelectFace(MeshModel * pmesh)
{
    if (!pmesh || pmesh->cm.sfn < 1)
        return false;
    m_pcm = pmesh;
    clearNonManifoldFace();     //先去掉非流形边
    if (!getAllVertexIndexFromSelectFace())
        return false;
    std::vector<T> tripletList;
    vector<Eigen::VectorXd> RHSPos;//超静定方程组右边  
    if (!bulidLeftRightMarix(tripletList, RHSPos))
        return false;

    //最小二乘解超静定方程组  
    SparseMatrix<double> Ls(m_vecIndex.size(), m_vecIndex.size());
    Ls.setFromTriplets(tripletList.begin(), tripletList.end());
    SparseMatrix<double> ls_transpose = Ls.transpose();
    SparseMatrix<double> LsLs = ls_transpose* Ls;
    Eigen::SimplicialCholesky<SparseMatrix<double>>MatricesCholesky(LsLs);

    int nrelIndex = 0;
    for (int i = 0;i<3;i++)
    {
        Eigen::VectorXd xyzRHS = ls_transpose*RHSPos[i];
        Eigen::VectorXd xyz = MatricesCholesky.solve(xyzRHS);
        for (int j = 0;j<m_vecIndex.size(); ++j)
        {
            nrelIndex = m_vecIndex[j];
            m_pcm->cm.vert[nrelIndex].P()[i] = xyz[j];
        }
    }

    return true;
}

//构建矩阵
bool Coptimize::bulidLeftRightMarix(vector<Triplet<double>>& vectorL, vector<Eigen::VectorXd>& vectorR)
{
    bool bhasBorder = false;
    vectorL.clear();
    int nall = 0;
    set<int> setAround;
    vectorR.clear();
    vectorR.resize(3);
    for (size_t i = 0; i < 3; i++)
    {
        vectorR[i].resize(m_vecIndex.size());
        vectorR[i].setZero();
    }
    for (size_t i = 0; i < m_vecIndex.size(); ++i)
    {
        CVertexO* pcur = &(m_pcm->cm.vert[m_vecIndex[i]]);
        nall += 1;
        if (getAllAroundVertex(pcur, setAround))
        {
            nall += setAround.size();
        }
        else
            bhasBorder = true;

    }
    //填充
    vectorL.reserve(nall);
    int nborder = 0;
    for (size_t i = 0; i < m_vecIndex.size(); ++i)
    {
        vectorL.push_back(T(i, i, 1));
        if (!bhasBorder && i % 10 == 0)
        {
            Point3f pt = m_pcm->cm.vert[m_vecIndex[i]].P();
            vectorR[0][i] = pt.X();
            vectorR[1][i] = pt.Y();
            vectorR[2][i] = pt.Z();
            ++nborder;
            continue;
        }
        CVertexO* pcur = &(m_pcm->cm.vert[m_vecIndex[i]]);

        if (getAllAroundVertex(pcur, setAround))
        {
            float w = -1.0f / setAround.size();
            for (auto iter = setAround.begin(); iter != setAround.end(); ++iter)
            {
                int n = m_mapIndex[(*iter)];
                vectorL.push_back(T(i, n, w));
            }
        }
        else
        {
            Point3f pt = m_pcm->cm.vert[m_vecIndex[i]].P();
            vectorR[0][i] = pt.X();
            vectorR[1][i] = pt.Y();
            vectorR[2][i] = pt.Z();
            ++nborder;
        }
    }
    return nborder > 0;
}

//获取顶点的所有邻接点
bool Coptimize::getAllAroundVertex(CVertexO* pv, set<int>& setAround)
{
    setAround.clear();
    bool bbord = false;
    CFaceO* pf = pv->VFp();
    if (!pf)
        return false;   //孤立点
    int ni = -1;
    for (size_t i = 0; i < 3; ++i)
    {
        if (pf->V(i) == pv)
        {
            ni = i;
            break;
        }
    }
    face::Pos<CFaceO> nmf(pf, ni);
    face::Pos<CFaceO> pos(pf, ni);
    pos.FlipV();
    setAround.insert(pos.v->Index());
    CFaceO* pLast = pf;
    nmf.FlipF();
    nmf.FlipE();
    while (nmf.f != pLast && nmf.f != pf)
    {
        if (!nmf.f->IsS())
            bbord = true;
        pos = nmf;
        pos.FlipV();
        setAround.insert(pos.v->Index());

        nmf.FlipF();
        nmf.FlipE();
    }
    if (nmf.f != pf)
        bbord = true;
    if (setAround.size() < 3)
        bbord = true;
    return !bbord;
}

效果：

原始网格

选中部分最小二乘处理之后