UVA12305 Polishing a Extruded Polygon 多面体切割 [Rujia Liu's Presents, Dedicated to Geometry and CG Lovers]

不知道这算不算个神题了，AC的时候只有我和汝哥两人过这题……

想了一天，调了一天，竟然因为三角剖分时就差了一句叉积判断相邻边的凹凸性WA了好几天……

1、三角剖分：

所给多面体是非凸的，难以处理，剖分成一个个三棱柱就都是凸多面体了。最后才开始写陌生的三角剖分的，已经写+调了二百多行疲惫不堪的时候看到一计算几何书上好复杂的nlogn算法，竟然还要再搞平衡树，简直要崩溃。百度、谷歌都不给力了，期刊论文看着也晕乎，只好凭着自己的理解来暴力剖分了。

读入每个点建立双向链表（单向应该也没关系，灵活性差点），然后开始循环剖分，对每个相邻边，先判断凹凸性，然后枚举其他所有边判是否和如图所示虚线交叉，通过判断后可以确定这是个凸出的角，可以切掉，就得到了一个三角形。不停地对剩下的多边形进行这样的操作，每一时刻一定会有这样的凸出的角可以切，最后可以完成剖分。

2、模拟切割：

剖分的同时就建立了一个个三棱柱，对“外侧”的面打上标记，这样才能计算正确的表面积。我用链表来存多面体的各个面，用构成多面体的各个面来表示这个多面体。切割的时候就是切这些面，思考起来挺麻烦的，实现的时候思路还是比较清晰，切割凸多面体得到的任何面一定还是凸多边形，而切凸多边形只需要把切面一侧的点去掉，与切面相交的两点连起来，和剩下的点围成新的凸包。对于切掉的整个面，乃至整个凸多面体，就可以有效地利用链表轻松删除。

3、计算体积和面积：

面积就用多边形面积公式，和二维坐标计算的原理是一样的。由于是凸多面体，那么凸多面体内任选一点，就可以把多面体分割成一个个棱锥，棱锥的体积就是点到面的距离乘以面积除以3。这个点可以随意取，取多面体的一个顶点就非常方便了。

调试虽然花了不少时间，但是工作量并不算大，基本上调出样例就差不多了。

应该是没有三点共线的数据，我把处理这个的代码删了是可以AC的。

好长的代码……

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<math.h>
#include<vector>
#include<list>
#include<algorithm>
#include<iostream>
using namespace std;
const int maxn = 3011;
const double eps = 1e-8;
inline int dcmp(double x)
{
    return x > eps ? 1 : (x < -eps ? -1 : 0);
}
inline double pz(double x)
{
    return dcmp(x) ? x : 0;
}
/******************************************************************************/
//定义三维点
struct Point3
{
    double x, y, z;
    Point3()
    {
        x = y = z = 0;
    }
    Point3(double a, double b, double c)
    {
        x = a, y = b, z = c;
    }
    Point3 cross(Point3 p)
    {
        return Point3(y * p.z - p.y * z,
                      z * p.x - x * p.z,
                      x * p.y - y * p.x);
    }
    double dot(Point3 p)
    {
        return x * p.x + y * p.y + z * p.z;
    }
    Point3 operator-(const Point3 &p)const
    {
        return Point3(x - p.x, y - p.y, z - p.z);
    }
    Point3 operator-()const
    {
        return Point3(-x, -y, -z);
    }
    Point3 operator+(const Point3 &p)const
    {
        return Point3(x + p.x, y + p.y, z + p.z);
    }
    Point3 operator*(const double b)const
    {
        return Point3(x * b, y * b, z * b);
    }
    bool operator==(const Point3 &p)const
    {
        return !dcmp(x - p.x) && !dcmp(y - p.y) && !dcmp(z - p.z);
    }
    bool operator!=(const Point3 &p)const
    {
        return !(*this == p);
    }
    bool operator<(const Point3 &p)const
    {
        if(!dcmp(x - p.x))
        {
            if(!dcmp(y - p.y))
                return dcmp(z - p.z) < 0;
            return dcmp(y - p.y) < 0;
        }
        return dcmp(x - p.x) < 0;
    }
    bool operator>(const Point3 &p)const
    {
        return p < *this;
    }
    bool operator>=(const Point3 &p)const
    {
        return !(*this < p);
    }
    bool operator<=(const Point3 &p)const
    {
        return !(*this > p);
    }
    Point3 fxl(Point3 b, Point3 c)//平面法向量
    {
        return (*this - b).cross(b - c);
    }
    double Dis(Point3 b)
    {
        return sqrt((*this - b).dot(*this - b));
    }
    double vlen()
    {
        return sqrt(dot(*this));
    }
    bool PinLine(Point3 b, Point3 c)//三点共线
    {
        return fxl(b, c).vlen() < eps;
    }
    bool PonPlane(Point3 b, Point3 c, Point3 d)//四点共面
    {
        return !dcmp(fxl(b, c).dot(d - *this));
    }
    bool PonSeg(Point3 b, Point3 c)//点在线段上，包括端点
    {
        return !dcmp((*this - b).cross(*this - c).vlen()) &&
               (*this - b).dot(*this - c) <= 0;
    }
    bool PinSeg(Point3 b, Point3 c)//点在线段上，不包括端点
    {
        return PonSeg(b, c) && *this != b && *this != c;
    }
    double PtoLine(Point3 b, Point3 c)//点到直线距离
    {
        return (*this - b).cross(c - b).vlen() / b.Dis(c);
    }
    double PtoPlane(Point3 b, Point3 c, Point3 d)//点到平面距离
    {
        return fabs(b.fxl(c, d).dot(*this - b)) / b.fxl(c, d).vlen();
    }
};
/******************************************************************************/
//定义平面+空间平面凸包
struct Plane
{
    double a, b, c, d;
    bool outplane;//计入表面积的面
    vector<Point3> p;
    Plane()
    {
        a = b = c = d = 0, outplane = false;
        p.clear();
    }
    inline void init(double a_, double b_, double c_, double d_)
    {
        a = a_, b = b_, c = c_, d = d_;
        p.clear();
    }
    inline void init(Point3 pa, Point3 pb, Point3 pc)
    {
        Point3 t = (pa - pb).cross(pa - pc);
        a = t.x, b = t.y, c = t.z;
        d = -(pa.x * t.x + pa.y * t.y + pa.z * t.z);
        p.clear();
    }
    Plane(double a_, double b_, double c_, double d_)
    {
        init(a_, b_, c_, d_);
    }
    Plane(Point3 pa, Point3 pb, Point3 pc)
    {
        init(pa, pb, pc);
    }
    double PtoPlane(const Point3 &pa)const//点面距离
    {
        return fabs(Sub(pa)) / sqrt(a * a + b * b + c * c);
    }
    double Sub(const Point3 &p)const//点代入方程
    {
        return p.x * a + p.y * b + p.z * c + d;
    }
    Point3 PcrossPlane(Point3 a, Point3 b)
    {
        return a + (b - a) * (PtoPlane(a) / (PtoPlane(a) + PtoPlane(b)));
    }
    int Parallel(Plane pl)//面平行
    {
        if(!dcmp(a * pl.b - b * pl.a) && !dcmp(a * pl.c - c * pl.a) && !dcmp(b * pl.c - c * pl.b))
            return dcmp(pl.Sub(p[0])) > 0 ? 1 : -1;
        return 0;
    }
    bool Cut(Plane &pl)//平面切割
    {
        switch(Parallel(pl))
        {
        case -1:
            return true;
        case 1:
            return false;
        }
        int i, j, k, n = p.size();
        bool flag1, flag2;
        Point3 p1, p2;
        vector<Point3> ret;
        for(i = flag1 = flag2 = 0; i < n; ++ i)
        {
            flag1 |= pl.Sub(p[i]) < 0;
            flag2 |= pl.Sub(p[i]) > 0;
        }
        if(flag1 != flag2) return flag1;
        if(!flag1) return true;
        for(i = 0; pl.Sub(p[i]) >= 0; ++ i);
        for(; pl.Sub(p[i]) < 0; i = (i + 1) % n);
        for(j = i; pl.Sub(p[j]) >= 0; j = (j + 1) % n);
        for(k = j; k != i; k = (k + 1) % n)
            ret.push_back(p[k]);
        p1 = pl.PcrossPlane(p[i], p[(i + n - 1) % n]);
        p2 = pl.PcrossPlane(p[j], p[(j + n - 1) % n]);
        ret.push_back(p1), ret.push_back(p2);
        p = ret;
        pl.p.push_back(p1), pl.p.push_back(p2);
        return true;
    }
    void Graham()//求空间平面凸包
    {
        int len, i, n, top = 1;
        vector<Point3> res;
        Point3 fx(a, b, c);
        sort(p.begin(), p.end());
        vector<Point3>::iterator iter = unique(p.begin(), p.end());
        p.erase(iter, p.end());
        n = p.size();
        res.push_back(p[0]), res.push_back(p[1]);
        for(i = 2; i < n; ++ i)
        {
            while(top && dcmp((res[top] - res[top - 1]).cross(p[i] - res[top]).dot(fx)) <= 0)
                res.pop_back(), -- top;
            res.push_back(p[i]), ++ top;
        }
        len = top;
        res.push_back(p[n - 2]), ++ top;
        for(i = n - 3; i >= 0; -- i)
        {
            while(top != len && dcmp((res[top] - res[top - 1]).cross(p[i] - res[top]).dot(fx)) <= 0)
                res.pop_back(), -- top;
            res.push_back(p[i]), ++ top;
        }
        res.pop_back();
        p = res;
    }
    double PolygonArea()//空间平面凸包面积
    {
        int n = p.size();
        double ret = 0;
        for(int i = 2; i < n; ++ i)
            ret += (p[i - 1] - p[0]).cross(p[i] - p[0]).vlen();
        return ret * 0.5;
    }
};
/******************************************************************************/
//定义变量
struct Polyhedron//立方体面集合
{
    list<Plane> pl;
};
list<Polyhedron> pol;
int n, h, m;//点个数，高度，切割次数
struct PointChain
{
    Point3 p;
    int pre, nex;
    bool outside;//标记发向nex的边为外部边，抬起的面计入总面积
} PC[maxn];
/******************************************************************************/
//判相交，辅助判断三角剖分合法性
inline double det(double x1, double y1, double x2, double y2)
{
    return x1 * y2 - x2 * y1;
}
inline double cross2(Point3 a, Point3 b, Point3 c)//平面叉积
{
    return det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}
inline bool SegCross(Point3 u1, Point3 u2, Point3 v1, Point3 v2)
{
    return ((dcmp(cross2(u1, u2, v1)) ^ dcmp(cross2(u1, u2, v2))) == -2 &&
            (dcmp(cross2(v1, v2, u1)) ^ dcmp(cross2(v1, v2, u2))) == -2) ||
           ((v1.PinSeg(u1, u2) || v2.PinSeg(u1, u2)) &&
            (u1.PonSeg(v1, v2) || u2.PonSeg(v1, v2))) ||
           ((v1.PonSeg(u1, u2) || v2.PonSeg(u1, u2)) &&
            (u1.PinSeg(v1, v2) || u2.PinSeg(v1, v2)));
}
/******************************************************************************/
bool CanDo(int s, int e)
{
    int tp = PC[s].nex;
    while(tp != s)
    {
        if(SegCross(PC[s].p, PC[e].p, PC[tp].p, PC[PC[tp].pre].p))
            return false;
        tp = PC[tp].nex;
    }
    return true;
}
void AddPolyhedron(int A, int B, int C)
{
    Polyhedron t;
    Plane pl;
    Point3 a = PC[A].p, b = PC[B].p, c = PC[C].p;
    Point3 a_ = Point3(a.x, a.y, h), b_ = Point3(b.x, b.y, h), c_ = Point3(c.x, c.y, h);
    pl.init(a, b, c);
    pl.p.push_back(a);
    pl.p.push_back(b);
    pl.p.push_back(c);
    pl.outplane = true;
    t.pl.push_front(pl);
    pl.init(a_, b_, c_);
    pl.p.push_back(a_);
    pl.p.push_back(b_);
    pl.p.push_back(c_);
    t.pl.push_front(pl);
    pl.init(a, b, b_);
    pl.p.push_back(a);
    pl.p.push_back(b);
    pl.p.push_back(b_);
    pl.p.push_back(a_);
    pl.outplane = PC[A].outside, PC[A].outside = false;
    t.pl.push_front(pl);
    pl.init(b, c, c_);
    pl.p.push_back(b);
    pl.p.push_back(c);
    pl.p.push_back(c_);
    pl.p.push_back(b_);
    pl.outplane = PC[B].outside, PC[B].outside = false;
    t.pl.push_front(pl);
    pl.init(c, a, a_);
    pl.p.push_back(c);
    pl.p.push_back(a);
    pl.p.push_back(a_);
    pl.p.push_back(c_);
    if(PC[C].nex == A && PC[C].outside)
        pl.outplane = true, PC[C].outside = false;
    else pl.outplane = false;
    t.pl.push_front(pl);
    pol.push_front(t);
}
void Triangulation()//循环枚举相邻三点划分三角形
{
    int tp = 0;
    while(PC[PC[tp].nex].nex != tp)
    {
        while(dcmp(cross2(PC[tp].p, PC[PC[tp].nex].p, PC[PC[PC[tp].nex].nex].p)) <= 0 || 
            !CanDo(tp, PC[PC[tp].nex].nex)) tp = PC[tp].nex;
        AddPolyhedron(tp, PC[tp].nex, PC[PC[tp].nex].nex);
        PC[PC[PC[tp].nex].nex].pre = tp;
        PC[tp].nex = PC[PC[tp].nex].nex;
    }
}
void init()
{
    double x, y;
    int i, tp;
    pol.clear();
    for(i = 0; i < n; ++ i)
    {
        scanf("%lf%lf", &x, &y);
        PC[i].p = Point3(x, y, 0);
        PC[i].pre = i - 1;
        PC[i].nex = i + 1;
        PC[i].outside = true;
    }
    PC[0].pre = n - 1, PC[n - 1].nex = 0;
    Triangulation();
}
/******************************************************************************/
list<Polyhedron>::iterator iti;
list<Plane>::iterator itj;
double CalSumVol()//求总体积
{
    double ans = 0;
    Point3 tmp;
    for(iti = pol.begin(); iti != pol.end(); ++ iti)
    {
        for(itj = (*iti).pl.begin(), tmp = (*itj).p[0]; itj != (*iti).pl.end(); ++ itj)
            ans += (*itj).PolygonArea() * (*itj).PtoPlane(tmp);
    }
    return ans / 3;
}
double CalSumArea()//求总面积
{
    double ans = 0;
    for(iti = pol.begin(); iti != pol.end(); ++ iti)
        for(itj = (*iti).pl.begin(); itj != (*iti).pl.end(); ++ itj)
            if((*itj).outplane) ans += (*itj).PolygonArea();
    return ans;
}

void MakeAns()
{
    double a, b, c, d;
    bool flag;
    while(m --)
    {
        scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
        for(iti = pol.begin(); iti != pol.end();)
        {
            Plane tmp(a, b, c, d);
            tmp.outplane = true;
            for(itj = (*iti).pl.begin(), flag = false; itj != (*iti).pl.end();)
            {
                if(!(*itj).Cut(tmp))
                {
                    list<Plane>::iterator tmpj = itj;
                    ++ itj;
                    (*iti).pl.erase(tmpj);
                }
                else flag = true, ++ itj;
            }
            if(!flag)
            {
                list<Polyhedron>::iterator tmpi = iti;
                ++ iti;
                pol.erase(tmpi);
            }
            else if(tmp.p.size())
            {
                tmp.Graham();
                (*iti).pl.push_front(tmp);
                ++ iti;
            }
            else ++ iti;
        }
        printf("%.3f %.3f\n", CalSumVol(), CalSumArea());
    }
}
int main()
{
    while(scanf("%d%d%d", &n, &h, &m), n | h | m)
    {
        init();
        printf("%.3f %.3f\n", CalSumVol(), CalSumArea());
        MakeAns();
    }
    return 0;
}