题目地址:LA4795点这里
总算把白书的例题学完了,再开始坑习题吧
额 先wa了,找不出错误,把eps从1-10 改成1e-8 就ac了 囧...
首先求出质心。
然后枚举每一个平面,最好的简短代码长度的方式就是把点放在数组里用for语句吧~虽然暴力了一点。
然后,一定要细分和桌面的接触面是什么,如果有四点共面,那么我们要考虑的就是这四个点的凸包(其实只用调用PointIntri就好了)
没有四点共面时就看投影是否离边距离>=0.2
我的代码:
#include<iostream>
#include<cmath>
#include<cstdio>
#include<iomanip>
#include<ctime>
#include<cstdlib>
#include<vector>
using namespace std;
const double eps=1e-8;
// 1e-8?
const double PI=acos(-1.0);
int dcmp(double x)
{
return (x>eps)-(x<-eps);
}
struct Point3
{
double x,y,z;
Point3(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
};
typedef Point3 Vector3;
Vector3 operator+(Point3 A,Point3 B)
{
return Point3(A.x+B.x,A.y+B.y,A.z+B.z);
}
Vector3 operator-(Point3 A,Point3 B)
{
return Point3(A.x-B.x,A.y-B.y,A.z-B.z);
}
Vector3 operator*(Point3 A,double p)
{
return Point3(A.x*p,A.y*p,A.z*p);
}
Vector3 operator/(Point3 A,double p)
{
return Point3(A.x/p,A.y/p,A.z/p);
}
bool operator<(Point3 &A,Point3 &B)
{
return A.x<B.x||(A.x==B.x&&A.y<B.y)||(A.x==B.x&&A.y==B.y&&A.z<B.z);
}
bool operator==(Point3 &A,Point3 &B)
{
return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0&&dcmp(A.z-B.z)==0;
}
ostream & operator<<(ostream & out,Vector3 A)
{
out<<A.x<<' '<<A.y<<' '<<A.z<<endl;
return out;
}
Point3 read_point3()
{
Point3 p;
scanf("%lf%lf%lf",&p.x,&p.y,&p.z);
return p;
}
double Dot(Vector3 A,Vector3 B)
{
return A.x*B.x+A.y*B.y+A.z*B.z;
}
double Length(Vector3 A)
{
return sqrt(Dot(A,A));
}
double Angle(Vector3 A,Vector3 B)
{
return acos(Dot(A,B)/Length(A)/Length(B));
}
Vector3 Cross(Vector3 A,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(Point3 A,Point3 B,Point3 C)
{
return Length(Cross(B-A,C-A));
}
double DistanceToPlane(const Point3 &p,const Point3 & p0,const Vector3 &n)
{
return fabs(Dot(p-p0,n))/Length(n);
}
// n不用管方向
Point3 GetPlaneProjection(Point3 p,Point3 p0,Vector3 n)
{
return p-n*Dot(p-p0,n)/Dot(n,n);
}
// 三点式
Point3 GetPlaneProjection(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
Vector3 n=Cross(p1-p0,p2-p0);
return p-n*Dot(p-p0,n)/Dot(n,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Vector3 n)
{
double t=Dot(n,p0-p1)/Dot(n,p2-p1);
return p1+(p2-p1)*t;
}
// 点在三角形内部
bool PointInTri(Point3 p,Point3 p0,Point3 p1,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))==0;
}
// 线段与三角形相交 不考虑共面情形
bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 A,Point3 B,Point3 & p)
{
Vector3 n=Cross(p1-p0, p2-p0);
if(dcmp(Dot(n,B-A))==0) return false;
else{
double t=Dot(n,p0-A)/Dot(n,B-A);
if(dcmp(t)<0||dcmp(t-1)>0) return false;
else
{
p=A+(B-A)*t;
return PointInTri(p,p0,p1,p2);
}
}
}
double DistanceToLine(Point3 P,Point3 A,Point3 B)
{
Vector3 v1=B-A;
Vector3 v2=P-A;
return Length(Cross(v1, v2))/Length(v1);
}
double DistanceToSegment(Point3 P,Point3 A,Point3 B)
{
if(A==B) return Length(P-A);
Vector3 v1=B-A,v2=P-A,v3=P-B;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
else if(dcmp(Dot(v1,v3))>0) return Length(v3);
else return Length(Cross(v1,v2))/Length(v1);
}
double volume6(Point3 A,Point3 B,Point3 C,Point3 D)
{
return Dot(D-A,Cross(B-A,C-A));
}
// 0-1随机数
double rand01() { return rand()/(double)RAND_MAX;}
double randeps() {
//srand((unsigned) time(NULL));
return (rand01()-0.5)*eps;}
Point3 add_noise(Point3 P)
{
return Point3(P.x+randeps(),P.y+randeps(),P.z+randeps());
}
// 三维凸包 增量法
struct Face {
int v[3];
Face(int a, int b, int c) { v[0] = a; v[1] = b; v[2] = c; }
Vector3 Normal(const vector<Point3>& P) const {
return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
}
// f是否能看见P[i]
int CanSee(const vector<Point3>& P, int i) const {
return Dot(P[i]-P[v[0]], Normal(P)) > 0;
}
};
// 增量法求三维凸包
// 注意:没有考虑各种特殊情况(如四点共面)。实践中,请在调用前对输入点进行微小扰动
vector<Face> CH3D(const vector<Point3>& P) {
int n = P.size();
vector<vector<int> > vis(n);
for(int i = 0; i < n; i++) vis[i].resize(n);
vector<Face> cur;
cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动,前三个点不共线
cur.push_back(Face(2, 1, 0));
for(int i = 3; i < n; i++) {
vector<Face> next;
// 计算每条边的“左面”的可见性
for(int j = 0; j < cur.size(); j++) {
Face& f = cur[j];
int res = f.CanSee(P, i);
if(!res) next.push_back(f);
for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
}
for(int j = 0; j < cur.size(); j++)
for(int k = 0; k < 3; k++) {
int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线,左边对P[i]可见
next.push_back(Face(a, b, i));
}
cur = next;
}
return cur;
}
double DistanceToPlane( const Point3 &p,const Point3 & p0,const Point3 & p1,const Point3 & p2)
{
return fabs(volume6(p,p0,p1,p2)/Area2(p0, p1, p2));
}
// 四面体质心 就是四个点的平均
Point3 Centroid(const Point3 & A,const Point3 & B,const Point3 & C,const Point3 & D)
{
return (A+B+C+D)/4.0;
}
struct ConvexPolyhedron
{
int n; // 总顶点数
vector<Point3> P,P2; // P是真实顶点, P1是扰乱后的顶点
vector<Face> faces;
void init(int n)
{
this->n=n;
P.resize(n);
P2.resize(n);
for(int i=0;i<n;i++)
{
P[i]=read_point3();
P2[i]=add_noise(P[i]);
}
faces=CH3D(P2);
}
// 本题需要,写在类里面
double mindist(Point3 C)
{
double ans=1e30;
for(int i=0;i<faces.size();i++)
{
ans=min(ans,DistanceToPlane(C,P[faces[i].v[0]], P[faces[i].v[1]], P[faces[i].v[2]])); // 都是P这个数组
}
return ans;
}
Point3 centroid()
{
Point3 C=Point3(0,0,0);
Point3 tot=Point3(0,0,0);
double totv=0;
for(int i=0;i<faces.size();i++)
{
Point3 p0=P[faces[i].v[0]];
Point3 p1=P[faces[i].v[1]];
Point3 p2=P[faces[i].v[2]];
double v=volume6(p0, p1, p2, C);
tot=tot+Centroid(p0,p1,p2,C)*v;
totv+=v;
}
return tot/totv;
}
};
// self
// p0,p1在面A,B,C 的同侧
bool SameSide(Point3 p0,Point3 p1,Point3 A,Point3 B,Point3 C)
{
Vector3 n=Cross(B-A, C-A);
Vector3 v1=p0-A;
Vector3 v2=p1-A;
return Dot(v1, n)*Dot(v2,n)>=0;
}
// 改写读点函数
bool read_point3(Point3 & P)
{
if(scanf("%lf%lf%lf",&P.x,&P.y,&P.z)==3)
{
return 1;
}
else return 0;
}
Point3 P[6];
Point3 F;
// 求质心
Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{
return (A+B+C+D)/4.0;
}
Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D,Point3 &E) // 注意A,B,C 为四面体公共平面
{
double v1=volume6(A, B, C, D);
double v2=volume6(A, B, C, E);
v1=fabs(v1);
v2=fabs(v2);
Point3 C1=Centroid(A, B, C, D);
Point3 C2=Centroid(A, B, C, E);
return (C1*v1+C2*v2)/(v1+v2);
}
bool ok(Point3 &P,Point3 &A,Point3 & B,Point3 &C) // P为重心
{
Point3 prj=GetPlaneProjection(P, A, B,C);
if(PointInTri(prj, A, B, C))
{
double d1=DistanceToLine(prj, A, B);
double d2=DistanceToLine(prj, B, C);
double d3=DistanceToLine(prj, C, A);
return dcmp(d1-0.2)>=0&&dcmp(d2-0.2)>=0&&dcmp(d3-0.2)>=0;
}
else return 0;
}
bool ok(Point3 &P,Point3 &A,Point3 & B,Point3 &C,Point3 &D) // P为重心 A,B 为公共边
{
Point3 prj=GetPlaneProjection(P, A, B,C);
if(PointInTri(prj, A, B, C)||PointInTri(prj, A, B, D))
{
double d1=DistanceToLine(prj, C, A);
double d2=DistanceToLine(prj, C, B);
double d3=DistanceToLine(prj, D, A);
double d4=DistanceToLine(prj, D, B);
return dcmp(d1-0.2)>=0&&dcmp(d2-0.2)>=0&&dcmp(d3-0.2)>=0&&dcmp(d4-0.2)>=0;
}
else return 0;
}
bool OnSamePlane(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{
// Vector3 n1=Cross(B-A, C-A);
// Vector3 n2=Cross(B-A, D-A);
//
// Vector3 n1n2=Cross(n1, n2);
// Point3 Zero=Point3(0,0,0);
//
// if(n1n2==Zero) return 1;
// else return 0;
return dcmp(volume6(A, B, C, D))==0;
}
void check4(Point3 A,Point3 B,Point3 D,Point3 E,vector<double> &ans,Point3 &Center) // A,B 是A,B,C 的棱
{
if(OnSamePlane(A, B, D, E))
{
if(PointInTri(A, B, D, E))
{
if(ok(Center,B,D,E))
{
ans.push_back(DistanceToPlane(F, B, D, E));
}
}
if(PointInTri(B, A, D, E))
{
if(ok(Center,A,D,E))
{
ans.push_back(DistanceToPlane(F, A, D, E));
}
}
else
{
if(ok(Center,A,B,D,E))
{
ans.push_back(DistanceToPlane(F, A, B, D));
}
}
}
}
int main()
{
Point3 P_cin;
int cnt=0;
int cas=0;
while(read_point3(P_cin))
{
P[(cnt++)%6]=P_cin;
if(cnt%6==0)
{
F=P[5];
cas++;
double d1=1e30,d2=0;
//
Point3 Center=Centroid(P[0], P[1], P[2], P[3],P[4]);
vector<double> ans;
for(int i=0;i<5;i++)
for(int j=i+1;j<5;j++)
for(int k=j+1;k<5;k++)
{
Point3 A=P[i],B=P[j],C=P[k];
bool fine=1;
for(int a=0;a<5;a++)
for(int b=a+1;b<5;b++)
{
if(!SameSide(P[a], P[b],A, B, C))
{
fine=0;
break;
}
}
// 不允许共面
for(int i=0;i<5;i++)
{
if(P[i]==A||P[i]==B||P[i]==C) continue;
if(OnSamePlane(P[i], A, B, C))
{
fine=0;
break;
}
}
if(!fine) continue;
else
{
if(ok(Center,A,B,C))
{
ans.push_back(DistanceToPlane(F, A, B, C));
}
}
}
check4(P[0],P[1],P[3],P[4],ans, Center);
check4(P[1],P[2],P[3],P[4],ans, Center);
check4(P[2],P[0],P[3],P[4],ans, Center);
// for(int i=0;i<ans.size();i++)
// {
// cout<<"ans"<<ans[i]<<endl;
// }
for(int i=0;i<ans.size();i++)
{
if(ans[i]<d1)
{
d1=ans[i];
}
if(ans[i]>d2)
{
d2=ans[i];
}
}
printf("Case %d: %.5f %.5f\n",cas,d1,d2);
}
}
}
LiuRujia的代码:
// LA4795 Paperweight
// Rujia Liu
#include<cstdio>
#include<cmath>
using namespace std;
const double eps = 1e-8;
int dcmp(double x) {
if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}
struct Point3 {
double x, y, z;
Point3(double x=0, double y=0, double z=0):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);
}
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)); }
double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D) { return Dot(D-A, Cross(B-A, C-A)); }
bool read_point3(Point3& p) {
if(scanf("%lf%lf%lf", &p.x, &p.y, &p.z) != 3) return false;
return true;
}
// 点p到平面p0-n的距离。n必须为单位向量
double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n) {
return fabs(Dot(p-p0, n)); // 如果不取绝对值,得到的是有向距离
}
// 点p在平面p0-n上的投影。n必须为单位向量
Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n) {
return p-n*Dot(p-p0, n);
}
// 点P到直线AB的距离
double DistanceToLine(const Point3& P, const Point3& A, const Point3& B) {
Vector3 v1 = B - A, v2 = P - A;
return Length(Cross(v1, v2)) / Length(v1);
}
// p1和p2是否在线段a-b的同侧
bool SameSide(const Point3& p1, const Point3& p2, const Point3& a, const Point3& b) {
return dcmp(Dot(Cross(b-a, p1-a), Cross(b-a, p2-a))) >= 0;
}
// 点在三角形P0, P1, P2中
bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2) {
return SameSide(P, P0, P1, P2) && SameSide(P, P1, P0, P2) && SameSide(P, P2, P0, P1);
}
// 四面体的重心
Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D) {
return (A + B + C + D)/4.0;
}
#include<algorithm>
using namespace std;
// 判断P是否在三角形A, B, C中,并且到三条边的距离都至少为mindist。保证P, A, B, C共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, double mindist) {
if(!PointInTri(P, A, B, C)) return false;
if(DistanceToLine(P, A, B) < mindist) return false;
if(DistanceToLine(P, B, C) < mindist) return false;
if(DistanceToLine(P, C, A) < mindist) return false;
return true;
}
// 判断P是否在凸四边形ABCD(顺时针或逆时针)中,并且到四条边的距离都至少为mindist。保证P, A, B, C, D共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, const Point3& D, double mindist) {
if(!PointInTri(P, A, B, C)) return false;
if(!PointInTri(P, C, D, A)) return false;
if(DistanceToLine(P, A, B) < mindist) return false;
if(DistanceToLine(P, B, C) < mindist) return false;
if(DistanceToLine(P, C, D) < mindist) return false;
if(DistanceToLine(P, D, A) < mindist) return false;
return true;
}
int main() {
for(int kase = 1; ; kase++) {
Point3 P[5], F;
for(int i = 0; i < 5; i++)
if(!read_point3(P[i])) return 0;
read_point3(F);
// 求重心坐标
Point3 c1 = Centroid(P[0], P[1], P[2], P[3]);
Point3 c2 = Centroid(P[0], P[1], P[2], P[4]);
double vol1 = fabs(Volume6(P[0], P[1], P[2], P[3])) / 6.0;
double vol2 = fabs(Volume6(P[0], P[1], P[2], P[4])) / 6.0;
Point3 centroid = (c1 * vol1 + c2 * vol2) / (vol1 + vol2);
// 枚举放置方案
double mindist = 1e9, maxdist = -1e9;
for(int i = 0; i < 5; i++)
for(int j = i+1; j < 5; j++)
for(int k = j+1; k < 5; k++) {
// 找出另外两个点的下标a和b
int vis[5] = {0};
vis[i] = vis[j] = vis[k] = 1;
int a, b;
for(a = 0; a < 5; a++) if(!vis[a]) { b = 10-i-j-k-a; break; }
// 判断a和b是否在平面i-j-k的异侧
int d1 = dcmp(Volume6(P[i], P[j], P[k], P[a]));
int d2 = dcmp(Volume6(P[i], P[j], P[k], P[b]));
if(d1 * d2 < 0) continue; // 是,则放置方案不合法
Vector3 n = Cross(P[j]-P[i], P[k]-P[i]); // 法向量
n = n / Length(n); // 单位化
Point3 proj = GetPlaneProjection(centroid, P[i], n); // 重心在平面i-j-k上的投影
bool ok = InsideWithMinDistance(proj, P[i], P[j], P[k], 0.2);
if(!ok) {
if(d1 == 0) { // i-j-k-a四点共面。i和j一定为ABC三个顶点之一,k和a是D或者E
if(!InsideWithMinDistance(proj, P[i], P[k], P[j], P[a], 0.2)) continue;
} else if(d2 == 0) { // i-j-k-b四点共面。i和j一定为ABC三个顶点之一,k和b是D或者E
if(!InsideWithMinDistance(proj, P[i], P[k], P[j], P[b], 0.2)) continue;
} else
continue;
}
// 更新答案
double dist = DistanceToPlane(F, P[i], n);
mindist = min(mindist, dist);
maxdist = max(maxdist, dist);
}
printf("Case %d: %.5lf %.5lf\n", kase, mindist, maxdist);
}
return 0;
}