計算幾何

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}$的大小

  • 不能太小，應相等的數判成不相等
  • 不能太大，不應相等的數判成相等


一般視題目而定，可大到 $10^{-2}$ 或小到 $10^{-14}$。
多落在 $10^{-6}$ 到 $10^{-12}$ 之間

估算$\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)$