2017 山东一轮集训 Day2 Shadow (三维凸包点在面上投影)

在三维坐标中，给定一个点光源，一个凸多面体，以及一个平面作为地面。

求该凸多面体在地面上阴影的面积。

这三个点共同确定了一个平面，这个平面就是地面。保证这三个点坐标互异且不共线。前三行每行三个实数，每行表示一个点。这三个点共同确定了一个平面，这个平面就是地面。保证这三个点坐标互异且不共线。
接下来一行三个实数，表示一个点。这个点就是点光源。
之后一个整数n，表示凸多面体顶点的数量。
之后n行，每行三个实数，表示凸多面体的一个顶点。

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cmath>
  4 #include <cstring>
  5 #include <vector>
  6 #include <algorithm>
  7 #include <queue>
  8 #include <stack>
  9 #include <map>
 10 #include <set>
 11 #include <utility>
 12 
 13 #ifdef DEBUG
 14     const int MAXN = 110;
 15 #else
 16     const int MAXN = 110;
 17 #endif
 18 
 19 const double eps = 1e-6;
 20 
 21 using namespace std;
 22 
 23 struct Vector {
 24     double x, y, z;
 25     Vector(double _x = 0, double _y = 0, double _z = 0): x(_x), y(_y), z(_z) {}
 26     Vector operator+(const Vector &rhs) const {
 27         return Vector(x + rhs.x, y + rhs.y, z + rhs.z);
 28     }
 29     Vector operator-(const Vector &rhs) const {
 30         return Vector(x - rhs.x, y - rhs.y, z - rhs.z);
 31     }
 32     Vector operator*(double k) {
 33         return Vector(x * k, y * k, z * k);
 34     }
 35     double len() {
 36         return sqrt(x * x + y * y + z * z);
 37     }
 38 };
 39 typedef Vector Point;
 40 
 41 inline double dot(const Vector &a, const Vector &b) {
 42     return a.x * b.x + a.y * b.y + a.z * b.z;
 43 }
 44 
 45 inline Vector cross(const Vector &a, const Vector &b) {
 46     return Vector(a.y * b.z - b.y * a.z, a.z * b.x - b.z * a.x, a.x * b.y - a.y * b.x);
 47 }
 48 
 49 bool flag;
 50 inline Vector readv() {
 51     double x, y, z;
 52     scanf("%lf%lf%lf", &x, &y, &z);
 53     if (flag) swap(y, z);
 54     return Vector(x, y, z);
 55 }
 56 
 57 inline void printv(const Vector &a) {
 58     printf("(%lf, %lf, %lf)
", a.x, a.y, a.z);
 59 }
 60 
 61 struct Line {
 62     Point s;
 63     Vector d;
 64     Line() {}
 65     Line(Point _s, Point _t): s(_s), d(_t - _s) {}
 66 };
 67 
 68 struct Plain {
 69     Point s;
 70     Vector d;
 71     Plain() {}
 72     Plain(Point a, Point b, Point c): s(a), d(cross(b - a, c - a)) {}
 73 } p;
 74 Point s;
 75 int n;
 76 Point pol[MAXN];
 77 Point sd[MAXN];
 78 
 79 inline int dcmp(double x) {
 80     return x < -eps ? -1 : x > eps ? 1 : 0;
 81 }
 82 
 83 inline Point itsct(Line l, Plain p) {
 84     double s1 = dot(p.d, l.s - p.s);
 85     double s2 = dot(p.d, l.s + l.d - p.s);
 86     if (!dcmp(s1 - s2)) {
 87         return itsct(Line(l.s, l.s + p.d), p);
 88     }
 89     double k = s1 / (s1 - s2);
 90     return l.s + l.d * k;
 91 }
 92 
 93 bool cmp(const Point &a, const Point &b) {
 94     return a.x < b.x;
 95 }
 96 
 97 Point conv[MAXN];
 98 inline int make_convex() {
 99     int cnt = 0;
100     sort(sd, sd + n, cmp);
101     conv[cnt++] = sd[0];
102     for (int i = 1; i < n; i++) {
103         while (cnt > 1 && dcmp(cross(sd[i] - conv[cnt - 1], conv[cnt - 1] - conv[cnt - 2]).z) < 0) cnt--;
104         if (dcmp((conv[cnt - 1] - sd[i]).len())) conv[cnt++] = sd[i];
105     }
106     for (int i = n - 2; i >= 0; i--) {
107         while (cnt > 1 && dcmp(cross(sd[i] - conv[cnt - 1], conv[cnt - 1] - conv[cnt - 2]).z) < 0) cnt--;
108         if (dcmp((conv[cnt - 1] - sd[i]).len())) conv[cnt++] = sd[i];
109     }
110     return cnt;
111 }
112 
113 inline void open_file() {
114     freopen("shadow.in", "r", stdin);
115     freopen("shadow.out", "w", stdout);
116 }
117 
118 int main() {
119     //open_file();
120     Point a = readv(), b = readv(), c = readv();
121     if (!dcmp(cross(b - a, c - a).z)) {
122         swap(a.y, a.z), swap(b.y, b.z), swap(c.y, c.z);
123         flag = true;
124     }
125     p = Plain(a, b, c);
126     s = readv();
127     scanf("%d", &n);
128     for (int i = 0; i < n; i++) pol[i] = readv();
129     for (int i = 0; i < n; i++) sd[i] = itsct(Line(s, pol[i]), p);
130     int cnt = make_convex();
131     double ans = 0;
132     for (int i = 1; i < cnt - 2; i++) {
133         ans += cross(conv[i] - conv[0], conv[i + 1] - conv[0]).len();
134     }
135     printf("%.2lf
", ans / 2);
136     
137     return 0;
138 }

直接给你平面方程

  1 #include <algorithm>
  2 #include <iostream>
  3 #include <cstdlib>
  4 #include <cstdio>
  5 #include <cmath>
  6 
  7 using namespace std;
  8 
  9 typedef struct p1node
 10 {
 11     double a,b,c,d;
 12 }plane;
 13 plane Pl;
 14 
 15 typedef struct p2node
 16 {
 17     double x,y,z;
 18 }point;
 19 point temp;
 20 point P[ 105 ];
 21 point S[ 105 ];
 22 
 23 //计算点在平面上的投影 
 24 int shadow( plane p, point s, int n )
 25 {
 26     //求出过s平行于plane的平面 ax+by+c=D
 27     double D = p.a*s.x+p.b*s.y+p.c*s.z;
 28     if ( D-p.d < 0 ) {//调整方向 
 29         p.a *= -1;p.b *= -1;p.c *= -1;
 30         p.d *= -1;
 31         D *= -1; 
 32     }
 33     
 34     //判断点与面的关系 
 35     int count = 0;
 36     for ( int i = 0 ; i < n ; ++ i ) {
 37         double det = p.a*P[i].x+p.b*P[i].y+p.c*P[i].z-D;
 38         if ( det < 0 ) count ++;
 39     }
 40     if ( count == 0 ) return 0;
 41     if ( count != n ) return n+1;
 42         
 43     for ( int i = 0 ; i < n ; ++ i ) {
 44         //直线方程: (Sx,Sy,Sz) + t(dx,dy,dz) 
 45         double dx = P[i].x - s.x;
 46         double dy = P[i].y - s.y;
 47         double dz = P[i].z - s.z;
 48         double t = (p.d-p.a*s.x-p.b*s.y-p.c*s.z)/(p.a*dx+p.b*dy+p.c*dz);
 49         
 50         P[i].x = s.x + t*dx;
 51         P[i].y = s.y + t*dy;
 52         P[i].z = s.z + t*dz;
 53     }
 54     return n;
 55 }
 56 
 57 //坐标系旋转，到x-y平面 
 58 void change( plane p, int n )
 59 {
 60     //平行于x-y平面的不用计算，也不能计算（分母为0）
 61     if  ( p.a*p.a + p.b*p.b == 0 ) return;
 62     for ( int i = 0 ; i < n ; ++ i ) {
 63         //绕z轴旋转 
 64         double cosC = p.b/sqrt(p.a*p.a+p.b*p.b);
 65         double sinC = p.a/sqrt(p.a*p.a+p.b*p.b);
 66         temp.x = P[i].x*cosC-P[i].y*sinC;
 67         temp.y = P[i].x*sinC+P[i].y*cosC;
 68         temp.z = P[i].z;
 69         P[i] = temp;
 70         //绕x轴旋转 
 71         double cosA = p.c/sqrt(p.a*p.a+p.b*p.b+p.c*p.c);
 72         double sinA = sqrt(p.a*p.a+p.b*p.b)/sqrt(p.a*p.a+p.b*p.b+p.c*p.c);
 73         temp.x = P[i].x;
 74         temp.y = P[i].y*cosA-P[i].z*sinA;
 75         temp.z = P[i].y*sinA+P[i].z*cosA;
 76         P[i] = temp; 
 77     }
 78 }
 79 
 80 //计算二维凸包面积 
 81 double crossproduct( point a, point b, point c )
 82 {
 83     return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
 84 }
 85  
 86 bool cmp1( point a, point b )
 87 {
 88     if ( a.x == b.x )
 89         return a.y < b.y;
 90     return a.x < b.x;
 91 }
 92 
 93 bool cmp2( point a, point b )
 94 {
 95     return crossproduct( P[0], a, b )>0;
 96 }
 97 
 98 void Graham( int n )
 99 {
100     sort( P+0, P+n, cmp1 );
101     sort( P+1, P+n, cmp2 );
102     
103     int top = -1;
104     if ( n > 0 ) S[++ top] = P[0];
105     if ( n > 1 ) S[++ top] = P[1];
106     if ( n > 2 ) {
107         for ( int i = 2 ; i < n ; ++ i ) {
108             while ( crossproduct( S[top-1], S[top], P[i] ) < 0 ) -- top;
109             S[++ top] = P[i];
110         }
111     }
112     
113     double area = 0.0;
114     for ( int i = 2 ; i <= top ; ++ i )
115         area += crossproduct( S[0], S[i-1], S[i] );
116     
117     printf("%.2lf
",area*0.5);
118 }
119 
120 int main()
121 {
122     int n;
123     while ( cin >> Pl.a >> Pl.b >> Pl.c >> Pl.d ) {
124         if ( Pl.a == 0 && Pl.b == 0 && Pl.c == 0 ) break;
125         cin >> n;
126         for ( int i = 0 ; i <= n ; ++ i )
127             cin >> P[i].x >> P[i].y >> P[i].z;
128         
129         int s_count = shadow( Pl, P[n], n );
130         if ( !s_count )
131             cout << "0.00" << endl;
132         else if ( s_count > n ) 
133             cout << "Infi" << endl;
134         else {
135             change( Pl, s_count );
136             Graham( s_count );
137         }
138     }
139 }