From c67ae08b62306d6dea19a83bcc32c572903c51e6 Mon Sep 17 00:00:00 2001 From: Grant Sanderson Date: Mon, 9 Jan 2023 16:46:38 -0800 Subject: [PATCH] Have stroke geometry shader pass to a coordinate system where the curve is some segment of y = x^2 --- .../quadratic_bezier_geometry_functions.glsl | 73 ++++++++++++++++--- .../shaders/quadratic_bezier_stroke/frag.glsl | 64 +++++++++++----- .../shaders/quadratic_bezier_stroke/geom.glsl | 52 ++++++------- 3 files changed, 131 insertions(+), 58 deletions(-) diff --git a/manimlib/shaders/inserts/quadratic_bezier_geometry_functions.glsl b/manimlib/shaders/inserts/quadratic_bezier_geometry_functions.glsl index 8cb0ce13..898ef399 100644 --- a/manimlib/shaders/inserts/quadratic_bezier_geometry_functions.glsl +++ b/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){ - mat3 shift = mat3( +vec2 complex_div(vec2 v, vec2 w){ + 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); - vec2 I = (b1 - b0) / sf; - vec2 J = vec2(-I.y, I.x); + // Compute complex division dest_vect / src_vect to determine rotation + vec2 complex_rot = complex_div(dest1 - dest0, src1 - src0); mat3 rotate = mat3( - I.x, J.x, 0.0, - I.y, J.y, 0.0, + complex_rot.x, complex_rot.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] + ); } diff --git a/manimlib/shaders/quadratic_bezier_stroke/frag.glsl b/manimlib/shaders/quadratic_bezier_stroke/frag.glsl index 9dc2f7e0..7022dc38 100644 --- a/manimlib/shaders/quadratic_bezier_stroke/frag.glsl +++ b/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; - -#INSERT quadratic_bezier_distance.glsl +const float QUICK_DIST_WIDTH = 0.1; -float dist_to_curve(){ - float dist = min_dist_to_curve(uv_coords, uv_b2, bezier_degree); - 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); +float cube_root(float x){ + return sign(x) * pow(abs(x), 1.0 / 3.0); +} + + +// Distance from (x0, y0) to the curve y = x^2 +float dist_to_curve(float x0, float y0){ + if(bezier_degree == 1.0){ + return y0; } - 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)); - float perp_dist = dot(uv_coords - uv_b2, v12); - if (perp_dist > 0) return max(perp_dist + buff, 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); } - return dist; + // 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); } \ No newline at end of file diff --git a/manimlib/shaders/quadratic_bezier_stroke/geom.glsl b/manimlib/shaders/quadratic_bezier_stroke/geom.glsl index 383ddb07..f4ff7f12 100644 --- a/manimlib/shaders/quadratic_bezier_stroke/geom.glsl +++ b/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_next; +// out float has_prev; +// 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 c1 = p0 + buff0 * p0_perp + aaw0 * v10; - vec2 c2 = p2 + buff2 * p2_perp + aaw2 * v12; - vec2 c3 = p2 - buff2 * p2_perp + aaw2 * v12; + vec2 c0 = p0 - buff0 * p0_perp; + vec2 c1 = p0 + buff0 * p0_perp; + vec2 c2 = p2 + buff2 * p2_perp; + 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); - if(has_next > 0) create_joint(angle_to_next, v21, buff2, c3, c3, c2, c2); + create_joint(angle_from_prev, v01, buff0, c0, c0, c1, c1); + 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 - mat3 xy_to_uv = get_xy_to_uv(flat_controls[0], flat_controls[1]); - float scale_factor = length(flat_controls[1] - flat_controls[0]); - uv_anti_alias_width = anti_alias_width * frame_shape.y / pixel_shape.y / scale_factor; - uv_b2 = (xy_to_uv * vec3(flat_controls[2], 1.0)).xy; + // Find uv conversion + mat3 xy_to_uv = get_xy_to_uv(flat_controls, bezier_degree); + float scale_factor = length(xy_to_uv[0].xy); + uv_anti_alias_width = scale_factor * anti_alias_width * (frame_shape.y / pixel_shape.y); // 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(