平面几何相关算法整理

自己在项目中用到的一些平面几何问题，主要都是参考其他博主的文章写的，会附上原文链接，本文仅作简要介绍和附上实现代码，具体原理可由参考链接中查看。

一、数据结构说明

二、判断线段是否相交

三、判断直线与线段是否相交

四、两线段交点

五、两直线交点

六、判断点在多边形内部

七、凸包排序

八、多边形顶点排序

九、判断多边形顶点是顺时针还是逆时针排序

十、多边形的有向面积

十一、点到线段的最短距离

十二、点关于直线的对称点

十三、曲线简化

一、数据结构说明

点类：vec2有一个成员变量为浮点型数组，v[0]表示x坐标，v[1]表示y坐标

class vec2
{
public:float v[2];vec2(){}vec2(float x,float y){v[0] = x; v[1] = y;}const float &operator [] (int i) const{return v[i];}float &operator [] (int i){return v[i];}
};

类外重载点类的部分运算符+、-、*以及求距离函数dist

static inline const vec2 operator + (const vec2 &v1, const vec2 &v2)
{vec2 result;for (int i = 0; i < 2; i++)result[i] = v1[i] + v2[i];return result;
}
static inline const vec2 operator * (const float &x, const vec2 &v)
{vec2 result;for (int i = 0; i < 2; i++)result[i] = x * v[i];return result;
}
static inline const vec2 operator - (const vec2 &v1, const vec2 &v2)
{vec2 result;for (int i = 0; i < 3; i++)result[i] = v1[i] - v2[i];return result;
}
static inline const float dist2(const vec2 &v1, const vec2 &v2)
{float d2 = (v2[0] - v1[0])*(v2[0] - v1[0]);for (int i = 1; i < 2; i++)d2 += (v2[i] - v1[i])*(v2[i] - v1[i]);return d2;
}
static inline const float dist(const vec2 &v1, const vec2 &v2)
{return sqrt(dist2(v1, v2));
}

二、判断线段是否相交

跨立互斥实验，将共线、重合，在端点处相交都视为相交。传入线段ab，cd

计算几何-判断线段是否相交 - 勿忘初心0924 - 博客园

//判断线段相交。包括共线、重合，端点处相交等特殊情况（跨立互斥实验）
static inline bool lineintersect(vec2 &a, vec2 &b, vec2 &c, vec2 &d)
{if (!(min(a[0], b[0]) <= max(c[0], d[0]) && min(c[1], d[1]) <= max(a[1], b[1]) && min(c[0], d[0]) <= max(a[0], b[0]) && min(a[1], b[1]) <= max(c[1], d[1])))return false;else{double u, v, w, z;//分别记录两个向量u = (c[0] - a[0])*(b[1] - a[1]) - (b[0] - a[0])*(c[1] - a[1]);v = (d[0] - a[0])*(b[1] - a[1]) - (b[0] - a[0])*(d[1] - a[1]);w = (a[0] - c[0])*(d[1] - c[1]) - (d[0] - c[0])*(a[1] - c[1]);z = (b[0] - c[0])*(d[1] - c[1]) - (d[0] - c[0])*(b[1] - c[1]);return (u*v <= 0.00000001 && w*z <= 0.00000001);}
}

三、判断直线与线段是否相交

传入线段2点，直线方向向量，直线上一点，直线的方向向量可由直线上2点相减得到。

static inline bool seg_line_is_intersect(vec2 segPt1, vec2 segPt2, vec2 dir, vec2 linep)
{bool result = false;double temp1 = dir[0] * (segPt1[1] - linep[1]) - dir[1] * (segPt1[0] - linep[0]);double temp2 = dir[0] * (segPt2[1] - linep[1]) - dir[1] * (segPt2[0] - linep[0]);if ((temp1 >= 0 & temp2 <= 0) | (temp1 <= 0 & temp2 >= 0))result = true;return result;
}

四、两线段交点

前提是线段已经相交，顺序传入线段line1,line2的点，返回交点坐标

线段与线段的交点_肘子zhouzi的博客-CSDN博客_线段与线段的交点

static inline vec2 Get_Intersect_Vertex(vec2 &line1p1, vec2 &line1p2, vec2 &line2p1, vec2 &line2p2)
{vec2 base = line2p2 - line2p1;vec2 hypo = line1p1 - line2p1;vec2 cad2 = line1p2 - line2p1;double S1 = fabs(base[0] * hypo[1] - hypo[0] * base[1]);double S2 = fabs(base[0] * cad2[1] - cad2[0] * base[1]);double M = dist(line2p2, line2p1);double d1 = S1 / M;   //线段1的两个端点到base向量（线段2）的距离double d2 = S2 / M;if (M != 0){double t = d1 / (d1 + d2);vec2 target;target[0] = line1p1[0] + (line1p2[0] - line1p1[0])*t;target[1] = line1p1[1] + (line1p2[1] - line1p1[1])*t;return target;}else  //用于检测bug，分母M为0出错{vec2 target(0, 0);return target;}
}

五、两直线交点

采用向量方法，直线表示为线上一点及方向向量，依次传入2条直线，返回交点

static inline vec2 intersection(vec2 &p, vec2 &v, vec2 &q, vec2 &w)
{//注意v,w是方向向量，p,q是直线上的一点//若已知直线1的方向向量为（a,b）,则在2维平面内，垂直于直线1的直线的方向向量为（b,-a）vec2 u = p - q;float a = w[0] * u[1] - w[1] * u[0];float b = v[0] * w[1] - v[1] * w[0];float t = a / b;return (p + t * v);
}

六、判断点在多边形内部

射线法，取一点，做该点的射线与多边形相交，看该射线与多边形的交点个数，如为奇数，点在多边形内部，否则在外部。

【计算几何】判断一个点是否在多边形内部_lazy-sheep的博客-CSDN博客_计算点在多边形内

static bool PointInPolygon(vec2 p, vector<vec2>&V)
{int   i, j = V.size() - 1;bool  oddNodes = false;for (i = 0; i < V.size(); i++){if ((V[i][1] < p[1] && V[j][1] >= p[1] || V[j][1] < p[1] && V[i][1] >= p[1]) && (V[i][0] <= p[0] || V[j][0] <= p[0])){if (V[i][0] + (p[1] - V[i][1]) / (V[j][1] - V[i][1])*(V[j][0] - V[i][0]) < p[0]){oddNodes = !oddNodes;}}j = i;}return oddNodes;
}

七、凸包排序

1.卷包裹法，输入PointSet,n为PointSet点数，len为输出凸包的点数，对一点集进行排序，寻找它的最大凸包点集合。

static inline vector<vec2> ConvexClosure(vector<vec2> &PointSet)
{int n = PointSet.size();vector<vec2>ch; ch.clear();int top = 0, i, index, first;double curmax, curcos, curdis;vec2 tmp,p1, p2, q1, q2;//声明线段p1p2和q1q2bool use[5000];//定义顶点数量最大值为5000tmp = PointSet[0];index = 0;// 选取y最小点，如果多于一个，则选取最左点for (i = 1; i < n; i++){if (PointSet[i][1] < tmp[1] || PointSet[i][1] == tmp[1]&&PointSet[i][0] < tmp[0]){index = i;}use[i] = false;}tmp = PointSet[index];first = index;use[index] = true;index = -1;ch[top++] = tmp;tmp[0] -= 100;p1 = tmp;p2 = ch[0];q1 = ch[0];while (index != first){curmax = -100;curdis = 0;// 选取与最后一条确定边夹角最小的点，即余弦值最大者for (i = 0; i<n; i++){if (use[i])continue;q2 = PointSet[i];// 根据cos值求夹角余弦，范围在（-1 -- 1 ）curcos = ((p2[0] - p1[0])*(q2[0] - q1[0]) + (p2[1] - p1[1])*(q2[1] - q1[1])) / (dist(p2, p1)*dist(q2, q1)); if (curcos>curmax || fabs(curcos - curmax)<1e-6 && dist(q1,q2)>curdis){curmax = curcos;index = i;curdis = dist(q1, q2);}}use[first] = false; //清空第first个顶点标志，使最后能形成封闭的hulluse[index] = true;ch[top++] = PointSet[index];p1 = ch[top - 2];p2 = ch[top - 1];q1 = ch[top - 1];}
}

八、多边形顶点排序

1.叉乘排序，针对凸多边形完成逆时针排序，两个向量叉乘得垂直于这两个向量所在平面的法向量，带有正负，根据正负可以判断两点的顺逆问题。

【计算几何】多边形点集排序 - 一点心青 - 博客园

static void vertexsort_1(vector<vec2>&pps, vector<int>&ppsindex)
{if (ppsindex.empty())return;for (int i = 0; i < ppsindex.size() - 1; i++){for (int j = i + 1; j < ppsindex.size(); j++){double x = (pps[ppsindex[j]][0] - pps[ppsindex[0]][0])*(pps[ppsindex[i]][1] - pps[ppsindex[0]][1]) -(pps[ppsindex[i]][0] - pps[ppsindex[0]][0])*(pps[ppsindex[j]][1] - pps[ppsindex[0]][1]);if (x >= 0.00000001){int tempp;tempp = ppsindex[j];ppsindex[j] = ppsindex[i];ppsindex[i] = tempp;}}}
}

2.按角度排序，先计算多边形的中点，即所有顶点的平均，以该点做向量起始点，各个点为向量终止点，计算该向量与x轴或y轴所成角度，根据角度排序。

//绕中点，按角度排序，中点O，顶点A，计算角AOx的大小，依次排序
static void vertexsort_2(vector<vec2>&pps, vector<int>&ppsindex, vec2 &cc)
{if (ppsindex.empty())return;vec2 B;B[0] = cc[0] + 10;B[1] = cc[1];for (int i = 0; i < ppsindex.size() - 1; i++){for (int j = 0; j < ppsindex.size() - i - 1; j++){if (Vvalule(pps[ppsindex[j]], B, cc) < Vvalule(pps[ppsindex[j + 1]], B, cc))  //0-180°{int tempp;tempp = ppsindex[j];ppsindex[j] = ppsindex[j + 1];ppsindex[j + 1] = tempp;}}}
}//vertexsort_2的子函数，求角P1cP2的sin值，只返回该值的正负，用于判断角是否大于180度
static bool sinv(const vec2 &P1, const vec2 &P2, const vec2 &c)
{double x = (P1[0] - c[0])*(P2[1] - c[1]) - (P2[0] - c[0])*(P1[1] - c[1]);if (x >= 0.00000001)return true;elsereturn false;
}
//vertexsort_2的子函数，计算角度大小，以弧度表示
static double Vvalule(const vec2 &P1, const vec2 &P2, const vec2 &c)
{double x = acos(cosv(P1, P2, c));double Pi = 3.1415926;if (sinv(P1, P2, c))  //0-180°{return x;}else{return 2 * Pi - x;}
}

九、判断多边形顶点是顺时针还是逆时针排序

鞋带公式，判断多边形是逆时针排序返回true,否则返回false

【Green公式】Hunter’s Apprentice（判断多边形为顺时针或逆时针）--鞋带公式_Charon_HN的博客-CSDN博客_格林公式顺时针和逆时针的区别

static bool Lacesformula(vector<vec2>&PS)
{double S = 0;for (int i = 0; i < PS.size(); i++){int e = (i + 1) % PS.size();S += -0.5*(PS[e][1] + PS[i][1])*(PS[e][0] - PS[i][0]);}if (S>0)return true;else return false;
}

十、多边形的有向面积

计算返回多边形有向面积，逆时针为正

static float DirectedArea(vector<vec2>&PS)
{double S = 0;for (int i = 0; i < PS.size(); i++){int e = (i + 1) % PS.size();S += -0.5*(PS[e][1] + PS[i][1])*(PS[e][0] - PS[i][0]);}return S;
}

十一、点到线段的最短距离

点pt到线段pq的距离,如果投影不在pq上，距离即为与端点的连线

static inline float DistancePointtoSegment(vec2 &pt, vec2 &p, vec2 &q)
{float pqx = q[0] - p[0];float pqy = q[1] - p[1];float dx = pt[0] - p[0];float dy = pt[1] - p[1];float d = pqx*pqx + pqy*pqy;float t = pqx*dx + pqy*dy;if (d > 0)t /= d;if (t < 0)t = 0;else if (t > 1)t = 1;// t = 0，计算 pt点 和 p点的距离; t = 1, 计算 pt点 和 q点 的距离; 否则计算 pt点 和 投影点 的距离。dx = p[0] + t*pqx - pt[0];dy = p[1] + t*pqy - pt[1];return dx*dx + dy*dy;
}

十二、点关于直线的对称点

static vec2 mirrorvec2(vec2 &p, vec2 &a, vec2 &b)
{vec2 mirp;//直线一般式float A = b[1] - a[1];float B = a[0] - b[0];float C = b[0] * a[1] - a[0] * b[1];mirp[0] = (B*B*p[0] - A*B*p[1] - A*C) / (A*A + B*B);mirp[1] = (-A*B*p[0] + A*A*p[1] - B*C) / (A*A + B*B);return interpolation(p, mirp, 2);
}

十三、曲线简化

道格拉斯算法，简化曲线，闭合曲线也可以做

曲线（笔迹）简化算法_Ironyho的博客-CSDN博客_曲线算法

//简化曲线-道格拉斯算法
static void CurveSimplification(const vector<vec2> &pointList, double epsilon, vector<vec2> &out)
{if (pointList.size() < 2)throw invalid_argument("Not enough points to simplify");// 找到起点和终点之间距离直线最大的点double dmax = 0.0;size_t index = 0;size_t end = pointList.size() - 1;for (size_t i = 1; i < end; i++){double d = DistPtoLine(pointList[i], pointList[0], pointList[end]);if (d > dmax){index = i;dmax = d;}}// 如果距离大于给定阈值epsilonif (dmax > epsilon){// 递归调用该函数vector<vec2> recResults1;vector<vec2> recResults2;vector<vec2> firstLine(pointList.begin(), pointList.begin() + index + 1);vector<vec2> lastLine(pointList.begin() + index, pointList.end());CurveSimplification(firstLine, epsilon, recResults1);CurveSimplification(lastLine, epsilon, recResults2);// 建立结果列表out.assign(recResults1.begin(), recResults1.end() - 1);out.insert(out.end(), recResults2.begin(), recResults2.end());if (out.size() < 2)throw runtime_error("Problem assembling output");}else{//Just return start and end pointsout.clear();out.push_back(pointList[0]);out.push_back(pointList[end]);}
}//闭合曲线简化
static vector<vec2> CurveSimplification(const vector<vec2> &pointList, double epsilon)
{//1.先对闭合曲线分成2段vector<vec2>pointlist1; pointlist1.clear();vector<vec2>pointlist2; pointlist2.clear();int indexi, indexj;//最远2点的位置下标double dmax = 0.0f;for (int i = 0; i < pointList.size(); i++){for (int j = i + 1; j < pointList.size() - 1; j++){double d = dist(pointList[i], pointList[j]);if (d > dmax){indexi = i;indexj = j;dmax = d;}}}pointlist1.insert(pointlist1.begin(), pointList.begin() + indexj, pointList.end());//indexj到end的元素赋予list1pointlist1.insert(pointlist1.end(), pointList.begin(), pointList.begin() + indexi + 1);//begin到indexi的元素赋予list1pointlist2.insert(pointlist2.begin(), pointList.begin() + indexi, pointList.begin() + indexj + 1);//indexi到indexj的元素赋予list2//2.分别对两段曲线进行道格拉斯压缩点集vector<vec2>out1; out1.clear();vector<vec2>out2; out2.clear();CurveSimplification(pointlist1, epsilon, out1);CurveSimplification(pointlist2, epsilon, out2);//3.融合2个压缩点集out1.insert(out1.end(), out2.begin() + 1, out2.end() - 1);//去掉首尾点，这2点是重复的？return out1;
}