1. 程式人生 > >【筆記篇】最良心的計算幾何學習筆記(一)

【筆記篇】最良心的計算幾何學習筆記(一)

變量類型 其他 條件 parallel node ons put 是否 通過

世界以痛吻我, 我卻報之以歌。

開新坑...
雖然不知道這坑要填多久...
文章同步上傳到github...
有想看的可以去看看→_→

*溫馨提示:

  1. 看本文之前請務必學習或回顧數學-必修2解析幾何數學-必修4平面向量有關內容...
  2. 本文中的代碼基本都是口胡的, 不過還是用了IDE, 只保證能過編譯, 不保證正確OvO...

    一 浮點數相關...

眾所周知, 浮點數的運算是有誤差的...所以有時候會出現一些很蛋疼的事情.. 比如兩個該相等的東西會得出不相等的結果... 為了避免這種事情的發生, 我們要允許一些誤差... 我們定義

const double eps=1e-9;

當然這裏如果要求精度的話也可以開long double... 而且這個1e-9也是不一定的...

有些題目開1e-7 1e-8 1e-10都是不一定的... (詛咒所有卡精度的出題人吃方便面只有調料包)
這個eps有什麽用呢...
這樣我們就可以自定義一下實數的比較...
然而double作為基本變量類型是不能重載==運算符的, 這樣就要寫個函數, 而且像作者這麽zz的人很可能會忘記調用OvO

inline int dcmp(const double &a){
    if(fabs(a)<eps) return 0;   //絕對值不超過eps即視為0(卡精度的來源)
    return a<0?-1:1;            //a<0 返回-1 a>0 返回1
}

這樣我們關於double的東西就基本寫完了..

二 平面向量相關

數學-必修4了沒? 沒有快去看...
根據平面向量基本定理, 每個向量\(\vec a\)都可以寫成\(x*\vec{e_1}+y*\vec{e_2}\)的形式...
這裏的\(\vec{e_1},\vec{e_2}\)可以是任意的, 但是為了方便我們通常取正交基底(就是\(\vec{e_1}\perp\vec{e_2}\)啦~)
每個向量都可以由數對\((x,y)\)唯一確定...
這樣我們就可以這麽定義一個向量:

struct vec{
    double x,y;
};

很顯然, 這個也可以表示一個點.
所以我們可以加一句這個, 當然也可以不加..

typedef
vec point;

在此基礎上, 我們就來定義各類運算.
1.為了方便, 我們可以在vec結構體裏加一個構造函數...

struct vec{
    double x,y;
    vec(int p,int q):x(p),y(q){}
    vec(){x=y=0.0;} //默認構造函數還是要有的..
};

2.向量的模?
表示為\(|\vec v|\), 就是向量的長度啦~

inline double len(const vec &a){
    return sqrt(a.x*a.x+a.y*a.y);
}

3.比較基礎的, 向量的四則運算.
這個由定義顯然

vec operator +(const vec &a,const vec &b){
    return vec(a.x+b.x,a.y+b.y);
}
vec operator -(const vec &a,const vec &b){
    return vec(a.x-b.x,a.y-b.y);
}
vec operator *(const vec &a,const double &b){
    return vec(a.x*b,a.y*b);
}
vec operator *(const double &b,const vec &a){ //變量類型不同, 不能直接用交換律..
    return vec(a.x*b,a.y*b);
}
vec operator /(const vec &a,const double &b){
    return vec(a.x/b,a.y/b);
}

4.然後就是點積(數量積 內積 標積)和叉積(向量積 外積 矢積)了..
不過似乎並不考慮叉積的方向..(大概是因為二維計算幾何的緣故吧..)
(不過好像叉積的定義是"矢量叉積定義為由(0,0)、p1、p2和p1+p2所組成的平行四邊形的帶符號的面積"?)

double operator ^(const vec &a,const vec &b){ //點積就用^就行了 反正應該是用不到異或的..
    return a.x*b.x+a.y*b.y;
}
double operator *(const vec &a,const vec &b){
    return a.x*b.y-a.y*b.x;
}

點積的話根據必修四, 可以判斷垂直.. 若\(\vec p\cdot\vec q=0\), 則\(\vec p\perp\vec q\).

叉積的一個重要的作用就是判斷向量的順逆時針關系..

  • \(\vec p\times\vec q>0\), 則\(\vec p\)\(\vec q\)的順時針方向;
  • \(\vec p\times\vec q<0\), 則\(\vec p\)\(\vec q\)的逆時針方向;
  • \(\vec p\times\vec q=0\), 則\(\vec p\parallel\vec q\), 但可能同向也可能反向.
    有一些小小的註意事項(雖然不一定有用)..
  • 據我猜測計算幾何中沒什麽用的(但數學考試有用啊←_←)
    \(\vec a\cdot\vec b=|a||b|cos\theta, \vec a\times\vec b=|a||b|sin\theta\)
  • 點積是滿足交換律的(顯然), 而叉積是不滿足交換律的(也顯然), 不過叉積滿足反交換律,
    \(\vec a\times\vec b=-\vec b\times\vec a\)
  • 點積和叉積都滿足加法的分配律.
  • 顯然不能連續點積, 叉積不滿足結合律.
  • 其他的似乎更沒啥用而且我也看不懂 想研究的可以點上方傳送門去baidu看..

三 [跟多邊形沒關系的]各種常見的簡單的判斷

1.兩點之間距離公式 先放在這裏

inline double dis(const vec &a,const vec &b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

其實寫多了之後可能會更喜歡寫len(a-b)...

2.兩個向量的順逆時針關系 見↑↑↑

3.折線段的拐向判斷
其實跟上面沒有太大區別..
兩條有公共頂點線段\(p_0p_1\)\(p_1p_2\)(因為是折線啊),
很顯然我們可以表示兩個向量\(\vec{v_1}=p_2-p_0, \vec{v_2}=p_1-p_0\)然後求\(\vec{v_1}\times\vec{v_2}\)...

  • \(\vec{v_1}\times\vec{v_2}>0\), 則折線往右拐;
  • \(\vec{v_1}\times\vec{v_2}<0\), 則折線往左拐;
  • \(\vec{v_1}\times\vec{v_2}=0\), 則\(p_0,p_1,p_2\)三點共線.

4.判斷點\(q\)是否在以\(p_0,p_1\)為對角頂點的矩形(不加說明則默認矩形邊與坐標軸平行)內
這個簡單啊 只要分別判斷\(q\)的橫縱坐標比\(p_0,p_1\)的最小值大, 最大值小(當然是可以相等的)即可..

inline bool ptInRect(const vec &q,const vec &p0,const vec &p1){
    return min(p0.x,p1.x)<=q.x&&q.x<=max(p0.y,p1.y)&&min(p0.y,p1.y)<=q.y&&q.y<=max(p0.y,p1.y);
}

5.判斷點\(q\)是否在線\(p0p1\)
根據2. , 很顯然的可以判斷\((q-p_0)\times(p_1-p_0)==0\)
但這樣就夠了嗎? 並不是... 因為\(q\)可能在延長線or反向延長線上...
所以我們要保證\(q\)在以\(p_0,p_1\)為對角頂點的矩形內...(你以為我為什麽要寫3. ?)

inline bool ptOnSeg(const vec &q,const vec &p0,const vec &p1){
    return (p1-p0)*q==0&&ptInRect(q,p0,p1);
}

這樣就可以咯..

6.判斷線段\(p_0p_1\)和直線\(q_0q_1\)是否相交
如果線段與直線相交, 那麽線段跨立直線..
技術分享圖片
首先 如圖,我們可以得出\(\vec{v_1}=p_0-q_0,\vec{v_2}=q_1-q_0,\vec{v_3}=p_1-q_0\).
由圖可以看出, 我們要讓\(p_0\)\(p_1\)放到直線\(q_0q_1\)的兩側, 就要讓\(\vec{v_1}\times\vec{v_2}\)\(\vec{v_2}\times\vec{v_3}\)同號(或者有一方為0)
所以有\((\vec{v_1}\times\vec{v_2})*(\vec{v_2}\times\vec{v_3})\geqslant0\), 也就是\(((p_0-q_0)\times(q_1-q_0))*((q_1-q_0)\times(p_1-q_0))\geqslant0\).

inline bool segCutLine(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    return ((p0-q0)*(q1-q0))*((q1-q0)*(p1-q0))>=0;
}

7.判斷兩條線段\(p_0p_1,q_0q_1\)是否相交

這個我們可以分為兩步

  • 快速排斥實驗
    兩個線段的端點為對角頂點的兩個矩形如果不相交, 那麽這兩條線段顯然不相交.

  • 跨立實驗
    明確一點: 如果兩線段相交, 則兩線段必然相互跨立對方.
    我們就可以根據上面的跨立方式得出要判斷\(((p_0-q_0)\times(q_1-q_0))*((q_1-q_0)\times(p_1-q_0))>0\ (p_0p_1跨立q_0q_1)\)\(((q_0-p_0)\times(p_1-q_0))*((p_1-q_0)\times(q_1-p_0))>0\ (q_0q_1跨立p_0p_1)\)
    而當前面這一坨等於0的時候呢? 我們可以得知共線, 而由於已經通過了快速排斥實驗, 所以端點一定在另一條線段上. 所以我們只要讓\(((p_0-q_0)\times(q_1-q_0))*((q_1-q_0)\times(p_1-q_0))\geqslant0\)\(((q_0-p_0)\times(p_1-q_0))*((p_1-q_0)\times(q_1-p_0))\geqslant0\)就可以了 (其實好啰嗦啊)

inline bool segCutSeg(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    if(min(p0.x,p1.x)>max(q0.x,q1.x)||min(q0.x,q0.y)>max(p0.x,p1.x)
        ||min(p0.y,p1.y)>max(q0.y,q1.y)||min(q0.y,q1.y)>max(p0.y,p1.y)) return false; //快速排斥實驗
    return ((p0-q0)*(q1-q0))*((q1-q0)*(p1-q0))>=0&&((q0-p0)*(p1-p0))*((p1-p0)*(q1-p0))>=0; //跨立實驗
}

8.判斷矩形是否在矩形中
只需要簡單的比較邊界就行了.

inline bool rectInRect(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    return min(p0.x,p1.x)<min(q0.x,q1.x)&&max(p0.x,p1.x)>max(q0.x,q1.x)&&min(p0.y,p1.y)<min(q0.y,q1.y)&&max(p0.y,p1.y)>max(q0.y,q1.y);
}

9.判斷圓是否在矩形中
圓在矩形中的充要條件是圓心在矩形內而且圓的半徑小於到矩形四邊距離的最小值.
話說圓要怎麽表示? 一個點表示圓心 一個double表示半徑嘛= =

inline bool cirInRect(const vec &o,const double &r,const vec &p0,const vec &p1){
    return ptInRect(o,p0,p1)&&min(
        min(abs(min(p0.x,p1.x)-o.x),abs(max(p0.x,p1.x)-o.x)),
        min(abs(min(p0.y,p1.y)-o.y),abs(max(p0.y,p1.y)-o.y)))<r; //這麽寫好像清楚一點?
}

10.判斷點是否在圓中
與圓心的距離小於半徑啊~

inline bool ptInCir(const vec &p,const vec &o,const double &r){
    return dis(p,o)<r;
}

11.判斷線段、折線、多邊形是否在圓中
圓是個凸集, 所以我們依次判斷每個端點是否在圓內即可..
代碼就不給了OvO

12.判斷圓\(O_1\)是否在圓\(O_2\)
首先嘛 肯定要\(r_1<r_2\)... 其次呢 要有\(|O_1O_2|<=r_2-r_1\)(其實好像可以和上一條一起判.. 因為負數顯然不會大於距離..) (等號成立時兩圓內切(突然不知道算不算在圓內了..先算是吧...(大不了到時候摳掉嘛)))

inline bool cirInCir(const vec &o1,const double &r1,const vec &o2,const double &r2){
    return dis(o1,o2)<=r2-r1;
}

13.點\(p\)到直線\(q0q1\)的距離,利用叉積求出平行四邊形的面積, 然後算出底邊長. 一比得到高即可.

inline bool ptDisLine(const vec &p,const vec &q0,const vec &q1){
    return (p-q0)*(q1-q0)/len(q1-q0);
}

14.點\(p\)到線段\(q0q1\)的距離,先判斷垂足是否在線段上.
是則同上, 不是則找距離兩個端點中近的那個
如何判斷呢?
技術分享圖片
如圖所示, 如果垂足在直線外, 那麽夾角\(\theta\)一定是鈍角, 則\(cos\theta<0\)
那麽\(\vec{e_1}\cdot\vec{e_2}=|\vec{e_1}||\vec{e_2}|cos\theta<0\),
而且這樣求出的\(\theta\)的頂點就是較近的點..
兩邊一求就好了...

double ptDisSeg(const point &p,const point &a,const point &b){
    if(dcmp((p-a)^(b-a))<0) return dist(p,a);
    if(dcmp((p-b)^(a-b))<0) return dist(p,b);
    return fabs((p-a)*(b-a))/dist(b,a);
}

15.兩條相交直線(線段)的交點. (所以要先判斷是否相交...)
好像兩種方法?
一種是利用解析幾何推導.
對於直線\(l:Ax+By+C=0\)和直線上的兩點\((x_0,y_0),(x_1,y_1)\)
斜率\(k=-\frac{A}{B}=\frac{y_0-y_1}{x_0-x_1}\),
那我們不妨設\(A=y_0-y_1,B=x_1-x_0\), 代入方程得到\(C=x_0y_1-x_1y_0\)
然後回到問題, 設兩條直線分別是
\[ l_1:a_1x+b_1y+c_1=0,\ l_2:a_1x+b_1y+c_2=0. \]
這樣(x,y)就是聯立得到方程組的解.
根據相關知識, 我們可以得到
\[ \left\{\begin{matrix} x=\frac{b_2c_1-b_1c_2}{b_1a_2-b_2a_1}\\ y=\frac{a_1c_2-a_2c_1}{b_1a_2-b_2a_1} \end{matrix}\right. \]
而這些量都可以用叉積的形式表示出來..
所以就可以這麽寫:

inline vec lineCutLineNode(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    double a1,b1,c1,a2,b2,c2,d;
    a1=p1.y-p0.y; b1=p0.x-p1.x; c1=p0*p1;
    a2=q1.y-q0.y; b2=q0.x-q1.x; c2=q0*q1;
    d=a1*b2-a2*b1;
    return vec((b2*c1-b1*c2)/d,(a1*c2-a2*c1)/d);
}

這種方法略微麻煩, 但是精度比較好.
另一種方法是利用叉積的比值.
利用參數方程的知識(數學 選修2-1 今天剛學噠XD), 我們可以用一個定點和一個方向向量確定一條直線.
令兩條直線
\[ l1=P+t\vec v, \ l2=Q+t\vec w \]
? 然後因為是交點, 所以同時滿足兩條直線的方程, 令交點在兩條直線的參數分別為\(\lambda\)\(\mu\), 則有
\[ P+\lambda\vec v=Q+\mu\vec w \]
? 兩邊分別叉乘\(\vec w\)
\[ (P+\lambda\vec v)\times\vec w=(Q+\mu\vec w)\times w \]
? 去括號,
\[ P\times\vec w+\lambda\vec v\times\vec w=Q\times\vec w+\mu\vec w\times\vec w \]
? 又\(\vec w\)一定與自己共線, 所以\(\vec w\times\vec w=0\), 這樣就沒有\mu了, 移項, 合並同類項,
\[ (Q-P)\times\vec w=\lambda\vec v\times\vec w \]
? 系數除到另一邊去, 得
\[ \lambda=\frac{(Q-P)\times\vec w}{\vec v\times\vec w} \]

? 那麽這樣就都是已知量了, 代入進去即可.
? 最後要求的交點就是\(P+\lambda\vec v\)啦~
?
留一道例題:poj1269

完整的代碼在github裏面喲~

【筆記篇】最良心的計算幾何學習筆記(一)