圖形處理（十二）拉普拉斯網格優化、最小二乘網格模型光順

hjimce發表於2015-06-15
優化模型

拉普拉斯網格優化、最小二乘網格模型光順

原文地址：http://blog.csdn.net/hjimce/article/details/46505863 

作者：hjimce

看這篇博文前，請先參考我的另外一篇博文《圖形處理（三）簡單拉普拉斯網格變形-Siggraph 2004》學習拉普拉斯座標的相關理論知識。這裡要分享的paper，是通過拉普拉斯的方法實現三角網格模型的優化。如果你已經非常熟悉三角網格曲面的拉普拉斯相關理論，實現這篇paper也就非常容易了。網格曲面的拉普拉斯座標不但可以用於變形、光順，還可以用於優化，總之好處多多，你只要學會了這一招，那麼就可以學會這些演算法了。


一、優化原理


利用Laplacian 座標重建的方法進行網格光順，原理很簡單，最簡單的就是隻要把源網格模型Laplacian 座標δ的模長縮小，方向不變，就可以然後進行Laplacian 網格重建，就可以實現簡單的光順效果。

然而如果要進行網格優化呢？怎麼實現？大牛們告訴我們一個比較規則的網格模型一個特點：當網格曲面上任意頂點的區域性片中包含的所有三角面片都為等腰三角形時，該頂點的同一Laplacian 座標和餘切Laplacian 座標相等。當將上述結論由某一個三角面片擴充套件到整個模型表面時可以發現：如果所有的三角面片都接近於正三角形，所有頂點的同一Laplacian 座標和餘切Laplacian 座標接近相等。

ok，上面的原理便是paper的思想，只要你懂得了抓住這個思想，那麼演算法實現起來就容易了。

在三角網格曲面上，頂點vi的拉普拉斯座標定義為，vi點與其一鄰接頂點加權組合的差：


當然權重wij需要滿足歸一化


權值的選取很多，但比較常用的權值是均勻權、餘切權，其計算公式如下：



二、網格優化演算法實現

演算法原理主要是：

1、通過先求解網格曲面餘切權計算得到的Laplacian 座標δ

2、構建均勻權下的拉普拉斯矩陣A

3、然後求解AX=δ.

void CMeshOptimize::OptimizeMesh(TriMesh*Tmesh)
{
	Tmesh->need_neighbors();
	int vn=Tmesh->vertices.size();
	//拉普拉斯矩陣設定
	int count0=0;
	vector<int>begin_N(vn);
	for (int i=0;i<vn;i++)
	{   
		begin_N[i]=count0;
		count0+=Tmesh->neighbors[i].size()+1;
	}
	typedef Eigen::Triplet<double> T;
	std::vector<T> tripletList(count0+vn);
	for(int i=0;i<vn;i++)
	{
		int nei_n=Tmesh->neighbors[i].size();
		tripletList[begin_N[i]]=T(i,i,-1.0*nei_n);
		for (int k = 0;k<nei_n;++k) 
		{
			tripletList[begin_N[i]+k+1]=T(i,Tmesh->neighbors[i][k],1);
		}
	}
	//約束矩陣設定
	m_boundary_wieght=100*m_contrain_weight;
	for (int i=0;i<vn;i++)
	{
		if(Tmesh->is_bdy(i))tripletList[count0+i]=T(vn+i,i,m_boundary_wieght);
		else tripletList[count0+i]=T(vn+i,i,m_contrain_weight);
	}

	SparseMatrixType Ls(2*vn,vn); 
	Ls.setFromTriplets(tripletList.begin(), tripletList.end());

	//最小二乘解超靜定方程組
	SparseMatrixType ls_transpose=Ls.transpose();
	SparseMatrixType LsLs =ls_transpose* Ls;
	vector<Eigen::VectorXd> RHSPos;//超靜定方程組右邊
	Compute_RHS(RHSPos);
	Eigen::SimplicialCholesky<SparseMatrixType>MatricesCholesky(LsLs);
	#pragma omp parallel for
	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<vn;j++)
		{
			Tmesh->vertices[j][i]=xyz[j];
		}
	}
	
}
//設定方程組右邊項
void CMeshOptimize::Compute_RHS(vector<Eigen::VectorXd> &RHSPos)
{
	int vn=m_OptimizeMesh->vertices.size();
	m_OptimizeMesh->need_neighbors();
	m_OptimizeMesh->need_adjacentfaces();
	RHSPos.clear();
	RHSPos.resize(3);
	for (int i=0;i<3;i++)
	{
		RHSPos[i].resize(2*vn);
		RHSPos[i].setZero();
	}

	int fn=m_OptimizeMesh->faces.size();
	m_OptimizeMesh->need_adjacentedges();
	#pragma omp parallel for
	for (int i=0;i<vn;i++)
	{   
		vector<float>CotWeight;
		float SumWeight;
		CotangentWeights(m_OptimizeMesh,i,CotWeight,SumWeight);
		int nei_n=m_OptimizeMesh->neighbors[i].size();
		//歸一化
		vector<int>&a=m_OptimizeMesh->neighbors[i];
		vec ls;
		for (int j=0;j<nei_n;j++)
		{	
			ls=ls+CotWeight[j]*m_OptimizeMesh->vertices[a[j]];
		}
		ls=ls-SumWeight*m_OptimizeMesh->vertices[i];
		for (int j=0;j<3;j++)
		{
			RHSPos[j][i]=ls[j];
		}
	}
	for (int i=vn;i<2*vn;i++)
	{
		for (int j=0;j<3;j++)
		{   
			if(m_OptimizeMesh->is_bdy(i-vn))RHSPos[j][i]=m_OptimizeMesh->vertices[i-vn][j]*m_boundary_wieght;
			else RHSPos[j][i]=m_OptimizeMesh->vertices[i-vn][j]*m_contrain_weight;
		}
	}


}
//計算一階鄰近點的各自cottan權重
void CMeshOptimize::CotangentWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum)
{   
	int NeighborNumber=TMesh->neighbors[vIndex].size();
	vweight.resize(NeighborNumber);
	WeightSum=0;
	vector<int>&NeiV=TMesh->neighbors[vIndex];
	for (int i=0;i<NeighborNumber;i++)
	{
		int j_nei=NeiV[i];
		vector<int>tempnei;
		Co_neighbor(TMesh,vIndex,j_nei,tempnei);
		float cotsum=0.0;
		for (int j=0;j<tempnei.size();j++)
		{
			vec vivo=TMesh->vertices[vIndex]-TMesh->vertices[tempnei[j]];
			vec vjvo=TMesh->vertices[j_nei]-TMesh->vertices[tempnei[j]];
			float dotvector=vivo DOT vjvo;
			dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);
			cotsum+=dotvector;
		}
		vweight[i]=cotsum/2.0;
		WeightSum+=vweight[i];
	}
	for (int k=0;k<NeighborNumber;++k)
	{
		vweight[k]=NeighborNumber*vweight[k]/WeightSum;
	}
	WeightSum=NeighborNumber;

}
void CMeshOptimize::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)
{
	Tmesh->need_adjacentedges();
	vector<int>&u_id_ae=Tmesh->adjancetedge[u_id]; 
	int en=u_id_ae.size();
	Tedge Co_Edge;
	for (int i=0;i<en;i++)
	{
		Tedge &ae=Tmesh->m_edges[u_id_ae[i]];
		int opsi=ae.opposite_vertex(u_id);
		if (opsi==v_id)
		{
			Co_Edge=ae;
			break;
		}
	}
	for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)
	{
		TriMesh::Face af=Tmesh->faces[Co_Edge.m_adjacent_faces[i]];
		for (int j=0;j<3;j++)
		{
			if((af[j]!=u_id)&&(af[j]!=v_id))
			{
				co_neiv.push_back(af[j]);
			}
		}
	}
}

至此三角網格的優化可以說演算法講完了。因為演算法比較簡答，本篇文章篇幅比較小，所以在這篇博文中順便講一下最小二乘網格的相關概念。

參考文獻：Laplacian Mesh Optimization

三、最小二乘網格相關理論

這邊順便講一下最小二乘網格，最小二乘網格最初的概念來源於paper《Least-squares Meshes》，我最初看到最小二乘網格這個概念是在paper《基於最小二乘網格的模型變形演算法》中看到的，因為之前學習微分域的網格變形演算法的時候，基本上對每種變形演算法都有看過相關的paper，用最小二乘網格的方法實現網格變形，其實跟基於多解析度的網格變形演算法是一樣的，都是對低頻空間中的網格模型進行變形操作。

最小二乘網格模型又稱為全域性的拉普拉斯光順網格，就是通過把網格模型的拉普拉斯座標設定為0，然後求解拉普拉斯方程：


即：


其中權重wij為：



這樣重建求解的網格即為最小二乘網格，也是一個光順後的網格模型，因為該模型的拉普拉斯座標全部為0。

當然求解上面的方程還需要控制頂點，控制頂點的個數對效果的影響還是蠻大的，可以看一下下面這個圖：


總之就是控制頂點的個數越少，越是光順。

在paper《Least-squares Meshes》中還演示了通過給定的拓撲連結關係，進行網格補洞，與原網格模型的區別。


*****************作者：hjimce     聯絡qq：1393852684   更多資源請關注我的部落格：http://blog.csdn.net/hjimce                原創文章，轉載請保留本行資訊。*****************

參考文獻：

1、《Least-squares Meshes》

2、《Laplacian Mesh Optimization》

3、《基於最小二乘網格的模型變形演算法》

相關文章