2D SDF(4): Regular Polygon - OhBonsai/art GitHub Wiki

Why learn SDF

Initially, I wanna to derive the SDF because I came arcoss the following code online, which i couldn't understand at all, but yet this code could magically draw a regular ploygon.

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;
}

Until now, I still don't understand this piece of code... .But through derivation, I've written a more elgant and understanable SDF function : ).

Implementation

The key to implementing a generic SDF for regular ploygons is using **Polar corrdinates. **Polar coordinates are a two-dimensional coordinate system that specifies the position of a point on a plane based on the angle and the distance from a fixed central point known as the pole (often corresponding to the origin in Cartesian coordinates). Unlike Cartesian coordinates, which use a grid of vertical and horizontal lines to specify a point in space with x (horizontal) and y (vertical) coordinates, polar coordinates represent a point with the following two values:

  1. Radius (r): The distance from the pole to the point. It is a non-negative value.
  2. Angle (θ): The angle measured from a fixed direction (usually the positive x-axis in Cartesian coordinates), going anti-clockwise to the line segment that connects the point with the pole.

Given a regular poloygon with N sides, centered at point $O(0, 0)$ with radiusr，the circle can be divided into N pieces , Each pieces angle is $\delta$ . Any point $P(r1, \theta)$ can be allocated to one of these pieces based on $\theta$ . As illustrated, the starting angle of a circular segment is 2PI / N * 0, and the ending angle is 2PI / N * 1. Assuming the starting point as $A$ , the ending point as $A^\prime$ . A perpendicular line from $O$ to $AA^\prime$ intersects at point $D$ . SDF is the length of projection of $PO$ on $PD$ minus the length of $PD$ . Hence we can write the following 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) ;
}

In fact, we only need to know to lengths. one is the projection length of $PO$ on $PD$ . antoher is the length of $PD$ . The length of $PD$ can be directly calculated using $r \times cos(\delta / 2)$ . The projection length can be determined by the angle betwwen $PO$ and . Then code can be simplified as follows

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;
}

Finally, we get the following graphic.

Angle between vectors

Implement a generic function float angle(vec2 iAxis, vec2 v), which taks two vectors as inputs, calculates the angle between the twe vectors, and normalizes the result to the range[0, 2 * PI]

Step1: Calculate the Angle through Dot Product $$cos(\theta) = \dfrac{iAxis \bigodot v} {|iAxis| \times |v|}$$
Step2: Convert the Angle Range The caculated $\theta$ will be within the $[0, \pi]$ range, To map the angle to the $[0, 2\pi]$ range, We can add brancing conditions based on the direction of the vectors, This direction judgment is done by comparing the cross product of the vectors and zero: $$xAxis \bigotimes v >= 0.0$$
And the code is below

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;
    }
}

Fix Bug

Although we can obtain a triangle with the above method, there is no arc on the exterior of the triangle's corners, and actually, the distance field is incorrect. As shown in the illustration below, the calculated result is the red line, but the correct equidistant line should be the blue one.

Draw perpendicular lines $p$ and $q$ through points $A$ , $A^\prime$ , respectively, with the starting edge as $s$ and the ending edge as $t$ .

There are several ways to determine the position of $P$ , For example, we can compare the lengths of the two blue sgements show in the diagram below

Or we can compare the sizes of the two angles show in the diagram below

Here, we opt for the second method by using function that calculartes vector angles, $Angle(PA, OA)$ and $Angle(PA^\prime, OA^\prime)$ . Ultimately, we can derive the following code.

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 :)

Symmetrical Point

Implement a generic function vec2 reflectPoint(vec2 p, vec2 v), which takes a vector v and a point p, and calculates the reflection point pp of p relative to the vector. There are many methods of implementation.

Choose a more intuitive method, with the point positions as illustrated above, the derivation is as follows. $$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;
}

In the above SDF (Signed Distance Function), we need to calculate two areas, area1 and area2, separately. However, in reality, area1 and area2 are symmetrical along Line(OD). If for any point $P$ , we can obtain its reflection point $PP$ relative to Line(OD), we only need to calculate the area for area2.

Thus, the GLSL code can be simplified as

  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);
      }
  }