2D SDF推导4: 通用正多边形

• Shadertoy: www.shadertoy.com/view/Mc3GWs
• Geogebra: www.geogebra.org/classic/mcd…

缘起

最开始想要看一看SDF，是因为我在网上找到下面一段代码，我发现我根本看不懂，但是挺好用的

float sdf_polygon(vec2 pt, float radius, int N) {
    float a = atan(pt.x, pt.y);
    float r =  PI * 2.0 / float(N);
    if (a < 0.) {
        a += PI * 2.0;
    }
    float d = length(pt) * (1. + f4 -  radius );
    float a2 = floor(f1 * .5 + a /r)*r;
    return cos(a2 - a) *d - f2 / 2.0;
}

到现在我还是看不懂这段代码....但是通过推导，我写出了上面代码更优雅，更好理解的Shader :)

实现

实现通用正多边形的SDF的秘诀，就在于坐标系的转换，通过极坐标可以很好的处理正多边形的每一条边。极坐标有以下性质

• 基于一个固定点（极点）和从这个点发出的射线来定义点的位置。
• 任何平面上的点都通过一对坐标(r, θ)来表示：r是点到极点的距离（半径），θ是从参考方向（通常是x轴正方向）到点的连线与参考方向之间的角度。

设正N边形中心为点 $O(0, 0)$ 半径为r，根据N的大小，将圆等分成N块。每块占据的角度为 $\delta$ 任意点 $P(r1, \theta)$ 根据 $\theta$ 都可以归属到某一块中。如下图，圆角块开始边角度为2PI / N * 0，结束边角度为2PI / N * 1. 假设开始点为 $A$ ,结束点为 $A^\prime$ . 过 $O$ 做 $AA^\prime$ 的垂线有交点 $D$ . 点 $P$ 到边的距离就是 $PO$ 在 $PD$ 上的投影长度减去 $PD$ 的长度。于是我们可以写出以下的shader

float sdf_polygon1(vec2 P, float radius, int sides) {
    float angle = atan(P.y, P.x);
    // Translate the value of angle into the range (0, 2 * PI).
    angle = P.y > 0. ? angle: angle + PI * 2.;
    
    float delta = 2. * PI / float(sides);
    float areaNumber = ceil(angle / delta);
    
    
    float theta1 = delta * areaNumber;
    float theta2 = delta * (areaNumber + 1.0);
   
    
    vec2 PA       = vec2(radius * cos(theta1), radius * sin(theta1) );
    vec2 PA_Prime = vec2(radius * cos(theta2), radius * sin(theta2) );
    vec2 PD       = (PA - PA_Prime)/2.0;
        
    float projectLength = abs(dot(PD, P)) / length(PD);
    return projectLength - length(PD) ;
}

实际上我们只需要知道两个长度即可， $PO$ 在 $PD$ 上的投影长度ProjectLength 和 $PD$ 的长度。PD的长度可以直接通过 $r \times cos(\delta / 2)$ 直接计算出。 投影长度可以通过 $PO$ 与 $PD$ 的夹角求得。 整体代码可以简化为

float sdf_polygon3(vec2 P, float radius, int sides) {
    float angle = atan(P.y, P.x);
    float delta = 2. * PI / float(sides);
    float theta = mod(angle, delta) - delta / 2.0;
  
    float lengthPD         = radius * cos(delta / 2.0);
    float lengthProjection = length(P) * cos(theta)

    return lengthProjection - lengthPD;
}

最后我们就可以得到以下漂亮的图形了

向量夹角

实现一个通用函数 float angle(vec2 iAxis, vec2 v), 输入两个向量，求两个向量的夹角，并转化到[0, 2 * PI]的区间

步骤 1: 通过点积计算角度 $cos(\theta) = \dfrac{iAxis \bigodot v} {\|iAxis\| \times \|v\|}$
步骤2: 转换角度范围 算出来的 $\theta$ 将会是在 ([-1, 1]) 范围内的弧度值，这意味着角度将会在 $[0, \pi]$ 范围内。需要将角度映射到 $[0, 2\pi]$ 范围，可以根据向量的方向判断使用下面的公式添加分支条件, 方向的判断是用向量叉积与0的大小比较 $xAxis \bigotimes v >= 0.0$
最后可以得到以下代码

float angle(vec2 iAxis, vec2 v) {
    float theta = acos(dot(iAxis, v) / length(iAxis) / length(v));
    float crossV = iAxis.x * v.y - iAxis.y * v.x;
    if (crossV < 0.0) {
      return PI * 2.0 - theta;
    } else {
      return theta;
    }
}

缺陷修复

以上的图形，虽然我们可以得到一个三角形，但是在三角形的角的外部，并没有得到一个弧线，实际上这个距离场是错误的。如下图所示，我们计算的结果是红色的这条线，正确等距应该是蓝色这条线。

过 $A$ , $A^\prime$ 两个点，做 $Line(A, A^\prime)$ 的垂线 $p$ $q$ ，另外有开始边和结束边为 $s$ $t$ . 直观的讲便是，如果点 $P$ 位于 $p$ $s$ 两条线的夹角内，距离为点 $P$ 到点 $A$ 的距离

确定A的位置有很多种方法 例如比较下图中两条蓝色线段的长短

或者是比较两个角度大小

这里采用第二种方法，利用向量夹角的函数求的 $Angle(PA, OA)$ 和 $Angle(PA^\prime, OA^\prime)$ 的大小. 最后可以得到以下代码

float sdf_polygon4(vec2 P, float radius, int sides) {
    // get polar angle
    float angle = atan(P.y, P.x);
    // make angle to range [0, 2*PI]
    angle = P.y > 0. ? angle: angle + PI * 2.;

    // get each piece angle
    float delta = 2. * PI / float(sides);
    // how many pieces?
    float areaNumber = floor(angle / delta);
    
    // start angle of current piece
    float theta1 = delta * areaNumber;
    // end angle of current piece
    float theta2 = delta * (areaNumber + 1.0);
   
    
    // start point on current piece
    vec2 POINT_A       = vec2(radius * cos(theta1), radius * sin(theta1) );
    // end point on current piece
    vec2 POINT_A_Prime = vec2(radius * cos(theta2), radius * sin(theta2) );
    // the middle of startPoint and endPoint
    vec2 POINT_D       = (POINT_A + POINT_A_Prime)/2.0;
    // corrdinate center
    vec2 POINT_O       = vec2(0.0); 
    
    // start axis of current piece
    vec2 iAxis1   =  vec2(-POINT_O+POINT_A);
    // end axis of current piece
    vec2 iAxis2   =  vec2(-POINT_O+POINT_A_Prime);


    // area 1
    vec2 vector1  = vec2(P - POINT_A);
    float a1 = vector_angle(iAxis1, vector1);
    if (a1 <  (delta / 2.0)) {
        return length(vector1);
    }
    
    // area 2
    vec2 vector2  = vec2(P - POINT_A_Prime);
    float a2 = vector_angle(iAxis2, vector2);
    if ((PI2 - a2) < (delta / 2.0) ) {
        return length(vector2);
    }
    
    // area 3 
    float theta = mod(angle, delta) - delta / 2.0;
    float l1 =  length(P) * cos(theta) - length(POINT_D);
    return l1;
}

Finally we got this :)

对称点

实现一个通用函数 vec2 reflectPoint(vec2 p, vec2 v), 输入向量 $v$ 和点 $p$ , 求 $p$ 相对于向量的对称点 $pp$ . 实现的方法有很多种

选择一个较为直观的方法，点位置如上图，推导如下 $P^\prime = P + \overrightarrow{P P^\prime} = P + 2 \times \overrightarrow{P E}$
$\overrightarrow{P E} = \overrightarrow{A P} - \overrightarrow{A E}$
$\overrightarrow{AE} = dot(\overrightarrow{AP}, norm(\overrightarrow{AD})) * norm(\overrightarrow{AD})$
转化为代码为

vec2 symmetrical_point(vec2 P, vec2 V) {
    // Ensure V is a unit vector
    vec3 V_normalized = normalize(V);
    
    // Calculate the projection point H of P onto the line, essentially the projection along vector V
    // Here we use the dot product and vector normalization to find the correct length along V
    vec3 H = dot(P, V_normalized) * V_normalized;
    
    // Calculate the vector from P to H, if P is already on the line then PH_vector will be (0,0,0)
    vec3 PH_vector = H - P;
    
    // Calculate the reflective point PP
    vec3 PP = P + 2.0 * PH_vector;
    
    return PP;
}

在上面SDF函数中，我们需要分别计算area1和 area2两个区域，但实际上area1 和 area2是沿着 Line(OD) 对称的。 如果针对任一点 $P$ 可以获取到相对于Line(OD)的对成点 $PP$ 。我们只需要计算area2这个区域。

于是GLSL的代码可以简化为

  for (int i=0; i<2; i++) {
      PP = P;
      if (i==1) {     // symmetrical for area2
          PP = symmetrical_point(P, POINT_D);
      }
      
      vec2 vector1  = vec2(PP - POINT_A);
      float a1 = vector_angle(iAxis1, vector1);
      if (a1 <  (delta / 2.0)) {
          return length(vector1);
      }
  }