diff --git a/Compute_Geometry.h b/Compute_Geometry.h index 66da5a8..7f1f6ab 100644 --- a/Compute_Geometry.h +++ b/Compute_Geometry.h @@ -1,1296 +1,1296 @@ -// 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 -#include -#include -#include -#include -#include -#include - -using namespace std; - -namespace Compute_Geometry{ -/*** 常用常量 ***/ -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); - -// 1.15、在（-1，-1）到（1，1）随机生成n个点 -vector randomGenPoints(int num); - -// 二、线 -// 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 firstOrderBezier(const Point& s, const Point& e, int inter_num); - -// 2.4、二阶贝塞尔曲线插值 -vector secondOrderBezier(const Point& s, const Point& e, const Point& p, int inter_num); - -// 2.5、三阶贝塞尔曲线插值 -vector thirdOrderBezier(const Point& s, const Point& e, const Point& p1, const Point& p2, int inter_num); - -// 2.6、在（-1，-1）到（1，1）随机生成n条线 -vector randomGenLines(int 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& polygon, vector& flags); - -// 4.2、判断多边形是否为凸多边形 -bool isConvex(const vector& polygon); - -// 4.3、求多边形围成的面积 -double areaOfPolygon(const vector& polygon); - -// 4.4、判断多边形是否按照逆时针排列 -bool isConterClock(const vector& polygon); - -// 4.5、判断点是否在多边形内部 -bool isPointInPolygon(const vector& polygon, const Point& p); - -// 4.6、判断线段是否在多边形内部 -bool isSegInPolygon(const vector& polygon, const Line& l); - -// 4.7、判断圆是否在多边形内部 -bool isCircleInPolygon(const vector& polygon, const Point& p, double radius); - -// 4.8、寻找点集凸包算法（graham算法） -vector findConvexGraham(const vector& points); - -// 4.9、寻找点集凸包算法（上下凸包法）时间复杂度O(nlogn) -vector findConvex(const vector& points); - -// 4.10、求简单多边形重心 -Point centerOfPolygon(const vector& polygon); - -// 4.11、求肯定在多边形内部的一个点 -Point pointInPolygon(const vector& polygon); - -// 4.12、计算多边形的范围 -void boxOfPolygon(const vector& 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); -} - -// 1.15、在（-1，-1）到（1，1）随机生成n个点 -// -// 参数： down_left : 区域左下角 up_right: 区域右上角 num : 生成点的数量 -// -vector randomGenPoints(int num) -{ - vector result; - - std::uniform_real_distribution dist(-0.9, 0.9); - std::mt19937 rng; - rng.seed(std::random_device{}()); - for (int i = 0; i < num; ++i) - { - double rand_x = dist(rng); - double rand_y = dist(rng); - - result.push_back(Point(rand_x, rand_y)); - } - return result; -} - -// 二、线 - -// 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 firstOrderBezier(const Point& s, const Point& e, int inter_num) -{ - vector 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 secondOrderBezier(const Point& s, const Point& e, const Point& p, int inter_num) -{ - vector 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 thirdOrderBezier(const Point& s, const Point& e, const Point& p1, const Point& p2, int inter_num) -{ - vector 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; -} - -// 2.6、在（-1，-1）到（1，1）随机生成n条线 -// -// 参数 num : 需要生成的线段的数量 -// -vector randomGenLines(int num) -{ - vector result; - - std::uniform_real_distribution dist(-0.9, 0.9); - std::mt19937 rng; - rng.seed(std::random_device{}()); - for (int i = 0; i < num; ++i) - { - double rand_sx = dist(rng); - double rand_sy = dist(rng); - Point p1(rand_sx, rand_sy); - double rand_ex = dist(rng); - double rand_ey = dist(rng); - Point p2(rand_ex, rand_ey); - - result.push_back(Line(p1,p2,true)); - } - return result; -} - -// 三、三角形 - -// 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& polygon, vector& 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& polygon) -{ - vector flags; - // 判断每个点的凹凸性 - checkConvex(polygon, flags); - // 如果有一个点不是凸的，则此多边形为非凸 - for (auto c : flags) - if (!c) - return false; - return true; -} - -// 4.3、求多边形围成的面积 -// -// 参数： polygon : 多边形 -// -double areaOfPolygon(const vector& 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& polygon) -{ - return areaOfPolygon(polygon) > 0; -} - -// 4.5、判断点是否在多边形内部 -// 判断从点引出一条线段与多边形的交点的个数 -// 奇数个则相交、偶数个则不相交 -// -// 参数： polygon : 多边形 p : 需要判断的点 -// -bool isPointInPolygon(const vector& 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& 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& 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 findConvexGraham(const vector& points) -{ - vector 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 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 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 findConvex(const vector& points) -//{ -// vector result; -// if (points.size() < 3) -// return result; -// -// vector tmp_points = points; -// // 首先将所有点按照字典序排序 -// sort(tmp_points.begin(), tmp_points.end(), cmp); -// -// // 上凸包 -// vector 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 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& 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& 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& 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; -} - -} +// 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 +#include +#include +#include +#include +#include +#include + +using namespace std; + +namespace Compute_Geometry{ +/*** 常用常量 ***/ +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); + +// 1.15、在（-1，-1）到（1，1）随机生成n个点 +vector randomGenPoints(int num); + +// 二、线 +// 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 firstOrderBezier(const Point& s, const Point& e, int inter_num); + +// 2.4、二阶贝塞尔曲线插值 +vector secondOrderBezier(const Point& s, const Point& e, const Point& p, int inter_num); + +// 2.5、三阶贝塞尔曲线插值 +vector thirdOrderBezier(const Point& s, const Point& e, const Point& p1, const Point& p2, int inter_num); + +// 2.6、在（-1，-1）到（1，1）随机生成n条线 +vector randomGenLines(int 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& polygon, vector& flags); + +// 4.2、判断多边形是否为凸多边形 +bool isConvex(const vector& polygon); + +// 4.3、求多边形围成的面积 +double areaOfPolygon(const vector& polygon); + +// 4.4、判断多边形是否按照逆时针排列 +bool isConterClock(const vector& polygon); + +// 4.5、判断点是否在多边形内部 +bool isPointInPolygon(const vector& polygon, const Point& p); + +// 4.6、判断线段是否在多边形内部 +bool isSegInPolygon(const vector& polygon, const Line& l); + +// 4.7、判断圆是否在多边形内部 +bool isCircleInPolygon(const vector& polygon, const Point& p, double radius); + +// 4.8、寻找点集凸包算法（graham算法） +vector findConvexGraham(const vector& points); + +// 4.9、寻找点集凸包算法（上下凸包法）时间复杂度O(nlogn) +vector findConvex(const vector& points); + +// 4.10、求简单多边形重心 +Point centerOfPolygon(const vector& polygon); + +// 4.11、求肯定在多边形内部的一个点 +Point pointInPolygon(const vector& polygon); + +// 4.12、计算多边形的范围 +void boxOfPolygon(const vector& 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); +} + +// 1.15、在（-1，-1）到（1，1）随机生成n个点 +// +// 参数： down_left : 区域左下角 up_right: 区域右上角 num : 生成点的数量 +// +vector randomGenPoints(int num) +{ + vector result; + + std::uniform_real_distribution dist(-0.9, 0.9); + std::mt19937 rng; + rng.seed(std::random_device{}()); + for (int i = 0; i < num; ++i) + { + double rand_x = dist(rng); + double rand_y = dist(rng); + + result.push_back(Point(rand_x, rand_y)); + } + return result; +} + +// 二、线 + +// 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 firstOrderBezier(const Point& s, const Point& e, int inter_num) +{ + vector 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 secondOrderBezier(const Point& s, const Point& e, const Point& p, int inter_num) +{ + vector 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 thirdOrderBezier(const Point& s, const Point& e, const Point& p1, const Point& p2, int inter_num) +{ + vector 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; +} + +// 2.6、在（-1，-1）到（1，1）随机生成n条线 +// +// 参数 num : 需要生成的线段的数量 +// +vector randomGenLines(int num) +{ + vector result; + + std::uniform_real_distribution dist(-0.9, 0.9); + std::mt19937 rng; + rng.seed(std::random_device{}()); + for (int i = 0; i < num; ++i) + { + double rand_sx = dist(rng); + double rand_sy = dist(rng); + Point p1(rand_sx, rand_sy); + double rand_ex = dist(rng); + double rand_ey = dist(rng); + Point p2(rand_ex, rand_ey); + + result.push_back(Line(p1,p2,true)); + } + return result; +} + +// 三、三角形 + +// 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& polygon, vector& 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& polygon) +{ + vector flags; + // 判断每个点的凹凸性 + checkConvex(polygon, flags); + // 如果有一个点不是凸的，则此多边形为非凸 + for (auto c : flags) + if (!c) + return false; + return true; +} + +// 4.3、求多边形围成的面积 +// +// 参数： polygon : 多边形 +// +double areaOfPolygon(const vector& 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& polygon) +{ + return areaOfPolygon(polygon) > 0; +} + +// 4.5、判断点是否在多边形内部 +// 判断从点引出一条线段与多边形的交点的个数 +// 奇数个则相交、偶数个则不相交 +// +// 参数： polygon : 多边形 p : 需要判断的点 +// +bool isPointInPolygon(const vector& 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& 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& 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 findConvexGraham(const vector& points) +{ + vector 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 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 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 findConvex(const vector& points) +//{ +// vector result; +// if (points.size() < 3) +// return result; +// +// vector tmp_points = points; +// // 首先将所有点按照字典序排序 +// sort(tmp_points.begin(), tmp_points.end(), cmp); +// +// // 上凸包 +// vector 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 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& 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& 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& 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; +} + +} diff --git a/ExportQBGISCtrlClass.cpp b/ExportQBGISCtrlClass.cpp index ab991cb..ac77651 100644 --- a/ExportQBGISCtrlClass.cpp +++ b/ExportQBGISCtrlClass.cpp @@ -407,7 +407,7 @@ extern "C" void WINAPI UpdateDataInfo(int index,const char* key,double value ) //功能：目标追踪测试接口 //输入：经纬度，点的像素大小，默认20 -extern "C" void WINAPI OnShowTargetPoint(const double lon, const double lat, const int pixelSize) +extern "C" void WINAPI OnShowTargetPoint(double lon, double lat, int pixelSize) { AFX_MANAGE_STATE(AfxGetAppModuleState()); if (g_bCreateMap) diff --git a/ExportQBGISCtrlClass.h b/ExportQBGISCtrlClass.h index c0f117f..848e684 100644 --- a/ExportQBGISCtrlClass.h +++ b/ExportQBGISCtrlClass.h @@ -108,7 +108,7 @@ extern "C" _declspec(dllexport) void DrawCallBackPoint(const BYTE callbackMode, //功能：目标追踪测试接口 //输入：经纬度，点的像素大小，默认20 -extern "C" _declspec(dllexport) void OnShowTargetPoint(const double lon, const double lat, const int pixelSize=20); +extern "C" _declspec(dllexport) void OnShowTargetPoint(double lon, double lat, int pixelSize); diff --git a/FlyLineDesign.cpp b/FlyLineDesign.cpp index c6212fb..8bc37dc 100644 --- a/FlyLineDesign.cpp +++ b/FlyLineDesign.cpp @@ -346,7 +346,7 @@ void CFlyLineDesign::OnSaveFlyLine() memset(&pt, 0, sizeof(PtStruct)); //原点数据写入 - fprintf(fp, "%d, %d, %lf, %lf, %.2lf, %d, %02X, %02X\n", designLineID, 0, g_gcsLon, g_gcsLat, 0.0, 0, 0, 2); + fprintf(fp, "%d, %d, %.7f, %.7f, %.2lf, %d, %02X, %02X\n", designLineID, 0, g_gcsLon, g_gcsLat, 0.0, 0, 0, 2); //遍历航点集合数据 for ( int i = 0; i < m_linePtNums; i++ ) @@ -356,11 +356,11 @@ void CFlyLineDesign::OnSaveFlyLine() //将除最后一个点的数据写入文件 if ( i < (m_linePtNums-1)) { - fprintf(fp, "%d, %d, %lf, %lf, %.2lf, %d, %02X, %02X\n", designLineID, pt.nPt, pt.dX, pt.dY, pt.nH, pt.nV, pt.ch1, pt.ch2); + fprintf(fp, "%d, %d, %.7f, %.7f, %.2lf, %d, %02X, %02X\n", designLineID, pt.nPt, pt.dX, pt.dY, pt.nH, pt.nV, pt.ch1, pt.ch2); } else//将最后一个点的航段特征数据写入文件 { - fprintf(fp, "%d, %d, %lf, %lf, %.2lf, %d, %02X, %02X\n", designLineID, pt.nPt, pt.dX, pt.dY, pt.nH, pt.nV, pt.ch1, lineProperty); + fprintf(fp, "%d, %d, %.7f, %.7f, %.2lf, %d, %02X, %02X\n", designLineID, pt.nPt, pt.dX, pt.dY, pt.nH, pt.nV, pt.ch1, lineProperty); } } diff --git a/GISControlDlg.rc b/GISControlDlg.rc index 7266140..a6b5c7c 100644 --- a/GISControlDlg.rc +++ b/GISControlDlg.rc @@ -679,6 +679,7 @@ BEGIN MENUITEM "开始编辑", ID_EDIT_NODE MENUITEM "修改航点", ID_EDIT_LINE MENUITEM "保存编辑", ID_EDIT_SAVE + MENUITEM "批量航线", ID_COPE_LINE END POPUP "限制区" BEGIN diff --git a/GISDlg.cpp b/GISDlg.cpp index 4d0629a..bd4ab35 100644 --- a/GISDlg.cpp +++ b/GISDlg.cpp @@ -188,6 +188,17 @@ CGISDlg::CGISDlg(CWnd* pParent /*=NULL*/) m_CircleGuideLayer = -1; //盘旋引导图层 + //多机子航线 + for (int i=0;i<20;++i) + { + for (int j=0;j<20;++j) + { + m_subLineLayerID[i][j]=-1; + m_subLinePtLayerID[i][j]=-1; + m_pHaveDrawCopyLineFlag[i][j] = -1; + } + } + m_bShowCopyLine = false; /*******************************测绘航线********************************/ surveyRegionLayerID = -1; tmpSurveyRegionLayerID = -1; @@ -583,6 +594,7 @@ BEGIN_MESSAGE_MAP(CGISDlg, CBCGPDialog) ON_COMMAND(ID_EDIT_NODE, OnEditLine) //编辑航线 ON_COMMAND(ID_EDIT_LINE, ShowModifyPointDlg) //修改航点 ON_COMMAND(ID_EDIT_SAVE, ShowEditSaveDlg) //保存编辑 + ON_COMMAND(ID_COPE_LINE, ShowMultiRouteSetting) //批量航线 /* ON_COMMAND(ID_MODIFY_LINEPOINT, OnMapSelect); //修改航点航线 ON_COMMAND(ID_EDIT_SAVE, OnMapSelect); //保存编辑*/ @@ -3717,7 +3729,7 @@ void CGISDlg::AddTargetPoint(const int nPt, const double dX, const double dY) //说明：创建专门用于绘制无人机的矢量面图层，每次先删除上次标绘的Shp，再重新标绘 void CGISDlg::DrawUAV(int uavid, const double dX, const double dY, const double yaw,bool control) { - if (uavid<0||uavid>1000) + if (uavid<0||uavid>2) { return; } @@ -4627,7 +4639,13 @@ void CGISDlg::DrawFlyLine(const DrawLineDataStruct lineData) } else { - strText.Format(_T("%d-%d"), lineData.lineID, lineData.pts[i].nPt); +/* + if (g_b981AMulti && lineData.pts[0].ch1 > 1) + { + strText.Format(_T("%d(%d)-%d"), lineData.lineID, lineData.pts[0].ch1,lineData.pts[i].nPt); + } + else*/ + strText.Format(_T("%d-%d"), lineData.lineID, lineData.pts[i].nPt); } //增加Labels @@ -4639,6 +4657,95 @@ void CGISDlg::DrawFlyLine(const DrawLineDataStruct lineData) m_map.Redraw(); } +//功能：标绘飞行航线(批量航线) +//输入：航线数据lineData +void CGISDlg::DrawCopyFlyLine(const DrawLineDataStruct lineData) +{ + //field索引值 + long fieldIndex = 0; + + //航线索引号 + int lineID = lineData.lineID-1; + //子航线索引号 + int sublineID = lineData.pts[0].ch1-1; +/* + map subLineLayerIDs; + map subLinePtLayerIDs; + if (m_subLineLayerID.find(lineData.lineID)!=m_subLineLayerID.end()) + { + subLineLayerIDs = m_subLineLayerID[lineData.lineID]; + } + if (m_subLinePtLayerID.find(lineData.lineID)!=m_subLinePtLayerID.end()) + { + subLinePtLayerIDs = m_subLinePtLayerID[lineData.lineID]; + }*/ + + //没有目标标绘SHP图层，创建该图层 + if (m_subLineLayerID[lineID][sublineID] ==-1) + { + //创建线图层 + CreateEmptyShapfile(m_subLineLayerID[lineID][sublineID], 1, /*RGB(0,255,0)*/LineClr[lineID]); + //创建线图层成功 + m_pHaveDrawCopyLineFlag[lineID][sublineID] = true; + } + + + //没有目标标绘SHP图层，创建该图层 + if (m_subLinePtLayerID[lineID][sublineID] ==-1) + { + //创建点图层 + CreateEmptyShapfile(m_subLinePtLayerID[lineID][sublineID], 0, /*RGB(0,255,0)*/LineClr[lineID]); + + //创建线图层成功 + m_pHaveDrawCopyLineFlag[lineID][sublineID] = true; + } + + //向点图层加入航点数据 + AddPoints2PointShapfile(m_subLinePtLayerID[lineID][sublineID], lineData.pts, lineData.pointNum); + + + //向线图层加入航点数据 + AddPoints2LineShapfile(m_subLineLayerID[lineID][sublineID], lineData.linePts, lineData.linePointNum); + + /////////////////////往SHP图层中加入目标点标注信息/////////////////////////////// + + //Label集合 + CLabels labesPtr; + labesPtr = (m_map.GetShapefile(m_subLinePtLayerID[lineID][sublineID])).GetLabels(); + + labesPtr.SetFontColor(/*RGB(255,0,0)*/LineClr[lineID]); + labesPtr.SetAlignment(1); + labesPtr.SetFontBold(true); + labesPtr.SetFontName(LPCTSTR("黑体")); + labesPtr.SetFontSize(16); + labesPtr.put_FontSize2(16); + labesPtr.SetFrameVisible(false); + labesPtr.SetAvoidCollisions(FALSE); + + + CString strText = _T(""); + + //标绘航点的名称 + for (long i=0; i14 /*&& lineDataGroup.lineID<10)*/ /*|| lineDataGroup.lineID==14 *//*|| lineDataGroup.lineID>13*/) + if(lineDataGroup.lineID<1 || lineDataGroup.lineID>14 ) { BCGPMessageBox("航线号不正确！"); return; } - + + m_lineFullPathName[lineDataGroup.lineID] = strLineFileName; + if (lineDataGroup.pts[0].ch1>1) + { + //清除标绘的航线 + ClearDrawedCopyLine(lineDataGroup.lineID-1,lineDataGroup.pts[0].ch1); + + //标绘航线 + DrawCopyFlyLine(lineDataGroup); + return; + } + int ch = lineDataGroup.linePts[lineDataGroup.pointNum-1].ch1; /***************************测绘区域处理************************************/ if (ch == 0x03) @@ -5635,7 +5802,7 @@ void CGISDlg::OnEnddesign() // if (IDOK == dlgSave.DoModal()) { strFlyLineName = dlgSave.GetPathName(); - pathName = ExtractFileName(strFlyLineName); + pathName = ExtractFileName(strFlyLineName,false); ////写出航线数据 FILE *fp = fopen(strFlyLineName, "w" ); @@ -5684,7 +5851,8 @@ void CGISDlg::OnEnddesign() // { vector> resultLines; GetMultiRouteLine(azmuth,lineInterval,numLine,m_pDesignLineStruct,m_designLinePointNum,resultLines); - SaveMultiRouteLine(pathName,heightInterval,resultLines); + vectorsavePathNameArr; + SaveMultiRouteLine(pathName,heightInterval,savePathNameArr,resultLines); } memset(m_pDesignLineStruct, 0, sizeof(PtStruct)*m_lineMaxPointNum); //初始化新建一条航线 @@ -8444,7 +8612,7 @@ void CGISDlg::GetMultiRouteLine(double azimuth,double lineInterval,int lineNumbe } } -void CGISDlg::SaveMultiRouteLine(CString pathDirName,double heightInterval,const vector>& resultLines) +void CGISDlg::SaveMultiRouteLine(CString pathDirName,double heightInterval,vector& savePathNameArr,const vector>& resultLines) { CString saveFileName; for (int i=0;ipathNameArr; + if (g_b981AMulti) + { + SetMultiRouteDlg multiRouteDlg; + if (multiRouteDlg.DoModal() == IDOK ) + { + azmuth = multiRouteDlg.azmuth; + heightInterval = multiRouteDlg.heightInterval; + lineInterval = multiRouteDlg.lineInterval; + numLine = multiRouteDlg.numLine; + bSaveMultiLine = true; + + CString selectLinePath = m_lineFullPathName[m_lineSelectedID]; + if (ExtractFileType(selectLinePath) == ".route") return; // 已装订的航线不允许批量生成 + DrawLineDataStruct selectLineData = m_ShowedLineDataList[m_lineSelectedID]; + CString dirPath = ExtractDirPath(selectLinePath); + CString fileName = ExtractFileName(selectLinePath,false); + pathName = dirPath + "\\" + fileName; + //批量航线保存 + vector> resultLines; + GetMultiRouteLine(azmuth,lineInterval,numLine,selectLineData.pts,selectLineData.pointNum,resultLines); + SaveMultiRouteLine(pathName,heightInterval,pathNameArr,resultLines); + } + else + { + bSaveMultiLine = false; + } + } + + ClearHighLightLine(); + m_bSelectFeatureFlag = false; + m_lineSelectedID==-1; + m_map.SetCursorMode(mapWindow::tkCursorMode::cmNone); + for (int i=0;i +#include #include "MapElevation.h" #include "../Include/8BMapDLL_type.h" #include "SaveZoneDlg.h" @@ -247,6 +248,13 @@ public: void ZoomToLocation(double lon,double lat); private: // by Wu + //map> m_subLineLayerID; //子航线线图层号 + //map> m_subLinePtLayerID; //子航线点图层号 + long m_subLineLayerID[20][20]; //子航线线图层号 + long m_subLinePtLayerID[20][20]; //子航线点图层号 + bool m_pHaveDrawCopyLineFlag[20][20]; + bool m_bShowCopyLine; + map m_lineFullPathName; // void ShowModifyPointDlg(int selectedPointID=0); //编辑航线时,显示航点设置对话框 void ShowEditSaveDlg(); //显示保存编辑对话框 @@ -319,14 +327,14 @@ private: bool m_bHaveAddMap; //航迹标绘图层ID - long m_flyTrackLayerID[1000]; + long m_flyTrackLayerID[3]; //UAV航迹标绘图层ID - long m_UAVFlyTrackLayerID[1000]; + long m_UAVFlyTrackLayerID[3]; long m_UAVFlyTrackLayerID_ADS[MAX_PLANE_NUM_ADS]; //是否已经标绘了无人机的飞行轨迹 - bool m_bFirstDrawUAVTrack[1000]; + bool m_bFirstDrawUAVTrack[3]; bool m_bFirstDrawUAVTrack_ADS[MAX_PLANE_NUM_ADS]; //地图显示区宽度、高度倒数 @@ -713,6 +721,9 @@ private: //功能：标绘飞行航线 //输入：航线数据lineData void DrawFlyLine(const DrawLineDataStruct lineData); + //功能：标绘飞行航线(批量航线) + //输入：航线数据lineData + void DrawCopyFlyLine(const DrawLineDataStruct lineData); afx_msg void OnUpdateLineDisplay1(CCmdUI* pCmdUI); afx_msg void OnUpdateLineDisplay2(CCmdUI* pCmdUI); @@ -727,6 +738,9 @@ private: //功能：清除已经标绘的航线 //输入：航线编号lineID，从0开始计数 void ClearDrawedLine(const int lineID); + //功能：清除已经标绘的航线(多机子航线) + //输入：航线编号lineID，从0开始计数 + void ClearDrawedCopyLine(const int lineID,const int sublineID); //功能：清除点图层标注信息 @@ -984,11 +998,12 @@ public: void OnShowDesignSurveyLineDlg(); /****************************多机航线一键生成************************************/ void GetMultiRouteLine(double azimuth,double lineInterval,int lineNumber,PtStruct* lineStruct,int nLinePts,vector>& resultLines); - void SaveMultiRouteLine(CString pathDirName,double heightInterval,const vector>& resultLines); + void SaveMultiRouteLine(CString pathDirName,double heightInterval,vector& savePathNameArr,const vector>& resultLines); + void ShowMultiRouteSetting(); /****************************吊舱目标追踪测试************************************/ long testLayer; - void OnShowTargetPoint(const double lon,const double lat,int pixelSize=20); + void OnShowTargetPoint(double lon,double lat,int pixelSize); }; //{{AFX_INSERT_LOCATION}} diff --git a/Globe.cpp b/Globe.cpp index ced09da..e22551c 100644 --- a/Globe.cpp +++ b/Globe.cpp @@ -97,6 +97,7 @@ extern bool CreateDirectory(const CString &strDir) } //提取没有后缀名的文件名 +/* extern CString ExtractFileName(CString fileName) { // 查找最后一个点的位置 @@ -114,7 +115,7 @@ extern CString ExtractFileName(CString fileName) // 如果没有找到点，说明没有后缀，直接显示文件名 return fileName; } -} +}*/ //功能：一个字节分成8位，分别存储到数组中 @@ -189,6 +190,58 @@ CString GetFlyLineName( const int lineID ) return (GetSoftwareCurrentDirectory() + _T( "\\route" ) + str + _T( ".txt" )); } +CString GetFlyLineDirName( const int lineID ) +{ + CString str; + if (lineID<=5) + { + str.Format( "航线%d", lineID ); + } + else if (lineID == 11) + { + str = "应急航线"; + } + else if( lineID == 12 ) + { + str = "电子围栏"; + } + else if (lineID == 14) + { + str = "回收航线"; + } + + return (GetSoftwareCurrentDirectory() + _T( "\\Route\\" ) + str ); +} + +//获取文件路径字符串 +CString ExtractDirPath(CString strFullPath) +{ + int n=strFullPath.ReverseFind('\\')+1; + return strFullPath.Left(n); +} +//获取文件名字字符串 +CString ExtractFileName(CString strFullPath,bool bIncludeType) +{ + CString strFile; + int n=strFullPath.GetLength()-strFullPath.ReverseFind('\\')-1; + strFile=strFullPath.Right(n); + if (bIncludeType==false) //去除文件名后缀 + { + n = strFile.ReverseFind('.'); + return strFile.Left(n); + } + else + { + return strFile; + } +} +//获取文件类型 +CString ExtractFileType(CString strFullPath) +{ + int n=strFullPath.GetLength()-strFullPath.ReverseFind('.')-1; + return strFullPath.Right(n); +} + //可视域分析相关变量 /************************************************************************/ diff --git a/Globe.h b/Globe.h index 1e91c1e..61e18ad 100644 --- a/Globe.h +++ b/Globe.h @@ -173,7 +173,7 @@ extern bool SearchDirectory(const CString &strDir); extern bool CreateDirectory(const CString &strDir); //提取没有后缀名的文件名 -extern CString ExtractFileName(CString fileName); +//extern CString ExtractFileName(CString fileName); //功能：一个字节分成8位，分别存储到数组中 //输入：字节oneByte @@ -193,6 +193,15 @@ extern CString GPSUTCTime2BeijingTime(const double utcTime); //获得航线文件名，输入航线号（1、2、3、4、5） extern CString GetFlyLineName( const int lineID ); +//获取航线文件夹名字 +extern CString GetFlyLineDirName( const int lineID ); + +//获取文件路径字符串 +extern CString ExtractDirPath(CString strFullPath); +//获取文件名字字符串 +extern CString ExtractFileName(CString strFullPath,bool bIncludeType=true); +//获取文件类型 +extern CString ExtractFileType(CString strFullPath); const int LINEPTNUM = 256; diff --git a/res/GISControlDlg.rc2 b/res/GISControlDlg.rc2 index 63186b0..23dc028 100644 --- a/res/GISControlDlg.rc2 +++ b/res/GISControlDlg.rc2 @@ -1,13 +1,13 @@ -// -// GISCONTROLDLG.RC2 - resources Microsoft Visual C++ does not edit directly -// - -#ifdef APSTUDIO_INVOKED - #error this file is not editable by Microsoft Visual C++ -#endif //APSTUDIO_INVOKED - - -///////////////////////////////////////////////////////////////////////////// -// Add manually edited resources here... - -///////////////////////////////////////////////////////////////////////////// +// +// GISCONTROLDLG.RC2 - resources Microsoft Visual C++ does not edit directly +// + +#ifdef APSTUDIO_INVOKED + #error this file is not editable by Microsoft Visual C++ +#endif //APSTUDIO_INVOKED + + +///////////////////////////////////////////////////////////////////////////// +// Add manually edited resources here... + +///////////////////////////////////////////////////////////////////////////// diff --git a/resource.h b/resource.h index b725614..2f68f82 100644 --- a/resource.h +++ b/resource.h @@ -426,13 +426,15 @@ #define ID__32940 32940 #define ID__DESIGNSURVEYLINE 32941 #define ID_DESIGNSURVEYLINE 32942 +#define ID_32943 32943 +#define ID_COPE_LINE 32944 // Next default values for new objects // #ifdef APSTUDIO_INVOKED #ifndef APSTUDIO_READONLY_SYMBOLS #define _APS_NEXT_RESOURCE_VALUE 1174 -#define _APS_NEXT_COMMAND_VALUE 32943 +#define _APS_NEXT_COMMAND_VALUE 32945 #define _APS_NEXT_CONTROL_VALUE 1056 #define _APS_NEXT_SYMED_VALUE 1000 #endif