计算几何系列 —— 计算几何多边形求交&求并&求交并比以及一道难题（含AC代码）

已于 2023-12-28 13:44:51 修改
阅读量1.3k 收藏 22
点赞数 35
分类专栏： 计算几何 文章标签： 算法 几何学 c++ c语言 青少年编程
于 2023-12-18 13:00:09 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Zz_0913/article/details/135049800
版权
       本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念，使读者快速入门，故不追求难度和深度，仅起到抛砖引玉的作用。     

  好久没有写计算几何题目了~（悲），最近想找题目好好练一练，因此撰写这份博客🤔

        这道题目便是   2023 ICPC 国际大学生程序设计竞赛亚洲区域赛（南京站）的     Problem B，当时的比赛赛场上呢是无人解出，这道题的难度属于夺冠水平~

      首先我们来看一下这道题目的题意，简化后就是：

            1.给定一个矩形（OBB），要求找到一个坐标轴平行矩形（AABB） ；

            2.两个矩形的相交面积比上面积并就是所谓的交并比（IOU）；

            3.我们需要求出这个矩形使IOU最大。

           上述是对一些交并比求解的可以在草稿纸上模拟的几种情况（读者们可以自己试试），总结出来重合部分多边形的边数在4-8，我们在求解的时候只需要多边形类求交就行了~😄

   我们一般怎么求交并比呢🤔？

               这个就是我们之前所学习的基本计算几何内容了：

     多边形面积交&多边形面积并,由于思路较为基础，基础部分已经做介绍，这里我们直接上源码可做调用：

      1.我们一般采用vector存储多边形，点类  or 线类（不含极角排序）：

   

        2.多边形求交：

namespace PolyIntersectArea{
	//ConvexPolygonIntersectArea
	db CPIA(Polygon a,Polygon b)
	{
		Polygon p,tmp;
		p = b;
		for(int i = 0;i < a.sze()n&& b.sze() > 2;i++)
		{
			int sflag = sgn((a[i + 1] - a[i]) ^ (p[0] - a[i])),eflag;
			tmp = Polygon(0);
			for(int j = 0;j < b.sze();j++,sflag = eflag)
			{
				if(sflag >= 0) tmp.add(p[j]);
				eflag = sgn((a[i + 1] - a[i]) ^ (p[j + 1] - a[i]));
				if((sflag ^ eflag) == -2)
				    tmp.add(Line(a[i],a[i + 1]).crossPoint(Line(p[j],p[j + 1])));
			}
			p = tmp;
			b.p.resize(tmp.sze());
		}
		if(b.sze() < 3) return 0.0;
		return p.getArea();
	}
	//simplePolygonIntersectArea
	db SPIA(Polygon a,Polygon b)
	{
		Polygon t1(3),t2(3);
		db res = 0,num1,num2;
		t1[0] = a[0],t2[0] = b[0];
		for(int i = 2;i < a.sze();i++)
		{
			t1[1] = a[i - 1],t1[2] = a[i];
			num1 = sgn((t1[1] - t1[0]) ^ (t1[2] - t1[0]));
			if(num1 < 0) swap(t1[1],t1[2]);
			for(int j = 2;j < b.sze();j++)
			{
				t2[1] = b[j - 1],t2[2] = b[j];
				num2 = sgn((t2[1] - t2[0]) ^ (t2[2] - t2[0]));
				if(num2 < 0) swap(t2[1],t2[2]);
				res = res + CPIA(t1,t2) * num1 * num2;
			}
		}
		return res;
	}
}

         3.多边形求并：

namespace PolyUnion{
	Point dir(Line in)
	{
		return in.p - in.v;
	}
	db pos(Point p,Line in)
	{
		return ((p - in.p) * dir(in)) / dir(in).len2();
	}
	db gao(vector <Polygon> po)
	{
		int n = po.size();
		db res = 0;
		for(int i = 0;i < n;i++){
		for(int ii = 0;ii < po[i].sze();ii++)
		{
			Point A = po[i][ii],B = po[i][(ii + 1) % p[i].sze()];
			Line AB = Line(A,B);
			vector<pair<db , int>>c;
			for(int j = 0;j < n;j++)
			if(i != j)
			{
				for(int jj = 0;jj < po[j].sze();jj++)
				{
					Point C = po[j][jj],D = po[j][(jj + 1) % po[j].sze()];
					Line CD = Line(C,D);
					int f1 = sgn((B - A) ^ (C - A));
					int f2 = sgn((B - A) ^ (D - A));
					if(! f1 && ! f2)
					{
						if(i < j && sgn((dir(AB) * dir(CD))) > 0)
						{
							c.push_back(make_pair(pos(C,AB),1));
							c.push_back(make_pair(pos(D,AB),-1));
						}
						continue;
					}
					db s1 = (D - C) ^ (A - C);
					db s2 = (D - C) ^ (B - C);
					db t = s1 / (s1 - s2);
					if(f1 >= 0 && f2 < 0) c.push_back(make_pair(t,1));
					if(f1 < 0 && f2 >= 0) c.push_back(make_pair(t,-1));
				}
			}
			c.push_back(make_pair(0.0,0));
			c.push_back(make_pair(1.0,0));
			sort(c.begin(),c.end());
			db s = 0.5 * (A ^ B),z = min(max((db)c[0].second,(db)0.0),(db)1.0);
			for(int j = 1,k = c[0].second;j < (int)c.size();j++)
			{
				db w = min(max(c[j].first,(db)0.0),(db)1.0);
				if(k == 0) res = res + s * (w - z);
				k = k + c[j].second;
				z = w;
			}
		}
	  }
	  return fabs(res);
	}
}

以上就是多边形面积交并的分别求法，(⊙﹏⊙)其实我们可以明显发现，讲这个归讲这个，在这道题目中呢我们其实不需要这么麻烦，因为我们已经知道两个都是矩形了，到时候直接写一段大模拟求交点和多边形就OK了😄！

好我们来剖析一下这道题的类型：

    这是一道函数求极值的问题~

  我们一般对于这种问题会有一下考虑：

                    1. 是否单调/单峰/凸/ 考虑 ——> 二分？三分？
                    2. 是否可导/可微？=> 解方程（组）/梯度下降/牛顿迭代
                    3. 局部最小陷阱？=> 爬山/模拟退火
                    4. 有没有约束？=> 拉格朗日乘子法/KKT/线性规划...

   解题步骤以及相关技巧：

        1.利用对称性可以简化问题。答案显然只和 OBB 的长宽比以 及角度有关，我们对 OBB 做平移，镜像都不影响答案。
        2.不妨把这个 OBB 平移，使得它的中心位于原点。

        3.AABB 参数化为 (x, y, w, h)，其中心位于 (x, y), 上下左右边界分别为y + (h / 2),y - (h / 2),x - (w / 2),x + (w / 2)。

      所以有 IOU ( 0 , 0 , w , h )  >= IOU ( x , y , w , h ),AABB的中心放在原点，也就是与OBB的中心重合的位置我们视为最优点。

        4.设 g(w, h) = IOU(0, 0,w, h)，u(h) 是面积并关于 h 的导数， 那 g(w, h) 可看成关于 h 的函数表示如下:

       5.注意到 u(h) 是关于 h 单调不递减且满足 0 ≤ u(h) ≤ w, 这 里的分子也是单调不递增的，而分母恒大于 0，故 g 关于 h是单峰的。对称地，g 关于 w 也是单峰的。

       6.g(w, h) 整体是单峰函数。可以利用三分法，爬山法，梯度下降法等求解，都可以收敛到全局最优解。接下来介绍公式方法:

         为了方便讨论，可以通过镜像以及等比缩放等操作（不影响答案），使得 OBB 长边为 cos ϕ，短边 sin ϕ，长边和 x-轴的夹角为 θ。(0 <= θ,ϕ <= 45°）

         进行分类讨论，最优解时，交集可能是矩形（θ = 0），六边形（OBB 细长且倾斜度较大时）和八边形。我们也可以对每种情况列出函数，然后求解梯度等于 0 的方程来获得公式解：        解法雀实非常的明确且抽象,代码如下·,可以好好理解一下：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double LD;
typedef pair <int, int> pii;


#define cp const point &
const LD eps = 0;
int sgn (LD x) { return x > eps ? 1 : (x < -eps ? -1 : 0); }
LD sqr (LD x) { return x * x; }
struct point {
	LD x, y;
	point () {}
	point (LD xx, LD yy) { x = xx, y = yy; }
	point operator + (cp a) const { return {x + a.x, y + a.y}; }
	point operator - (cp a) const { return {x - a.x, y - a.y}; }
	point operator * (LD a) const { return {x * a, y * a}; }
	point operator / (LD a) const { return {x / a, y / a}; }
	point rot90() const { return {-y, x}; }
	bool operator == (cp a) const {
		return x == a.x && y == a.y;	
	}
	void read() {
		int xx, yy;
		cin >> xx >> yy;
		x = xx, y = yy;
	}
};

LD det (cp a, cp b) { return a.x * b.y - b.x * a.y; }
LD dot (cp a, cp b) { return a.x * b.x + a.y * b.y; }
LD dis (cp a, cp b) { return sqrt (sqr (a.x - b.x) + sqr(a.y - b.y)); }

bool turn_left(cp a, cp b, cp c) {
	return sgn (det (b - a, c - a)) > 0;
}

#define cl const line &
struct line {
	point s, t;
	line () {}
	line (point ss, point tt) { s = ss, t = tt; }
	bool operator == (cl a) const {
		return s == a.s && t == a.t;
	}
};

point line_inter (cl a, cl b) {
	LD s1 = det (a.t - a.s, b.s - a.s);
	LD s2 = det (a.t - a.s, b.t - a.s);
	return (b.s * s2 - b.t * s1) / (s2 - s1);
}
bool turn_left (cl l, cp p) { return turn_left(l.s, l.t, p); }

line h[8];
LD hpi_nosort() {
	line q[8]; int l = 0, r = -1;
	point ret[8];
	q[0] = h[0]; q[1] = h[1]; r = 1;
	ret[1] = line_inter(q[0], q[1]);
	for (int t = 2; t < 8; t++) {
		auto &i = h[t];
		while (l < r && !turn_left(i, ret[r]))
			-- r;
		while (l < r && !turn_left(i, ret[l + 1]))
			++ l;
		++ r; q[r] = i;
		if (l != r) ret[r] = line_inter(q[r - 1], q[r]);
	}
	ret[l] = line_inter(q[r], q[l]);
	LD area = 0;
	for (int i = l; i <= r; i++) area += det(ret[i], ret[i == r ? l : i + 1]);
	return area;
}


LD ans, sq;
LD check (LD x, LD y) {
	h[1] = {{-1, -y}, {1, -y}};
	h[3] = {{x, -1}, {x, 1}};
	h[5] = {{1, y}, {-1, y}};
	h[7] = {{0, 1}, {0, -1}};
	LD a = hpi_nosort();
	LD f = x * y * 4; 
	return a / (f + sq - a);
}

const LD Rp = (sqrt(5) - 1) / 2, Lp = 1 - Rp;
LD check (LD x) {
	LD l = 0, r = 1, lv = -1, rv = -1, v = 0;
	for (int t = 42; t; t--) {
		LD b = (r - l) * Rp;
		LD lmid = r - b, rmid = l + b;
		if (lv == -1) lv = check(x, lmid);
		if (rv == -1) rv = check(x, rmid);
		if (rv < lv) rv = lv, lv = -1, r = rmid;
		else lv = rv, rv = -1, l = lmid;
	}
	ans = max(ans, v = max(lv, rv));
	return v;
}


void work() {
	vector <point> p;
	for (int i = 0; i < 4; i++) {
		point u; u.read();
		p.push_back(u);
	}
	if (!turn_left(p[0], p[1], p[2])) {
		reverse(p.begin(), p.end());
	}
	if (sgn(det(p[1] - p[0], {1, 0})) == 0 || sgn(det(p[2] - p[1], {1, 0})) == 0) {
		cout << "1\n";
		return;
	}
	auto center = (p[0] + p[2]) / 2;
	{
		int k = 0;
		for (int i = 0; i < 4; i++) {
			if (p[i].x < p[k].x) k = i;
		}
		vector <point> q;
		for (int i = 0; i < 4; i++) q.push_back(p[(i + k) % 4] - center);
		p = move(q);
	}
	sq = abs(det(p[1] - p[0], p[2] - p[1]));
	LD Lx = p[2].x;
	LD Ly = p[3].y;
	for (auto &[x, y] : p) {
		x /= Lx;
		y /= Ly;
	}
	sq /= Lx * Ly;
	h[0] = {p[0], p[1]};
	h[2] = {p[1], p[2]};
	h[4] = {p[2], p[3]};
	h[6] = {p[3], p[0]};
	LD xl = 0, xr = 1, lv = -1, rv = -1;
	ans = check(1 - 1e-9);
	for (int t = 38; t; t--) {
		LD b = (xr - xl) * Rp;
		LD lmid = xr - b, rmid = xl + b;
		if (lv == -1) lv = check(lmid);
		if (rv == -1) rv = check(rmid);
		if (rv < lv) rv = lv, lv = -1, xr = rmid;
		else lv = rv, rv = -1, xl = lmid;
	}
	check((xl + xr) / 2);
	cout << (double)ans << '\n'; 
}

int main() {
	ios::sync_with_stdio(false); cin.tie(0);
	int T = 1;
	cin >> T;
	cout << fixed << setprecision(10);
	for (int ca = 1; ca <= T; ca ++) {
		work();
	}
}

          

         

Zz_0913
关注
  • 35
    点赞
  • 踩
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
C语言-计算几何常用算法源码——多边形基本算法
06-07
1. 判断多边形是否简单多边形 8 2. 检查多边形顶点的凸凹性 9 3. 判断多边形是否凸多边形 9 4. 求多边形面积 9 5. 判断多边形顶点的排列方向，方法一 10 6. 判断多边形顶点的排列方向，方法二 10 7. 射线法判断点是否在多边形内 10 8. 判断点是否在凸多边形内 11 9. 寻找点集的graham算法 12 10.寻找点集凸包的卷包裹法 13 11.判断线段是否在多边形内 14 12.求简单多边形的重心 15 13.求凸多边形的重心 17 14.求肯定在给定多边形内的一个点 17 15.求从多边形外一点出发到该多边形的切线 18 16.判断多边形的核是否存在 19
2023 ICPC 亚洲区域赛南京站计算几何 —— Problem B.交并比(源代码出炉+解析）
Zz_0913的博客
01-03 1175
2.两个矩形的相交面积比上面积并就是所谓的交并比（IOU）；1.给定一个矩形（OBB），要求找到一个坐标轴平行矩形（AABB）；3.我们需要求出这个矩形使IOU最大。博客👆中，我们通过非常复杂的代数方法证明关于IOU的这个函数是。三分法，爬山法，梯度下降法。
《玩不够的数学：算术与几何的妙趣》：第一章　平面上的几何艺术
热门推荐
GarfieldEr007的专栏
02-28 1万+
第一章　平面上的几何艺术 人们往往从悖论中获得思维的乐趣，而几何学的悖论就是不可能图形。如今我们已创造出数千种这样的二维图像，不断挑战我们的眼睛和思维。三角形、披萨饼、七巧板也蕴藏着无穷的变化和巧妙的发现。 不可能！你确信吗？ 人们从透视错觉得来灵感，创造了神秘的“不可能图形”。人类的视觉系统让我们觉得这样的图形很奇怪。然而这些图形确实是可行的，并为我们带来双重乐趣——先是惊奇，
【转载】ACM计算几何题目推荐
chenwenxiaocom的专栏
09-01 4493
一。基础题目  1.1 有固定算法的题目  A， 最近点对问题 最近点对问题的算法基于扫描线算法。 ZOJ     2107    Quoit Design    典型最近点对问题 POJ    3714    Raid    变种最近点对问题 B，最小包围
科学-几何：几何
weixin_30398227的博客
03-07 834
ylbtech-科学-几何：几何 几何，就是研究空间结构及性质的一门学科。它是数学中最基本的研究内容之一，与分析、代数等等具有同样重要的地位，并且关系极为密切。几何学发展历史悠长，内容丰富。它和代数、分析、数论等等关系极其密切。几何思想是数学中最重要的一类思想。暂时的数学各分支发展都有几何化趋向，即用几何观点及思想方法去探讨各数学理论。常见定理有勾股定理，欧拉定理，斯图尔...
计算几何线段求交 (包含软件使用说明书)
08-14
本程序支持绘制任意线段，支持捕捉线段端点(前提是光标距某一端点必须小于一定的阈值，程序设定为15个像素)，可以绘制任意多边形，求交点完毕后还可以在此基础上继续绘制再次求交。示例效果如下图：
计算几何–线段求交
01-07
线段求交 问题描述：已知两条线段P1P2P_1P_2P1​P2​和Q1Q2Q_1Q_2Q1​Q2​，判断P1P2P_1P_2P1​P2​和Q1Q2Q_1Q_2Q1​Q2​是否相交，若相交，求出交点。 两条线段的位置关系可以分为三类：有重合部分、无重合部分但有...
计算几何线段求交
03-13
线段求交针对的是平面空间中线段的求交，在所有线段都不是水平的情况下，任意的两条线段至多交于一点，且没有三条或多条线段交于一点。
（HDUACM201803版_14）计算几何基础 - 副本_不规则多边形_计算几何基础pptx_
09-30
ACM算法 计算几何基础 用于计算不规则多边形，凹多边形和凸多边形
【算法】二分算法——x的平方根
2302_79031646的博客
05-17 789
用二分边界算法来求解——x的平方根，从暴力求解到二分的优化细说
力扣62 不同路径 Java版本
最新发布
m0_47066863的博客
05-20 186
机器人每次只能向下或者向右移动一步。机器人试图达到网格的右下角（在下图中标记为 “Finish” ）。一个机器人位于一个 m x n 网格的左上角 （起始点在下图中标记为 “Start” ）。从左上角开始，总共有 3 条路径可以到达右下角。题目数据保证答案小于等于 2 * 109。输入：m = 3, n = 7。输入：m = 3, n = 2。输入：m = 7, n = 3。输入：m = 3, n = 3。问总共有多少条不同的路径？
MT2057 门票
l141930402的博客
05-15 761
前缀和+归并排序
【算法】【二叉树，DFS，哈希集合，分类讨论】力扣1110. 删点成林
未来就在脚下。
05-20 496
本题通过DFS遍历二叉树，结合分类讨论的方法，逐步删除指定节点并生成新的森林。该算法有效地处理了节点删除后树结构的调整问题，并通过哈希集合优化删除判断的时间复杂度，最终实现了高效的解决方案。
技艺高超的魔法师：Java运算符
epitomizelu的专栏
05-14 791
Java运算符是构建复杂程序的基础，掌握它们的使用和优先级对于写出高效、清晰的代码至关重要。通过本文的介绍和实例，希望能帮助读者更好地理解和运用这些运算符，让代码表达更加精准和高效。记住，实践是学习的最好方式，不断编写代码并探索不同的场景，你会越来越熟练地掌握它们。
力扣：209. 长度最小的子数组
icbbm的博客
05-17 388
如果不存在符合条件的子数组，返回。是该条件下的长度最小的子数组。找出该数组中满足其总和大于等于。个正整数的数组和一个正整数。
22. 括号生成
编程SHARE的博客
05-20 214
数字 n 代表生成括号的对数，请你设计一个函数，用于能够生成所有可能的并且 有效的 括号组合。
算法模块方法总结（数组）
Cobrander2_0的博客
05-17 152
二分查找（Binary Search）是一种在有序数组（必须是不能重复的，因为如果有重复的，会从几个符合的选项中输出一个，可能不合题意）中查找某一特定元素的搜索算法。搜索过程从数组的中间元素开始，如果中间元素正好是要查找的元素，则搜索过程结束；如果某一特定元素大于或者小于中间元素，则在数组大于或小于中间元素的那一半中查找，而且跟开始一样从中间元素开始比较。如果在某一步骤数组为空，则代表找不到。二分查找的思想是便于理解的，但是边界条件是很容易写错的。比如下面的。
常见加解密算法03 - RC4逆向认识
a5right的博客
05-14 962
各位聪明绝顶，才高八斗的读者们你们好！今天我们主要讨论编译之后的RC4算法识别。题外话，之前看到一个蛋疼的小知识，说“势”这个字最好不好查词典释义。我是很好奇的，果然后来无法直视势不可挡这个成语。言归正传，我们将上一篇的Java与C的代码都编译一下，分别反编译看看。
代码求两个几何多边形的交点坐标
05-27
求两个几何多边形的交点坐标可以用线段求交的思路来解决。具体步骤如下： 1. 对于每个多边形，遍历其所有边界线段，求出每条线段的参数方程。 2. 对于第一个多边形的每条线段，与第二个多边形的所有线段进行求交。 3. 对于每一组相交的线段，求出其交点坐标。 4. 将所有交点坐标存入列表中，并去除重复的点。 下面是一个简单的 Python 代码实现： ```python def intersection(poly1, poly2): # poly1, poly2 分别为两个多边形的顶点列表 def line_intersection(line1, line2): # 求两条线段的交点坐标 x1, y1, x2, y2 = line1 x3, y3, x4, y4 = line2 d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4) if d == 0: return None x = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d y = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d return (x, y) def segment_intersection(seg1, seg2): # 求两个线段是否相交，相交则返回交点坐标 line1 = (poly1[seg1[0]][0], poly1[seg1[0]][1], poly1[seg1[1]][0], poly1[seg1[1]][1]) line2 = (poly2[seg2[0]][0], poly2[seg2[0]][1], poly2[seg2[1]][0], poly2[seg2[1]][1]) intersection_point = line_intersection(line1, line2) if intersection_point is None: return None if (intersection_point[0] < min(poly1[seg1[0]][0], poly1[seg1[1]][0]) or intersection_point[0] > max(poly1[seg1[0]][0], poly1[seg1[1]][0]) or intersection_point[1] < min(poly1[seg1[0]][1], poly1[seg1[1]][1]) or intersection_point[1] > max(poly1[seg1[0]][1], poly1[seg1[1]][1]) or intersection_point[0] < min(poly2[seg2[0]][0], poly2[seg2[1]][0]) or intersection_point[0] > max(poly2[seg2[0]][0], poly2[seg2[1]][0]) or intersection_point[1] < min(poly2[seg2[0]][1], poly2[seg2[1]][1]) or intersection_point[1] > max(poly2[seg2[0]][1], poly2[seg2[1]][1])): return None return intersection_point intersections = [] for i in range(len(poly1)): for j in range(len(poly2)): intersection_point = segment_intersection((i, (i+1)%len(poly1)), (j, (j+1)%len(poly2))) if intersection_point is not None: intersections.append(intersection_point) # 去除重复的点 intersections = list(set(intersections)) return intersections ``` 其中，`poly1` 和 `poly2` 分别为两个多边形的顶点列表，每个顶点为一个二元组表示坐标。该函数返回两个多边形的交点坐标列表。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值