From 86e2541bb1fc44a675bbdffecde4ed2c8b7bc05d Mon Sep 17 00:00:00 2001 From: cbwu <504-wuchengbo@htsdfp.com> Date: Mon, 8 Jul 2024 10:57:15 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=B5=8B=E7=BB=98=E8=88=AA=E7=BA=BF?= =?UTF-8?q?=E7=94=9F=E6=88=90=E7=AE=97=E6=B3=95=E4=BC=98=E5=8C=96=EF=BC=8C?= =?UTF-8?q?=E8=88=AA=E6=B5=8B=E5=8C=BA=E5=9F=9F=E7=A6=81=E6=AD=A2=E7=BB=98?= =?UTF-8?q?=E5=88=B6=E5=87=B9=E5=A4=9A=E8=BE=B9=E5=BD=A2=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Compute_Geometry.h | 1296 +++++++++++++++++++++++++++++++++++++++ GISControlDlg.rc | 22 +- GISDlg.cpp | 20 +- TopologicalAnalysis.cpp | 140 ++++- TopologicalAnalysis.h | 236 ++++++- designsurveylinedlg.cpp | 105 +++- designsurveylinedlg.h | 2 +- resource.h | 3 +- 8 files changed, 1792 insertions(+), 32 deletions(-) create mode 100644 Compute_Geometry.h diff --git a/Compute_Geometry.h b/Compute_Geometry.h new file mode 100644 index 0000000..66da5a8 --- /dev/null +++ b/Compute_Geometry.h @@ -0,0 +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; +} + +} diff --git a/GISControlDlg.rc b/GISControlDlg.rc index 108578a..7266140 100644 --- a/GISControlDlg.rc +++ b/GISControlDlg.rc @@ -397,7 +397,7 @@ BEGIN COMBOBOX IDC_COMBO_DIRECTION,67,49,73,63,CBS_DROPDOWNLIST | CBS_SORT | WS_VSCROLL | WS_TABSTOP END -IDD_DLG_DESIGNSURVEYLINE DIALOGEX 0, 0, 211, 119 +IDD_DLG_DESIGNSURVEYLINE DIALOGEX 0, 0, 211, 138 STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU CAPTION "航测航线设计" FONT 8, "MS Shell Dlg", 400, 0, 0x1 @@ -406,18 +406,20 @@ BEGIN EDITTEXT IDC_EDIT_HEIGHT,70,11,41,15,ES_AUTOHSCROLL,WS_EX_RIGHT LTEXT "航线间隔(m):",IDC_STATIC,10,32,57,10 EDITTEXT IDC_EDIT_LINEINTERVAL,70,29,41,15,ES_AUTOHSCROLL,WS_EX_RIGHT - PUSHBUTTON "导入航测区",IDC_BTN_INPUTREGION,10,68,55,18 - PUSHBUTTON "绘制航测区",IDC_BTN_DRAWREGION,69,68,55,18 - PUSHBUTTON "保存航测区",IDC_BTN_SAVEREGION,128,68,55,18 - PUSHBUTTON "装订航线",IDC_BTN_BINDLINE,10,90,55,18 - PUSHBUTTON "生成航线",IDC_BTN_CALCULATELINE,69,90,55,18 + PUSHBUTTON "导入航测区",IDC_BTN_INPUTREGION,10,90,55,18 + PUSHBUTTON "绘制航测区",IDC_BTN_DRAWREGION,69,90,55,18 + PUSHBUTTON "保存航测区",IDC_BTN_SAVEREGION,128,90,55,18 + PUSHBUTTON "装订航线",IDC_BTN_BINDLINE,10,112,55,18 + PUSHBUTTON "生成航线",IDC_BTN_CALCULATELINE,69,112,55,18 LTEXT "航程:",IDC_STATIC,114,14,44,10 LTEXT "0 m",IDC_TXT_ROUTELENGTH,140,14,50,10 LTEXT "航时:",IDC_STATIC,115,32,44,10 - LTEXT "飞行速度(m/s):",IDC_STATIC,6,50,66,10 - EDITTEXT IDC_EDIT_FLYSPEED,70,47,41,15,ES_AUTOHSCROLL,WS_EX_RIGHT - PUSHBUTTON "航时计算",IDC_BTN_CALCULATETIME,123,46,52,16 + LTEXT "飞行速度(m/s):",IDC_STATIC,6,68,66,10 + EDITTEXT IDC_EDIT_FLYSPEED,70,65,41,15,ES_AUTOHSCROLL,WS_EX_RIGHT + PUSHBUTTON "航时计算",IDC_BTN_CALCULATETIME,123,64,52,16 LTEXT "0时0分0秒",IDC_TXT_FLYTIME,141,32,59,10 + LTEXT "航线外扩(m):",IDC_STATIC,10,49,57,10 + EDITTEXT IDC_EDIT_EXTERNALLENGTH,70,46,41,15,ES_AUTOHSCROLL,WS_EX_RIGHT END IDD_DLG_SETMULTIROUTE DIALOGEX 0, 0, 210, 159 @@ -574,7 +576,7 @@ BEGIN LEFTMARGIN, 3 RIGHTMARGIN, 207 TOPMARGIN, 3 - BOTTOMMARGIN, 114 + BOTTOMMARGIN, 134 END IDD_DLG_SETMULTIROUTE, DIALOG diff --git a/GISDlg.cpp b/GISDlg.cpp index 2a286e2..13df20b 100644 --- a/GISDlg.cpp +++ b/GISDlg.cpp @@ -1932,8 +1932,24 @@ void CGISDlg::MouseDownMap1(short Button, short Shift, long x, long y) //像素坐标转换地理坐标 m_map.PixelToProj(x,y, &dX, &dY); - surveyRegionLons.push_back(dX); - surveyRegionLats.push_back(dY); + if (surveyRegionLons.size()>2) + { + TopologicalAnalysis analysis; + Point2D p; + p.x = dX; + p.y = dY; + if(!analysis.isPointInPolygon(p,surveyRegionLons,surveyRegionLats)) //禁止绘制凹多边形 + { + surveyRegionLons.push_back(dX); + surveyRegionLats.push_back(dY); + } + } + else + { + surveyRegionLons.push_back(dX); + surveyRegionLats.push_back(dY); + } + } if (m_bHaveShowDistanceDlg && (Button == 1)) diff --git a/TopologicalAnalysis.cpp b/TopologicalAnalysis.cpp index 12d3c35..b31894b 100644 --- a/TopologicalAnalysis.cpp +++ b/TopologicalAnalysis.cpp @@ -263,11 +263,13 @@ bool TopologicalAnalysis::GetBoundingBoxVertices(const vector& polygonX, } //计算直线与多边形的交点 -void TopologicalAnalysis::linePolygonIntersections(const Point2D& linePt1,const Point2D& linePt2,const vector& polygonX,const vector& polygonY, vector& result) +int TopologicalAnalysis::linePolygonIntersections(const Point2D& linePt1,const Point2D& linePt2,const vector& polygonX,const vector& polygonY, vector& result) { + Line2D line1; line1.p1 = linePt1; line1.p2 = linePt2; + line1.is_seg = false; Point2D intersectionPt; vector resPoints; for (int i=0;i0) + if (resPoints.size()>1) { //获取最左或者最右的两个点,忽略凹多边形中间交点 double minX = resPoints[0].x; @@ -304,6 +325,8 @@ void TopologicalAnalysis::linePolygonIntersections(const Point2D& linePt1,const result.push_back(resPoints[maxIndex]); } + return resPoints.size(); + } // 判断两条线段是否相交 @@ -332,8 +355,121 @@ int TopologicalAnalysis::isLineIntersecting(const Line2D& l1, const Line2D& l2, intersection.x = l1.p1.x + ua*(l1.p2.x - l1.p1.x); intersection.y = l1.p1.y + ua*(l1.p2.y - l1.p1.y); +/* + double a1 = l1.p2.y - l1.p1.y; + double b1 = l1.p1.x - l1.p2.x; + double c1 = l1.p2.x * l1.p1.y - l1.p1.x * l1.p2.y; + + double a2 = l2.p2.y - l2.p1.y; + double b2 = l2.p1.x - l2.p2.x; + double c2 = l2.p2.x * l2.p1.y - l2.p1.x * l2.p2.y; + double det = a1 * b2 - a2 * b1; + + intersection.x = (b2 * c1 - b1 * c2) / det; + intersection.y = (a1 * c2 - a2 * c1) / det;*/ + return 1; //相交 } return 0; //不相交 } + +//计算多边形最小宽度 +int TopologicalAnalysis::getLineWithPolygonMinWidth(const vector& polygonX,const vector& polygonY) +{ + int n = polygonX.size(); + double minWidth = INT64_MAX; + int lineID = 0; + + for (int i=0;i width) width = dist; + } + if(width < minWidth) + { + minWidth = width; + lineID = startID; + } + } + + return lineID; +} + +//判断点与直线的位置关系 +int TopologicalAnalysis::pointPosition(const Point2D& pt,const Line2D& line) +{ + double cross = crossProduct(line.p1, line.p2, pt); + if (cross > 0) { + return 1; //左侧 + } else if (cross < 0) { + return -1; //右侧 + } else { + return 0; //线上 + } +} + + +// 判断点P是否在多边形polygon内 +bool TopologicalAnalysis::isPointInPolygon(const Point2D& P, vector& regionLons,vector& regionLats) { + int n = regionLons.size(); + bool inside = false; + double xints; + for(int i=0, j=n-1; i P.y) != (yj > P.y)) && // 边的两端点在射线两侧 + (P.x < (xj - xi) * (P.y - yi) / (yj - yi) + xi)) { // 射线与边相交的条件 + inside = !inside; // 交点数加一，改变inside状态 + } + } + return inside; +} + +//直线与线段的交点 +bool TopologicalAnalysis::findIntersection(double x1, double y1, double x2, double y2, // 直线1的两个点 + double x3, double y3, double x4, double y4,Point2D& intersection) { // 线段的两个端点 + // 计算直线1的方向向量 + double dx1 = x2 - x1; + double dy1 = y2 - y1; + + // 计算线段的方向向量 + double dx2 = x4 - x3; + double dy2 = y4 - y3; + + // 计算向量积，判断是否平行（即无交点） + double cross = -dx1 * dy2 + dy1 * dx2; + if (cross == 0) return false; + + // 计算参数t1（直线1上的点）和t2（线段上的点） + double t1 = (dx2 * (y3 - y1) - dy2 * (x3 - x1)) / cross; + double t2 = (dx1 * (y3 - y1) - dy1 * (x3 - x1)) / cross; + + // 判断t2是否在线段上 + if (t2 >= 0 && t2 <= 1) { + // 计算并返回交点坐标 + intersection.x = x1 + t1 * dx1; + intersection.y = y1 + t1 * dy1; + return true; + } else { + // 如果t2不在[0, 1]区间内，说明不相交于线段上 + return false; + } +} \ No newline at end of file diff --git a/TopologicalAnalysis.h b/TopologicalAnalysis.h index 235b06d..9626cb8 100644 --- a/TopologicalAnalysis.h +++ b/TopologicalAnalysis.h @@ -5,14 +5,26 @@ #include #include #include "geocompute.h" +#include using namespace std; struct Point2D { - double x, y; + //double x, y; + double x; // x坐标 + double y; // y坐标 + double z; // z坐标（默认为0，如果需要三维点则给z赋值） + + Point2D(double a = 0, double b = 0, double c = 0) { x = a; y = b; z = c; } // 构造函数 }; struct Line2D { - Point2D p1, p2; + //Point2D p1, p2; + Point2D p1; // 起点 + Point2D p2; // 终点 + bool is_seg; // 是否是线段 + + Line2D() {}; // 默认构造函数 + Line2D(Point2D a, Point2D b, bool _is_seg = true) { p1 = a; p2 = b; is_seg = _is_seg; } // 构造函数(默认是线段) }; struct Rectangle2D { @@ -55,12 +67,230 @@ public: bool GetBoundingBoxVertices(const vector& polygonX,const vector& polygonY,vector& rectangleX,vector& rectangleY); //计算直线与多边形的交点 - void linePolygonIntersections(const Point2D& linePt1,const Point2D& linePt2,const vector& polygonX,const std::vector& polygonY,std::vector& result); + int linePolygonIntersections(const Point2D& linePt1,const Point2D& linePt2,const vector& polygonX,const std::vector& polygonY,std::vector& result); + + /** 计算多边形最小宽度对应的边 + * param: polygonX,polygonY,多边形顶点坐标(首尾坐标相同) + * return: 返回线段起始点ID号 + */ + int getLineWithPolygonMinWidth(const vector& polygonX,const vector& polygonY); + //判断点与直线的位置关系 + int pointPosition(const Point2D& pt,const Line2D& line); + // 判断点P是否在多边形polygon内 + bool isPointInPolygon(const Point2D& P, vector& regionLons,vector& regionLats); private: // 判断两条线段是否相交 int isLineIntersecting(const Line2D& l1, const Line2D& l2, Point2D& intersection); + + //数学计算 +private: + //计算两点（向量）的点积 + double dotProduct(const Point2D& p1, const Point2D& p2) { + return p1.x * p2.x + p1.y * p2.y; + } + //计算两个点（向量）之间的差向量 + Point2D subtract(const Point2D& p1, const Point2D& p2) { + Point2D point; + point.x = p1.x - p2.x; + point.y = p1.y - p2.y; + return point; + } + + //计算向量叉积 + double crossProduct(const Point2D& A, const Point2D& B, const Point2D& P) { + return (B.x - A.x) * (P.y - A.y) - (B.y - A.y) * (P.x - A.x); + } + +/* + double length(const Point2D& p) { + return std::sqrt(p.x * p.x + p.y * p.y); + }*/ + +/* + Point2D normalize(const Point2D& p) { + double len = length(p); + Point2D pt; + pt.x = p.x / len; + pt.y = p.y / len; + return pt; + }*/ + +private: + // 点的加法 + Point2D add(const Point2D& lhs, const Point2D& rhs) + { + Point2D res; + + res.x = lhs.x + rhs.x; + res.y = lhs.y + rhs.y; + res.z = lhs.z + rhs.z; + + return res; + } + // 点的减法 + Point2D sub(const Point2D& lhs, const Point2D& rhs) + { + Point2D res; + + res.x = lhs.x - rhs.x; + res.y = lhs.y - rhs.y; + res.z = lhs.z - rhs.z; + + return res; + } + // 向量的乘法 + Point2D mul(const Point2D& p, double ratio) + { + Point2D res; + + res.x = p.x * ratio; + res.y = p.y * ratio; + res.z = p.z * ratio; + + return res; + } + // 向量的除法 + Point2D div(const Point2D& p, double ratio) + { + Point2D res; + + res.x = p.x / ratio; + res.y = p.y / ratio; + res.z = p.z / ratio; + + return res; + } + + // 点判断相等 + bool equal(const Point2D& lhs, const Point2D& rhs) + { + return(lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z); + } + + // 1.3、矢量标准化（矢量的长度规约到1） + // + // 参数： vec ： 矢量 + // + Point2D normalize(const Point2D& vec) + { + Point2D res; + + res = div(vec, length(vec)); + + return res; + } + + // 1.9、点是否在线上 + // 线分为直线和线段，直线表示的是直线是否经过点 + // + // 参数：p : 点 l : 线段或者线 + // + bool isponl(const Point2D& p, const Line2D& l) + { + Point2D line_vec = sub(l.p2, l.p1); + Point2D point_vec1 = sub(p, l.p1); + Point2D point_vec2 = sub(p, l.p2); + + Point2D mul_vec = multiply(line_vec, point_vec1); + double dot = dotMultiply(point_vec1, point_vec2); + // 点是否在线段上 + if (l.is_seg) + { + if (equal(p, l.p1) || equal(p, l.p2)) + return true; + return (0.0 == length(mul_vec) && dot < 0.0); + + } + // 点是否在直线上 + else + { + return (0.0 == length(mul_vec)); + } + } + + // 1.4、矢量点乘 + // + // 参数：(p1-op)为矢量1，（p2-op）为矢量2 + // + double dotMultiply(const Point2D& op, const Point2D& p1, const Point2D& 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 Point2D& vec1, const Point2D& vec2) + { + return(vec1.x * vec2.x + vec1.y * vec2.y + vec1.z * vec2.z); + } + + // 1.5、矢量叉乘 + // + // 参数：(p1-op)为矢量1，（p2-op）为矢量2 + // + Point2D multiply(const Point2D& op, const Point2D& p1, const Point2D& p2) + { + Point2D 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 + // + Point2D multiply(const Point2D& vec1, const Point2D& vec2) + { + Point2D 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; + } + + // 参数： vec1 矢量1 vec2 矢量2 + // + double Cos(const Point2D& vec1, const Point2D& vec2) + { + Point2D unit_vec1 = normalize(vec1); + Point2D unit_vec2 = normalize(vec2); + + return dotMultiply(unit_vec1, unit_vec2); + } + + // 1.2、矢量的长度 + // + // 参数： vec 矢量 + // + double length(const Point2D& vec) + { + return (sqrt(pow(vec.x, 2) + pow(vec.y, 2) + pow(vec.z, 2))); + } + + + + + // 计算两点之间的向量 + Point2D mathVector(const Point2D& p1, const Point2D& p2) { + Point2D pt; + pt.x = p2.x - p1.x; + pt.y = p2.y - p1.y; + return pt; + } + // 计算两点之间的叉积 + double crossProduct(const Point2D& p1, const Point2D& p2) { + return p1.x * p2.y - p1.y * p2.x; + } + + //直线与线段的交点 + bool findIntersection(double x1, double y1, double x2, double y2, // 直线1的两个点 + double x3, double y3, double x4, double y4,Point2D& intersection); + }; diff --git a/designsurveylinedlg.cpp b/designsurveylinedlg.cpp index b3a4955..1f4dc2b 100644 --- a/designsurveylinedlg.cpp +++ b/designsurveylinedlg.cpp @@ -48,8 +48,8 @@ BOOL DesignSurveyLineDlg::OnInitDialog() GetDlgItem( IDC_BTN_CALCULATELINE )->EnableWindow( false ); GetDlgItem( IDC_BTN_CALCULATETIME )->EnableWindow( false ); SetDlgItemText(IDC_EDIT_FLYSPEED,"10"); + SetDlgItemText(IDC_EDIT_EXTERNALLENGTH,"0"); - return TRUE; } @@ -76,6 +76,8 @@ CString DesignSurveyLineDlg::extractFileName(CString fileName) //设置航测区域 void DesignSurveyLineDlg::SetSurveyRegion(const vector& lons,const vector& lats) { + surveyRegionLons.clear(); + surveyRegionLats.clear(); for (int i=0;i& lons,const vecto } //计算测绘航线 -void DesignSurveyLineDlg::calculateSurveyLine(double lineInterval,vector&surveyLineLons,vector& surveyLineLats) +void DesignSurveyLineDlg::calculateSurveyLine(double lineInterval,vector&surveyLineLons,vector& surveyLineLats,double outLength) { vector recLons; vector recLats; TopologicalAnalysis analysis; GeoCompute geocompute; - analysis.GetBoundingBoxVertices(surveyRegionLons,surveyRegionLats,recLons,recLats); + //计算外包矩形 + //analysis.GetBoundingBoxVertices(surveyRegionLons,surveyRegionLats,recLons,recLats); + + int startPt1 = analysis.getLineWithPolygonMinWidth(surveyRegionLons,surveyRegionLats); + recLons.push_back(surveyRegionLons[startPt1]); + recLons.push_back(surveyRegionLons[startPt1+1]); + recLats.push_back(surveyRegionLats[startPt1]); + recLats.push_back(surveyRegionLats[startPt1+1]); + double flyAngle;//飞行航向 + CalculateTwoPtsAzimuth(flyAngle,surveyRegionLons[startPt1],surveyRegionLats[startPt1],surveyRegionLons[startPt1+1],surveyRegionLats[startPt1+1],3); + //计算偏移航向 + Line2D line; + line.p1.x = surveyRegionLons[startPt1]; + line.p1.y = surveyRegionLats[startPt1]; + line.p2.x = surveyRegionLons[startPt1 + 1]; + line.p2.y = surveyRegionLats[startPt1 + 1]; + Point2D p; //直线的下一个顶点 + if (startPt1+1 == (surveyRegionLons.size()-1)) //直线为多边形最后一条边 + { + p.x = surveyRegionLons[1]; + p.y = surveyRegionLats[1]; + } + else + { + p.x = surveyRegionLons[startPt1 + 2]; + p.y = surveyRegionLats[startPt1 + 2]; + } + int ptPos = analysis.pointPosition(p,line); + if (ptPos == 1) //左侧 + { + int i = (flyAngle + 270)/360; + flyAngle = (flyAngle + 270) - 360*i; + } + else if(ptPos == -1) //右侧 + { + int i = (flyAngle + 90)/360; + flyAngle = (flyAngle + 90) - 360*i; + } + double lon1,lat1,lon2,lat2; //第一条切割线为航线间隔一半 - geocompute.computeOffsetGeoPosition(recLons[0],recLats[0],180,lineInterval/2000,lon1,lat1); - geocompute.computeOffsetGeoPosition(recLons[1],recLats[1],180,lineInterval/2000,lon2,lat2); + geocompute.computeOffsetGeoPosition(recLons[0],recLats[0],flyAngle,lineInterval/2000,lon1,lat1); + geocompute.computeOffsetGeoPosition(recLons[1],recLats[1],flyAngle,lineInterval/2000,lon2,lat2); Point2D pt1,pt2; pt1.x = lon1; pt1.y = lat1; pt2.x = lon2; pt2.y = lat2; vector result; - analysis.linePolygonIntersections(pt1,pt2,surveyRegionLons,surveyRegionLats,result); + //交点个数 + int nCrossPt = analysis.linePolygonIntersections(pt1,pt2,surveyRegionLons,surveyRegionLats,result); - while(lat1>recLats[2]){ + while(nCrossPt>1){ //平移切割线 - geocompute.computeOffsetGeoPosition(lon1,lat1,180,lineInterval/1000,lon1,lat1); - geocompute.computeOffsetGeoPosition(lon2,lat2,180,lineInterval/1000,lon2,lat2); + geocompute.computeOffsetGeoPosition(lon1,lat1,flyAngle,lineInterval/1000,lon1,lat1); + geocompute.computeOffsetGeoPosition(lon2,lat2,flyAngle,lineInterval/1000,lon2,lat2); pt1.x = lon1; pt1.y = lat1; pt2.x = lon2; pt2.y = lat2; //计算交点 - analysis.linePolygonIntersections(pt1,pt2,surveyRegionLons,surveyRegionLats,result); + nCrossPt = analysis.linePolygonIntersections(pt1,pt2,surveyRegionLons,surveyRegionLats,result); + + } + + int nlinePts = result.size(); + if (nlinePts%2!=0) + { + nlinePts--; + } + //计算外扩航点 + if (nlinePts>1 && outLength>0) + { + double angle; + CalculateTwoPtsAzimuth(angle,result[0].x,result[0].y,result[1].x,result[1].y,3); + //计算外扩航向 + double outAngle1,outAngle2; + outAngle2 = angle; + int n = (angle + 180)/360; + outAngle1 = (angle + 180) - 360*n; + geocompute.computeOffsetGeoPosition(result[1].x,result[1].y,outAngle2,outLength/1000,lon2,lat2); + result[1].x = lon2; + result[1].y = lat2; + for (int i=2;i surveyLineLons; vector surveyLineLats; //计算航测航线 - calculateSurveyLine(lineInterval,surveyLineLons,surveyLineLats); + calculateSurveyLine(lineInterval,surveyLineLons,surveyLineLats,outLength); //保存航测航线 saveSurveyLine(surveyLineLons,surveyLineLats,height); //发送消息到主程序，进行标绘 diff --git a/designsurveylinedlg.h b/designsurveylinedlg.h index 44801cd..321d80f 100644 --- a/designsurveylinedlg.h +++ b/designsurveylinedlg.h @@ -28,7 +28,7 @@ public: enum { IDD = IDD_DLG_DESIGNSURVEYLINE }; private: - void calculateSurveyLine(double lineInterval,vector&surveyLineLons,vector& surveyLineLats);//计算测绘航线 + void calculateSurveyLine(double lineInterval,vector&surveyLineLons,vector& surveyLineLats,double outLength=0);//计算测绘航线 void saveSurveyLine(const vector&surveyLineLons,const vector& surveyLineLats,double height); //保存测绘航线 bool saveSurveyRegion(const vector&surveyRegionLons,const vector& surveyRegionLats); //保存航测区域 double calculateRouteLength(const vector&surveyLineLons,const vector& surveyLineLats); //计算航程 diff --git a/resource.h b/resource.h index 01239ba..b725614 100644 --- a/resource.h +++ b/resource.h @@ -166,9 +166,10 @@ #define IDC_EDIT_FLYSPEED 1056 #define IDB_BITMAP2 1057 #define IDC_EDIT_MARKER_NAME 1057 -#define IDC_TXT_ROUTELENGTH2 1057 #define IDC_TXT_FLYTIME 1057 #define IDC_TREE_MARKERS 1058 +#define IDC_EDIT_ 1058 +#define IDC_EDIT_EXTERNALLENGTH 1058 #define IDD_DIALOG_MARKER 1063 #define IDR_MENU3 1066 #define IDR_LIST_POPUP 1066 -- 2.37.1.windows.1