Have stroke geometry shader pass to a coordinate system where the curve is some segment of y = x^2

manimlib/shaders/inserts/quadratic_bezier_geometry_functions.glsl

@@ -3,22 +3,77 @@ float cross2d(vec2 v, vec2 w){
 }
 
 
-mat3 get_xy_to_uv(vec2 b0, vec2 b1){
+vec2 complex_div(vec2 v, vec2 w){
-    mat3 shift = mat3(
+    return vec2(dot(v, w), cross2d(w, v)) / dot(w, w);
+}
+
+
+vec2 xs_on_clean_parabola(vec2 controls[3]){
+    /*
+    Given three control points for a quadratic bezier,
+    this returns the two values (x0, x2) such that the
+    section of the parabola y = x^2 between those values
+    is isometric to the given quadratic bezier.
+
+    Adapated from https://github.com/raphlinus/raphlinus.github.io/blob/master/_posts/2019-12-23-flatten-quadbez.md
+    */
+    vec2 b0 = controls[0];
+    vec2 b1 = controls[1];
+    vec2 b2 = controls[2];
+
+    vec2 dd = normalize(2 * b1 - b0 - b2);
+
+    float u0 = dot(b1 - b0, dd);
+    float u2 = dot(b2 - b1, dd);
+    float cp = cross2d(b2 - b0, dd);
+
+    return vec2(u0 / cp, u2 / cp);
+}
+
+
+mat3 map_point_pairs(vec2 src0, vec2 src1, vec2 dest0, vec2 dest1){
+    /*
+    Returns an orthogonal matrix which will map
+    src0 onto dest0 and src1 onto dest1.
+    */
+    mat3 shift1 = mat3(
         1.0, 0.0, 0.0,
         0.0, 1.0, 0.0,
-        -b0.x, -b0.y, 1.0
+        -src0.x, -src0.y, 1.0
+    );
+    mat3 shift2 = mat3(
+        1.0, 0.0, 0.0,
+        0.0, 1.0, 0.0,
+        dest0.x, dest0.y, 1.0
     );
 
-    float sf = length(b1 - b0);
+    // Compute complex division dest_vect / src_vect to determine rotation
-    vec2 I = (b1 - b0) / sf;
+    vec2 complex_rot = complex_div(dest1 - dest0, src1 - src0);
-    vec2 J = vec2(-I.y, I.x);
     mat3 rotate = mat3(
-        I.x, J.x, 0.0,
+        complex_rot.x, complex_rot.y, 0.0,
-        I.y, J.y, 0.0,
+        -complex_rot.y, complex_rot.x, 0.0,
         0.0, 0.0, 1.0
     );
-    return (1.0 / sf) * rotate * shift;
+
+    return shift2 * rotate * shift1;
+}
+
+
+mat3 get_xy_to_uv(vec2 controls[3], float bezier_degree){
+    vec2[2] dest;
+    if (bezier_degree == 1.0){
+        dest[0] = vec2(0, 0);
+        dest[1] = vec2(1, 0);
+    }else{
+        vec2 xs = xs_on_clean_parabola(controls);
+        float x0 = xs.x;
+        float x2 = xs.y;
+        dest[0] = vec2(x0, x0 * x0);
+        dest[1] = vec2(x2, x2 * x2);
+    }
+    return map_point_pairs(
+        controls[0], controls[2], dest[0], dest[1]
+    );
 }
 
 

manimlib/shaders/quadratic_bezier_stroke/frag.glsl

@@ -3,44 +3,72 @@
 #INSERT camera_uniform_declarations.glsl
 
 in vec2 uv_coords;
-in vec2 uv_b2;
 
 in float uv_stroke_width;
-in vec4 color;
 in float uv_anti_alias_width;
+in vec4 color;
 
-in float has_prev;
-in float has_next;
 in float bezier_degree;
 
 out vec4 frag_color;
 
-
+const float QUICK_DIST_WIDTH = 0.1;
-#INSERT quadratic_bezier_distance.glsl
 
 
-float dist_to_curve(){
+float cube_root(float x){
-    float dist = min_dist_to_curve(uv_coords, uv_b2, bezier_degree);
+    return sign(x) * pow(abs(x), 1.0 / 3.0);
-    if(has_prev == 0 && uv_coords.x < 0){
-        float buff = 0.5 * uv_stroke_width - uv_anti_alias_width;
-        return max(-uv_coords.x + buff, dist);
 }
-    if(has_next == 0 && uv_coords.x > uv_b2.x){
+
-        float buff = 0.5 * uv_stroke_width - uv_anti_alias_width;
+
-        vec2 v12 = normalize(uv_b2 - vec2(1, 0));
+// Distance from (x0, y0) to the curve y = x^2
-        float perp_dist = dot(uv_coords - uv_b2, v12);
+float dist_to_curve(float x0, float y0){
-        if (perp_dist > 0) return max(perp_dist + buff, dist);
+    if(bezier_degree == 1.0){
+        return y0;
     }
-    return dist;
+    if(false && uv_stroke_width < QUICK_DIST_WIDTH){
+        // This is a quick approximation for computing
+        // the distance to the curve.
+        // Evaluate F(x, y) = y - x^2
+        // divide by its gradient's magnitude
+        return (y0 - x0 * x0) / sqrt(1 + 4 * x0 * x0);
+    }
+    // Otherwise, explicit solve for the minmal distance using the cubic formula
+    //
+    // The distance squared between (x0, y0) and a point (x, x^2) looks like
+    //
+    // (x0 - x)^2 + (y0 - x^2)^2 = x^4 + (1 - 2y0)x^2 - 2x0 * x + (x0^2 + y0^2)
+    //
+    // Setting the derivative equal to zero (and rescaling) looks like
+    // 
+    // x^3 + (0.5 - y0) * x - 0.5 * x0 = 0
+    //
+    // float p = 0.5 - y0;
+    // float mhq = 0.25 * x0;  // negative half of q
+    // float sqrt_disc = sqrt(mhq * mhq + p * p * p / 27.0);
+    // float x = cube_root(mhq + sqrt_disc) + cube_root(mhq - sqrt_disc);
+    // return distance(uv_coords, vec2(x, x * x));
+
+    float x = x0;
+    float p = (0.5 - y0);
+    float q = -0.5 * x0;
+    for(int i = 0; i < 2; i++){
+        float fx = x * x * x + p * x + q;
+        float dfx = 3 * x * x + p;
+        x = x - fx / dfx;
+    }
+    return distance(uv_coords, vec2(x, x * x));
 }
 
 
 void main() {
     if (uv_stroke_width == 0) discard;
 
+    float x0 = uv_coords.x;
+    float y0 = uv_coords.y;
     // An sdf for the region around the curve we wish to color.
-    float signed_dist = abs(dist_to_curve()) - 0.5 * uv_stroke_width;
+    float signed_dist = abs(dist_to_curve(x0, y0)) - 0.5 * uv_stroke_width;
 
     frag_color = color;
+    // if(uv_stroke_width > QUICK_DIST_WIDTH) frag_color = vec4(1, 0, 0, 1);
     frag_color.a *= smoothstep(0.5, -0.5, signed_dist / uv_anti_alias_width);
 }

manimlib/shaders/quadratic_bezier_stroke/geom.glsl

@@ -31,12 +31,11 @@ out vec4 color;
 out float uv_stroke_width;
 out float uv_anti_alias_width;
 
-out float has_prev;
+// out float has_prev;
-out float has_next;
+// out float has_next;
 out float bezier_degree;
 
 out vec2 uv_coords;
-out vec2 uv_b2;
 
 // Codes for joint types
 const float AUTO_JOINT = 0;
@@ -103,17 +102,15 @@ int get_corners(
     float aaw = anti_alias_width * frame_shape.y / pixel_shape.y;
     float buff0 = 0.5 * stroke_widths[0] + aaw;
     float buff2 = 0.5 * stroke_widths[2] + aaw;
-    float aaw0 = (1 - has_prev) * aaw;
-    float aaw2 = (1 - has_next) * aaw;
 
-    vec2 c0 = p0 - buff0 * p0_perp + aaw0 * v10;
+    vec2 c0 = p0 - buff0 * p0_perp;
-    vec2 c1 = p0 + buff0 * p0_perp + aaw0 * v10;
+    vec2 c1 = p0 + buff0 * p0_perp;
-    vec2 c2 = p2 + buff2 * p2_perp + aaw2 * v12;
+    vec2 c2 = p2 + buff2 * p2_perp;
-    vec2 c3 = p2 - buff2 * p2_perp + aaw2 * v12;
+    vec2 c3 = p2 - buff2 * p2_perp;
 
     // Account for previous and next control points
-    if(has_prev > 0) create_joint(angle_from_prev, v01, buff0, c0, c0, c1, c1);
+    create_joint(angle_from_prev, v01, buff0, c0, c0, c1, c1);
-    if(has_next > 0) create_joint(angle_to_next, v21, buff2, c3, c3, c2, c2);
+    create_joint(angle_to_next, v21, buff2, c3, c3, c2, c2);
 
     // Linear case is the simplest
     if(degree == 1){
@@ -137,10 +134,15 @@ void main() {
     bezier_degree = (abs(v_joint_angle[1]) < ANGLE_THRESHOLD) ? 1.0 : 2.0;
     vec3 unit_normal = camera_rotation * vec3(0.0, 0.0, 1.0); // TODO, track true unit normal globally
 
-    // Adjust stroke width based on distance from the camera
+    // Control points are projected to the xy plane before drawing, which in turn
+    // gets tranlated to a uv plane.  The z-coordinate information will be remembered
+    // by what's sent out to gl_Position, and by how it affects the lighting and stroke width
+    vec2 flat_controls[3];
     float scaled_strokes[3];
     for(int i = 0; i < 3; i++){
         float sf = perspective_scale_factor(verts[i].z, focal_distance);
+        flat_controls[i] = sf * verts[i].xy;
+
         if(bool(flat_stroke)){
             vec3 to_cam = normalize(vec3(0.0, 0.0, focal_distance) - verts[i]);
             sf *= abs(dot(unit_normal, to_cam));
@@ -148,15 +150,6 @@ void main() {
         scaled_strokes[i] = v_stroke_width[i] * sf;
     }
 
-    // Control points are projected to the xy plane before drawing, which in turn
-    // gets tranlated to a uv plane.  The z-coordinate information will be remembered
-    // by what's sent out to gl_Position, and by how it affects the lighting and stroke width
-    vec2 flat_controls[3];
-    for(int i = 0; i < 3; i++){
-        float sf = perspective_scale_factor(verts[i].z, focal_distance);
-        flat_controls[i] = sf * verts[i].xy;
-    }
-
     // If the curve is flat, put the middle control in the midpoint
     if (bezier_degree == 1.0){
         flat_controls[1] = 0.5 * (flat_controls[0] + flat_controls[2]);
@@ -165,15 +158,13 @@ void main() {
     // Set joint information
     float angle_from_prev = v_joint_angle[0];
     float angle_to_next = v_joint_angle[2];
-    has_prev = 1.0;
-    has_next = 1.0;
     if(angle_from_prev == DISJOINT_CONST){
+        // TODO, mark the fact that there is no previous
         angle_from_prev = 0.0;
-        has_prev = 0.0;
     }
     if(angle_to_next == DISJOINT_CONST){
+        // TODO, mark the fact that there is no next
         angle_to_next = 0.0;
-        has_next = 0.0;
     }
 
     // Corners of a bounding region around curve
@@ -187,16 +178,15 @@ void main() {
     int index_map[5] = int[5](0, 0, 1, 2, 2);
     if(n_corners == 4) index_map[2] = 2;
 
-    // Find uv conversion matrix
+    // Find uv conversion
-    mat3 xy_to_uv = get_xy_to_uv(flat_controls[0], flat_controls[1]);
+    mat3 xy_to_uv = get_xy_to_uv(flat_controls, bezier_degree);
-    float scale_factor = length(flat_controls[1] - flat_controls[0]);
+    float scale_factor = length(xy_to_uv[0].xy);
-    uv_anti_alias_width = anti_alias_width * frame_shape.y / pixel_shape.y / scale_factor;
+    uv_anti_alias_width = scale_factor * anti_alias_width * (frame_shape.y / pixel_shape.y);
-    uv_b2 = (xy_to_uv * vec3(flat_controls[2], 1.0)).xy;
 
     // Emit each corner
     for(int i = 0; i < n_corners; i++){
         uv_coords = (xy_to_uv * vec3(corners[i], 1.0)).xy;
-        uv_stroke_width = scaled_strokes[index_map[i]] / scale_factor;
+        uv_stroke_width = scale_factor * scaled_strokes[index_map[i]];
         // Apply some lighting to the color before sending out.
         vec3 xyz_coords = vec3(corners[i], verts[index_map[i]].z);
         color = finalize_color(