三角網格上高斯曲率和平均曲率

ref:https://www.cnblogs.com/eat-too-much/p/12595574.html

https://computergraphics.stackexchange.com/questions/1718/what-is-the-simplest-way-to-compute-principal-curvature-for-a-mesh-triangle
同學的CSDN部落格 https://blog.csdn.net/qq_38517015/article/details/105185241

中國科大傅孝明老師的框架：框架下載及配置執行

code

https://github.com/lishaohsuai/digital_geo

/**************************************************
@brief   : 輸出弧度制的角度 使用兩個向量
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
double MeshAlgorithm::vector3fAngle(vector3d a, vector3d b) {
	double cosTheta = a.dot(b) / (a.norm() * b.norm());
#ifndef debug
	std::cout << "[DEBUG] theta is " << acos(cosTheta) << std::endl;
	std::cout << "[DEBUG] vec a" << a.x << "," << a.y << "," << a.z << std::endl;
	std::cout << "[DEBUG] vec b" << b.x << "," << b.y << "," << b.z << std::endl;
#endif
	return acos(cosTheta);
}


/**************************************************
@brief   : 應用海倫公式計算面積
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
double MeshAlgorithm::areaUseThreePoints(Point3d a, Point3d b, Point3d c) {
	//應用海倫公式   S=1/4sqrt[(a+b+c)(a+b-c)(a+c-b)(b+c-a)]
	double lenA = sqrt(pow(b.x - c.x, 2) + pow(b.y - c.y, 2) + pow(b.z - c.z, 2));// b - c 兩點的座標
	double lenB = sqrt(pow(a.x - c.x, 2) + pow(a.y - c.y, 2) + pow(a.z - c.z, 2));// a - c 兩點的座標
	double lenC = sqrt(pow(b.x - a.x, 2) + pow(b.y - a.y, 2) + pow(b.z - a.z, 2));// a - b 兩點的座標
	double Area = 1.0 / 4.0 * sqrt((lenA + lenB + lenC) * (lenA + lenB - lenC) * (lenA + lenC - lenB) * (lenB + lenC - lenA));
#ifndef debug
	std::cout << "[DEBUG] area " << Area << std::endl;
	std::cout << "[DEBUG] a" << a.x << " " << a.y << " " << a.z << std::endl;
	std::cout << "[DEBUG] b" << b.x << " " << b.y << " " << b.z << std::endl;
	std::cout << "[DEBUG] c" << c.x << " " << c.y << " " << c.z << std::endl;
#endif
	return Area;
}


/**************************************************
@brief   : 計算三角形的面積通過兩個向量
@author  : lee
@input   : none
@output  : none
@time    : none
***************************************************/
double MeshAlgorithm::areaUseVector(vector3d a, vector3d b) {
	double area = 0.0;
	area = (a.crossp(b)).norm() * 0.5f;
	return fabs(area) > std::numeric_limits<double>::min() ? area : 1e-8;
}


/**************************************************
@brief   : 通過三個點計算其外接圓的圓心的座標
           參考 https://www.zhihu.com/question/40422123  李玉昆
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
Point3d MeshAlgorithm::calPointCircumcircle(Point3d a, Point3d b, Point3d c) {

	double a1, b1, c1, d1;
	double a2, b2, c2, d2;
	double a3, b3, c3, d3;

	double x1 = a.x, y1 = a.y, z1 = a.z;
	double x2 = b.x, y2 = b.y, z2 = b.z;
	double x3 = c.x, y3 = c.y, z3 = c.z;

	a1 = (y1*z2 - y2 * z1 - y1 * z3 + y3 * z1 + y2 * z3 - y3 * z2);
	b1 = -(x1*z2 - x2 * z1 - x1 * z3 + x3 * z1 + x2 * z3 - x3 * z2);
	c1 = (x1*y2 - x2 * y1 - x1 * y3 + x3 * y1 + x2 * y3 - x3 * y2);
	d1 = -(x1*y2*z3 - x1 * y3*z2 - x2 * y1*z3 + x2 * y3*z1 + x3 * y1*z2 - x3 * y2*z1);

	a2 = 2 * (x2 - x1);
	b2 = 2 * (y2 - y1);
	c2 = 2 * (z2 - z1);
	d2 = x1 * x1 + y1 * y1 + z1 * z1 - x2 * x2 - y2 * y2 - z2 * z2;

	a3 = 2 * (x3 - x1);
	b3 = 2 * (y3 - y1);
	c3 = 2 * (z3 - z1);
	d3 = x1 * x1 + y1 * y1 + z1 * z1 - x3 * x3 - y3 * y3 - z3 * z3;
	Point3d rlt;
	rlt.x = -(b1*c2*d3 - b1 * c3*d2 - b2 * c1*d3 + b2 * c3*d1 + b3 * c1*d2 - b3 * c2*d1)
		/ (a1*b2*c3 - a1 * b3*c2 - a2 * b1*c3 + a2 * b3*c1 + a3 * b1*c2 - a3 * b2*c1);
	rlt.y = (a1*c2*d3 - a1 * c3*d2 - a2 * c1*d3 + a2 * c3*d1 + a3 * c1*d2 - a3 * c2*d1)
		/ (a1*b2*c3 - a1 * b3*c2 - a2 * b1*c3 + a2 * b3*c1 + a3 * b1*c2 - a3 * b2*c1);
	rlt.z = -(a1*b2*d3 - a1 * b3*d2 - a2 * b1*d3 + a2 * b3*d1 + a3 * b1*d2 - a3 * b2*d1)
		/ (a1*b2*c3 - a1 * b3*c2 - a2 * b1*c3 + a2 * b3*c1 + a3 * b1*c2 - a3 * b2*c1);

	return rlt;
}



/**************************************************
@brief   : 根據外接圓的圓心是否要修正，做合理的修正
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
Point3d MeshAlgorithm::calPointVoronoiMixed(Point3d a, Point3d b, Point3d c) {
	Point3d p = calPointCircumcircle(a, b, c);//計算外接圓的圓心
	if (!isInTriangle(a, b, c, p)) {//不在三角形內部的時候 返回bc的中點
		return { (b.x + c.x) / 2.0, (b.y + c.y) / 2.0, (b.z + c.z) / 2.0 };
	}
	return p;
}

/**************************************************
@brief   : 判斷兩個向量是否同向
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
bool sameSide(Point3d a, Point3d b, Point3d c, Point3d p) {
	vector3d AB(b.x - a.x, b.y - a.y, b.z - a.z);
	vector3d AC(c.x - a.x, c.y - a.y, c.z - a.z);
	vector3d AP(p.x - a.x, p.y - a.y, p.z - a.z);

	vector3d v1 = AB.crossp(AC);
	vector3d v2 = AB.crossp(AP);

	// v1 and v2 should point to the same direction
	return v1.dot(v2) >= 0;
}

/**************************************************
@brief   : 判斷點是否在三角形內  參考連結 https://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
bool MeshAlgorithm::isInTriangle(Point3d a, Point3d b, Point3d c, Point3d p) {
	return sameSide(a, b, c, p) &&
		sameSide(b, c, a, p) &&
		sameSide(c, a, b, p);
}



/**************************************************
@brief   : 計算兩個三角形的對標 cotalpha cotbelta
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
void MeshAlgorithm::calCotAlphaCotBeta(Point3d p, Point3d a, Point3d b, Point3d c, double &cotAlpha, double &cotBeta) {
	vector3d AP(p.x - a.x, p.y - a.y, p.z - a.z);
	vector3d AB(b.x - a.x, b.y - a.y, b.z - a.z);
	double cosAlpha = cos(vector3fAngle(AP, AB));
	cotAlpha = cosAlpha / sqrt(1 - cosAlpha * cosAlpha);
	vector3d CP(p.x - c.x, p.y - c.y, p.z - c.z);
	vector3d CB(b.x - c.x, b.y - c.y, b.z - c.z);
	double cosBeta = cos(vector3fAngle(CP, CB));
	cotBeta = cosBeta / sqrt(1 - cosBeta * cosBeta);
}


/**************************************************
@brief   : 顏色對映歸一化  
@author  : 王丹丹
@input   : none
@output  : none
@time    : none
**************************************************/
void ImageAlgorithm::normalize(std::vector<double> & val) {
	//找到最大值和最小值，然後對映到[0,1]
	double max = -10000.0, min = 10000.0;
	int n = val.size();
	for (int i = 0; i < n; i++){
		if (val[i] > max) max = val[i];
		if (val[i] < min) min = val[i];
	}
	double t = max - min;
	//需要討論一下相等的情況
	for (int i = 0; i < n; i++){
		val[i] = (val[i] - min) / t;
	}
}


/**************************************************
@brief   : 對映高斯曲率到 r g b 參照libigl中的相關程式碼
@author  : 王丹丹
@input   : none
@output  : none
@time    : none
**************************************************/
void ImageAlgorithm::colorMap(double gaussCur, double &r, double &g, double &b)
{
	const double rone = 0.8;
	const double gone = 1.0;
	const double bone = 1.0;
	double  x = gaussCur;
	x = (gaussCur < 0 ? 0 : (x > 1 ? 1 : x));
	//可以簡單地理解：紅色的曲率最大，藍色的最小
	if (x < 1. / 8.){
		r = 0;
		g = 0;
		b = bone * (0.5 + (x) / (1. / 8.)*0.5);
	}
	else if (x < 3. / 8.){
		r = 0;
		g = gone * (x - 1. / 8.) / (3. / 8. - 1. / 8.);
		b = bone;
	}
	else if (x < 5. / 8.){
		r = rone * (x - 3. / 8.) / (5. / 8. - 3. / 8.);
		g = gone;
		b = (bone - (x - 3. / 8.) / (5. / 8. - 3. / 8.));
	}
	else if (x < 7. / 8.){
		r = rone;
		g = (gone - (x - 5. / 8.) / (7. / 8. - 5. / 8.));
		b = 0;
	}
	else{
		r = (rone - (x - 7. / 8.) / (1. - 7. / 8.)*0.5);
		g = 0;
		b = 0;
	}
}

=======================================另一個檔案========================================================
/**************************************************
@brief   : 計算高斯曲率  參考 https://www.cnblogs.com/VVingerfly/p/4428722.html 中的離散公式
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
void MeshViewerWidget::GaussianCurvatureProcess(void) {
	// 遍歷所有點計算每個點的高斯曲率

	//std::ofstream cout("case1.txt");
	std::vector<double> gauss;
	for (Mesh::VertexIter v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it) {
		double value = calGaussianCurvature(v_it);
		gauss.push_back(value);
		//std::cout << "Gaussian " <<  value << std::endl;
	}
	myImageAlgorithm.normalize(gauss);
	int k = 0;
	for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it)
	{
		double r, g, b;
		myImageAlgorithm.colorMap(gauss[k++], r, g, b);
		mesh.set_color(*v_it, OpenMesh::Vec3uc(int(r * 255), int(g * 255), int(b * 255)));
	}
}


/**************************************************
@brief   : 計算平均曲率   參考    https://computergraphics.stackexchange.com/questions/1718/
											what-is-the-simplest-way-to-compute-principal-curvature-for-a-mesh-triangle
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
void MeshViewerWidget::MeanCurvatureProcess(void) {
	// 便利所有點計算每個點的平均曲率
	for (Mesh::VertexIter v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it) {
		std::cout << "MeanCurvature " << calMeanCurvature(v_it) << std::endl;
	}
}


/**************************************************
@brief   : 計算每個點的平均曲率
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
double MeshViewerWidget::calMeanCurvature(Mesh::VertexIter vertexIndex) {
	OpenMesh::Vec3d P = (mesh).point(*vertexIndex);//中心點座標
	Point3d vecPoint(P[0], P[1], P[2]);
	std::vector<Point3d> neighborPoints;// N_1 rings points
	for (Mesh::VertexOHalfedgeIter vo_it = mesh.voh_begin(*vertexIndex); vo_it != mesh.voh_end(*vertexIndex); ++vo_it) {//這個頂點所帶有的半邊迭代器
		OpenMesh::ArrayKernel::VertexHandle to_v = mesh.to_vertex_handle(*vo_it);
		OpenMesh::Vec3d toPoint = mesh.point(to_v);// 鄰接點
		neighborPoints.push_back({ toPoint[0],toPoint[1],toPoint[2] });
	}

	// 計算 voronoi 點集
	std::vector<Point3d> voronoiPoints;
	double angle = 0;
	vector3d sumVector(0, 0, 0);
	for (int i = 0; i < neighborPoints.size(); i++) {
		Point3d voronoiPoint = myMeshAlgorithm.calPointVoronoiMixed(vecPoint, neighborPoints[i], (i == (neighborPoints.size() - 1) ? neighborPoints[0] : neighborPoints[i + 1]));
		voronoiPoints.push_back(voronoiPoint);
		// 計算∑(cotα + cotβ)(Qi - P)
		double cotAlpha;
		double cotBeta;
		myMeshAlgorithm.calCotAlphaCotBeta(vecPoint, (i == 0 ? neighborPoints[neighborPoints.size() - 1] : neighborPoints[i - 1]), neighborPoints[i],
			(i == (neighborPoints.size() - 1) ? neighborPoints[0] : neighborPoints[i + 1]), cotAlpha, cotBeta);
		sumVector.x += (cotAlpha + cotBeta)*(neighborPoints[i].x - vecPoint.x);
		sumVector.y += (cotAlpha + cotBeta)*(neighborPoints[i].y - vecPoint.y);
		sumVector.z += (cotAlpha + cotBeta)*(neighborPoints[i].z - vecPoint.z);
	}
	// 採用海倫公式計算voronoi三角形的面積和
	double sumArea = 0;
	for (int i = 0; i < voronoiPoints.size(); i++) {
		sumArea += myMeshAlgorithm.areaUseThreePoints(vecPoint, voronoiPoints[i], (i == voronoiPoints.size() - 1) ? voronoiPoints[0] : voronoiPoints[i + 1]);
	}
	double lengthSumVector = sqrt(sumVector.x * sumVector.x + sumVector.y * sumVector.y + sumVector.z * sumVector.z);
	return 0.5f * (1.0 / (sumArea * 2) * lengthSumVector);
}


/**************************************************
@brief   : 計算一個點所對應的高斯曲率
@author  : lee
@input   : none
@output  : none
@time    : none
**************************************************/
double MeshViewerWidget::calGaussianCurvature(Mesh::VertexIter vertexIndex) {
	OpenMesh::Vec3d P = (mesh).point(*vertexIndex);//中心點座標
	Point3d vecPoint(P[0], P[1], P[2]);
	std::vector<Point3d> neighborPoints;// N_1 rings points
	for (Mesh::VertexOHalfedgeIter vo_it = mesh.voh_begin(*vertexIndex); vo_it != mesh.voh_end(*vertexIndex); ++vo_it) {//這個頂點所帶有的半邊迭代器
		OpenMesh::ArrayKernel::VertexHandle to_v = mesh.to_vertex_handle(*vo_it);
		OpenMesh::Vec3d toPoint = mesh.point(to_v);// 鄰接點
		neighborPoints.push_back({ toPoint[0],toPoint[1],toPoint[2] });
	}
	
	// 計算 voronoi 點集
	std::vector<Point3d> voronoiPoints;
	double angle = 0;
	for (int i = 0; i < neighborPoints.size(); i++) {
		Point3d voronoiPoint = myMeshAlgorithm.calPointVoronoiMixed(vecPoint, neighborPoints[i], (i == (neighborPoints.size() - 1) ? neighborPoints[0] : neighborPoints[i + 1]));
		voronoiPoints.push_back(voronoiPoint);
		vector3d P1(neighborPoints[i].x - vecPoint.x, neighborPoints[i].y - vecPoint.y, neighborPoints[i].z - vecPoint.z);
		Point3d tmp = (i == (neighborPoints.size() - 1) ? neighborPoints[0] : neighborPoints[i + 1]);
		vector3d P2(tmp.x - vecPoint.x, tmp.y - vecPoint.y, tmp.z - vecPoint.z);
		angle += myMeshAlgorithm.vector3fAngle(P1, P2);
	}
	// 採用海倫公式計算voronoi三角形的面積和
	double sumArea = 0;
	for (int i = 0; i < voronoiPoints.size(); i++) {
		sumArea += myMeshAlgorithm.areaUseThreePoints(vecPoint, voronoiPoints[i], (i == voronoiPoints.size() - 1) ? voronoiPoints[0] : voronoiPoints[i + 1]);
	}
	// 計算選中點1鄰域角度的集合

	return 1.0 / sumArea * (2 * pi - angle);
}