梁友棟-柏世奇演算法

梁友棟—柏世奇演算法（以梁友棟和柏世奇的名字命名）是電腦圖形學中的一個線段裁剪演算法。梁友棟—柏世奇演算法使用直線的參數方程和不等式組來描述線段和裁剪窗口的交集。求解出的交集將被用於獲知線的哪些部分是應當繪製在螢幕上的。這一演算法比科恩－蘇澤蘭演算法（Cohen-Sutherland algorithm）要更加高效，梁友棟—柏世奇演算法的基本思想是：在計算線段與裁剪窗交集之前做儘可能多的判斷。

## 演算法描述

考慮直線的參數方程：

$$x = x_0 + t (x_1 - x_0) = x_0 + t \Delta x,$$

$$y = y_0 + t (y_1 - y_0) = y_0 + t \Delta y.$$

點在裁剪窗內，若

$$x_\text{min} \le x_0 + t \Delta x \le x_\text{max}$$

且

$$y_\text{min} \le y_0 + t \Delta y \le y_\text{max},$$

其可用4個不等式表達：

$$t p_i \le q_i, \quad i = 1, 2, 3, 4,$$

其中

$$\begin{align} <pre> p_1 &= -\Delta x, & q_1 &= x_0 - x_\text{min}, & &\text{( 左 )} \\ p_2 &= \Delta x, & q_2 &= x_\text{max} - x_0, & &\text{( 右 )} \\ p_3 &= -\Delta y, & q_3 &= y_0 - y_\text{min}, & &\text{( 下 )} \\ p_4 &= \Delta y, & q_4 &= y_\text{max} - y_0. & &\text{( 上 )} </pre> \end{align}$$

計算最終線段：

- 與裁剪窗平行的直線在平行的邊界上有 $p_i = 0$
- 若對於這樣的 $i$ $q_i < 0$，則線段全部在裁剪窗的外面，可以被消除
- 當 $p_i < 0$時，線從裁剪窗外向內走；$p_i > 0$
- 對非零的 $p_k$, $u = q_i / p_i$
- 對每條線，計算 $u_1$ 和 $u_2$。對 $u_1$檢查 $p_i < 0$ 的邊界（即從外向內）。令 $u_1$ 為 $\{0, q_i / p_i\}$ $u_2$檢查 $p_i > 0$ 的邊界（即從內向外）。令 $u_2$ 為 $\{1, q_i / p_i\}$ $u_1 > u_2$

## 演算法講解

Liang-Barsky 演算法（梁友棟演算法）是圖形幾何中一種非常有效的找到一個線段存在於一個 AABB 範圍內的那一段，也就是一個「裁剪」演算法。

### 線段表達式

想要理解這個演算法之前，我們先需要了解一種表達線段的方式：`P(t) = s + d * t`。

其中，t 表示一個線上的點是整條線的距離的多少部分（百分之多少），t ∈ 1, 0，因此 t 就可以理解為歸一化（normalized）的線長度比例。

所有關於 t 的計算都是從線段的一個端點：s 出發。因此整個線段也可以理解為射線。

s 是線段的一個端點，在這裡被算作是起點；d 是整個線段起點到終點的差距，也就是和這個線段全等的向量：dx, dy。

整個用來表達線段的方程：`P(t) = s + d * t` 就是在描述線上的一點在指定 t 的情況在是在什麼位置。

### 算法方程

```text
p * t <= q
```

整個方程可以被理解為在一個給定的 t 值時，這個線滿足一個 constraint。

這個 constraint，就是整個演算法的核心：

有一個在線上的點在 AABB 範圍之內。

而這個整體的 constraint，可以被 break down into 4 個小 constrain：上下左右邊的 constraint：

有一個在線上的點在某邊的一側。例如：有一點在左側邊的右邊；有一點在上側邊的下方；有一點在右側邊的左邊；有一點在下側邊的上方。

當一條線段全部滿足這四個 constraint 的時候，它就 guarantee 是在這個 AABB 空間內的。而本演算法就是要求出它在 AABB 範圍內的起點和終點位置。

#### 過程詳解

在計算線段在 AABB 範圍內的起點和終點位置的過程中，本質上的原理就是以 AABB 空間的四條邊作為「刀」，切除空間外線段的部分。

在這過程中，我們會修改兩個 t 的值：他們代表裁剪後的線段的 begin 和 end 的 t：tStart = 0，tEnd = 1。

在修改這兩個 t 的時候，each t can only go forward, which means going towards the direction of the other t and cannot go back. 這也就意味著 tStart 只能加，tEnd 只能減。

當發現唯一滿足 constraint 的 t 值是需要往後走的，或者需要走到超過另外一個 t 值，那就證明這個 constraint 無法被滿足，也就證明線段不在 AABB 空間裡。

#### 方程詳解

這個方程是抽象後的結果，因此一下看不懂很正常。以下所有解釋將以左側邊作為例子，其他邊同理。

在解釋所有東西之前，我們先假設開始點：start 在 AABB 範圍內，也就是假設 xmin < sx < xmax && ymin < sy < ymax。

單個這樣的方程表示四個 constraint 中的一個。每一個方程處理的是單個維度（either x or y value），單條線（of the AABB range）。

一個這樣的方程所構成的約束就是是否存在一個 t 值，這個 t 值在線上的位置滿足約束。

而這個「滿足約束」指的就是這個點在邊界的朝向 AABB 範圍的一邊。

##### p

p 總體上代表的是整個線段在對應維度，in this case 是 x，的長度，也就是線段起點到線段終點的 x 差距（dx）。

這個 p 值在這裡可以被理解為速度，在經過 t 時間後就達到了某個特定的 x 值位置。

但是 p 值並不直接等於 dx，但是 p 的絕對值等於 dx。

p 的真正含義是：線段在特定維度，相對於邊界方向的長度。

記住這個定理：在假設 start 在約束內時，p 為正時：線段向 end 的方向前進就是在向邊界的方向前進；p 為負時：線段向 end 方向前進就是在向遠離邊界的方向前進。

而當假設 start 不在約束內時，那麼 p 為正時的行為就和之前 p 為負的行為一樣，反之亦然。

##### q

q 代表開始點（sx）與邊界的距離。

但是，這個距離也是相對的，這代表它也是有符號的。

可以這麼理解這個距離：

將任何正在解的邊界看做 y = 0 的那條線（x 軸），將開始點的位置看做是在以邊界作為 x 軸的坐標系中的一個坐標，q 就是該坐標的 y 值。

所有 q 為正的值就代表該點滿足約束：在正在解的邊界的靠近 AABB 範圍的一邊

##### 方程推導

對於整個演算法嘗試解決的問題，具象化後可以表達為這樣一個約束：

```text
x >= xmin && x <= xmax && y >= ymin && y <= ymax
```

這個約束實際上包含四個約束，分別對應 AABB 範圍的四條邊。例如這是左側邊：

```text
x >= xmin
```

而在這裡，x 是實際世界坐標，我們想要的是它對於線段的坐標，確切地說是對於線段 start，也就是 sx 的坐標：

由於 x 是線上的一個點，所以我們可以把 x 換成：

```text
sx + t * dx >= xmin
```

然後，為了所有的方程最終都可以以同樣的格式表示，我們在這裡進行一些 rearrangement：

```text
t * dx >= xmin - sx
-dx * t <= sx - xmin
```

而在這裡，-dx 就是 p，sx - xmin 就是 q。

而其他三個約束也可以轉化為這樣的格式：p * t <= q。這一步是利於我們通過一個函數解決所有約束而不需要寫 8 個函數。

再舉一個例子：x <= xmax

```text
sx + t * dx <= xmax
t * dx <= xmax - sx
```

在這裡 dx 就是 p，xmax - sx 就是 q。

### 判斷邏輯

我們需要判斷的就是能夠和邊界相交的那個點的 t 和我目前有的 t 是什麼關係。

由於我們知道只要兩條直線不平行，他們就絕對有一個相交點。

因此第一件事就是處理平行的情況：

```cpp
if (p == 0.0f) return q >= 0.0f;
```

然後，就是處理所有正常情況，也就是找到這個相交點的 t。

由於我們已經把數據進行了抽象，所以我們並不需要 solve Linear Equations，而是直接用已有數據計算就可以：

```cpp
float tMeet = q / p;
```

再算出交點 t 了之後，我們就可以正式進行判斷了。

### 最終單約束代碼（C++）

```cpp
bool liang_solveOneConstraint(float p, float q, float& tStart, float& tEnd) {
    if (p == 0.0f) return q >= 0.0f;

    float tMeet = q / p;

    if (p > 0.0f) {
        if (tMeet < tStart) return false;
        if (tMeet < tEnd) tEnd = tMeet;
    } else {
        if (tMeet > tEnd) return false;
        if (tMeet > tStart) tStart = tMeet;
    }
}
```

## 最終代碼（C++）

```cpp
bool liang_solveOneConstraint(float p, float q, float& tStart, float& tEnd) {
    if (p == 0.0f) return q >= 0.0f;

    float tMeet = q / p;

    if (p > 0.0f) {
        if (tMeet < tStart) return false;
        if (tMeet < tEnd) tEnd = tMeet;
    } else {
        if (tMeet > tEnd) return false;
        if (tMeet > tStart) tStart = tMeet;
    }
}

tx::vec2 liang_findInRangeStartEndT(const tx::vec2& start, const tx::vec2 end,
                                   const tx::vec2& rangeBottomLeft,
                                   const tx::vec2& rangeTopRight) {
    tx::vec2 dVec = end - start;
    float dx = dVec.x();
    float dy = dVec.y();

    const tx::vec2& min = rangeBottomLeft;
    const tx::vec2& max = rangeTopRight;

    float tStart = 0.0f;
    float tEnd = 1.0f;

    if (!liang_solveOneConstraint(-dx, start.x() - min.x(), tStart, tEnd)) return tx::InvalidVec;
    if (!liang_solveOneConstraint(-dy, start.y() - min.y(), tStart, tEnd)) return tx::InvalidVec;
    if (!liang_solveOneConstraint( dx, max.x() - start.x(), tStart, tEnd)) return tx::InvalidVec;
    if (!liang_solveOneConstraint( dy, max.y() - start.y(), tStart, tEnd)) return tx::InvalidVec;

    return tx::vec2{ tStart, tEnd };
}
```

## 示例代碼

```c++
// Liang--Barsky line-clipping algorithm
#include<iostream>
#include<graphics.h>
#include<math.h>
using namespace std;
// this function gives the maximum
float maxi(float arr[],int n) {
 float m = 0;
 for (int i = 0; i < n; ++i)
   if (m < arr[i])
     m = arr[i];
 return m;
}
// this function gives the minimum
float mini(float arr[], int n) {
 float m = 1;
 for (int i = 0; i < n; ++i)
   if (m > arr[i])
     m = arr[i];
 return m;
}
void liang_barsky_clipper(float xmin, float ymin, float xmax, float ymax,
                         float x1, float y1, float x2, float y2) {
 // defining variables
 float p1 = -(x2 - x1);
 float p2 = -p1;
 float p3 = -(y2 - y1);
 float p4 = -p3;
 float q1 = x1 - xmin;
 float q2 = xmax - x1;
 float q3 = y1 - ymin;
 float q4 = ymax - y1;
 float posarr[5], negarr[5];
 int posind = 1, negind = 1;
 posarr[0] = 1;
 negarr[0] = 0;
 rectangle(xmin, 467 - ymin, xmax, 467 - ymax); // drawing the clipping window
 if ((p1 == 0 && q1 < 0) || (p3 == 0 && q3 < 0)) {
     outtextxy(80, 80, "Line is parallel to clipping window!");
     return;
 }
 if (p1 != 0) {
   float r1 = q1 / p1;
   float r2 = q2 / p2;
   if (p1 < 0) {
     negarr[negind++] = r1; // for negative p1, add it to negative array
     posarr[posind++] = r2; // and add p2 to positive array
   } else {
     negarr[negind++] = r2;
     posarr[posind++] = r1;
   }
 }
 if (p3 != 0) {
   float r3 = q3 / p3;
   float r4 = q4 / p4;
   if (p3 < 0) {
     negarr[negind++] = r3;
     posarr[posind++] = r4;
   } else {
     negarr[negind++] = r4;
     posarr[posind++] = r3;
   }
 }
 float xn1, yn1, xn2, yn2;
 float rn1, rn2;
 rn1 = maxi(negarr, negind); // maximum of negative array
 rn2 = mini(posarr, posind); // minimum of positive array

  if (rn1 > rn2)  { // reject
    outtextxy(80, 80, "Line is outside the clipping window!");
    return;
  }

 xn1 = x1 + p2 * rn1;
 yn1 = y1 + p4 * rn1; // computing new points
 xn2 = x1 + p2 * rn2;
 yn2 = y1 + p4 * rn2;
 setcolor(CYAN);
 line(xn1, 467 - yn1, xn2, 467 - yn2); // the drawing the new line
 setlinestyle(1, 1, 0);
 line(x1, 467 - y1, xn1, 467 - yn1);
 line(x2, 467 - y2, xn2, 467 - yn2);
}
int main() {
 cout << "\nLiang-barsky line clipping";
 cout << "\nThe system window outlay is: (0,0) at bottom left and (631, 467) at top right";
 cout << "\nEnter the co-ordinates of the window(wxmin, wxmax, wymin, wymax):";
 float xmin, xmax, ymin, ymax;
 cin >> xmin >> ymin >> xmax >> ymax;
 cout << "\nEnter the end points of the line (x1, y1) and (x2, y2):";
 float x1, y1, x2, y2;
 cin >> x1 >> y1 >> x2 >> y2;
 int gd = DETECT, gm;
 // using the winbgim library for C++, initializing the graphics mode
 initgraph(&gd, &gm, "");
 liang_barsky_clipper(xmin, ymin, xmax, ymax, x1, y1, x2, y2);
 getch();
 closegraph();
}
```

## 外部連結

- http://hinjang.com/articles/04.html#eight
- Skytopia: The Liang-Barsky line clipping algorithm in a nutshell!