polygon - 演算法筆記

half-plane

一條直線把二維平面劃分為兩半，其中一半就是「半平面」。半平面可以包含直線，也可以不包含直線。

半平面的一些圖示方式：

實作時，通常以兩個點來記錄半平面的直線、以兩個點的順序來記錄半平面的方向。

half-plane intersection

「半平面交集」就是許多個半平面的交集區域。半平面交集的結果可能是：一個開放區域、一個凸多邊形、一條線、一條線段、一個點、空集合。

演算法（incremental method）

一、預先使用四個半平面，設定一個極大的正方形邊界，讓半平面交集擁有邊界。
二、逐一加入每個半平面，求出當下的半平面交集（凸多邊形）。

online演算法，隨時維護一個半平面交集。每次更新需時O(N)，總時間複雜度O(N²)，N是半平面數目。

typedef complex<float> Point;
typedef vector<Point> Polygon;
typedef pair<Point,Point> Line;
#define x real()
#define y imag()

// 兩向量叉積
float cross(Point a, Point b)
{
	return a.x * b.y - a.y * b.x;
}

// 向量oa與向量ob進行叉積
float cross(Point o, Point a, Point b)
{
	return (a.x-o.x) * (b.y-o.y) - (a.y-o.y) * (b.x-o.x);
}

// 多邊形面積
float area(Polygon p)
{
	float a = 0;
	int n = p.size();
	for (int i=0; i<n; ++i)
		a += cross(p[i], p[(i+1)%n]);
	return fabs(a) / 2;
}

// 兩線交點
Point intersection(Point a1, Point a2, Point b1, Point b2)
{
	Point a = a2 - a1, b = b2 - b1, s = b1 - a1;
	return a1 + a * cross(b, s) / cross(b, a);
}

// 一個凸多邊形與一個半平面的交集
Polygon halfplane_intersection(Polygon p, Line line)
{
	Polygon q;
	Point p1 = line.first, p2 = line.second;

// 依序窮舉凸多邊形所有點，判斷是否在半平面上。
	// 如果凸多邊形與半平面分界線有相交，就求交點。
	int n = p.size();
	for (int i=0; i<n; ++i)
	{
		float c = cross(p1, p2, p[i]);
		float d = cross(p1, p2, p[(i+1)%n]);
		if (c >= 0) q.push_back(p[i]);
		if (c * d < 0) q.push_back(intersection(p1, p2, p[i], p[(i+1)%n]));
	}
	return q;
}

Line lines[10];	// 半平面

void demo()
{
	// 預先設定一個極大的正方形邊界
	Polygon p;
	p.push_back((Point){-1e9,-1e9});
	p.push_back((Point){-1e9,+1e9});
	p.push_back((Point){+1e9,-1e9});
	p.push_back((Point){+1e9,+1e9});

// 每一個半平面都與目前的半平面交集求交集
	for (int i=0; i<10; ++i)
	{
		p = halfplane_intersection(p, lines[i]);
		if (area(p) == 0) break;	// 退化或者空集合
	}
}

UVa 10084 10117 11265

演算法（randomized incremental method）

預先排序、預先隨機排列，最差、平均時間複雜度O(NlogN)。過程如同「convex hull: incremental method」，求切點改為求交點。

演算法（divide-and-conquer method）

時間複雜度O(NlogN)，N是半平面數目。

divide：把半平面分成兩堆。
conquer：分別遞迴求解。
combine：求兩個凸多邊形的交集。O(N)

兩個凸多邊形的交集，可以用旋轉卡尺求解，也可以用掃描線求解。

演算法【尚無正式名稱】

2006年由競賽選手南京外国语学校朱泽园《半平面交的新算法及其实用价值》提出。我不清楚是否已有正式學術論文，也不確定演算法是否正確。

半平面交集是一個凸多邊形，凸多邊形的邊很有順序的沿著外圍繞行一圈，越來越傾斜。援引Graham's scan的精神，所有半平面，預先按照角度（不是斜率）排序，逐步找出凸多邊形的邊。

一、所有半平面，按照角度排序。O(NlogN)
　　角度相同的半平面只需留下最內側的一個。
二、建立一個deque，加入前面兩個半平面。O(1)
三、從第三個半平面開始，依序將半平面加入deque。O(N)
　甲、deque右端持續彈出，直到deque右端的兩個半平面的交點，位於此半平面內。
　乙、deque左端持續彈出，直到deque左端的兩個半平面的交點，位於此半平面內。
　丙、deque右端加入此半平面。
四、刪除deque兩端多餘的半平面。O(N)
　甲、deque右端持續彈出，直到deque右端的兩個半平面的交點，位於deque左端的半平面內。

時間複雜度O(NlogN)，主要取決於排序的時間。N是半平面數目。

有點類似簡單多邊形的凸包演算法Melkman's algorithm，但是兩者並非點線對偶。

typedef complex<float> Point;
#define x real()
#define y imag()
typedef vector<Point> Polygon;
typedef pair<Point,Point> Line;
#define p1 first
#define p2 second

float cross(Point a, Point b)
{
	return a.x * b.y - a.y * b.x;
}

float cross(Point o, Point a, Point b)
{
	return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

float cross(Line l, Point p)
{
	return cross(l.p1, l.p2, p);
}

bool parallel(Line l1, Line l2)
{
	return cross(l1.p2 - l1.p1, l2.p2 - l2.p1) == 0;
}

Point intersection(Line l1, Line l2)
{
	Point& a1 = l1.p1, &a2 = l1.p2;
	Point& b1 = l2.p1, &b2 = l2.p2;
	Point a = a2 - a1, b = b2 - b1, s = b1 - a1;
	return a1 + a * cross(b, s) / cross(b, a);
}

bool cmp(Line l1, Line l2)
{
	return arg(l1.p2 - l1.p1) < arg(l2.p2 - l2.p1);
}

Polygon halfplane_intersection(vector<Line>& hp)
{
	// 所有半平面，按照角度排序。
	sort(hp.begin(), hp.end(), cmp);

// deque
	int L = 0, R = 0;
	vector<Line> l(hp.size());	// 凸多邊形的邊
	vector<Point> p(hp.size());	// 凸多邊形的頂點
	l[R] = hp[0];

for (int i=1; i<hp.size(); i++)
	{
		// deque右端持續彈出（若不允許退化，則添上等於）
		while (L < R && cross(hp[i], p[R-1]) < 0) R--;
		// deque左端持續彈出（若不允許退化，則添上等於）
		while (L < R && cross(hp[i], p[L])   < 0) L++;

// 如果已經大於等於180度，
		// deque為空，可以提早結束演算法。
		// 不好實作，故省略。

// 加入當前半平面
		l[++R] = hp[i];

// 兩個平行的半平面，刪除靠外側的半平面。
		// 預先加入四個半平面作為邊界，
		// 此處就不會接連出現斜率相等、方向相反的半平面。
		if (parallel(l[R-1], hp[i]) &&
			cross(l[--R], hp[i].p1) < 0) l[R] = hp[i];

// 計算頂點
		if (L < R) p[R-1] = intersection(l[R], l[R-1]);
	}
	// deque右端持續彈出（避免L型絲帶）（若不允許退化，則添上等於）
	while (L < R && cross(l[L], p[R-1]) < 0) R--;
	// 空集合
	if (R-L <= 1) return Polygon();
	// 計算頂點
	if (L < R) p[R] = intersection(l[L], l[R]);

// 重新儲存凸包頂點
	Polygon ch;
	for (int i=L; i<=R; i++) ch.push_back(p[i]);
	// 刪除重複的點。可能有多個半平面，交於同一點。
	ch.resize(unique(ch.begin(), ch.end()) - ch.begin());
	// 刪除頭尾重複的點。
	if (ch.size() > 1 && ch.front() == ch.back())
		ch.pop_back();
	return ch;	// 包含退化的情況
}

polygon intersection

兩個凸多邊形的交集

有三種演算法，時間複雜度都是O(N+M)。

半平面交集：半平面按照角度排序，採用merge sort。【尚待確認】

掃描線：凸多邊形從左右極點切開成上下兩條鏈。四鏈齊進。

交互前進：外側邊超車，直到擋住內側邊前方、或者撞到內側邊；內側邊超車，直到前方看不到外側邊、或者撞到外側邊。任選起點，先跑一圈找到任意交點；交點為起點，再跑一圈圍出交集。

UVa 137 10321

兩個簡單多邊形的交集、聯集、差集

掃描線。「segment intersection」略作修改。時間複雜度O(NlogN + K)，N是所有頂點數量，K是交點數量。

UVa 805 ICPC 7583

polygon clipping

用直線（半平面）裁切凸多邊形

凸多邊形劃一道直線，切割成兩個凸多邊形。

直線改成半平面，裁切成一個凸多邊形。

裁切前後，頂點先後順序一樣！可以循序掃描！

以點為主角、邊為配角，沿凸多邊形走一圈，點在半平面上則保存，邊與直線有交點則保存。時間複雜度O(N)。

UVa 11989

用直線（半平面）裁切簡單多邊形

我沒有研究。

UVa 10867 10974

用凸多邊形裁切簡單多邊形

polygon offsetting

簡單多邊形，所有邊一齊往垂直方向移動。往外則擴張、往內則收縮。

收縮時，頂點不斷合併，直到剩下一點。夾角決定速率，速率決定頂點合併順序。建立排序資料結構，時間複雜度O(NlogN)。

polygon offsetting

甲以各種姿態碰觸乙的邊緣，甲沿著乙的邊緣掃一圈，盡可能擴張（收縮）乙。結果包含了直線和弧線，厚度均勻。

甲求凸包、求直徑（最遠點對），甲可以簡化成線段。時間複雜度O(N)。

Minkowski sum / Minkowski difference

兩個簡單多邊形（內部所有點）的「pairwise sum」。

甲不偏不倚。甲內部任取一個基準點，甲的基準點沿著乙的邊緣走一圈，補厚（削薄）乙。基準點通常是頂點。

凸多邊形，有著高速演算法，時間複雜度O(N+M)。

簡單多邊形，切割成大量凸多邊形（三角剖分、凸分割），分頭計算，最後再聯集。時間複雜度O(N²M²)。

visibility of polygon

一個點的可見區域，被一個凸多邊形遮擋

點在凸多邊形內部。一定可以看見整個凸多邊形。回顧一下凸多邊形的定義：多邊形內部所有兩點連線，都在多邊形內部。令其中一點是給定點，令另外一點是可見點，兩點之間一定不會被凸多邊形遮擋。

點在凸多邊形上面。如果歸類於內部，可見區域就是整個凸多邊形；如果歸類於外部，可見區域的邊界，就是該點碰觸的邊的延長線。

點在凸多邊形外部。上下兩條切線圍住的區域受到影響，可見區域變成切線與凹面圍住的區域。找切線很簡單：建立點到凸多邊形各頂點的向量，運用叉積找出轉至最右、轉至最左的向量。

時間複雜度O(N)。

UVa 746

一個點的可見區域，被一個簡單多邊形遮擋

點在多邊形內部。被一個簡單多邊形遮擋的可見區域，也是簡單多邊形。採用incremental method，從一個可見點開始，依照逆時針順序，一次讀入一條邊，隨時更新當下的可見區域。一共分成六種情況：

一、可見區域直接增加一條邊。
二、可見區域增加一條割線。
　　割線會撞到多邊形以後的邊，屆時就能求得割點。
三、可見區域切除看不見的部分，改為一條割線。
　　割線會撞到可見區域以前的邊，必須倒退找割點。
四、同三。
五、同一。
六、同二。

倒退找割點，頂多倒退N條邊。時間複雜度O(N)。

附帶一提，如果預先找到兩個以上可見點，就可以切斷成鏈，實施divide-and-conquer method。

點在多邊形上面、外部。原理相同，不再贅述。

【待補程式碼】

UVa 10907

一個點的可見區域，被數個凸多邊形遮擋

點在凸多邊形外部、凸多邊形們都不相交，才有討論價值。

旋轉的掃描線。以給定點為基準點，所有頂點按照角度排序。依序掃描每個點，當下掃中的點，如果被前一點的鄰邊遮擋，就馬上刪除此點。時間複雜度O(NlogN)。

可以加速到O(N+HlogH)。牽涉到直線與凸多邊形交點的演算法。

一個點的可見區域，被數個簡單多邊形遮擋

旋轉的掃描線。找到給定點到每個多邊形的兩條切線，以及到切點的距離。所有切線按照角度排序。依序掃描每條切線，用二元樹動態維護每個多邊形的可見順序。時間複雜度O(NlogH)，N是所有頂點的數目，H是所有多邊形的數目。

比較容易實作的方式。找到給定點到各頂點的射線，以及到頂點的距離。所有射線按照角度排序。依序處理每條射線，用二元樹動態維護每條邊的可見順序。時間複雜度O(NlogN)。

visibility graph

大量的簡單多邊形頂點連線，連線不與障礙物交錯。亦可添加起點、終點、中繼點。主要用途是尋找最短路徑，繞過障礙物。

有時候主角不是點，而是多邊形。當主角只位移、不旋轉，則以Minkowski sum補厚障礙物，並且相對地削薄主角、變成點。

UVa 10762 10376 ICPC 3306

kernel of polygon

凸多邊形的核

位於多邊形內部，可以看到整個多邊形內部的地方，稱作「核」。「核」是設置監視攝影機、地標、廣告看板的好地方。「核」可能是風水最旺的地方。

一個多邊形的核，即是所有邊的半平面交集。一個多邊形的核，可能是一個凸多邊形、一條線段、一個點、空集合。

凸多邊形的核，顯然就是原本的凸多邊形。

簡單多邊形的核

直接套用半平面交集演算法，時間複雜度O(NlogN)。借用visibility of polygon演算法的手法，時間複雜度得壓至O(N)。

ICPC 2512 UVa 588