2D基本图形的Sign Distance Function (SDF)详解（上）

前言

符号距离函数（sign distancefunction），简称SDF，又可以称为定向距离函数（oriented distance function），在空间中的一个有限区域上确定一个点到区域边界的距离并同时对距离的符号进行定义：点在区域边界内部为正，外部为负，位于边界上时为0。

SDF在光线追踪领域有很重要的应用。实际上写这篇博客是受到了 闫令琪 大神的课程B站视频-GAMES101-现代计算机图形学入门的影响，他在课程中分享了一个神奇的网站 Shadertoy ，并且将另一位大神 inigo quilez 的作品进行了展示，为我打开了一扇新世界的大门。实际上你能在iq大神自己写的文章：raymarching distance fields，中看到他利用raymarching技术生成的许多精美的图像。这里不得不再次表达对shadertoy网站感激之情，这个网站伟大的地方就在于你能从中完完整整地看到那些绚丽夺目的特效背后的代码，并且你可以对它们进行修改然后实时编译查看效果，这一切仅需要你有一点GLSL语言基础以及一定的数学能力（数学令人头大），好在主要是几何学的知识。
咳咳扯远了，iq大神也提到，SDF在raymarching发挥了重要作用，因此他自己也有两篇博客分别列出了2D以及3D基本图形的SDF，问题就在于iq大神这个人吧，一看他就是大神，代码写得都那么简洁凝练，追求高效。能不写if分支就绝不写，能一个语句处理五种复杂情况就绝不写哪怕一句多余的话。这样的代码确实看起来短小精悍优美，但是在缺乏注释的情况下对于阅读者学习者来说可就是一种折磨了。因此本文实际上就是对其中的2D篇中各个函数的注释和分析，代码不是我写的（我也写不出(lll￢ω￢)），光是看懂就耗费了很长的时间了。
另外，如果您对计算机图形学方向不感兴趣甚至根本就不是程序员出身，那本文抛开代码也会提及纯数学意义上的“如何计算平面上任意一点到一个给定图形之间的最短距离”这个问题，希望对您有所帮助O(∩_∩)O
废话不多说，让我们开始吧！

实践出真知

希望您在阅读本文时，可以同时打开这篇文章：2D distance functions, 并且打开iq大神在每个函数后面对应在shadertoy网站上面的实现，这样直观地感受到每一条语句它所发挥的实际作用，您可以动态调整某些数据以观察其改变，这样能帮您更好地理解该函数的实现思路。当然，手边的纸笔在理解数学意义上或许能带来更大的帮助。

基础2D图形的SDF

1. 圆形

代码：（注：代码中传入的参数p在每个函数中都表示需要计算最短距离的平面上的任意一点）

/**
* 圆形：  1. 原点位于中心点
*        2. r表示半径
*/
float sdCircle( vec2 p, float r )
{// 与圆心距离位r的点，在该圆上，SDF取值0return length(p) - r;
}

圆形几乎是最简单的2D图形了，它的定义就是与圆心的距离等于半径的所有点的集合。接下来我们会经常见到下面这种风格的图片，实际上他就是iq大神在展示自己的2D距离函数的效果，我们可以看到白色的线条连起来的就是和目标图形距离为0的点，也就是在图形上的点。而蓝白色就是图形内部点，与图形的最短距离为负数，外面黑黄的部分则是图形外部的点，距离为正。每一个封闭的圈都是一条等距离线。

2.线段

代码：

/**
* 线段：  1. a，b表示线段两个端点的坐标
*/
float sdSegment( in vec2 p, in vec2 a, in vec2 b )
{// pa表示a点指向p点的向量， ba表示a点指向b点的向量vec2 pa = p-a, ba = b-a;// h表示pa在ba上投影的长度占ba长度的比例，限定到[0,1]float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );// ba*h 是以a为起点，方向与ba相同，长度等于pa在ba方向上投影的长度的向量// pa视为斜边，ba*h是直角边，那么pa - ba*h则是另一条直角边，也就是从p点做垂线垂直于ab，显然该垂线的长度就是所求最短距离return length( pa - ba*h );
}

实际上虽然线段看起来简单，但他确实接下来大多数图形距离函数的基础，因为总有一些图形它拥有一条直边，那时计算最短距离时肯定要求p点到某一条线段的距离。那么我们现在来解释一下，对于任意一条线段，我们可以按照下面这样把它所在的平面划分为三个区域，边界就是两条过两个端点垂直于线段的直线，划分的标准是什么呢？就是平面内任意点P在区域1内时，线段AB上与它最短的距离点（暂且称为C点）始终是A点，最短距离的长度是|AP|，或者说，长度可以写成|AP - AC|。区域3也是同样，点C始终是B点，长度是|AP - AB|或者写成|AP - AC|。而在区域2内很显然，是c点，也就是点p在ab方向上的投影，最短距离是|AP - AC|。行了，我们将三个区域的表示方法统一了（实际上这是iq大神的惯用手段）。
确定AC是什么很关键，从下图看很容易知道其实AC就是AB乘一个系数得到的，p点在在区域1时该系数始终是0，区域3内时该系数是1，区域2内该系数就是AP在AB上的投影占AB的比例，这其实就是clamp函数的作用了。剩下的其实就不难理解，我觉得关于线段的SDF这一个函数应该算是iq大神写的比较容易理解的代码，其他地方他会直接把第二和第三句合在一起写（心累）。之所以这里要啰嗦这么多，完全是为了给下面其他的几何图形打基础，因为这个逻辑会经常出现，clamp函数也会经常出场的。

最终结果如下图，这里又得啰嗦几句了，大家观察下图白线外面的圈像什么？实际上就是胶囊体。线段的SDF结果只要再减去一个常数就可以得到一个胶囊体，这个神奇的特性下面会分析。

3. 长方形

代码：

/**
* 长方形 box：  1. 原点位于长方形的中心点，形状是轴对称的
*               2. b表示长方形右上角顶点的坐标
*/
float sdBox( in vec2 p, in vec2 b )
{// abs(p)是常用技巧，由于该图形四个象限都是相同的，因此都映射到第一象限即可// 现在的d表示长方体右上角顶点直线p点的向量vec2 d = abs(p)-b;// p点在外部：length(max(d,0.0)), 在内部则是min(max(d.x,d.y),0.0), 这两项总至少有一项为0return length(max(d,0.0)) + min(max(d.x,d.y),0.0);
}

同样的思路，我们将依据图形上与P点最短距离位置点选择的不同，将第一象限划分（其他象限都映射到第一象限即可）四个区域，其中三个属于长方形外部，一个属于内部。我们先来看看外部的SDF怎么计算，观察下图我们可以看到代码中的d其实就是四条红色的向量（区域1内的和绿色向量重叠了），落在区域2内的点距离函数取d.y即可，因为d.x是负数，落在区域4内的点取d.x，因为d.y是负数，负数就表示你在那一个方向上处于长方形的“内部”，若两个分量都是正数，那就取|d|即可。三种情况统一到一句代码length(max(d, 0.0))里面了。
看长方形内部的点，也就是区域3，此时的d两个分量都是负数，那只要选择绝对值小的那个分量即可（绿色的向量），即max(d.x, d.y))。那加个min()是什么意思呢？还不是因为iq大神想要把这四个区域的判定都挤压在一条语句内完成嘛~

再次注意观察下面的图形，如果我们想要得到一个圆角的长方形的SDF该怎么做呢？没错，只需要将普通长方形的SDF结果减去一个常数，也就是圆角的半径即可。你看白圈外面的等高线不就是一个带圆角的长方形吗？

4. 菱形

代码：

/**
* 菱形：  1. 原点在菱形的中心点，四个顶点都在坐标轴上
*         2. b.x = 与x轴正半轴交点， b.y = 与y轴正半轴交点
*/
float ndot(vec2 a, vec2 b ) { return a.x*b.x - a.y*b.y; }
float sdRhombus( in vec2 p, in vec2 b)
{vec2 q = abs(p);  // 坐标轴对称// 计算的是线段b的中点到p点的向量，在b上的投影限制到[-1, 1]// 负数表示向量偏向于y轴， 正数表示偏向x轴float h = clamp((-2.0*ndot(q,b)+ndot(b,b))/dot(b,b),-1.0, 1.0);/* 实际上h从-1到1的滑动过程，[0.5*b*vec2(1.0-h,1.0+h)] 表示一条从原点出发，终点在向量b上由左上到右下的滑动的向量 */float d = length( q - 0.5*b*vec2(1.0-h,1.0+h) );// 符号：可计算(b.x, -b.y)和(p.x, p.y-b.y)的叉积，得两向量的相对位置return d * sign( q.x*b.y + q.y*b.x - b.x*b.y );
}

这个函数可能乍一看会让人懵逼，实际上我也觉得有点过于复杂且没必要了。它的核心思想其实还是和上面的线段的SDF一模一样，但是在线段的SDF中，我们是选取了线段的一个端点和目标点P进行连线然后投影，但是在这里iq大神不选择端点了，选择了一个中点，然后继续投影，那我们知道现在就不应该再限制在[0, 1]，而是[-1, 1]，为了处理选择中点带来的影响，下一句计算d时也得进行一步看上去比较绕的操作，不过核心思路没变，那就是当P点在区域1是，P点的投影点应该固定为A点是吧，我们来看看0.5*b*vec2(1.0-h,1.0+h)，当h取-1时，这句表达式确实得到A点坐标，那就行了。区域3同理。至于符号的计算，可以计算P和某一端点（代码里选择了A点）的连线，然后看OP在AP的左边还是右边即可，使用叉积的正负可以用来判断。最终化简结果就像return后面代码写的那样。
可能有些同学一开始会疑惑定义的ndot函数的几何意义，起码在这里它没表现出什么特殊的几何意义，如果你手算公式，会发现那仅仅是一个简化后的运算结果而已。真要说它有什么几何意义其实也有，下文中“一般三角形”的部分会做讨论。

我自己使用前面的线段SDF思路写了一个菱形的简单实现，思路可能更清晰一点，结果是一样的：

float sdRhombus( in vec2 p, in vec2 b)
{// 参考线段的SDFvec2 a = vec2(0, b.y); // 与y轴交点vec2 c = vec2(b.x, 0);vec2 pa = q-a, ba = c-a;float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );float d = length( pa - ba*h );return d * sign( q.x*b.y + q.y*b.x - b.x*b.y );
}

5. 等边三角形

代码：

/*** 等边三角形：  1. 原点在三角形中心点*              2. r表示边长的一半（不是中心点到顶点的距离）
*/
float sdEquilateralTriangle(in vec2 p, in float r )
{const float k = sqrt(3.0);p.x = abs(p.x);// 由于三角形可以看成是由三个部分绕中心点三个角度生成的，x负半轴不用考虑// 所以只需要考虑在区域 1的那部分，这一部分可以通过关于直线y=-1/sqrt(3)*x对称映射到下方的三分之一// 映射关系如下：if( p.x+k*p.y>0.0 ) p=vec2(p.x-k*p.y,-k*p.x-p.y)/2.0;// 一切都归一到只需要看区域 2-3-4即可// x在减去r后，若其值-2*r<x<0，说明原来范围是-r<0<r,即知道点p在x轴上处于三角形内部，只需看y轴长度即可// 否则，需要原x轴坐标减去左/右下角坐标x值p.x -= r;p.y += r/k;p.x -= clamp( p.x, -2.0*r, 0.0 );return -length(p)*sign(p.y);
}

首先提一句，对于这些正多边形，后面会给出比较通用的一个函数来处理。这里单把正三角拿出来可能是为了图一乐（