一个三维模型，可以看做是一个无向图。顶点是图的点，边可以作为图的边，那么可以得到一个邻接矩阵，该矩阵叫做拉普拉斯矩阵，计算这个邻接矩阵的特征向量就是所谓的谱分析。

那么问题来了。如何定义边的权重？额，当然有不同的方法，看你需要分析的特性是什么。一般来说是边的两个对角的cotangle之和。另外，矩阵对角线元素是该行其它元素之和再取负。下一步就是计算该矩阵的特征向量了。如果计算了？我前面写过两篇文章讲述如何使用MATLAB计算稀疏矩阵的特征向量，现在我就是这么做的。首先，拉普拉斯矩阵是一个对称矩阵，而且是一个大型的稀疏矩阵。因此，只能使用稀疏矩阵来存储数据了。而且对称矩阵能保证一些特殊的性质，比如第一个特征向量是k*(1,1,1,1…)。第二个特征向量分布得比较有规律，但是通过实验发现，也许第二个也没有规律，可能是第三或者第四个，因为要利用特征向量的分布，需要综合多个特征向量。

下面是我计算拉普拉斯矩阵的代码，需要考虑角度大于90的情况，以及除0等等。

//a和b夹角的cot
	double CSEMesh::GetCotangent(double a, double b, double c)
	{
		if (a*b < 1e-8)
		{
			return 0;
		}
		else
		{
			double cosc = (a*a+b*b-c*c)/(2*a*b);

			if (cosc < -1)
				cosc = -1;
			if (cosc > 1)
				cosc = 1;

			double angle = acos(cosc);

			double tan_angle = tan(angle);

			if (fabs(tan_angle) < 1e-8)
				return 1 / 1e-8;
			else
				return 0.5 / tan_angle;
		}
	}


	void CSEMesh::ComputeLaplaceWeight()
	{
		for (Enriched_Model::Edge_iterator pEdge = m_pModel->edges_begin();
			pEdge != m_pModel->edges_end();
			pEdge++)
		{
			double cotA = 1.0, cotB = 1.0;
			double a, b, c;

			Enriched_Model::Edge_iterator pEdgeTmp = pEdge;
			if (!pEdgeTmp->is_border())
			{
				c = pEdgeTmp->length();
				b = pEdgeTmp->next()->length();
				a = pEdgeTmp->next()->next()->length();
				cotA = GetCotangent(a, b, c);
			}

			pEdgeTmp = pEdge->opposite();
			if (!pEdgeTmp->is_border())
			{
				c = pEdgeTmp->length();
				b = pEdgeTmp->next()->length();
				a = pEdgeTmp->next()->next()->length();
				cotB = GetCotangent(a, b, c);
			}

			if(cotA + cotB < 1e-7)
			{
				pEdge->weight() = 1e-7;
				pEdge->opposite()->weight() = 1e-7;
			}
			else
			{
				pEdge->weight() = cotA + cotB;
				pEdge->opposite()->weight() = cotA + cotB;
			}
		}
	}

接下来是计算特征向量的代码，因为混合着cgal等，不是很好分离，但是可以用来参考下。

void CSEMesh::ComputeEigen(int type)
{
		static vector siVec;

		int nRow = m_pModel->size_of_vertices();
		int nColum = nRow;
		m_eigenValues.resize(EIGEN_VECTOR_NUM);
		for (int i = 0; i < EIGEN_VECTOR_NUM; ++i)
		{
			m_eigenVectors[type][i].resize(nRow);
		}
		
		int nNum = m_pModel->size_of_halfedges();
		siVec.resize(nNum + nRow);

		int i = 0;
		for (Enriched_Model::Edge_iterator pEdge = m_pModel->edges_begin();
			pEdge != m_pModel->edges_end();
			pEdge++)
		{
			int beg = pEdge->opposite()->vertex()->tag();
			int end = pEdge->vertex()->tag();

			siVec[i].r = beg;
			siVec[i].c = end;
			siVec[i].v = -pEdge->weight();
			++i;

			siVec[i].r = end;
			siVec[i].c = beg;
			siVec[i].v = -pEdge->weight();
			++i;
		}
		assert(i == nNum);

		for(Enriched_Model::Vertex_iterator pVertex = m_pModel->vertices_begin(); 
			pVertex != m_pModel->vertices_end(); pVertex++)
		{
			Enriched_Model::Vertex::Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
			Enriched_Model::Vertex::Halfedge_around_vertex_circulator d = pHalfEdge;
			double sum = 0.0;

			CGAL_For_all(pHalfEdge, d)
			{
				sum += pHalfEdge->weight();
			}

			if (fabs(sum) < 1e-7)
			{
				sum = 0.0;
			}

			siVec[i].r = pVertex->tag();
			siVec[i].c = pVertex->tag();
			siVec[i].v = sum;
			++i;
		}
		assert(i == nNum + nRow);
		std::sort(siVec.begin(), siVec.end());

		CMatlabSparseMatrix sm(&siVec;[0], nNum + nRow, nRow, nColum);
		int nEigenNum = std::min(nRow, (int)EIGEN_VECTOR_NUM);
		sm.GetEigens(m_eigenValues, m_eigenVectors[type], nEigenNum);
}

关于如何实现MATLAB稀疏矩阵类，可以参考我其它的两篇文章。下面是第二个特征向量的分布效果，能够观察到是一个非常有规律的标量场。

您可能也喜欢：
使用Matlab计算稀疏矩阵的特征向量	使用C++调用MATLAB引擎计算稀疏矩阵的特征向量	使用CGAL计算SDF分割模型	使用CGAL实现平面分割三维模型	使用CGAL实现简单的模型分割
无觅