几何计算与图形处理-优快云博客

“`
// sgn返回x经过eps处理的符号，负数返回-1，正数返回1，x的绝对值如果足够小，就返回0。
const double eps = 1e-8;
int sgn(double x) { return x < -eps ? -1 : x > eps ? 1 : 0; }

```
double mysqrt(double x) { return max(0.0, sqrt(x)); }

// Pt是Point的缩写
//int版
struct Pt {
    int x, y;
    Pt() { }
    Pt(int x,int y) : x(x), y(y) { }
};
Pt operator - (Pt a, Pt b) { return Pt(a.x - b.x, a.y - b.y); }
Pt operator + (Pt a, Pt b) { return Pt(a.x + b.x, a.y + b.y); }
Pt operator * (int A, Pt p) { return Pt(p.x*A, p.y*A); }
Pt operator * (Pt p, int A) { return Pt(p.x*A, p.y*A); }
int sqr(int x){return x*x;}
int norm(Pt p) { return p.x*p.x + p.y*p.y; }
double dist (Pt a, Pt b) { return sqrt(norm(a-b)*1.0); }
double len(Pt p) { return sqrt(sqr(p.x)+sqr(p.y)); }
void print(Pt p) { printf("(%f, %f)", p.x, p.y); }
int dot(Pt a, Pt b) { return a.x * b.x + a.y * b.y; }
int det(Pt a, Pt b) { return a.x * b.y - a.y * b.x; }
int det(Pt a,Pt b,Pt o) { return (a.x-o.x)*(b.y-o.y) - (b.x-o.x)*(a.y-o.y);}
double triArea(Pt a, Pt b, Pt c) { return abs(det(a,b,c))/2.0; }

//double版
struct Pt {
    double x, y;
    Pt() { }
    Pt(double x, double y) : x(x), y(y) { }
};
double sqr(double x){return x*x;}
Pt operator - (Pt a, Pt b) { return Pt(a.x - b.x, a.y - b.y); }
Pt operator + (Pt a, Pt b) { return Pt(a.x + b.x, a.y + b.y); }
Pt operator * (double A, Pt p) { return Pt(p.x*A, p.y*A); }
Pt operator * (Pt p, double A) { return Pt(p.x*A, p.y*A); }
double norm(Pt p) { return sqrt(p.x*p.x + p.y*p.y); }
double dist (Pt a, Pt b) { return norm(a-b); }
void print(Pt p) { printf("(%f, %f)", p.x, p.y); }
double len(Pt p) { return sqrt(sqr(p.x)+sqr(p.y)); }
double dot(Pt a, Pt b) { return a.x * b.x + a.y * b.y; }
double det(Pt a, Pt b) { return a.x * b.y - a.y * b.x; }
double det(Pt a,Pt b,Pt o) { return (a.x-o.x)*(b.y-o.y) - (b.x-o.x)*(a.y-o.y);}
double triArea(Pt a, Pt b, Pt c) { return fabs(det(a,b,c))/2.0; }
Pt rotate(Pt p, double a) {
    return Pt(p.x*cos(a) - p.y*sin(a), p.x*sin(a) + p.y*cos(a));
}

struct Sg {
    Pt s, t;
    Sg() { }
    Sg(Pt s, Pt t) : s(s), t(t) { }
    Sg(double a, double b, double c, double d) : s(a, b), t(c, d) { }
};
bool PtOnSegment(Pt s, Pt t, Pt a) {
    return !det(a-s, a-t) && min(s.x, t.x) <= a.x && a.x <= max(s.x, t.x) &&
        min(s.y, t.y) <= a.y && a.y <= max(s.y, t.y);
}
bool PtOnSegment(Pt p, Pt a, Pt b) {
    return !sgn(det(p-a, b-a)) && sgn(dot(p-a, p-b)) <= 0;
}
bool PtOnLine(Pt p, Pt s, Pt t) {
    return !sgn(det(p-a, b-a));
}
bool SgOnLine(Sg seg, Pt ls, Pt lt) {
    return sgn(det(ls,lt,seg.s)*det(ls,lt,seg.t))>0;
}
Pt PtLineProj(Pt s, Pt t, Pt p) {
    double r = dot(p-s, t-s) / (t - s).norm();
    return s + (t - s) * r;
}
bool parallel(Pt a, Pt b, Pt s, Pt t) {
    return !sgn(det(a-b, s-t));
}
double PtSegmentDist(Pt a, Pt b, Pt p) {
    if (sgn(dot(p-a, b-a)) <= 0) return (p-a).norm();
    if (sgn(dot(p-b, a-b)) <= 0) return (p-b).norm();
    return fabs(det(a-p, b-p)) / (a-b).norm();
}
bool SegCro(Sg u,Sg v)   /** 判断线段是否相交 */  
{  
    return((max(u.s.x,u.t.x)>=min(v.s.x,v.t.x))&&  
        (max(v.s.x,v.t.x)>=min(u.s.x,u.t.x))&&  
        (max(u.s.y,u.t.y)>=min(v.s.y,v.t.y))&&  
        (max(v.s.y,v.t.y)>=min(u.s.y,u.t.y))&&  
        (sgn(det(v.s,u.t,u.s)*det(u.t,v.t,u.s))>=0)&&  
        (sgn(det(u.s,v.t,v.s)*det(v.t,u.t,v.s))>=0));  
}

//由三个点求平面方程
//平面方程A*x+B*y+C*z+D=0;
int A = ((a[2].y-a[1].y)*(a[3].z-a[1].z)-(a[2].z-a[1].z)*(a[3].y-a[1].y));
int B = ((a[2].z-a[1].z)*(a[3].x-a[1].x)-(a[2].x-a[1].x)*(a[3].z-a[1].z));
int C = ((a[2].x-a[1].x)*(a[3].y-a[1].y)-(a[2].y-a[1].y)*(a[3].x-a[1].x));
int D = -(A * a[1].x + B * a[1].y + C * a[1].z);

//重心
Pt triangleMassCenter(Pt a, Pt b, Pt c) {
    return (a+b+c) / 3.0;
}
//外心
Pt circumCenter(Pt a, Pt b, Pt c) {
    double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1) / 2.0;
    double a2 = b.x - a.x, b2 = b.y - a.y, c2 = (a1*a1 + b1*b1) / 2.0;
    double d = a1*b2 -a2*b1;
    return a + Pt(c1*b2 - c2*b1, a1*c2 - a2*c1) / d;
}
//垂心
Pt Orthocenter(Pt a, Pt b, Pt c) {
    return triangleMassCenter(a, b, c) * 3.0 - circumCenter(a, b, c) * 2.0;
}

#define nxt(x) ((x+1)%n)
struct Polygon{
    vector<Pt> poi;
    Polygon(){}
    Polygon(vector<Pt> _poi):poi(_poi){}
    int size(){
        return poi.size();  
    }   
};
int PtInPolygon(Pt p, Polygon &a) {
    int num = 0, d1, d2, k, n = a.size();
    for (int i = 0; i < n; ++i) {
        if (PtOnSegment(p, a.poi[i], a.poi[nxt(i)])) {
            return 2;
        }
        k = det(a.poi[nxt(i)]-a.poi[i], p-a.poi[i]);//double时要加sgn
        d1 = a.poi[i].y - p.y;
        d2 = a.poi[nxt(i)].y - p.y;
        if (k > 0 && d1 <= 0 && d2 > 0) num++;
        if (k < 0 && d2 <= 0 && d1 > 0) num--;
    }
    return num != 0;
}

typedef vector<Pt> Convex;
typedef vector<Pt> Polygon;
// 排序比较函数，水平序
bool comp_less(Pt a, Pt b) {
    return a.y<b.y || a.y==b.y && a.x<b.x;
}
// 返回a中点计算出的凸包，结果存在res中
void convex_hull(Convex &res, vector<Pt> &a) {
    res.resize(2 * a.size() + 5);
    sort(a.begin(), a.end(), comp_less);
    a.erase(unique(a.begin(), a.end()), a.end());
    int m = 0;
    for (int i = 0; i < int(a.size()); ++i) {
        while (m>1 && sgn(det(res[m-1] - res[m-2], a[i] - res[m-2])) <= 0)
            --m;
        res[m++] = a[i];
    }
    int k = m;
    for (int i = int(a.size()) - 2; i >= 0; --i) {
        while (m>k && sgn(det(res[m-1] - res[m-2], a[i] - res[m-2])) <= 0)
            --m;
        res[m++] = a[i];
    }
    res.resize(m);
    if (a.size() > 1) res.resize(m-1);
}

//求多边形面积
double S(Polygon &p)  
{  
    double res = 0;  
    p.push_back(p[0]);  
    for(int i=1;i<p.size()-1;i++)  
       res += det(p[0],p[i],p[i+1]);  
    if(res < 0) res = -res;  
    return res / 2.0;  
}

//求凸包的最远距离
double convex_diameter(const Convex &a, int &first, int &second) {
    int n = a.size();
    double ans = 0.0;
    first = second = 0;
    if (n == 1) return ans;
    for (int i = 0, j = 1; i < n; ++i) {
        while (sgn(det(a[nxt(i)]-a[i], a[j]-a[i]) - det(a[nxt(i)]-a[i], a[nxt(j)]-a[i])) < 0)
            j = nxt(j);
        double d = max((a[i]-a[j]).norm(), (a[nxt(i)]-a[nxt(j)]).norm());
        if (d > ans) ans=d, first=i, second=j;
    }
    return ans;
}

计算几何模板