HDU5733 tetrahedron（计算几何）

HDU5733 tetrahedron（计算几何）

（本题解由 ljw 提供）

题面

题意：给四个点，问是不是四面体，如果是四面体则给出内切球的球心坐标和球的半径。

范围：$x,y$ 坐标 $\in [-1e6,1e6]$

分析：判断四个点是否共面即可判断是不是四面体，这里可以利用行列式的值计算（不会的去复习高数下的前几章和线代的相关内容）。

我们令$S=S_{△ABC}+S_{△ABD}+S_{△ACD}+S_{△BCD}$

则可以得到内切球球心的表达式为：

$x=(\frac{S_{△ABC}\ast D.x+S_{△ABD}\ast C.x+S_{△ACD}\ast B.x+S_{△BCD}\ast A.x)}{S}$

$y=(\frac{S_{△ABC}\ast D.y+S_{△ABD}\ast C.y+S_{△ACD}\ast B.y+S_{△BCD}\ast A.y)}{S}$

$z=(\frac{S_{△ABC}\ast D.z+S_{△ABD}\ast C.z+S_{△ACD}\ast B.z+S_{△BCD}\ast A.z)}{S}$

（具体证明可以去找高中竞赛教辅）

之后即可利用内切球心到四面体的距离相等，利用四个面的面积和整个四面体的体积即可得到球的半径。($r\ast S=S_{△ABC}\ast h$，$h$为$D$到$S_{△ABC}$的距离)

Notice：不存在的时候输出为四个O，不是0

Code：

#include <cstdio>
#include <cmath>
#include <vector>
using namespace std;
struct Point{
	double x, y, z;
	Point(double x,double y,double z):x(x),y(y),z(z){}
	Point(){
	}
}point[4];
// 判断是否共面
int judge(Point point[])
{
	double x[4], y[4], z[4];
	for(int i = 1; i < 4; i++){
		x[i] = point[i].x - point[0].x;
		y[i] = point[i].y - point[0].y;
		z[i] = point[i].z - point[0].z;
	}
	return (x[1] * y[2] * z[3]) + (x[2] * y[3] * z[1]) + (x[3] * y[1] *z[2]) - (x[3] * y[2] * z[1]) - (y[3] * z[2] * x[1]) - (z[3] * x[2] * y[1]);
}
double dis(Point a, Point b){
	return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2) + pow(a.z - b.z, 2));
}
// 海伦公式
double area(Point a, Point b, Point c){
	double ab = dis(a, b);
	double ac = dis(a, c);
	double bc = dis(b, c);
	double p = (ab + ac + bc) / 2;
	return sqrt(p * (p - ab) * (p - ac) * (p - bc));
}
// 球心坐标
vector<double> cal(Point point[])
{
	vector <double> ans;
	ans.resize(3);
	for(int i = 0; i < 3; i++){
		for(int j = 0; j < 4; j++){
			switch(i){
				case 0:
					ans[i] += area(point[(j + 1) % 4], point[(j + 2) % 4], point[(j + 3) % 4]) * point[j].x; break;
				case 1:
					ans[i] += area(point[(j + 1) % 4], point[(j + 2) % 4], point[(j + 3) % 4]) * point[j].y; break;
				case 2:
					ans[i] += area(point[(j + 1) % 4], point[(j + 2) % 4], point[(j + 3) % 4]) * point[j].z; break;
			}
		}
	}
	return ans;
}
Point Cross(Point a, Point b)
{    
	return Point((a.y * b.z - a.z * b.y), -(a.x * b.z - b.x * a.z), (a.x * b.y - a.y * b.x));
}
double Dot(Point A, Point B)
{    
	return (A.x * B.x + A.y * B.y + A.z * B.z);
} 
int main()
{
	while(~scanf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &point[0].x, &point[0].y, &point[0].z, &point[1].x, &point[1].y, &point[1].z, &point[2].x, &point[2].y, &point[2].z, &point[3].x, &point[3].y, &point[3].z)){
		if(judge(point)){
			vector<double> ans = cal(point);
			double S = 0;
			for(int j = 0; j < 4; j++){
				S += area(point[(j + 1) % 4], point[(j + 2) % 4], point[(j + 3) % 4]);
			}
			Point n = Cross(Point(point[0].x - point[1].x, point[0].y - point[1].y, point[0].z - point[1].z), Point(point[2].x - point[1].x, point[2].y - point[1].y, point[2].z - point[1].z));
			
			Point AD = Point(point[0].x - point[3].x, point[0].y - point[3].y, point[0].z - point[3].z);
			//printf("%lf %lf %lf\n", n.x, n.y, n.z);
			double h = fabs(Dot(AD, n) / sqrt(Dot(n, n)));
			double s = area(point[0], point[1], point[2]);
			//printf("%f %f %f\n", s, S, h);
			printf("%.4f %.4f %.4f %.4f\n", ans[0] / S, ans[1] / S, ans[2] / S, s * h / S);
		} else {
			puts("O O O O");
		}
	}
	return 0;
}

【END】感谢观看