三维三角形质量计算

已知三角形三个点，计算三角形质量：

（1）常见公式

#include <vcg/space/point3.h>	
template<class P3ScalarType>
	P3ScalarType GetTriQuality(vcg::Point3<P3ScalarType> const& p0, vcg::Point3<P3ScalarType> const& p1, vcg::Point3<P3ScalarType> const& p2)
	{
		P3ScalarType a = (p1 - p0).Norm();
		P3ScalarType b = (p2 - p0).Norm();
		P3ScalarType c = (p1 - p2).Norm();
		P3ScalarType sum = (a + b + c) * 0.5; //half perimeter
		P3ScalarType area2 = sum * (a + b - sum) * (a + c - sum) * (b + c - sum);//area
		if (area2 <= 0)
			return 0;
		P3ScalarType maxlength = std::max(std::max(a, b), c);

		return 6.0 / sqrt(3.0) * area2 / (sum * maxlength);//Quality
	}

（2）其他公式

	/// Compute a shape quality measure of the triangle composed by points p0,p1,p2
	/// It Returns mean ratio 2sqrt(a, b)/(a+b) where a+b are the eigenvalues of the M^tM of the
	/// transformation matrix into a regular simplex
	/// the range is range [0, 1]
	template<class P3ScalarType>
	P3ScalarType QualityMeanRatio(vcg::Point3<P3ScalarType> const& p0,vcg::Point3<P3ScalarType> const& p1,vcg::Point3<P3ScalarType> const& p2) 
	{

		P3ScalarType a = (p1 - p0).Norm();
		P3ScalarType b = (p2 - p0).Norm();
		P3ScalarType c = (p1 - p2).Norm();
		P3ScalarType sum = (a + b + c) * 0.5; //semiperimeter
		P3ScalarType area2 = sum * (a + b - sum) * (a + c - sum) * (b + c - sum);
		if (area2 <= 0) 
			return 0;
		return (4.0 * sqrt(3.0) * sqrt(area2)) / (a * a + b * b + c * c);
	}

	/// Compute a shape quality measure of the triangle composed by points p0,p1,p2
	/// It Returns inradius/circumradius
	/// the range is range [0, 1]
	/// e.g. Equilateral triangle 1, halfsquare: 0.81, ... up to a line that has zero quality.
	template<class P3ScalarType>
	P3ScalarType QualityRadii(vcg::Point3<P3ScalarType> const& p0,vcg::Point3<P3ScalarType> const& p1,vcg::Point3<P3ScalarType> const& p2)
	{

		P3ScalarType a = (p1 - p0).Norm();
		P3ScalarType b = (p2 - p0).Norm();
		P3ScalarType c = (p1 - p2).Norm();

		P3ScalarType sum = (a + b + c) * 0.5;
		P3ScalarType area2 = sum * (a + b - sum) * (a + c - sum) * (b + c - sum);
		if (area2 <= 0) return 0;
		//circumradius: (a*b*c)/(4*sqrt(area2))
		//inradius: (a*b*c)/(4*circumradius*sum) => sqrt(area2)/sum;
		return (8 * area2) / (a * b * c * sum);
	}

	/// Compute a shape quality measure of the triangle composed by points p0,p1,p2
	/// It Returns 2*AreaTri/(MaxEdge^2),
	/// the range is range [0.0, 0.866]
	/// e.g. Equilateral triangle sqrt(3)/2, halfsquare: 1/2, ... up to a line that has zero quality.
	template<class P3ScalarType>
	P3ScalarType Quality(vcg::Point3<P3ScalarType> const& p0, vcg::Point3<P3ScalarType> const& p1, vcg::Point3<P3ScalarType> const& p2)
	{
		Point3<P3ScalarType> d10 = p1 - p0;
		Point3<P3ScalarType> d20 = p2 - p0;
		Point3<P3ScalarType> d12 = p1 - p2;
		Point3<P3ScalarType> x = d10 ^ d20;

		P3ScalarType a = Norm(x);
		if (a == 0) return 0; // Area zero triangles have surely quality==0;
		P3ScalarType b = SquaredNorm(d10);
		if (b == 0) return 0; // Again: area zero triangles have surely quality==0;
		P3ScalarType t = b;
		t = SquaredNorm(d20); if (b < t) b = t;
		t = SquaredNorm(d12); if (b < t) b = t;
		return a / b;
	}