凸组合参数化

凸组合参数化

框架：

用的是中国科大傅孝明老师的框架：框架下载 安装教程

推导

最常用的一个方法就是Tutte’s embedding algorithm，将一个三角形网格面参数化至2D平面上。

首先将边界点固定在2D的正方形或者圆形上（凸多边形），可以根据不同的权值进行对内部点坐标的计算。

内点坐标的计算如下：

使用均匀权值：

均匀参数化是最早的固定边界参数化算法了。均匀参数化的基本思想是把网格的边界点映射到一个平面上的凸多边形,而内点坐标为其一环邻近点的凸组合，其权值取以它为中心点的一环的平均值。

根据坐标计算公式：

代码：

//找到第一个边界点放入boundry
void MeshViewerWidget::FindFirstBoundry() {
	for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); v_it++) {
		if (mesh.is_boundary(*v_it)) {
			boundry.push_back((*v_it).idx());
			break;
		}
	}
}



//将所有边界点放入boundry中
void MeshViewerWidget::FindAllBoundry() {
	int count = 0;
	//遍历第一个点周围的所有点，找到一个边界邻点
	for (auto vv_it = mesh.vv_begin(mesh.vertex_handle(boundry[0])); vv_it.is_valid(); vv_it++) {
		if (mesh.is_boundary(*vv_it)) {
			boundry.push_back((*vv_it).idx());
			++count;
			break;
		}
	}//此时boundry中有两个边界点的index
	//根据序号找到这两个点
	Mesh::VertexHandle v_hd = mesh.vertex_handle(boundry[1]);
	Mesh::VertexHandle v_begin = mesh.vertex_handle(boundry[0]);
    
	while (v_hd != v_begin)
	{
		//每次从数组的最后一个点开始找下一个点
		for (auto vv_it = mesh.vv_begin(mesh.vertex_handle(boundry[count])); vv_it.is_valid(); vv_it++) {
			//只要不与前一个相同就继续循环
			if (mesh.is_boundary(*vv_it) && (*vv_it).idx() != boundry[count - 1]) {
				boundry.push_back((*vv_it).idx());
				++count;
				break;
			}
		}
		v_hd = mesh.vertex_handle(boundry[count]);
	}
}


//设置边界
void MeshViewerWidget::SetBoundry(Eigen::MatrixXd &Boundry) {
	Boundry.resize(mesh.n_vertices(), 2);
	Boundry = Eigen::MatrixXd::Zero(mesh.n_vertices(), 2);
	double angle = 3.1415926 / (boundry.size()-1);  //根据点的数量将圆等分
	for (int i = 0; i < boundry.size() - 1; ++i) {   //固定边界点的坐标在圆周上
		//mesh2.add_vertex(Mesh::Point(sin(2 * i*angle), cos(2 * i*angle), 0));
		Boundry(boundry[i], 0) = cos(2 * i * angle);  //x
		Boundry(boundry[i], 1) = sin(2 * i * angle);  //y
		
	}
}


//计算A矩阵
void MeshViewerWidget::ComputeA(Eigen::MatrixXd &A) { 
	//重置A矩阵大小
	A.resize(mesh.n_vertices(), mesh.n_vertices());
	A = Eigen::MatrixXd::Zero(mesh.n_vertices(), mesh.n_vertices());
	//遍历所有点
	for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it) {
		A((*v_it).idx(), (*v_it).idx()) = 1;
		if (mesh.is_boundary(*v_it)) {
			continue;
		}
		int vvCount = 0;
		//与该点相邻的点的数量
		for (auto vv_it = mesh.vv_begin(*v_it); vv_it.is_valid(); ++vv_it) {
			++ vvCount;
		}
		double weight = -1.0 / double(vvCount);   //权重简单的认为是邻点数量的倒数
		for (auto vv_it = mesh.vv_begin(*v_it); vv_it.is_valid(); ++vv_it) {  //更新A矩阵的值
			A((*v_it).idx(), (*vv_it).idx()) = weight;
		}
	}
	//cout << "A:    " << A << endl;

}

//计算B矩阵，B矩阵就是2D平面上边界点的位置矩阵
void MeshViewerWidget::ComputeB(Eigen::MatrixXd &B, Eigen::MatrixXd &Boundry) {
	B.resize(mesh.n_vertices(), 2);
	B = Boundry;
	//cout << "B:    " << B << endl;
}


Mesh MeshViewerWidget::getParamMesh_Convex(void) {
	Eigen::MatrixXd A;
	Eigen::MatrixXd B;
	Eigen::MatrixXd Boundry;
	Eigen::MatrixXd NewCoordinate;
	Mesh temp;      //新的2D平面

	FindFirstBoundry();
	FindAllBoundry();
	SetBoundry(Boundry);
	ComputeA(A);
	ComputeB(B, Boundry);

	NewCoordinate.resize(mesh.n_vertices(), 2);

	//NewCoordinate = A.inverse() * B;  //矩阵逆求x
	//NewCoordinate = A.colPivHouseholderQr().solve(B);  //QR分解求x
	//NewCoordinate = A.llt().solve(B);   //cholesky求x
	Eigen::SparseMatrix<float> A_sparse(mesh.n_vertices(), mesh.n_vertices());
	Eigen::SparseMatrix<float> B_sparse(mesh.n_vertices(), 2);
	Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<float> > Solver_sparse;
	buildSparseMatrix(A_sparse, A, mesh.n_vertices(), mesh.n_vertices());
	buildSparseMatrix(B_sparse, B, mesh.n_vertices(), 2);
	// 设置迭代精度
	Solver_sparse.setTolerance(0.001);
	Solver_sparse.compute(A_sparse);
	//NewCoordinate 即为解
	NewCoordinate = Solver_sparse.solve(B_sparse);


	//std::cout << "NewCoordinate:" << NewCoordinate << std::endl;

	Mesh::VertexHandle *vhandle;    //点的迭代器指针
	vhandle = new Mesh::VertexHandle[mesh.n_vertices()];  //开辟的空间大小
	double my_pi = acos(-1);

	int i = 0;
	for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it) {
		vhandle[i++] = temp.add_vertex(Mesh::Point(NewCoordinate((*v_it).idx(), 0), NewCoordinate((*v_it).idx(), 1), 0));
	}

	for (auto f_it = mesh.faces_begin(); f_it != mesh.faces_end(); ++f_it) {
		std::vector<Mesh::VertexHandle>  face_vhandles;   //存放2D的每个三角形平面
		for (auto fv_it = mesh.fv_iter(*f_it); fv_it.is_valid(); ++fv_it) {  //遍历该面的所有点
			face_vhandles.push_back(vhandle[(*fv_it).idx()]);  //通过vhandle将3D的点转化到2D平面上
		}
		temp.add_face(face_vhandles);
	}
	return temp;
}

//构建稀疏矩阵
void buildSparseMatrix(Eigen::SparseMatrix<float> &A1_sparse, Eigen::MatrixXd A,int A_rows,int A_cols) {
	std::vector<Eigen::Triplet<float>> tripletlist;  //构建类型为三元组的vector
	for (int i = 0; i < A_rows; i++)
	{
		for (int j = 0; j < A_cols; j++)
		{
			if (std::fabs(A(i, j)) > 1e-6)
			{
				//按Triplet方式填充，速度快
				tripletlist.push_back(Eigen::Triplet<float>(i, j, A(i, j)));

				// 直接插入速度慢
				//A1_sparse.insert(i, j) = A(i, j);
			}
		}
	}
	A1_sparse.setFromTriplets(tripletlist.begin(), tripletlist.end());
	// 压缩优化矩阵
	A1_sparse.makeCompressed();
}

演示效果：