题意：三维空间中，给出两个三角形的左边，问是否相交。
面积法判断点在三角形内：
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<memory.h>
#include<cstdlib>
#include<vector>
#define clc(a,b) memset(a,b,sizeof(a))
#define LL long long int
#define up(i,x,y) for(i=x;i<=y;i++)
#define w(a) while(a)
const double inf=0x3f3f3f3f;
const double PI = acos(-1.0);
using namespace std; struct Point3
{
double x, y, z;
Point3(double x=, double y=, double z=):x(x),y(y),z(z) { }
}; typedef Point3 Vector3; Vector3 operator + (const Vector3& A, const Vector3& B)
{
return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);
}
Vector3 operator - (const Point3& A, const Point3& B)
{
return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);
}
Vector3 operator * (const Vector3& A, double p)
{
return Vector3(A.x*p, A.y*p, A.z*p);
}
Vector3 operator / (const Vector3& A, double p)
{
return Vector3(A.x/p, A.y/p, A.z/p);
} const double eps = 1e-;
int dcmp(double x)
{
if(fabs(x) < eps) return ;
else return x < ? - : ;
} double Dot(const Vector3& A, const Vector3& B)
{
return A.x*B.x + A.y*B.y + A.z*B.z;
}
double Length(const Vector3& A)
{
return sqrt(Dot(A, A));
}
double Angle(const Vector3& A, const Vector3& B)
{
return acos(Dot(A, B) / Length(A) / Length(B));
}
Vector3 Cross(const Vector3& A, const Vector3& B)
{
return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x);
}
double Area2(const Point3& A, const Point3& B, const Point3& C)
{
return Length(Cross(B-A, C-A));
} Point3 read_point3()
{
Point3 p;
scanf("%lf%lf%lf", &p.x, &p.y, &p.z);
return p;
} // 点在三角形P0, P1, P2中
bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2)
{
double area1 = Area2(P, P0, P1);
double area2 = Area2(P, P1, P2);
double area3 = Area2(P, P2, P0);
return dcmp(area1 + area2 + area3 - Area2(P0, P1, P2)) == ;
} // 三角形P0P1P2是否和线段AB相交
bool TriSegIntersection(const Point3& P0, const Point3& P1, const Point3& P2, const Point3& A, const Point3& B, Point3& P)
{
Vector3 n = Cross(P1-P0, P2-P0);
if(dcmp(Dot(n, B-A)) == ) return false; // 线段A-B和平面P0P1P2平行或共面
else // 平面A和直线P1-P2有惟一交点
{
double t = Dot(n, P0-A) / Dot(n, B-A);
if(dcmp(t) < || dcmp(t-) > ) return false; // 不在线段AB上
P = A + (B-A)*t; // 交点
return PointInTri(P, P0, P1, P2);
}
} bool TriTriIntersection(Point3* T1, Point3* T2)
{
Point3 P;
for(int i = ; i < ; i++)
{
if(TriSegIntersection(T1[], T1[], T1[], T2[i], T2[(i+)%], P)) return true;
if(TriSegIntersection(T2[], T2[], T2[], T1[i], T1[(i+)%], P)) return true;
}
return false;
} int main()
{
int T;
scanf("%d", &T);
while(T--)
{
Point3 T1[], T2[];
for(int i = ; i < ; i++) T1[i] = read_point3();
for(int i = ; i < ; i++) T2[i] = read_point3();
printf("%d\n", TriTriIntersection(T1, T2) ? : );
}
return ;
}
用Barycentric坐标法判断点在三角形内（重心法）
#include<cstdio>
#include<cmath>
using namespace std; struct Point3
{
double x, y, z;
Point3(double x=, double y=, double z=):x(x),y(y),z(z) { }
}; typedef Point3 Vector3; Vector3 operator + (const Vector3& A, const Vector3& B)
{
return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);
}
Vector3 operator - (const Point3& A, const Point3& B)
{
return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);
}
Vector3 operator * (const Vector3& A, double p)
{
return Vector3(A.x*p, A.y*p, A.z*p);
}
Vector3 operator / (const Vector3& A, double p)
{
return Vector3(A.x/p, A.y/p, A.z/p);
} const double eps = 1e-;
int dcmp(double x)
{
if(fabs(x) < eps) return ;
else return x < ? - : ;
} double Dot(const Vector3& A, const Vector3& B)
{
return A.x*B.x + A.y*B.y + A.z*B.z;
}
double Length(const Vector3& A)
{
return sqrt(Dot(A, A));
}
double Angle(const Vector3& A, const Vector3& B)
{
return acos(Dot(A, B) / Length(A) / Length(B));
}
Vector3 Cross(const Vector3& A, const Vector3& B)
{
return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x);
}
double Area2(const Point3& A, const Point3& B, const Point3& C)
{
return Length(Cross(B-A, C-A));
} Point3 read_point3()
{
Point3 p;
scanf("%lf%lf%lf", &p.x, &p.y, &p.z);
return p;
} // 点在三角形P0, P1, P2中
// http://www.blackpawn.com/texts/pointinpoly/default.html
bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2)
{
Vector3 v0 = P2 - P0;
Vector3 v1 = P1 - P0;
Vector3 v2 = P - P0; // Compute dot products
double dot00 = Dot(v0, v0);
double dot01 = Dot(v0, v1);
double dot02 = Dot(v0, v2);
double dot11 = Dot(v1, v1);
double dot12 = Dot(v1, v2); // Compute barycentric coordinates
double invDenom = / (dot00 * dot11 - dot01 * dot01);
double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
double v = (dot00 * dot12 - dot01 * dot02) * invDenom; // Check if point is in triangle
return (dcmp(u) >= ) && (dcmp(v) >= ) && (dcmp(u + v - ) <= );
} // 三角形P0P1P2是否和线段AB相交
bool TriSegIntersection(const Point3& P0, const Point3& P1, const Point3& P2, const Point3& A, const Point3& B, Point3& P)
{
Vector3 n = Cross(P1-P0, P2-P0);
if(dcmp(Dot(n, B-A)) == ) return false; // 线段A-B和平面P0P1P2平行或共面
else // 平面A和直线P1-P2有惟一交点
{
double t = Dot(n, P0-A) / Dot(n, B-A);
if(dcmp(t) < || dcmp(t-) > ) return false; // 不在线段AB上
P = A + (B-A)*t; // 交点
return PointInTri(P, P0, P1, P2);
}
} bool TriTriIntersection(Point3* T1, Point3* T2)
{
Point3 P;
for(int i = ; i < ; i++)
{
if(TriSegIntersection(T1[], T1[], T1[], T2[i], T2[(i+)%], P)) return true;
if(TriSegIntersection(T2[], T2[], T2[], T1[i], T1[(i+)%], P)) return true;
}
return false;
} int main()
{
int T;
scanf("%d", &T);
while(T--)
{
Point3 T1[], T2[];
for(int i = ; i < ; i++) T1[i] = read_point3();
for(int i = ; i < ; i++) T2[i] = read_point3();
printf("%d\n", TriTriIntersection(T1, T2) ? : );
}
return ;
}
证明可以参考http://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html
三角形的三个点在同一个平面上，如果选中其中一个点，其他两个点不过是相对该点的位移而已，比如选择点A作为起点，那么点B相当于在AB方向移动一段距离得到，而点C相当于在AC方向移动一段距离得到。
所以对于平面内任意一点，都可以由如下方程来表示
P = A + u * (C – A) + v * (B - A) // 方程1
如果系数u或v为负值，那么相当于朝相反的方向移动，即BA或CA方向。那么如果想让P位于三角形ABC内部，u和v必须满足什么条件呢？有如下三个条件
u >= 0
v >= 0
u + v <= 1
几个边界情况，当u = 0且v = 0时，就是点A，当u = 0,v = 1时，就是点B，而当u = 1, v = 0时，就是点C
整理方程1得到P – A = u(C - A) + v(B - A)
令v0 = C – A, v1 = B – A, v2 = P – A，则v2 = u * v0 + v * v1，现在是一个方程，两个未知数，无法解出u和v，将等式两边分别点乘v0和v1的到两个等式
(v2) • v0 = (u * v0 + v * v1) • v0
(v2) • v1 = (u * v0 + v * v1) • v1
注意到这里u和v是数，而v0，v1和v2是向量，所以可以将点积展开得到下面的式子。
v2 • v0 = u * (v0 • v0) + v * (v1 • v0) // 式1
v2 • v1 = u * (v0 • v1) + v * (v1• v1) // 式2
解这个方程得到
u = ((v1•v1)(v2•v0)-(v1•v0)(v2•v1)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0))
v = ((v0•v0)(v2•v1)-(v0•v1)(v2•v0)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0))