Milo的编程
首发于Milo的编程
用 C 语言画光(三):形状

用 C 语言画光(三):形状

再前两篇中,我们一直只用到圆形作为基础形状,本文再介绍一些常用形状的 SDF。但这些形状也只是作为例子,重要的是掌握如何推导出这些 SDF。

1. 圆形

我们先重温一下圆形的 SDF:

\phi_\text{circle}(\mathbf{x}) = \left \| \mathbf{x} - \mathbf{c} \right \|-r\\

当中 \left \| \mathbf{x} - \mathbf{c} \right \| 就是坐标 \mathbf{x} 与一个点 \mathbf{c} (圆心)的欧氏距离(Euclidean distance)。把这个距离减去 r ,实际上等同于把一点向外扩展 r 的范围。其实,对于任何形状的 SDF,减去一个常数 r ,都等同于把该形状它向外扩展 r 的范围。本文稍后也会看到更多应用例子。

在数学上,我们还可考虑这是任何形状与(圆心位于原点的)圆盘的闵可夫斯基和(Minkowski sum)。

2. 平面

除了圆形,最简单的形状应该是平面(plane)。在二维中这个形状应该称作直线,或者更精确地应称作二维的半空间(half-space),就是把二维空间用直线切割成两个无限大面积的部分。我们这里暂且称这个形状为平面吧。

任意维度空间的(超)平面都可以用平面上的一点 \mathbf{p} 与法线(normal vector) \mathbf{n} 来定义,由于平面与法线垂直,平面上的任意一个矢量 \mathbf{x} - \mathbf{p} 与法线的点积为零:

(\mathbf{x} - \mathbf{p}) \cdot \mathbf{n} = 0\\

如果法矢量为单位矢量 \hat{\mathbf{n}} ,空间中任意一点与平面的带符号距离为:

\phi_\text{plane}(\mathbf{x})=(\mathbf{x}-\mathbf{p})\cdot\hat{\mathbf{n}}\\

我们还可把点积展开成这种形式:

\phi_\text{plane}(\mathbf{x})=Ax + By + C\\

但为了方便使用,我们的实现用点和法线作为参数:

float planeSDF(float x, float y, float px, float py, float nx, float ny) {
    return (x - px) * nx + (y - py) * ny;
}

Result scene(float x, float y) {
    Result b = { planeSDF(x, y, 0.0f, 0.5f, 0.0f, 1.0f), 0.8f };
    return b;
}

这个例子中,上半部分是发射 0.8 自发光的形状:

我们可以用 CSG,例如把圆形与上面的平面相交,得出半圆:

Result scene(float x, float y) {
    Result a = { circleSDF(x, y, 0.5f, 0.5f, 0.2f), 1.0f };
    Result b = {  planeSDF(x, y, 0.0f, 0.5f, 0.0f, 1.0f), 0.8f };
    return intersectOp(a, b);
}

3. 胶囊形

在《用 C 语言画直线》已使用到胶囊形(capsule)的 SDF,在这里简单推导一下。

我们把胶囊形以两点 \mathbf{a}, \mathbf{b} 和半径 r 来表示:

有些场合可能每使用一点 \mathbf{a} 及矢量 \mathbf{u} = \mathbf{b} - \mathbf{a} 去表示,实际上我们在计算其 SDF 时也需要做这个减法。

计算胶囊形的 SDF,等价于计算一点与线段 ab 的距离,再扩大半径 r (即把结果减去 r )。而计算一点与线段的距离的方法,是把点 \mathbf x 投影在线段的延长线上,然后把投影点 \mathbf x' 控制在线段的范围内:

\begin{align} \mathbf{x}' &= \mathbf{a} + \max \left( \min \left ( \mathbf{v} \cdot \frac{\mathbf{u}}{\|\mathbf{u}\|}, \|\mathbf{u}\| \right), 0 \right) \frac{\mathbf{u}}{\|\mathbf{u}\|}\\ &= \mathbf{a} + \max \left( \min \left ( \mathbf{v} \cdot \frac{\mathbf{u}}{\|\mathbf{u}\|^2}, 1 \right), 0 \right) \mathbf{u}\\ &= \mathbf{a} + \max \left( \min \left ( \frac{ \mathbf{v} \cdot\mathbf{u}}{\mathbf{u}\cdot\mathbf{u}}, 1 \right), 0 \right) \mathbf{u} \end{align}\\

上式把最后的 1 / \|\mathbf{u}\| 因子放进括号里,减少开方运算。然后,所需的距离为:

d = \|\mathbf{p} - \mathbf{p}'\|\\

这里我们把线段和胶囊形的 SDF 分开为两个函数实现:

float segmentSDF(float x, float y, float ax, float ay, float bx, float by) {
    float vx = x - ax, vy = y - ay, ux = bx - ax, uy = by - ay;
    float t = fmaxf(fminf((vx * ux + vy * uy) / (ux * ux + uy * uy), 1.0f), 0.0f);
    float dx = vx - ux * t, dy = vy - uy * t;
    return sqrtf(dx * dx + dy * dy);
}

float capsuleSDF(float x, float y, float ax, float ay, float bx, float by, float r) {
    return segmentSDF(x, y, ax, ay, bx, by) - r;
}

Result scene(float x, float y) {
    Result c = {  capsuleSDF(x, y, 0.4f, 0.4f, 0.6f, 0.6f, 0.1f), 1.0f };
    return c;
}

结果:


4. 矩形

这里实现的矩形是二维定向矩形包围盒(oriented bounding box, OBB),可以由中心点 \mathbf{c} 、旋转角 \theta 和半长(half length) \mathbf{s} 来表示:

首先我们要把点 \mathbf{x} 变换至 OBB 的坐标系。我们可反过来思考,如果从 OBB 的坐标系的一点,转换至世界坐标系:

\mathbf{x} = \begin{bmatrix} \cos \theta & -\sin\theta\\ \sin \theta & \cos\theta \end{bmatrix}\mathbf{x}'+\mathbf{c}\\

那么,其逆变换就是(旋转矩阵是正交矩阵,有 Q^{-1}=Q^\mathrm{T} ):

\mathbf{x}' = \begin{bmatrix} \cos \theta & \sin\theta\\ -\sin \theta & \cos\theta \end{bmatrix}(\mathbf{x}-\mathbf{c})\\

因为矩形的对称性,我们只需考虑第一象限,可用绝对值把 \mathbf{x}' 移到第一象限。然后,要对内和外分计算距离,再取得最少值。

float boxSDF(float x, float y, float cx, float cy, float theta, float sx, float sy) {
    float costheta = cosf(theta), sintheta = sinf(theta);
    float dx = fabs((x - cx) * costheta + (y - cy) * sintheta) - sx;
    float dy = fabs((y - cy) * costheta - (x - cx) * sintheta) - sy;
    float ax = fmaxf(dx, 0.0f), ay = fmaxf(dy, 0.0f);
    return fminf(fmaxf(dx, dy), 0.0f) + sqrtf(ax * ax + ay * ay);
}

Result scene(float x, float y) {
    Result d = { boxSDF(x, y, 0.5f, 0.5f, TWO_PI / 16.0f, 0.3f, 0.1f), 1.0f };
    return d;
}

结果:

减去一个半径(如 0.1)就获得圆角长方形:

Result e = { boxSDF(x, y, 0.5f, 0.5f, TWO_PI / 16.0f, 0.3f, 0.1f) - 0.1f, 1.0f };
return e;

5. 三角形

三角形的 SDF 需计算点至三个线段的最短距离,然后再判断点在内还是在外。这里我们用缠绕顺序(winding order)来定义三角形的内外。详细的判断方法请参考 [1]。

float triangleSDF(float x, float y, float ax, float ay, float bx, float by, float cx, float cy) {
    float d = fminf(fminf(
        segmentSDF(x, y, ax, ay, bx, by),
        segmentSDF(x, y, bx, by, cx, cy)),
        segmentSDF(x, y, cx, cy, ax, ay));
    return (bx - ax) * (y - ay) > (by - ay) * (x - ax) && 
           (cx - bx) * (y - by) > (cy - by) * (x - bx) && 
           (ax - cx) * (y - cy) > (ay - cy) * (x - cx) ? -d : d;
}

Result scene(float x, float y) {
    Result f = { triangleSDF(x, y, 0.5f, 0.2f, 0.8f, 0.8f, 0.3f, 0.6f), 1.0f };
    return f;
}

结果:

减去一个半径就获得圆角三角形:

6. 结语

本文介绍了 5 种二维形状的 SDF,当中用到不同的技巧,例如运用形状的对称性质。长方形实际上也可用三角形的做法实现,然而不考虑对称性,所需的运算会更多。

读者可利用这些形状及 CSG 来绘画各种形状的光,也可以尝试推导其他形状的 SDF。

我们之后会再讨论用 SDF 的其他技巧,但接下来的两篇将先回到光学相关内容──反射折射

本文的源程序位于 shapes.cm.c

参考

[1] Dan Sunday, Inclusion of a Point in a Polygon.

编辑于 2017-12-11

文章被以下专栏收录