最小二乘网格优化

基于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;
}

效果：

原始网格



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