This file is implementation of Common Common Computational Geometry Algorithms.Please please pay attention to input according to the specified data type.

個人實現的一些計算幾何中常見的演算法，包括點，線，多邊形等；所有演算法只依賴於C++標準庫，不用包含任何其他第三方庫，包含此標頭檔案即可使用。使用時請注意按照規定的資料型別進行輸入，目前只使用C++來實現演算法，具體演算法原理會陸續在Github上更新。

目前實現的演算法包括點、向量、線段、直線、三角形、多邊形、圓等基本計算幾何模型。

// Copyright (C) Common Computational Geometry Algorithms e.U, ZutterHao
//
// This file is implementation of Common Common Computational Geometry Algorithms.
//
// Please please pay attention to input according to the specified data type.
//
// Author: ZutterHao .Nanjing University ,VISG
// Github: https://github.com/fanghao6666

// 個人實現的一些計算幾何中常見的演算法，包括點，線，多邊形等
// 所有演算法只依賴於C++標準庫，不用包含任何其他第三方庫
// 使用時請注意按照規定的資料型別進行輸入
// 目前只使用C++來實現演算法，具體演算法原理會陸續在Github上更新
// Github: https://github.com/fanghao6666

/*** 包含的標頭檔案 ***/
#include <iostream>
#include <vector>
#include <map>
#include <stack>
#include <algorithm>

using namespace std;

/*** 常用常量 ***/
const double PI = 3.14159265;

/******************************* 基本幾何資料型別 *******************************/
// 點,二維和三維,同時點也可以表示一個向量
struct Point
{
    double x;    // x座標
    double y;    // y座標
    double z;    // z座標（預設為0，如果需要三維點則給z賦值）

    Point(double a = 0, double b = 0, double c = 0) { x = a; y = b; z = c; } // 建構函式
};
// 點的加法
Point add(const Point& lhs, const Point& rhs)
{
    Point res;

    res.x = lhs.x + rhs.x;
    res.y = lhs.y + rhs.y;
    res.z = lhs.z + rhs.z;

    return res;
}
// 點的減法
Point sub(const Point& lhs, const Point& rhs)
{
    Point res;

    res.x = lhs.x - rhs.x;
    res.y = lhs.y - rhs.y;
    res.z = lhs.z - rhs.z;

    return res;
}
// 向量的乘法
Point mul(const Point& p, double ratio)
{
    Point res;

    res.x = p.x * ratio;
    res.y = p.y * ratio;
    res.z = p.z * ratio;

    return res;
}
// 向量的除法
Point div(const Point& p, double ratio)
{
    Point res;
    
    res.x = p.x / ratio;
    res.y = p.y / ratio;
    res.z = p.z / ratio;

    return res;
}
// 點判斷相等
bool equal(const Point& lhs, const Point& rhs)
{
    return(lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z);
}
// 線,包括線段和直線
struct Line
{
    Point s;    // 起點
    Point e;    // 終點
    bool is_seg; // 是否是線段

    Line() {};    // 預設建構函式
    Line(Point a, Point b, bool _is_seg = true) { s = a; e = b; is_seg = _is_seg; }    // 建構函式(預設是線段)
};
// 三角形平面
struct Triangle
{
    Point v0;
    Point v1;
    Point v2;
    bool is_plane;

    Triangle() {}; // 預設建構函式
    Triangle(Point a, Point b, Point c, bool _is_plane = false) { v0 = a; v1 = b; v2 = c; is_plane = _is_plane; }// 建構函式（預設是三角形）
};

/******************************* 計算幾何演算法目錄 *******************************/

// 一、點
// 1.1、兩點之間的距離
double distance(const Point& p1, const Point& p2);

// 1.2、向量長度
double length(const Point& vec);

// 1.3、向量標準化
Point normalize(const Point& vec);

// 1.4、向量點乘
double dotMultiply(const Point& op, const Point& p1, const Point& p2);
double dotMultiply(const Point& vec1, const Point& vec2);

// 1.5、向量叉乘
Point multiply(const Point& op, const Point& p1, const Point& p2);
Point multiply(const Point& vec1, const Point& vec2);

// 1.6、點到線的距離
double ptolDistance(const Point& p, const Line& l);

// 1.7、點到線的投影點
Point ptolProjection(const Point& p, const Line& l);

// 1.8、點關於線的對稱點
Point ptolSymmetry(const Point& p, const Line& l);

// 1.9、點是否在線上
bool isponl(const Point& p, const Line& l);

// 1.10、向量夾角正弦
double Sin(const Point& op, const Point& p1, const Point& p2);
double Sin(const Point& vec1, const Point& vec2);

// 1.11、向量夾角餘弦
double Cos(const Point& op, const Point& p1, const Point& p2);
double Cos(const Point& vec1, const Point& vec2);

// 1.12、向量夾角正切
double Tan(const Point& op, const Point& p1, const Point& p2);
double Tan(const Point& vec1, const Point& vec2);

// 1.13、向量夾角角度
double Angle(const Point& op, const Point& p1, const Point& p2,bool is_radian = true);
double Angle(const Point& vec1, const Point& vec,bool is_radian = true);

// 1.14、判斷三點是否共線
bool isPointsCollinear(const Point& p1, const Point& p2, const Point& p3);

// 二、線
// 2.1、線段是否相交
bool isSegIntersect(const Line& l1, const Line& l2,Point& inter_p);

// 2.2、求直線的夾角
double angleOfLines(const Line& l1, const Line& l2,bool is_radian = true);

// 2.3、一階貝塞爾曲線插值
vector<Point> firstOrderBezier(const Point& s, const Point& e, int inter_num);

// 2.4、二階貝塞爾曲線插值
vector<Point> secondOrderBezier(const Point& s, const Point& e, const Point& p,int inter_num);

// 2.5、三階貝塞爾曲線插值
vector<Point> thirdOrderBezier(const Point& s, const Point& e, const Point& p1,const Point& p2,int inter_num);

// 三、三角形
// 3.1、三角形三個點是否能夠構成三角形
bool isTriangle(const Triangle& t);

// 3.2、點是否在三角形內部
bool isPointInTriangle(const Triangle& t, const Point& p,double& u,double& v);

// 3.3、點到平面的投影點（距離最近的點）
Point ptotProjection(const Triangle& t, const Point& p);

// 3.4、點到平面的距離
double ptotDistance(const Triangle& t, const Point& p);

// 3.5、線段和平面的交點
Point ltotInterPoint(const Triangle& t, const Line& l);

// 3.6、計算平面的單位法向量
Point getUnitNormal(const Triangle& t);

// 3.7、計算三角形的面積
double areaOfTriangle(const Triangle& t);


// 四、多邊形
// 4.1、判斷多邊形頂點的凹凸性
void checkConvex(const vector<Point>& polygon, vector<bool>& flags);

// 4.2、判斷多邊形是否為凸多邊形
bool isConvex(const vector<Point>& polygon);

// 4.3、求多邊形圍成的面積
double areaOfPolygon(const vector<Point>& polygon);

// 4.4、判斷多邊形是否按照逆時針排列
bool isConterClock(const vector<Point>& polygon);

// 4.5、判斷點是否在多邊形內部
bool isPointInPolygon(const vector<Point>& polygon, const Point& p);

// 4.6、判斷線段是否在多邊形內部
bool isSegInPolygon(const vector<Point>& polygon, const Line& l);

// 4.7、判斷圓是否在多邊形內部
bool isCircleInPolygon(const vector<Point>& polygon, const Point& p, double radius);

// 4.8、尋找點集凸包演算法（graham演算法）
vector<Point> findConvexGraham(const vector<Point>& points);

// 4.9、尋找點集凸包演算法（上下凸包法）時間複雜度O(nlogn)
vector<Point> findConvex(const vector<Point>& points);

// 4.10、求簡單多邊形重心
Point centerOfPolygon(const vector<Point>& polygon);

// 4.11、求肯定在多邊形內部的一個點
Point pointInPolygon(const vector<Point>& polygon);

// 4.12、計算多邊形的範圍
void boxOfPolygon(const vector<Point>& polygon, Point& down_left, Point& up_right);

// 五、圓
// 5.1、點和圓的關係
int pointToCircle(const Point& c, double radius, const Point& p);

// 5.2、直線和圓的關係
int lineToCircle(const Point& c, double radius, const Line& l);

// 5.3、線段和圓的關係
int segToCircle(const Point& c, double radius, const Line& l);

// 5.4、兩圓之間的關係
int circleToCircle(const Point& c1, double raduis1, const Point& c2, double radius2);



/******************************* 計算幾何演算法實現 *******************************/
//一、點

// 1.1、兩點之間的距離
//
// 引數：p1 : 第一個點 p2: 第二個點
//  
double distance(const Point& p1, const Point& p2)
{
    return(sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2)));
}

// 1.2、向量的長度
//
// 引數： vec 向量
//
double length(const Point& vec)
{
    return (sqrt(pow(vec.x, 2) + pow(vec.y, 2) + pow(vec.z, 2)));
}

// 1.3、向量標準化（向量的長度規約到1）
//
// 引數： vec ： 向量
//
Point normalize(const Point& vec)
{
    Point res;

    res = div(vec, length(vec));

    return res;
}

// 1.4、向量點乘
//
// 引數：(p1-op)為向量1，（p2-op）為向量2
//
double dotMultiply(const Point& op, const Point& p1, const Point& p2)
{
    return ((p1.x - op.x) * (p2.x - op.x) + (p1.y - op.y) * (p2.y - op.y) + (p1.z - op.z) * (p2.z - op.z));
}
// 引數：vec1為向量1，vec2為向量2
//
double dotMultiply(const Point& vec1, const Point& vec2)
{
    return(vec1.x * vec2.x + vec1.y * vec2.y + vec1.z * vec2.z);
}

// 1.5、向量叉乘
//
// 引數：(p1-op)為向量1，（p2-op）為向量2
// 
Point multiply(const Point& op, const Point& p1, const Point& p2)
{
    Point result;

    result.x = (p1.y - op.y) * (p2.z - op.z) - (p2.y - op.y) * (p1.z - op.z);
    result.y = (p1.z - op.z) * (p2.x - op.x) - (p2.z - op.z) * (p1.x - op.x);
    result.z = (p1.x - op.x) * (p2.y - op.y) - (p2.x - op.x) * (p1.y - op.y);

    return result;
}
// 引數： vec1為向量1，vec2為向量2
//
Point multiply(const Point& vec1, const Point& vec2)
{
    Point result;

    result.x = vec1.y * vec2.z - vec2.y * vec1.z;
    result.y = vec1.z * vec2.x - vec2.z * vec1.x;
    result.z = vec1.x * vec2.y - vec2.x * vec1.y;

    return result;
}

// 1.6、點到線的距離
// 
// 引數： p : 點  l：直線
//
double ptolDistance(const Point& p, const Line& l)
{
    Point line_vec = sub(l.e,l.s);
    Point point_vec = sub(p, l.s);

    // 首先計算點線上段投影長度
    double project_len = dotMultiply(line_vec, point_vec) / length(line_vec);

    // 根據勾股定理計算點的距離
    double distance = sqrt(pow(length(line_vec), 2) - pow(project_len, 2));

    return distance;
}

// 1.7、點到線的投影點
//
// 引數：p : 點  l : 線
//
Point ptolProjection(const Point& p, const Line& l)
{
    Point line_vec = sub(l.e, l.s);
    Point point_vec = sub(p, l.s);
    Point unit_line_vec = normalize(line_vec);

    // 計算投影長度
    double project_len = dotMultiply(point_vec, unit_line_vec);

    // 投影點
    Point project_p = add(l.s,mul(unit_line_vec, project_len));

    return project_p;
}

// 1.8、點關於線的對稱點
//
// 引數： p : 點  l : 對稱線
//
Point ptolSymmetry(const Point& p, const Line& l)
{
    // 首先求出點在直線上的投影點
    Point project_p = ptolProjection(p, l);

    // 點到投影點的向量
    Point project_vec = sub(project_p, p);

    // 對稱點
    Point symmetry_p = add(p, mul(project_vec, 2));

    return symmetry_p;
}

// 1.9、點是否在線上
// 線分為直線和線段，直線表示的是直線是否經過點
//
// 引數：p : 點  l : 線段或者線
// 
bool isponl(const Point& p, const Line& l)
{
    Point line_vec = sub(l.e, l.s);
    Point point_vec1 = sub(p, l.s);
    Point point_vec2 = sub(p, l.e);

    Point mul_vec = multiply(line_vec, point_vec1);
    double dot = dotMultiply(point_vec1, point_vec2);
    // 點是否線上段上
    if (l.is_seg)
    {
        if (equal(p,l.s) || equal(p,l.e))
            return true;
        return (0.0 == length(mul_vec) && dot < 0.0);
        
    }
    // 點是否在直線上
    else
    {
        return (0.0 == length(mul_vec));
    }
}

// 1.10、向量夾角正弦
//
// 引數： op : 向量公共點 p1 : 向量1端點 p2 : 向量2端點
//
double Sin(const Point& op, const Point& p1, const Point& p2)
{
    Point vec1 = sub(p1, op);
    Point vec2 = sub(p2, op);

    return Sin(vec1, vec2);
}
// 引數： vec1 向量1  vec2 向量2
// 
double Sin(const Point& vec1, const Point& vec2)
{
    return sqrt(1.0 - pow(Cos(vec1, vec2), 2));
}

// 1.11、向量夾角餘弦
//
// 引數： op : 向量公共點 p1 : 向量1端點 p2 : 向量2端點
//
double Cos(const Point& op, const Point& p1, const Point& p2)
{
    Point vec1 = sub(p1, op);
    Point vec2 = sub(p2, op);

    return Cos(vec1, vec2);
}
// 引數： vec1 向量1  vec2 向量2
// 
double Cos(const Point& vec1, const Point& vec2)
{
    Point unit_vec1 = normalize(vec1);
    Point unit_vec2 = normalize(vec2);

    return dotMultiply(unit_vec1, unit_vec2);
}

// 1.12、向量夾角正切
//
// 引數： op : 向量公共點 p1 : 向量1端點 p2 : 向量2端點
//
double Tan(const Point& op, const Point& p1, const Point& p2)
{
    Point vec1 = sub(p1, op);
    Point vec2 = sub(p2, op);

    return Tan(vec1, vec2);
}
// 引數： vec1 向量1  vec2 向量2
// 
double Tan(const Point& vec1, const Point& vec2)
{
    double cos = Cos(vec1, vec2);
    double sin = Sin(vec1, vec2);

    // 分母不為零
    if (0.0 == cos)
        return -1;
    else
        return (sin / cos);
}

// 1.13、計算點的夾角角度
// 引數:  op : 向量公共點 p1 : 向量1端點 p2 : 向量2端點 is_radian : 預設為弧度制
//
double Angle(const Point& op, const Point& p1, const Point& p2, bool is_radian)
{
    double cos_value = Cos(op, p1, p2);

    if (is_radian)
    {
        return acos(cos_value);
    }
    else
    {
        return (acos(cos_value) / PI * 180.0);
    }
}
// 引數： vec1 : 向量1 vec2 : 向量2
// 
double Angle(const Point& vec1, const Point& vec2,bool is_radian)
{
    double cos_value = Cos(vec1, vec2);

    if (is_radian)
    {
        return acos(cos_value);
    }
    else
    {
        return (acos(cos_value) / PI * 180.0);
    }
}

// 1.14、判斷三點是否共線
bool isPointsCollinear(const Point& p1, const Point& p2, const Point& p3)
{
    Line line(p1, p2, false);

    // 判斷第三個點是否在前兩個點的線段上
    return isponl(p3, line);
}

// 二、線

// 2.1、線段是否相交
// 其中如果線段的端點重合或者某個線段端點在另一個線段上也算相交
// 線段判斷是否相交，如果是直線則相當於判斷是否平行
//
// 引數： l1 : 線段1  l2 : 線段2  inter_p : 如果相交返回交點
//
bool isSegIntersect(const Line& l1, const Line& l2,Point& inter_p)
{
    Point line1 = sub(l1.e, l1.s);
    Point line2 = sub(l2.e, l2.s);
    Point norm1 = normalize(line1);
    Point norm2 = normalize(line2);
    // 線段相交
    if (l1.is_seg)
    {
        // 端點線上段上
        if (isponl(l1.s, l2))
        {
            inter_p = l1.s;
            return true;
        }
        if (isponl(l1.e, l2))
        {
            inter_p = l1.e;
            return true;
        }
        if (isponl(l2.s, l1))
        {
            inter_p = l2.s;
            return true;
        }
        if (isponl(l2.e, l1))
        {
            inter_p = l2.e;
            return true;
        }
        // 判斷線段是否相互跨立
        double dot1 = dotMultiply(multiply(sub(l2.s, l1.s), line1), multiply(sub(l2.e, l1.s), line1));
        double dot2 = dotMultiply(multiply(sub(l1.s, l2.s), line2), multiply(sub(l1.e, l2.s), line2));
        if (dot1 < 0.0 && dot2 < 0.0)
        {
            double t1 = length(multiply(sub(l1.s, l2.s), norm2)) / length(multiply(norm2, norm1));
            double t2 = length(multiply(sub(l2.s, l1.s), norm1)) / length(multiply(norm1, norm2));

            inter_p = add(l1.s, mul(norm1, t1));
            return true;
        }
        else
        {
            return false;
        }

    }
    // 直線相交
    else
    {
        if (Cos(line1, line2) == 1.0)
            return false;

        double t1 = length(multiply(sub(l1.s, l2.s), norm2)) / length(multiply(norm2, norm1));
        double t2 = length(multiply(sub(l2.s, l1.s), norm1)) / length(multiply(norm1, norm2));

        inter_p = add(l1.s, mul(norm1, t1));
        return true;
    }
}

// 2.2、求直線的夾角
// 
// 
// 引數： l1 : 線段1 l2 : 線段2
//
double angleOfLines(const Line& l1, const Line& l2,bool is_radian)
{
    Point line1 = sub(l1.e, l1.s);
    Point line2 = sub(l2.e, l2.s);

    return Angle(line1, line2, is_radian);
}

// 2.3、一階貝塞爾曲線插值
// 
// 引數： s： 起點 e : 終點 inter_num：插值點數量（不包括起始點）
// 返回值包括起始點
//
vector<Point> firstOrderBezier(const Point& s, const Point& e, int inter_num)
{
    vector<Point> res;
    res.push_back(s);
    for (int i = 1; i <= inter_num; ++i)
    {
        double a1 = double(i) / double(inter_num + 1);
        double a2 = 1.0 - a1;
        res.push_back(add(mul(s, a2), mul(e, a1)));
    }
    res.push_back(e);

    return res;
}

// 2.4、二階貝塞爾曲線插值
// 
// 引數： s： 起點 e : 終點 p : 控制點 inter_num：插值點數量（不包括起始點）
// 返回值包括起始點
//
vector<Point> secondOrderBezier(const Point& s, const Point& e, const Point& p, int inter_num)
{
    vector<Point> res;
    res.push_back(s);
    for (int i = 1; i <= inter_num; ++i)
    {
        double a = double(i) / double(inter_num + 1);
        double a1 = pow(a,2);
        double a2 = 2 * a * (1.0 - a);
        double a3 = pow(1.0 - a, 2);
        res.push_back(add(add(mul(s, a3), mul(p, a2)), mul(e, a1)));
    }
    res.push_back(e);

    return res;
}

// 2.5、三階貝塞爾曲線插值
// 
// 引數： s： 起點 e : 終點 p1，p2 : 控制點 inter_num：插值點數量（不包括起始點）
// 返回值包括起始點
//
vector<Point> thirdOrderBezier(const Point& s, const Point& e, const Point& p1, const Point& p2, int inter_num)
{
    vector<Point> res;
    res.push_back(s);
    for (int i = 1; i <= inter_num; ++i)
    {
        double a = double(i) / double(inter_num + 1);
        double a1 = pow(a, 3);
        double a2 = 3 * pow(a, 2) * (1.0 - a);
        double a3 = 3 * pow(1.0 - a, 2) * a;
        double a4 = pow(1.0 - a, 3);
        res.push_back(add(add(add(mul(s, a4), mul(p1, a3)), mul(p2, a2)),mul(e,a1)));
    }
    res.push_back(e);

    return res;
}

// 三、三角形

// 3.1、三角形三個點是否能夠構成三角形
// 不共線的三個點組成一個三角形
// 
// 引數： t : 三角形
// 
bool isTriangle(const Triangle& t)
{
    return isPointsCollinear(t.v0, t.v1, t.v2);
}

// 3.2、點是否在三角形內部（重心法）
// 演算法方法連結： https://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html
//
// 引數： t : 三角形 p : 需要判斷的點 u,v分別為用於表示點在兩條邊上投影係數
//
bool isPointInTriangle(const Triangle& t, const Point& p, double& u, double& v)
{
    Point vec1 = sub(t.v1, t.v0);
    Point vec2 = sub(t.v2, t.v0);
    Point vec_p = sub(p, t.v0);

    double dot00 = dotMultiply(vec1, vec1);
    double dot01 = dotMultiply(vec1, vec2);
    double dot02 = dotMultiply(vec1, vec_p);
    double dot11 = dotMultiply(vec2, vec2);
    double dot12 = dotMultiply(vec2, vec_p);

    double inverDeno = 1 / (dot00 * dot11 - dot01 * dot01);

    u = (dot11 * dot02 - dot01 * dot12) * inverDeno;
    v = (dot00 * dot12 - dot01 * dot02) * inverDeno;

    if (u < 0 || u > 1) return false;
    if (v < 0 || v > 1) return false;
    if (u + v < 1)return true;
    else return false;
}

// 3.3、點到平面最近的點，即點到平面的投影
// 
// 引數：t : 三角形 p : 點
// 
Point ptotProjection(const Triangle& t, const Point& p)
{
    Point vec_p = sub(p, t.v0);
    Point unit_normal = getUnitNormal(t);

    double ratio = dotMultiply(vec_p, unit_normal);

    return sub(p, mul(unit_normal, ratio));
}

// 3.4、點到平面的距離
//
// 引數： t : 三角形所在的平面 p : 需要判斷的點
//
double ptotDistance(const Triangle& t, const Point& p)
{
    Point project_p = ptotProjection(t, p);

    return distance(p, project_p);
}

// 3.5、線段和平面的交點
// 
// 引數： t : 三角形所在平面 l : 直線
//
Point ltotInterPoint(const Triangle& t, const Line& l)
{
    Point line_vec = sub(l.e, l.s);
    Point point_vec = sub(t.v0, l.s);
    Point unit_plane_normal = getUnitNormal(t);

    double ratio = dotMultiply(point_vec, unit_plane_normal) / dotMultiply(unit_plane_normal, line_vec);

    return add(l.s, mul(line_vec, ratio));
}

// 3.6、計算平面的單位法向量
// 
// 引數： t : 三角形平面
// 
Point getUnitNormal(const Triangle& t)
{
    Point vec1 = sub(t.v1, t.v0);
    Point vec2 = sub(t.v2, t.v0);

    return normalize(multiply(vec1, vec2));
}

// 3.7、計算三角形的面積
//
// 引數： t : 三角形平面
//
double areaOfTriangle(const Triangle& t)
{
    return (0.5 * length(multiply(sub(t.v1,t.v0),sub(t.v2,t.v0))));
}

// 四、多邊形
// 如果沒有特別說明，則預設輸入多邊形頂點按照逆時針排列，點為二維點即z=0

// 4.1、判斷多邊形頂點的凹凸性
// 
// 引數： polygon : 多邊形點集合 flags : 標誌每個點是否是凸的
//
void checkConvex(const vector<Point>& polygon, vector<bool>& flags)
{
    flags.resize(polygon.size());

    // 找到第一個凸的頂點
    int index = 0;
    for (int i = 1; i < polygon.size(); ++i)
    {
        if (polygon[i].y < polygon[index].y ||
            (polygon[i].y == polygon[index].y && polygon[i].x < polygon[index].x))
        {
            index = i;
        }
    }
    /* 判斷每個點的凹凸性 
    *  通過判斷前後兩個點的向量叉乘判斷是否滿足逆時針
    */
    int size = polygon.size() - 1;
    flags[index] = true;
    while (size)
    {
        if (multiply(polygon[index], polygon[(index + 1) % size], polygon[(index + 2) % size]).z >= 0)
        {
            flags[(index + 1) % size] = true;
        }
        else
        {
            flags[(index + 1) % size] = false;
        }
        index++;
        size--;
    }
}

// 4.2、判斷多邊形是否為凸多邊形
//
// 引數 ： polygon : 輸入的多邊形頂點
//
bool isConvex(const vector<Point>& polygon)
{
    vector<bool> flags;
    // 判斷每個點的凹凸性
    checkConvex(polygon, flags);
    // 如果有一個點不是凸的，則此多邊形為非凸
    for (auto c : flags)
        if (!c)
            return false;
    return true;
}

// 4.3、求多邊形圍成的面積
//
// 引數： polygon : 多邊形
//
double areaOfPolygon(const vector<Point>& polygon)
{
    // 如果多邊形點的數量少於三個則非法
    int size = polygon.size();
    if (size < 3) return 0;

    double area(0.0);
    for (int i = 0; i < size; ++i)
    {
        area += polygon[i].y * (polygon[(i - 1 + size) % size].x - polygon[(i + 1) % size].x);
    }

    return (area / 2);
}

// 4.4、判斷多邊形是否按照逆時針排列
//
// 引數： polygon : 多邊形
//
bool isConterClock(const vector<Point>& polygon)
{
    return areaOfPolygon(polygon) > 0;
}

// 4.5、判斷點是否在多邊形內部
// 判斷從點引出一條線段與多邊形的交點的個數
// 奇數個則相交、偶數個則不相交
// 
// 引數： polygon : 多邊形 p : 需要判斷的點
//
bool isPointInPolygon(const vector<Point>& polygon, const Point& p)
{
    Point down_left, up_right;
    boxOfPolygon(polygon, down_left, up_right);
    
    // 位於多邊形外部一點
    Point out_p = sub(down_left, Point(10.0, 0.0));
    
    int cnt(0);
    Line p_line(p, out_p, true);
    for (int i = 0; i < polygon.size(); ++i)
    {
        Point s = polygon[i];
        Point e = polygon[(i + 1) % polygon.size()];
        Line seg(s, e, true);
        Point inter_p;
        if (isSegIntersect(p_line, seg, inter_p))
        {
            cnt++;
        }
    }

    return (cnt % 2 == 1);
}

// 4.6、判斷線段是否在多邊形內部
// 線段在多邊形內部的條件是兩個端點都在多邊形內且不與多邊形相交
// 
// 引數： polygon : 多邊形 l ： 線段
bool isSegInPolygon(const vector<Point>& polygon, const Line& l)
{
    // 首先判斷線段端點是否都在多邊形內部
    bool is_s_in = isPointInPolygon(polygon, l.s);
    bool is_e_in = isPointInPolygon(polygon, l.e);

    // 然後判斷線段是否相交
    if (is_s_in && is_e_in)
    {
        for (int i = 0; i < polygon.size(); ++i)
        {
            Point s = polygon[i];
            Point e = polygon[(i + 1) % polygon.size()];
            Line seg(s, e, true);
            Point inter_p;
            if (isSegIntersect(l, seg, inter_p))
            {
                return false;
            }
        }
        return true;
    }
    else
    {
        return false;
    }
}

// 4.7、判斷圓是否在多邊形內部
// 只有多邊形所有的邊都在圓的外部，圓才處於多邊形內部
//
//    引數： polygon : 多邊形 c : 圓心 radius ： 半徑
//
bool isCircleInPolygon(const vector<Point>& polygon, const Point& c, double radius)
{
    for (int i = 0; i < polygon.size(); ++i)
    {
        const Point& p1 = polygon[i];
        const Point& p2 = polygon[(i + 1) % polygon.size()];
        Line line(p1, p2, true);
        if (segToCircle(c, radius, line) != 2)
            return false;
    }
    return true;
}

// 4.8、尋找點集凸包演算法（graham演算法）
// 演算法連結：https://blog.csdn.net/acm_zl/article/details/9342631
//
// 引數： points ： 平面點集
//
vector<Point> findConvexGraham(const vector<Point>& points)
{
    vector<Point> result;

    // 點的數量必須大於三個
    if (points.size() < 3)
        return result;

    // 尋找最底部的點
    int index = 0;
    for (int i = 0; i < points.size(); ++i)
    {
        if (points[i].y < points[index].y)
        {
            index = i;
        }
    }
    Point convex_p = points[index];

    // 計算每個點的極角
    map<double, int> cos_map;
    Point x_vec(1.0, 0.0);
    for (int i = 0; i < points.size(); ++i)
    {
        if (i != index)
        {
            double cos_value = Cos(sub(points[i], convex_p), x_vec);
            // 如果有多個點有相同的極角，則取最遠的點
            if (cos_map.count(-cos_value) != 0)
            {
                if (length(points[i]) > length(points[cos_map[-cos_value]]))
                    cos_map[-cos_value] = i;
            }
            else
                cos_map[-cos_value] = i;
        }
    }
    
    // 儲存結果的棧
    stack<int> result_stack;
    // 存入開始的兩個點
    result_stack.push(index);
    result_stack.push(cos_map.begin()->second);

    for (auto iter = (++cos_map.begin()); iter != cos_map.end(); ++iter)
    {
        int first = result_stack.top();
        result_stack.pop();
        int second = result_stack.top();

        Point vec1 = sub(points[first], points[second]);
        Point vec2 = sub(points[iter->second], points[first]);
        if (multiply(vec1, vec2).z >= 0)
            result_stack.push(first);
        result_stack.push(iter->second);
    }

    // 將資料從棧中讀取
    while (!result_stack.empty())
    {
        result.push_back(points[result_stack.top()]);
        result_stack.pop();
    }

    std::reverse(result.begin(), result.end());

    return result;
}

// 4.9、尋找點集凸包演算法（上下凸包法）時間複雜度O(nlogn)
//
//    引數： points : 平面點集
//
// 點根據字典序的比較函式
bool cmp(Point a, Point b)
{
    if (a.x == b.x)
        return a.y < b.y;
    return a.x < b.x;
}
vector<Point> findConvex(const vector<Point>& points)
{
    vector<Point> result;
    if (points.size() < 3)
        return result;

    vector<Point> tmp_points = points;
    // 首先將所有點按照字典序排序
    sort(tmp_points.begin(), tmp_points.end(), cmp);

    // 上凸包
    vector<Point> upper_hull;
    // 存入第一個和第二個點
    upper_hull.push_back(tmp_points[0]);
    upper_hull.push_back(tmp_points[1]);
    for (int i = 2; i < tmp_points.size(); ++i)
    {
        upper_hull.push_back(tmp_points[i]);
        while (upper_hull.size() > 2 && multiply(sub(upper_hull[upper_hull.size() - 2], upper_hull[upper_hull.size() - 3]), sub(upper_hull[upper_hull.size() - 1], upper_hull[upper_hull.size() - 3])).z >= 0)
        {
            upper_hull.erase(upper_hull.end() - 2);
        }
    }
    // 下凸包
    vector<Point> lower_hull;
    // 存入倒數第一第二個點
    lower_hull.push_back(tmp_points[tmp_points.size() - 1]);
    lower_hull.push_back(tmp_points[tmp_points.size() - 2]);
    for (int i = tmp_points.size() - 3; i >= 0; --i)
    {
        lower_hull.push_back(tmp_points[i]);
        while (lower_hull.size() > 2 && multiply(sub(lower_hull[lower_hull.size() - 2], lower_hull[lower_hull.size() - 3]), sub(lower_hull[lower_hull.size() - 1], lower_hull[lower_hull.size() - 3])).z >= 0)
        {
            lower_hull.erase(lower_hull.end() - 1);
        }
    }
    // 刪除重複點
    lower_hull.erase(lower_hull.begin());
    lower_hull.erase(lower_hull.end() - 1);

    // 合併上下凸包
    upper_hull.insert(upper_hull.end(), lower_hull.begin(), lower_hull.end());

    result = upper_hull;

    return result;
}

// 4.10、求簡單多邊形重心
// 演算法原理連結：
//
// 引數： polygon ： 簡單多邊形
//
Point centerOfPolygon(const vector<Point>& polygon)
{
    double polygon_area(0.0);
    Point center;
    Point origin;

    for (int i = 0; i < polygon.size(); ++i)
    {
        Point curr_p = polygon[i];
        Point next_p = polygon[(i + 1) % polygon.size()];
        Triangle t(origin, curr_p, next_p);

        double curr_area = areaOfTriangle(t);
        polygon_area += curr_area;

        center = add(center, mul(div(add(curr_p, next_p), 3), curr_area));
    }

    center = div(center, polygon_area);

    return center;
}

// 4.11、求肯定在多邊形內部的一個點
// 定理1: 每個多邊形至少有一個凸頂點，
//          x座標最大、最小的點肯定是凸頂點，y座標最大、最小的點肯定是凸頂點
// 定理2：頂點數>= 4的簡單多邊形至少有一條對角線
//
// 引數： polygon ： 簡單多邊形
//
Point pointInPolygon(const vector<Point>& polygon)
{
    // 凸頂點和索引
    int index = 0;
    Point convex_p = polygon[0];
    // 尋找一個凸頂點
    for (int i = 0; i < polygon.size(); ++i)
    {
        if (polygon[i].y < convex_p.y)
        {
            index = i;
            convex_p = polygon[i];
        }
    }
    // 獲取凸頂點前後一個點
    int size = polygon.size();
    Point pre_p = polygon[(index - 1 + size) % size];
    Point next_p = polygon[(index + 1) % size];
    Triangle t(convex_p, pre_p, next_p);
    double min_d = double(INT_MAX);
    bool flag = false;
    Point min_p;
    for (int i = 0; i < polygon.size(); ++i)
    {
        if (i == index || i == ((index - 1 + size) % size) || i == ((index + 1) % size))
            continue;
        flag = true;
        if (distance(convex_p, polygon[i]) < min_d)
        {
            min_p = polygon[i];
            min_d = distance(convex_p, polygon[i]);
        }
    }
    // 如何沒有頂點在三角形內部，則返回前後點的中點
    if (!flag)
    {
        return div(add(pre_p, next_p), 2);
    }
    // 返回最近點和凸頂點的中點
    return div(add(convex_p, min_p), 2);
}

// 4.12、獲取多邊形的包圍輪廓
// 即多邊形的最小包圍盒，由左下和右上兩個點表示
//
// 引數： polygon : 多邊形 down_left : 左下點  up_right : 右上點
//
void boxOfPolygon(const vector<Point>& polygon, Point& down_left, Point& up_right)
{
    double max_x = double(INT_MIN), min_x = double(INT_MAX);
    double max_y = double(INT_MIN), min_y = double(INT_MAX);

    for (auto c : polygon)
    {
        max_x = (c.x > max_x) ? c.x : max_x;
        min_x = (c.x < min_x) ? c.x : min_x;
        max_y = (c.y > max_y) ? c.y : max_y;
        min_y = (c.y < min_y) ? c.y : min_y;
    }

    down_left = Point(min_x, min_y);
    up_right = Point(max_x, max_y);
}


// 五、圓
// 5.1、點和圓的關係
//
// 引數： c: 圓心 radiuns ： 圓的半徑  p : 判斷的點
// 返回值 ： 0 ： 圓內 1 ： 圓上 2： 圓外
//
int pointToCircle(const Point& c, double radius, const Point& p)
{
    double ptoc_d = distance(c, p);

    if (ptoc_d < radius)
        return 0;
    else if (ptoc_d == radius)
        return 1;
    else
        return 2;
}

// 5.2、直線和圓的關係
//
// 引數： c: 圓心 radiuns ： 圓的半徑  l : 判斷的直線
// 返回值 ： 0 ： 相交 1 ：相切  2： 相離
//
int lineToCircle(const Point& c, double radius, const Line& l)
{
    double ctol_d = ptolDistance(c, l);

    if (ctol_d < radius)
        return 0;
    else if (ctol_d == radius)
        return 1;
    else
        return 2;
}

// 5.3、線段和圓的關係
//
// 引數： c: 圓心 radiuns ： 圓的半徑  l : 判斷的線段
// 返回值 ： 0 ： 圓內 1 ： 與圓相交  2： 圓外
//
int segToCircle(const Point& c, double radius, const Line& l)
{
    double ctol_d = ptolDistance(c, l);

    if (ctol_d > radius)
        return 2;
    else if (ctol_d == radius)
        return 1;
    else
    {
        Point project_p = ptolProjection(c, l);
        if (isponl(project_p, l))
            return 1;
        else
            return 2;
    }
}

// 5.4、兩圓之間的關係
// 
// 引數： c1 : 圓1圓心，r1 圓1半徑 c2 : 圓2圓心，r2 圓2半徑
// 返回值：0 ：內含 1：內切 2：相交 3：外切 4：外離
//
int circleToCircle(const Point& c1, double r1, const Point& c2, double r2)
{
    double ctoc_d = distance(c1, c2);

    if (ctoc_d < abs(r1 - r2))
        return 0;
    else if (ctoc_d == abs(r1 - r2))
        return 1;
    else if (ctoc_d > abs(r1 - r2) && ctoc_d < (r1 + r2))
        return 2;
    else if (ctoc_d == (r1 + r2))
        return 3;
    else if (ctoc_d > (r1 + r2))
        return 4;
}

&n