點積與叉積

電腦實施運算，通常會有浮點數誤差。為了避免浮點數誤差，當使用電腦計算幾何問題，會採用不同於一般的數學公式和定理。

點積（dot product）、叉積（cross product）這兩個運算只有加減法和乘法，而不包括除法，能夠有效避免除法產生的浮點數誤差，另一方面也能夠減少計算時間。點積與叉積有著許多好用的特性，大部分的幾何問題，都可以運用點積與叉積來計算答案。

以下都是用二維空間當作範例。

資料結構

點積與叉積是向量運算，所以先設計一個向量的資料結構。

struct Vector {int x, y;}; // 二維向量的資料結構 // 點積運算 int dot(Vector v1, Vector v2) { // 沒有除法，儘量避免誤差。 return v1.x * v2.x + v1.y * v2.y; } // 叉積運算，回傳純量（除去方向） int cross(Vector v1, Vector v2) { // 沒有除法，儘量避免誤差。 return v1.x * v2.y - v1.y * v2.x; }

向量資料結構擁有一個座標，並且擁有一支點積函式與一支叉積函式。

兩個向量做點積的結果是一個純量。兩個向量做叉積的結果為一個向量。二維向量的情況下，叉積的結果只有第三個數值不是零。我們只會用到第三個數值，所以讓叉積函式的回傳值為純量。

點積、叉積跟長度的關係

點積的結果為垂直投影的某種量。這某種量取絕對值，再除以底向量長度，得到底。

叉積的結果為平行四邊形的面積量。面積量取絕對值，再除以底向量長度，得到高。

struct Point {float x, y;}; // 點的資料結構 typedef Point Vector; // 向量的資料結構，和點一樣 // 向量的長度 float length(Vector v) { return sqrt(v1.x * v1.x + v2.y * v2.y); // return sqrt(dot(v, v)); } void base_height(Point p, Point p1, Point p2) { Vector v1 = p1 - p; Vector v2 = p2 - p; float base = fabs(dot(v1, v2)) / length(v1); float height = fabs(cross(v1, v2)) / length(v1); }

點積、叉積跟角度的關係

void sin_cos_θ(Point p, Point p1, Point p2) { Point p, p1, p2; Vector v1 = p1 - p; Vector v2 = p2 - p; float l1 = length(v1); float l2 = length(v2); float cosθ = dot(v1, v2) / l1 / l2; float sinθ = cross(v1, v2) / l1 / l2; float θ = acos(cosθ); // [0, π] float θ = asin(sinθ); // [-π/2, π/2] }

注意到acos與asin的回傳值，回傳的結果是弳度量（radian）而非度度量（grade），而且回傳值的範圍也不同。一般都以點積與acos求得介於0˚到180˚之間的夾角大小。

點積與向量夾角

利用點積的性質，可以粗略判斷夾角大小：點積大於0時，兩向量夾角小於90˚；等於0時，夾角等於90˚；小於零時，夾角大於90˚且小於180˚。

// 向量oa與向量ob進行點積，判斷∠aob之大小。 float dot(Point o, Point a, Point b) { return (a.x-o.x) * (b.x-o.x) + (a.y-o.y) * (b.y-o.y); }

叉積與向量旋轉

利用叉積的性質，可以粗略判斷夾角方向：叉積大於0時，兩向量前後順序為逆時針順序（在180˚之內）；等於0時，兩向量平行，也就是指夾角等於0˚或180˚；小於0時，兩向量前後順序為順時針順序（在180˚之內）。

// 向量oa與向量ob進行叉積，判斷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); }

UVa 10445

Distance

以下簡單介紹二維座標平面上計算距離的方式。

UVa 152 10514 10709

點到原點距離

float length(Point p) { return sqrt(p.x * p.x + p.y * p.y); } // 也可以運用pow來計算平方。但是沒有多大好處。 float length(Point p) { return sqrt(pow(p.x, 2) + pow(p.y, 2)); } // 也可以運用dot來計算距離。但是沒有多大好處。 float length(Point p) { return sqrt(dot(p, p)); } // 也可以直接使用標準函式庫 float length(Point p) { return hypot(p.x, p.y); }

開根號相當耗費時間。有些幾何計算，可將數學式子簡化到不必開根號，以節省計算時間。比方說，兩個距離比大小。因此，設計不開根號的程式碼，有時候也是有用處的。

// 長度的平方倍，沒有開根號。 float length2(Point p) { return p.x * p.x + p.y * p.y; }

點到點距離

float distance(Point p1, Point p2) { return sqrt((p2.x - p1.x) * (p2.x - p1.x) + (p2.y - p1.y) * (p2.y - p1.y)); } float distance(Point p1, Point p2) { float dx = p2.x - p1.x, dy = p2.y - p1.y; return sqrt(dx * dx + dy * dy); } float distance(Point p1, Point p2) { return sqrt(pow(p2.x - p1.x, 2) + pow(p2.y - p1.y, 2)); } // 重複利用之前的程式碼 // 好處是程式碼結構較穩定、可讀性高 // 壞處是程式執行時要多呼叫一次函式，各有利弊 float distance(Point p1, Point p2) { return length(p2 - p1); }

點到線距離

數學常常以ax+by+c=0表示直線，運用abs(ax+by+c)/sqrt(a²+b²)公式計算點到線距離。計算學則是直接以相異兩點表示直線，運用叉積計算點到線距離。

叉積求出v1v2組成的平行四邊形面積，然後除以v2的長度，便是垂直距離。叉積可能會有負值，請記得取絕對值，才不會得到負的距離。

// 線的資料結構，注意p1和p2不能相等！ struct Line {Point p1, p2;}; float distance(Point p, Line l) { // 面積除以底得高 Vector v1 = p - l.p1, v2 = l.p2 - l.p1; return fabs(cross(v1, v2)) / length(v2); } // 便宜行事：直線的相異兩點不包成struct，參數不使用&。 float distance(Point p, Point p1, Point p2) { // 面積除以底得高 Vector v1 = p - p1, v2 = p2 - p1; return fabs(cross(v1, v2)) / length(v2); }

點到線段距離

點到線段的最短距離，有時候是點到線上的垂直距離，有時候卻是點到線段端點的距離。

點到線段的距離，和三點共線、點在線上這些因素無關，所以這裡將空間劃分為垂直距離區和端點距離區兩塊，用點積進行判斷。這只是一種劃分方式，各位也可以自行發明適合的劃分方式。

// 線段的資料結構，注意p1和p2不能相等！ struct Line {Point p1, p2;}; typedef Line Segment; float distance(Point p, Segment s) { // 建立向量求點積，以判斷點在哪一區。 Vector v = s.p2 - s.p1; Vector v1 = p - s.p1; Vector v2 = p - s.p2; // 端點距離區：點到點距離 // 等於=可以省略 if (dot(v, v1) <= 0) return length(v1); if (dot(v, v2) >= 0) return length(v2); // 垂直距離區：點到線距離 return fabs(cross(v, v1)) / length(v); }

線段到線段距離

兩線段相交，距離為零；兩線段不相交，窮舉所有的端點到線段距離，取最短者。兩線段相交請參考後面章節，點到線段距離請參考前面章節。

各位可依照上圖所列舉的各種情況，驗證此方法是有效的。

float distance(Segment s1, Segment s2) { // 兩線段相交，距離為零。 if (intersect(s1, s2)) return 0; // 兩線段不相交，兩線段平行或者共線。 return min( distance(s1.p1, s2), // 點到線段距離 distance(s1.p2, s2), distance(s2.p1, s1), distance(s2.p2, s1) ); }

線到線距離

兩線相交，距離為零；兩線平行，距離為l1上任取一點到l2的距離。用叉積判斷平行。

float distance(Line l1, Line l2) { Vector v1 = l1.p2 - l1.p1, Vector v2 = l2.p2 - l2.p1; // 兩線不平行、兩線相交，距離為零。 if (cross(v1, v2) != 0) return 0; // 兩線平行：點到線距離 return distance(l1.p1, l2); }

Intersection

人類比電腦擅長判斷相交。人類可以追著線條移動，快速找到交點；人類也有很強的空間感，能夠迅速劃分地理位置，看一眼就能區隔出成堆的線段。但是電腦卻做不到這些，電腦只會算數字、分條件。

判斷相交原本是極容易的事情，主角改為電腦之後，卻變成極複雜的事情了。下面介紹二維座標平面上判斷相交的方式、計算交點的方式。

UVa 191 273 378 527 866 10902

點與線段相交

利用點到點距離。開根號，有誤差。

bool intersect(Point p, Point p1, Point p2) { return distance(p, p1) + distance(p, p2) == distance(p1, p2); }

比較妥當的方式，是先用叉積判斷點與線段是否共線，再用點積判斷點是否位於線段端點之間。

bool intersect(Point p, Point p1, Point p2) { // 點與線段共線，而且點要在線段之間。 return cross(p, p1, p2) == 0 && dot(p, p1, p2) <= 0; }

兩共線線段相交

線段端點位於另一條線段上面，則相交。

// 點與線段已確定共線，判斷相交。 bool intersect(Point p, Point p1, Point p2) { return p.x >= min(p1.x, p2.x) && p.x <= max(p1.x, p2.x) && p.y >= min(p1.y, p2.y) && p.y <= max(p1.y, p2.y); } // 兩線段已確定共線，判斷相交 bool intersect(Point a1, Point a2, Point b1, Point b2) { return intersect(a1, a2, b1) || intersect(a1, a2, b2) || intersect(b1, b2, a1) || intersect(b1, b2, a2); }

兩線段相交，之一

兩線段相交，也就是一條線段被另一條線段分為兩邊。兩線段端點位於另一條線段的兩側，則相交；小心處理端點共線的情況。

bool intersect(Point a1, Point a2, Point b1, Point b2) { float c1 = cross(a1, a2, b1); float c2 = cross(a1, a2, b2); float c3 = cross(b1, b2, a1); float c4 = cross(b1, b2, a2); // 端點不共線 if (c1 * c2 < 0 && c3 * c4 < 0) return true; // 端點共線 if (c1 == 0 && intersect(a1, a2, b1)) return true; if (c2 == 0 && intersect(a1, a2, b2)) return true; if (c3 == 0 && intersect(b1, b2, a1)) return true; if (c4 == 0 && intersect(b1, b2, a2)) return true; return false; }

兩線段相交，之二

先判斷兩線段在X座標軸、Y座標軸上面的投影是否相交（數線上的線段相交），然後再用叉積判斷線段端點位於另一條線段的兩側。不必煩惱端點共線的情況。

bool intersect1D(float a1, float a2, float b1, float b2) { if (a1 > a2) swap(a1, a2); if (b1 > b2) swap(b1, b2); return max(a1, b1) <= min(a2, b2); } bool intersect(Point a1, Point a2, Point b1, Point b2) { return intersect1D(a1.x, a2.x, b1.x, b2.x) && intersect1D(a1.y, a2.y, b1.y, b2.y) && cross(a1, a2, b1) * cross(a1, a2, b2) <= 0 && cross(b1, b2, a1) * cross(b1, b2, a2) <= 0; } // 正號+1，零0，負號-1。 int sgn(float x) { // 正常做法 // return x == 0 ? 0 : (x > 0 ? +1 : -1); // 取巧做法 return (x > 0) - (x < 0); // 考慮誤差 // return fabs(x) < 1e-9 ? 0 : (x > 0 ? +1 : -1); // three-way comparison operator，但是無法轉換成int。 // return x <=> 0; } bool intersect(Point a1, Point a2, Point b1, Point b2) { return intersect1D(a1.x, a2.x, b1.x, b2.x) && intersect1D(a1.y, a2.y, b1.y, b2.y) && sgn(cross(a1, a2, b1)) * sgn(cross(a1, a2, b2)) <= 0 && sgn(cross(b1, b2, a1)) * sgn(cross(b1, b2, a2)) <= 0; }

兩線段相交，之三

接著介紹不切實際的方法，但是其觀念值得一提。

線段相交，可以想像成是兩條交錯的四邊形對角線。換句話說，就是將線段的端點安排成四邊形的頂點，讓四邊形的對角線成為原來的兩條線段。如此一來，只要用一個四邊形，便可代表這兩條線段。

凸四邊形的對角線，都會相交；凹四邊形、交叉四邊形的對角線，不會相交──判斷線段相交，化作判斷凸四邊形。要判斷凸多邊形，只要順著多邊形的外圍繞一圈，看看是否一直往同側轉彎即可。判斷轉彎得利用叉積：順時針轉彎叉積得正值，逆時針得負值。叉積等於零則表示線段端點產生三點以上共線。

// 判斷線段相交，但不能處理端點共線的情況 bool intersect(Point a1, Point a2, Point b1, Point b2) { Vector v1 = b1 - a1; Vector v2 = a2 - b1; Vector v3 = b2 - a2; Vector v4 = a1 - b2; float c1 = cross(v1, v2); float c2 = cross(v2, v3); float c3 = cross(v3, v4); float c4 = cross(v4, v1); // 朝同一方向轉彎。都是順時針，或者都是逆時針。 if (c1 > 0 && c2 > 0 && c3 > 0 && c4 > 0) return true; if (c1 < 0 && c2 < 0 && c3 < 0 && c4 < 0) return true; if (c1 == 0 && intersect(a1, a2, b1)) return true; if (c2 == 0 && intersect(a1, a2, b2)) return true; if (c3 == 0 && intersect(b1, b2, a1)) return true; if (c4 == 0 && intersect(b1, b2, a2)) return true; return false; }

兩線交點

數學公式解，聯立方程式。

Point intersection(Point a1, Point a2, Point b1, Point b2) { Point a = a2 - a1, b = b2 - b1, s = b1 - a1; // 兩線平行，交點不存在。 // 兩線重疊，交點無限多。 if (cross(a, b) == 0) return INF; // 計算交點 return a1 + a * cross(s, b) / cross(a, b); }

兩共線線段交點

兩線段重疊，交點無限多；兩線段端點相接觸，交點恰為一點；兩線段不相交，交點不存在。

如果兩共線線段恰有一個交點，那麼交點就一定是線段端點。

Point intersection(Point a1, Point a2, Point b1, Point b2) { // 確定交點後，剩下的線段端點必須位於交點的不同側。 if (a1 == b1 && dot(a1, a2, b2) <= 0) return a1; if (a1 == b2 && dot(a1, a2, b1) <= 0) return a1; if (a2 == b1 && dot(a2, a1, b2) <= 0) return a2; if (a2 == b2 && dot(a2, a1, b1) <= 0) return a2; // 交點無限多、交點不存在。 return INF; }

兩線段交點

Point intersection(Point a1, Point a2, Point b1, Point b2) { float c1 = cross(a1, a2, b1); float c2 = cross(a1, a2, b2); float c3 = cross(b1, b2, a1); float c4 = cross(b1, b2, a2); if (c1 * c2 < 0 && c3 * c4 < 0) return 兩線交點; if (c1 == 0 && c2 == 0) return 兩共線線段交點; // 兩線不平行、不共線、不交叉，有可能接觸於一點。 if (c1 == 0 && intersect(a1, a2, b1)) return b1; if (c2 == 0 && intersect(a1, a2, b2)) return b2; if (c3 == 0 && intersect(b1, b2, a1)) return a1; if (c4 == 0 && intersect(b1, b2, a2)) return a2; return INF; }

兩線交點，之二

介紹另一個方法，本質等同於數學公式解。請見下圖：

以a1為基準點，以b1b2為平行四邊形的底，利用兩個平行四邊形的高的比例，便能求出a1到a2與a1到交點的距離比例。

平行四邊形的面積可用叉積運算求出，所以這個方法相當方便。實作程式碼時，要注意叉積的順序，叉積的順序將導致正負號的差異。

Point intersection(Point a1, Point a2, Point b1, Point b2) { Point a = a2 - a1, b = b2 - b1, s = b1 - a1; float c1 = cross(a, b), c2 = cross(s, b); if (c1 == 0) if (c2 == 0) return INF; // 兩線重疊 else return INF; // 兩線平行但不重疊 else return a1 + (a1a2 * (c2 / c1)); }

兩線段交點，之二

接續上一段內容。小心判斷兩個平行四邊形的高的比例是不是0到1，當線段b1b2位於a1到a2的範圍，兩者相交。至於範圍是指a1和a2兩點，以平行b1b2的方向畫出的直線作為範圍，如下圖所示：

另外，除了以a1a2為範圍，還要再以b1b2為範圍再算一遍，才能確定兩線段到底有沒有相交。

Point intersection(Point a1, Point a2, Point b1, Point b2) { Point a = a2 - a1, b = b2 - b1, s = b1 - a1; float c1 = cross(a, b), c2 = cross(s, b); Point t = a1 - b1; float c3 = cross(b, a), c4 = cross(t, a); // 調整一下正負號，方便判斷 if (c1 < 0) c1 = -c1, c2 = -c2; if (c3 < 0) c3 = -c3, c4 = -c4; if (c1 == 0) if (c2 == 0) return INF; // 兩線段共線 else return INF; // 兩線段平行但不共線 else if (c2 >= 0 && c2 <= c1 && c4 >= 0 && c4 <= c3) return a1 + (a1a2 * (c2 / c1)); else return INF; // 不平行、不相交 }

由於s與t只差一負號，故可以用s代替t。另外c1與c3也只差一負號，故可以用c1代替c3。程式碼可再精簡。

Point intersection(Point a1, Point a2, Point b1, Point b2) { Point a = a2 - a1, b = b2 - b1, s = b1 - a1; float c1 = cross(a, b); float c2 = cross(s, b); float c4 = cross(s, a); // c4方向顛倒 if (c1 < 0) c1 = -c1, c2 = -c2, c4 = -c4; if (c1 != 0 && c2 >= 0 && c2 <= c1 && c4 >= 0 && c4 <= c1) return a1 + (a * (c2 / c1)); else return INF; // 兩線段共線 }

Transformation

以下介紹二維座標平面上的常見動作。

運用C++的複數函式庫，以複數表示二維座標，就能少寫很多程式碼。

// C與C++各有自己的複數函式庫，千萬不要混用了。 // 此處是C++的複數函式庫。 #include <complex> typedef complex<float> Point; #define x real() #define y imag()

UVa 10263 12303

Translate

「位移」是直直移動，大小不變。

運用Incremental Method的精神，圖形的位移，分解成線段的位移，分解成點的位移。

點的位移有兩種想法，第一種是座標相加的概念，位移量便是座標差；第二種是向量相加的概念，位移量便是向量差。

雖然第二種想法不太直覺，但是由於向量很強大，還是習慣一下向量吧！

Point a(2, 3); // 原本的點 Point t(4, 1); // 位移量 Point b = a + t; // a位移至b

Rotate

「旋轉」，先取一個旋轉中心點，旋轉整張圖。

點的旋轉，先把旋轉中心點位移至原點，就容易處理了。

複數相乘等於長度相乘、角度相加。製造一個複數，長度等於一，角度等於旋轉角度，就可以運用複數乘法，完成點的旋轉。

Point a(5, 7); // 原本的點 Point o(2, 5); // 旋轉中心點 // 製造一個複數，長度等於一，角度等於30度。 float θ = 30.0/360.0 * 2.0*pi; // 逆時針、正30度。 Point r(cos(θ), sin(θ)); Point r = polar(1, θ); // 同上 // a以o為中心旋轉30度，旋轉至b。 Point b = (a - o) * r + o;

Scale

「縮放」，先取一個縮放中心點，縮放整張圖。

點的縮放，先把縮放中心點位移至原點，就容易處理了。

Point a(5, 7); // 原本的點 Point o(2, 5); // 縮放中心點 // p1以o為中心縮放5倍，縮放至p2。 Point b = (a - o) * 5.0 + o;

Project

「投影」，先取一條投影線，整張圖垂直投射至線上。

project這個英文單字有「投影」和「計畫」兩種意義，此處講的是「投影」。

點的投影有兩種想法，第一種想法是直線交點的概念，首先求出投影線與其法線，再解聯立方程式；第二種想法是向量點積的概念，首先把投影線位移至原點，就容易處理了。

雖然第二種想法不太直覺，但是由於向量很強大，老話一句，還是習慣一下向量吧！

Point a(5, 7); // 原本的點 Point p(1, 2), q(3, 4); // 投影線 // p位移至原點 Point pa = p - a; Point pq = p - q; // a投影至直線(p,q)，得到b。 Point pb = pq * dot(pa, pq) / norm(pq); Point b = a + pb;

Reflect

「鏡射」，先取一條對稱線，整張圖沿線翻轉，正面變反面。

reflect這個英文單字有「鏡射」和「反射」兩種意義，此處講的是「鏡射」。

有兩種實作方式，後者的程式碼比較簡潔。

Point a(5, 7); // 原本的點 Point p(1, 2), q(3, 4); // 投影線 // p位移至原點 Point pa = p - a; Point pq = p - q; // a投影至直線(p,q)，得到b。 // a鏡射至直線(p,q)，得到c。 Point pb = pq * dot(pa, pq) / norm(pq); Point c = a + pa + ((pb - pa) * 2); Point c = a + (pb * 2 - pa);