几何相关计算

喵喵不爱吃鱼

已于 2024-07-17 19:32:08 修改

阅读量154

点赞数 3

分类专栏： GIS 文章标签：几何学 c++

于 2024-07-15 20:10:37 首次发布

本文链接：https://blog.csdn.net/weixin_42391513/article/details/140447877

版权

GIS 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、判断两个矩形是否相交

当矩形1的最大值比矩形2的最小值都小，那矩形1和矩形2一定不相交，其他同理。

struct Point {
	double x;
	double y;
};
struct Rec {
	Point min;
	Point max;
};
bool is_rec_overlap(const Rec& r1, const Rec& r2) {
	if(r1.max.x < r2.min.x
		|| r1.max.y < r2.min.y
		|| r2.max.x < r1.min.x
		|| r2.max.y < r1.min.y) {
		return false;
	} else {
		return true;
	}	
}

二、判断两条线段是否相交

排斥实验+跨立实验

#include<iostream>
#include<vector>

struct Point {
	double x;
	double y;
	Point(double _x, double _y) : x(_x),y(_y) {}
};

struct Edge {
    Point s;
    Point e;
    Edge(const Point& _s, const Point& _e) : s(_s),e(_e) {}
};

// 判断点在直线的左右侧：计算叉积(p_e_p_s) × (p-p_s)
// 若返回结果>0，则点p在线段ps_p_e的左侧；
// 若返回结果<0，则点p在线段ps_p_e的右侧；
// 若返回结果=0，则点p与ps、p_e共线；
double is_left(const Point& p, const Point& p_s, const Point& p_e) {
    return (p_e.x - p_s.x)*(p.y - p_s.y) - (p_e.y - p_s.y)*(p.x - p_s.x);
}

bool edges_intersect(const Edge& edge1, const Edge& edge2) {
    // 排斥实验：最小外包框不想交，edges不可能相交
    if(std::min(edge1.s.x, edge1.e.x) > std::max(edge2.s.x, edge2.e.x)
        || std::min(edge1.s.y, edge1.e.y) > std::max(edge2.s.y, edge2.e.y)
        || std::min(edge2.s.x, edge2.e.x) > std::max(edge1.s.x, edge1.e.x)
        || std::min(edge2.s.y, edge2.e.y) > std::max(edge1.s.y, edge1.e.y)) {
        return false;
    }

    // 跨立实验：edge1的两个端点在edge2的两侧，且edge2的两个端点在edge1的两侧，则edges相交。
    if(is_left(edge1.s, edge2.s, edge2.e) * is_left(edge1.e, edge2.s, edge2.e) <= 0
        && is_left(edge2.s, edge1.s, edge1.e) * is_left(edge2.e, edge1.s, edge1.e) <= 0) {
        return true;
    }

    return false;
}

参考：https://www.cnblogs.com/TangMoon/p/7611115.html

三、判断点是否在多边形内

简单版

// 判断点在多边形内部还是外部（射线法）：
// 从待判断的点出发，向某一方向（这里假设水平向右）发射一条射线;
// 遍历多边形的每一条边，检查射线是否相交
// 统计交点数量，若为偶数，则点在外部，否则，点在内部。
// 这里交点的判断简化为：edge的两个端点y值一个大于pt.y，一个小于pt.y，同时统计pt所在水平直线与edge交点的x值小于pt.x的个数。
bool is_in_polygon(const Point& pt, const std::vector<Point>& polygon) {
	int count = polygon.size();
	if(count < 3) {
		return false;
	}

    int intersecting_count = 0;
	for(int i = 0; i < polygon.size(); i++) {
		Point edge_s = polygon.at(i);
		Point edge_e = polygon.at((i+1)%count);
        if((pt.y <= edge_s.y && pt.y > edge_e.y 
            || pt.y <= edge_e.y && pt.y > edge_s.y)) {
            // 外层的if保证了edge_e.y != edge_s.y
            double pt_on_edge_x = (pt.y - edge_s.y) * (edge_e.x - edge_s.x) / (edge_e.y - edge_s.y) + edge_s.x;
            if(pt_on_edge_x < pt.x) {
                intersecting_count++;
            }
        }
	}
    return intersecting_count % 2 == 1;
}

参考：https://www.cnblogs.com/tgycoder/p/4901600.html

四、垂足计算

double cal_distance_square(const Point& pt1, const Point& pt2) {
    return (pt2.x - pt1.x) * (pt2.x - pt1.x) + (pt2.y - pt1.y) * (pt2.y - pt1.y);
}

// 计算(p_e-p_s)和(p-p_s)的点积
double dot_product(const Point& p, const Point& p_s, const Point& p_e) {
    return (p_e.x - p_s.x)*(p.x - p_s.x) + (p_e.y - p_s.y)*(p.y - p_s.y);
}

// 判断点pt到edge的投影点project_pt与edge的关系
/*  r = AP dot AB / ||AB||^2
    r=0 : project_pt = A
    r=1 : project_pt = B
    r<0 : project_pt is on backward extension of AB
    r>1 : project_pt is on forward extension of AB
   0<r<1: project_pt is on AB
*/
double relation(const Point& pt, const Edge& edge) {
    double len_square = cal_distance_square(edge.s, edge.e);
    // 特殊情况处理
    if(len_square < EP) {
        if(cal_distance_square(pt, edge.s) < EP) {
            return 0;
        } else {
            return -1; // 随便给个<0或>1的值
        }
    }
    return dot_product(pt, edge.s, edge.e) / len_square;
}

// 计算点pt到线段edge投影点.
// 若投影点在线段上，则project_pt赋值投影点，返回true
// 若投影点不在线段上，则project_pt赋值最近点，返回false。
bool calc_project_pt(const Point& pt, const Edge& edge, Point& project_pt) {
    double r = relation(pt, edge);
    if(r < 0) {
        project_pt = edge.s;
        return false;
    } else if (r > 1) {
        project_pt = edge.e;
        return false;
    }
    project_pt.x = edge.s.x + r * (edge.e.x - edge.s.x);
    project_pt.y = edge.s.y + r * (edge.e.y - edge.s.y);
    return true;
}

五、贝塞尔曲线

n阶贝塞尔曲线： $B(t)=\sum_{i=0}^{n}\binom{n}{i}P_i(1-t)^{n-i}t^i$

n阶贝塞尔曲线递归表达： $P_0^n=(1-t)P_0^{n-1}+tP_1^{n-1},t\in [0,1]$

二次贝塞尔曲线： $B(t)=(1-t)^2P_0+2t(1-t)P_1+t^2P_2,t\in [0,1]$

三次贝塞尔曲线： $B(t)=(1-t)^3P_0+3t(1-t)^2P_1+3t^2(1-t)P_2+t^3P_3,t\in [0,1]$

$\binom{n}{i}$ 表示从n个数中取i个数构成一个组合有多少种取法。

贝塞尔曲线（Bezier Curve）原理、公式推导及matlab代码实现-CSDN博客l

六、判断多边形顺时针还是逆时针

辛普森面积法

// 向量v1与向量v2的叉乘:
// 三维空间中为法向量，该向量垂直于向量v1与v2构成的平面。
// 二维空间中，有公式：v1×v2 = |v1||v2|sin(θ)，该角度有正负，是指从v1到v2的角度
double cross_product(const Point& v1, const Point& v2) {
    return v1.x * v2.y - v1.y * v2.x;
}

// 辛普森面积法判断多边形坐标点是顺时针还是逆时针存储
// 面积<0，则为顺时针；面积>0，则为逆时针
bool is_clockwise_polygon(const std::vector<Point>& polygon) {
    int n = polygon.size();
    if(polygon.size() < 3) {
        return false;
    }
    double area = 0.0;
    for(int i = 1; i < n-1; i++) {
        Point v1 = Point(polygon.at(i).x - polygon.at((i-1+n)%n).x,
            polygon.at(i).y - polygon.at((i-1+n)%n).y);
        Point v2 = Point(polygon.at((i+1)%n).x - polygon.at(i).x,
            polygon.at((i+1)%n).y - polygon.at(i).y);
        area += cross_product(v1, v2) / 2;
        std::cout << "i: " << i << ",area: " << area << std::endl;
        
    }
    return area < 0;
}

七、判断凹多边形

（1）角度法：判断每个顶点的角度是否大于180度，若存在这样的顶点，则多边形为凹多边形。

（2）叉乘法：假设多边形形点为逆时针存储，依次判断 $P_{i-1}P_i,P_iP_{i+1}$ 叉乘结果是否存在负值，若存在则为凹多边形。

// 向量v1与向量v2的叉乘:
// 三维空间中为法向量，该向量垂直于向量v1与v2构成的平面。
// 二维空间中，有公式：v1×v2 = |v1||v2|sin(θ)，该角度有正负，是指从v1到v2的角度
double cross_product(const Point& v1, const Point& v2) {
    return v1.x * v2.y - v1.y * v2.x;
}

// 向量法判断多边形凹凸性。
// 以下代码要求：polygon的坐标点须为逆时针存储。
// 若坐标点为顺时针存储，代码中判断条件需要改为：cp>0。
bool is_convex_polygon1(const std::vector<Point>& polygon) {
    int n = polygon.size();
    if(polygon.size() < 3) {
        return false;
    }
    for(int i = 0; i < n; i++) {
        Point v1 = Point(polygon.at(i).x - polygon.at((i-1+n)%n).x,
            polygon.at(i).y - polygon.at((i-1+n)%n).y);
        Point v2 = Point(polygon.at((i+1)%n).x - polygon.at(i).x,
            polygon.at((i+1)%n).y - polygon.at(i).y);
        double cp = cross_product(v1, v2);
        std::cout << "i: " << i << ",cp: " << cp << std::endl;
        if(cp < 0 ) {
            return false;
        }
    }
    return true;
}