梁友栋-柏世奇算法

梁友栋—柏世奇算法（以梁友栋和柏世奇的名字命名）是计算机图形学中的一个线段裁剪算法。梁友栋—柏世奇算法使用直线的参数方程和不等式组来描述线段和裁剪窗口的交集。求解出的交集将被用于获知线的哪些部分是应当绘制在屏幕上的。这一算法比科恩－苏泽兰算法（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!