計算幾何 – Computational Geometry – 課程大綱
計算幾何
Computational Geometry
eddy1021
課程大綱
- 座標與向量
- 有向面積
- 線段相交
- 誤差分析
何爲計算幾何
座標與向量
表示平面上的幾何圖形
- 長度
- 角度
- 座標
- 向量
實作上表示平面幾何
座標、向量
實作上表示平面幾何
#include <utility> // include pair
//typedef std::pair<int,int> Pt;
typedef std::pair<double,double> Pt;
#define X first
#define Y second
Pt point( double x , double y ){
return make_pair( x , y );
}
int main(){
Pt a = point( 4 , 3 );
printf( "%d %d\n" , a.X , a.Y );
}
向量(數學補充)
內積(Dot): $A \cdot B = |A| |B| cos \theta $ $A \cdot B = A_x B_x + A_y B_y$ 外積(cross): $A \times B = |A| |B| sin \theta $ $A \times B = A_x B_y - A_y B_x$
向量基本操作
typedef std::pair<double,double> Pt;
#define X first
#define Y second
Pt operator+( const Pt& p1 , const Pt& p2 ){
return Pt( p1.X + p2.X , p1.Y + p2.Y );
}
Pt operator-( const Pt& p1 , const Pt& p2 ){
return Pt( p1.X - p2.X , p1.Y - p2.Y );
}
double operator*( const Pt& p1 , const Pt& p2 ){
return p1.X * p2.X + p1.Y * p2.Y;
}
double operator^( const Pt& p1 , const Pt& p2 ){
return p1.X * p2.Y - p1.Y * p2.X;
}
向量基本操作
Pt operator*( const Pt& p1 , const double& k ){
return Pt( p1.X * k , p1.Y * k );
}
Pt operator/( const Pt& p1 , const double& k ){
return Pt( p1.X / k , p1.Y / k );
}
double abs( const Pt& p1 ){
return sqrt( p1 * p1 );
}
有向面積
兩向量形成三角形
兩向量可夾出三角形 其面積可由外積得出： $\Delta OAB = \frac{1}{2} \vec{A} \times \vec{B}$ 此面積我們稱爲$\textbf{有向面積}$ (外積有正有負)
有向面積的方向
$\vec{AB} \times \vec{AC}$ $=(5,0) \times (3,4)$ $= 5 \times 4 - 0 \times 3$ $=20 > 0$
有向面積的方向
$\vec{AC} \times \vec{AB}$ $=(3,4) \times (5,0)$ $= 3 \times 0 - 4 \times 5$ $= -20 < 0$
多邊形的有向面積
如同三角形 多邊形的有向面積： $\textbf{ 逆時針爲正}$ $\textbf{ 順時針爲負}$
多邊形的有向面積
多邊形的有向面積
任意在平面上選一點 $A$ 將所有點與 $A$ 點連線 相鄰兩點將與 $A$ 點形成三角形 依逆時針(或順時針)序依序加總各三角形有向面積 即爲該多邊形之有向面積
多邊形的有向面積
一般而言，會挑選原點作爲參考點 若一個多邊形的頂點依序爲： $P_0,P_1,\ldots,P_{N-1},P_N=P_0$ 則多邊形的有向面積公式爲： $\text{Area}=\frac{1}{2} \sum\limits_{i=0}^{N-1} \vec{P_i} \times \vec{P_{i+1}}$
線段相交
如何判斷？
直接求出交點，再判斷交點是否在線段內
如何判斷？
直接求出交點，再判斷交點是否在線段內
- 求交點 $\Rightarrow$ 沒有交點？無限多交點？
- 交點不好求 e.g. 垂直線
- 交點位置誤差？
更精準的求線段相交
利用外積定義方向函式：int ori( const Pt& o , const Pt& a , const Pt& b ){
double cross = ( a - o ) ^ ( b - o );
if( fabs( cross ) < eps ) return 0;
return cross > 0 ? 1 : -1;
}
更精準的求線段相交
若線段 $P_1P_2$ 與 $P_3P_4$ 相交。 則點 $P_3 與點 P_4$ 會在線段 $P_1P_2$ 異側。更精準的求線段相交
若線段 $P_1P_2$ 與 $P_3P_4$ 相交。 則點 $P_3 與點 P_4$ 會在線段 $P_1P_2$ 異側。 利用方向函式即代表： $\texttt{ori}(P_1,P_2,P_3)\times \texttt{ori}(P_1,P_2,P_4) < 0$這樣就考慮完了嗎？
兩條平行線也可能相交
兩條平行線也可能相交
共線才有可能相交！ 平行： $(P_2 - P_1) \text{^} (P_4 - P_3) = 0$ 共線： $\texttt{ori}(P_1,P_2,P_3) = 0$平行線相交交點位置
點 $P_1$ 在 線段 $P_2P_3$ 上： $(P_2 - P_1) \cdot (P_3 - P_1) < 0$誤差分析
爲什麼需要誤差分析
- 計算幾何中經常遇到浮點數
- 浮點數以二進位儲存必然會產生誤差
- $\frac{1}{3},\sqrt{2},\pi$
浮點數儲存誤差
puts( 0.1 + 0.2 == 0.3 ? "equal" : "not equal" );形態 大小 精度 float 4 B $10^{-7}$ double 8 B $10^{-16}$ long double 10 B $10^{-19}$
誤差容忍值($\texttt{eps}$)
將所有數字膨脹 eps 的大小 $x \Rightarrow (x-\texttt{eps},x+\texttt{eps})$ 兩數相差 $\texttt{eps}$ 內視爲相等容忍誤差的比較
bool equal( const double& a , const double& b ){
return b - eps < a && a < b + eps;
}
bool less( const double& a , const double& b ){
return a < b - eps;
}
bool lessOrEqual( const double& a , const double& b ){
return a < b + eps;
}
如何決定$\texttt{eps}$的大小
- 不能太小，應相等的數判成不相等
- 不能太大，不應相等的數判成相等
如何決定$\texttt{eps}$的大小
- 不能太小，應相等的數判成不相等
- 不能太大，不應相等的數判成相等
估算$\texttt{eps}的下界$
誤差會有疊加的現象！估算$\texttt{eps}的下界$
加減法 $\Rightarrow$ 絕對誤差相加 $(x+\Delta x) \pm (y+\Delta y) = (x \pm y) + (\Delta x \pm \Delta y)$估算$\texttt{eps}的下界$
乘除法 $\Rightarrow$ 相對誤差相加 $(x+\Delta x) \cdot (y+\Delta y) = xy( 1 + \frac{ \Delta x }{x} )(1+\frac{ \Delta y }{ y })$ $\approx xy( 1 + \frac{ \Delta x }{x} + \frac{ \Delta y }{ y } )$估算$\texttt{eps}的下界$
經過 $K$ 次運算，相對誤差不會超過 $K \epsilon$ 因此，數字範圍在 $V$ 內時，eps 至少要是 $V K \epsilon$估算$\texttt{eps}的上界$
要能分辨不同的數字！估算$\texttt{eps}的上界$
- 保證變動 $10^{-7}$ 不影響答案
- 答案輸出至小數點後第六位
- 答案與正確答案絕對或相對誤差小於 $10^{-6}$ 視爲正確
計算幾何避免誤差大法！
不到最後關頭，絕不用浮點數！更多數學
徑度
$\pi=3.14159265358979...=180^{\circ}$
三角函數
$cos(x)=\frac{b}{c}$ $sin(x)=\frac{a}{c}$ $tan(x)=\frac{a}{b}$
廣義三角函數
$cos(\theta)=\frac{x}{r}$ $sin(\theta)=\frac{y}{r}$ $tan(\theta)=\frac{y}{x}$ $atan2(y,x)=\theta$
凸包
凸包
Monotone Chain
極角排序
極角排序
Sort by angle
Pt o;
D angle( const Pt& x ){
return atan2( x.Y , x.X );
}
bool cmp( Pt a , Pt b ){
return angle( a - o ) < angle( b - o );
}
Sort by cross
Pt o;
bool cmp( Pt a , Pt b ){
return ( a - o ) ^ ( b - o ) > 0;
}
平面上 $N$ 個點，問一條直線最多通過幾個點
$O(N^3) \Rightarrow O(N^2\lg N)$