SE_Work2_交點個數

項目	內容
課程：北航-2020-春-軟件工程	博客園班級博客
要求：求交點個數	我的項目做業
班級：005	Sample
GitHub地址	IntersectProject

1、PSP估算

1. 在開始實現程序以前，在下述 PSP 表格記錄下你估計將在程序的各個模塊的開發上耗費的時間。（0.5'）
2. 在你實現完程序以後，在下述 PSP 表格記錄下你在程序的各個模塊上實際花費的時間。（0.5'）

PSP2.1	Personal Software Process Stages	預估耗時（分鐘）	實際耗時（分鐘）
Planning	計劃
· Estimate	· 估計這個任務須要多少時間	10	10
Development	開發
· Analysis	· 需求分析 (包括學習新技術)	20	60
· Design Spec	· 生成設計文檔	10	10
· Design Review	· 設計複審 (和同事審覈設計文檔)	5	0
· Coding Standard	· 代碼規範 (爲目前的開發制定合適的規範)	5	0
· Design	· 具體設計	20	60
· Coding	· 具體編碼	60	120
· Code Review	· 代碼複審	20	0
· Test	· 測試（自我測試，修改代碼，提交修改）	30	60
Reporting	報告
· Test Report	· 測試報告	20	60
· Size Measurement	· 計算工做量	20	10
· Postmortem & Process Improvement Plan	· 過後總結, 並提出過程改進計劃	20	60
	合計	220	450

如圖所示，第一次項目的估計狀況是嚴重與事實不符的，時間上花費了兩倍有餘，多是由於事先根本對這些程序沒有很是清晰的概念。

2、思路及設計

解題思路描述。即剛開始拿到題目後，如何思考，如何找資料的過程。（3'）

設計實現過程。設計包括代碼如何組織，好比會有幾個類，幾個函數，他們之間關係如何，關鍵函數是否須要畫出流程圖？單元測試是怎麼設計的？（4'）

代碼說明。展現出項目關鍵代碼，並解釋思路與註釋說明。（3'）

1. 解題思路

最開始最直接的想法就是暴力枚舉，每次取一條，和其餘全部圖形求交點，加入集合中，讓集合自動去重。複雜度爲\(O(N^2)\) ，對於\(1\leq N\leq 1000\) 徹底沒問題，可是通過試驗，\(N>10^5\) 時不能達到性能要求。

實際上每次取出一條，不用和其餘全部圖形都求一遍，以前取出過的圖形已經求過了不用再求，只用和沒有取出的圖形求一遍，可是這種思路複雜度沒有降低，仍然是一種「暴力求解法」。

轉而思考其餘方法，動態規劃就是一種極佳的辦法，也在網絡上搜到相似的題目：求直線交點數目，可是並不符合題意，由於他們有一個很重要的假設，任意三條直線不想交與一點。

咱們知道，若是只有直線相交，任意三條直線不交一點，且沒有直線平行，則總數是固定的\(\frac{n(n-1)}{2}\)，而每出現下面所示的狀況，總數都會相應減少：

• 聚合點（n條線交於一點）：總數減\(\frac{n(n-1)}{2}-1\)
• 平行線（n條線平行）：總數減\(\frac{n(n-1)}{2}\)

平行線很好說，創建集合列一個個加進來就好，然而真正要命的是聚合點怎麼去尋找。三條直線交於一點，不只須要考慮斜率，還得考慮截距，即ax+by+c=0 中a、b、c都須要考慮。n條直線相交於一點即如下方程組有解：

$$\begin{cases} a_1 x+b_1y=-c_1\\ a_2 x+b_2y=-c_2\\...\\ a_n x+b_ny=-c_n\\ \end{cases} $$

至關於一個維數大於秩的矩陣有解：

$$ \begin{pmatrix} a_{1} & b_{1}\\ \vdots & \vdots\\ a_{n} & b_{n} \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} = -\begin{pmatrix} c_1\\ \vdots\\c_n \end{pmatrix}$$

至關於：

$$ \begin{pmatrix} a_{1} & b_{1} & c_1\\ \vdots & \vdots &\vdots\\ a_{n} & b_{n} &c_n \end{pmatrix} \begin{pmatrix} x\\ y \\1 \end{pmatrix} = 0$$

其中，

$$\begin{pmatrix} a_{1} \ \vdots\ a_{n} \end{pmatrix}$$ $$\begin{pmatrix} b_{1} \ \vdots\ b_{n} \end{pmatrix}$$ $$\begin{pmatrix} c_{1} \ \vdots\ c_{n} \end{pmatrix}$$ 是在n維空間中的一個向量，而它們的線性組合要能爲0向量（其中c爲必要項），最終要保證的是三個向量共面！

然而就算你能以\(O(1)\)的複雜度判斷多條直線是否相交與一點，要遍歷全部的組合又比如一顆巨大的樹。最終只能望而卻步。

2. 設計實現過程

主函數：每次輸入一個圖形，與其餘全部圖形進行相交運算，將所得的結果存入交點集合

Dot類：繼承與\(Pair<float, float>\) ，表示點的座標或者一個向量，能進行加減乘除運算，同時有求模長、單位化方法

class Dot : public pair<float, float> {
public:
    Dot(float a, float b) {
        first = a;
        second = b;
    }

    inline Dot operator+(Dot dot) { return {first + dot.first, second + dot.second}; }
    inline Dot operator-(Dot dot) { return {first - dot.first, second - dot.second}; }
    inline Dot operator*(float t) { return {first * t, second * t}; }
    inline Dot operator/(float t) { return {first / t, second / t}; }

    inline float norm() { return first*first + second*second; }
    inline float abs() { return sqrt(norm()); }
    inline void unify() { *this=*this/abs(); }

};

Line類：有經過兩個點座標初始化方法和三個參數a,b,c初始化方法，有判斷是否與另外一條線平行，求與另外一條線交點的方法

class Line {
public:
    float a, b, c;

    Line(float x0, float y0, float x1, float y1);
    Line(float aa, float ab, float ac);

    bool parallel(Line l);
    void intersect(set<Dot> *intersections, Line l);
    Dot intersect(Line l);
}

Circle類：有經過圓心座標和半徑初始化方法，求到一條直線距離，求與另外一個圓交點，求與直線的交點方法。

class Circle {
public:
    float x, y, r;

    Circle(float ax, float ay, float ar) : x(ax), y(ay), r(ar) {}

    float distance(Line l);
    void intersect(set<Dot> *intersections, Line l);
    void intersect(set<Dot> *intersections, Circle c);
};

3. 單元測試

單元測試分別設計了直線與直線相交（包括平行和相交），圓與直線相交（包括相離、相切、相交），圓與圓相交（包括相離、相切、相交）三種狀況，輸出最後的交點座標。

如圖所示，在這三種測試下，結果均正確。

4. 關鍵代碼說明

由於交點個數多是0、一、2個，無法設定一個定長的返回值，故只好傳入一個指向點集的指針，有幾個交點就存幾個進去。

首先是Line的兩個intersect方法：

// 直接經過公式，返回交點，不想交則報錯
Dot Line::intersect(Line l) { 
    if(parallel(l)) throw exception();
    return Dot(b * l.c - l.b * c, l.a * c - a * l.c) / (a * l.b - l.a * b);
}

// 將交點加入集合中，忽略錯誤
void Line::intersect(set<Dot> *intersections, Line l) {
    try {
        intersections->insert(intersect(l));
    } catch (exception e) {}
}

其次是Circle的兩個intersect方法，和circle相交的算法直接經過公式計算，再也不贅述。而與直線相交須要先作一個到該直線的垂線，求出垂足：

void Circle::intersect(set<Dot> *intersections, Line l) {
    float d = distance(l);
    if (d > r) return;
    Line ll(l.b, -l.a, l.a * y - l.b * x);	// 垂線
    Dot dot = ll.intersect(l);	// 垂足
    if (d == r) {
        intersections->insert(dot);
        return;
    }
    float t = sqrt(r * r - d * d);	// 距垂足長度
    Dot direction(l.b, -l.a);	// 方向向量
    direction.unify();			// 單位化
    intersections->insert(dot + direction * t);
    intersections->insert(dot - direction * t);
}

3、性能及優化

記錄在改進程序性能上所花費的時間，描述你改進的思路，並展現一張性能分析圖（由 VS 2019 的性能分析工具自動生成），並展現你程序中消耗最大的函數。（3'）

N	時間(ms)
200	9
400	36
600	92
800	175
1000	303
2000	1429
3000	3458
4000	6463
5000	9820
10000	40873

如圖所示採用暴力解法，結果大體呈現平方級複雜度，\(t=9(\frac{N}{200})^2\)，當N取 500w，須要時間5625Mms =1562.5h是根本不可接受的。

經過性能分析圖的描述如上，我發現程序中消耗最大的函數是void intersect(set<Dot> *intersections, Circle c);方法：

不可避免地出現了大量不可簡化的計算。

void Circle::intersect(set<Dot> *intersections, Circle c) {
    float a1, b1, R1, a2, b2, R2;
    a1 = x;
    b1 = y;
    R1 = r;

    a2 = c.x;
    b2 = c.y;
    R2 = c.r;

    float R1R1 = R1 * R1;
    float a1a1 = a1 * a1;
    float b1b1 = b1 * b1;

    float a2a2 = a2 * a2;
    float b2b2 = b2 * b2;
    float R2R2 = R2 * R2;

    float subs1 = a1a1 - 2 * a1 * a2 + a2a2 + b1b1 - 2 * b1 * b2 + b2b2;
    if(subs1<=0) return;

    float subs2 = -R1R1 * a1 + R1R1 * a2 + R2R2 * a1 - R2R2 * a2 + a1a1 * a1 - a1a1 * a2 - a1 * a2a2 + a1 * b1b1 -
                  2 * a1 * b1 * b2 + a1 * b2b2 + a2a2 * a2 + a2 * b1b1 - 2 * a2 * b1 * b2 + a2 * b2b2;
    float subs3 = -R1R1 * b1 + R1R1 * b2 + R2R2 * b1 - R2R2 * b2 + a1a1 * b1 + a1a1 * b2 - 2 * a1 * a2 * b1 -
                  2 * a1 * a2 * b2 + a2a2 * b1 + a2a2 * b2 + b1b1 * b1 - b1b1 * b2 - b1 * b2b2 + b2b2 * b2;
    float sigma = sqrt((R1R1 + 2 * R1 * R2 + R2R2 - a1a1 + 2 * a1 * a2 - a2a2 - b1b1 + 2 * b1 * b2 - b2b2) *
                       (-R1R1 + 2 * R1 * R2 - R2R2 + subs1));

    Dot dot1(subs2, subs3), dot2(b1 - b2, a2 - a1);
    intersections->insert((dot1 + dot2 * sigma) / (2 * subs1));
    intersections->insert((dot1 - dot2 * sigma) / (2 * subs1));
}