Got stroke shaders working in 3d

manimlib/mobject/types/vectorized_mobject.py

@@ -1,7 +1,9 @@
 import itertools as it
+import operator as op
 import moderngl
 
 from colour import Color
+from functools import reduce
 
 from manimlib.constants import *
 from manimlib.mobject.mobject import Mobject
@@ -23,7 +25,7 @@ from manimlib.utils.space_ops import angle_between_vectors
 from manimlib.utils.space_ops import cross2d
 from manimlib.utils.space_ops import earclip_triangulation
 from manimlib.utils.space_ops import get_norm
-from manimlib.utils.space_ops import normalize
+from manimlib.utils.space_ops import get_unit_normal
 from manimlib.utils.space_ops import z_to_vector
 from manimlib.utils.shaders import get_shader_info
 
@@ -75,6 +77,7 @@ class VMobject(Mobject):
             ("point", np.float32, (3,)),
             ("prev_point", np.float32, (3,)),
             ("next_point", np.float32, (3,)),
+            ('unit_normal', np.float32, (3,)),
             ("stroke_width", np.float32, (1,)),
             ("color", np.float32, (4,)),
             ("joint_type", np.float32, (1,)),
@@ -668,6 +671,14 @@ class VMobject(Mobject):
             self.get_end_anchors(),
         ))))
 
+    def get_points_without_null_curves(self, atol=1e-9):
+        nppc = self.n_points_per_curve
+        distinct_curves = reduce(op.or_, [
+            (abs(self.points[i::nppc] - self.points[0::nppc]) > atol).any(1)
+            for i in range(1, nppc)
+        ])
+        return self.points[distinct_curves.repeat(nppc)]
+
     def get_arc_length(self, n_sample_points=None):
         if n_sample_points is None:
             n_sample_points = 4 * self.get_num_curves() + 1
@@ -679,6 +690,38 @@ class VMobject(Mobject):
         norms = np.array([get_norm(d) for d in diffs])
         return norms.sum()
 
+    def get_area_vector(self):
+        # Returns a vector whose length is the area bound by
+        # the polygon formed by the anchor points, pointing
+        # in a direction perpendicular to the polygon according
+        # to the right hand rule.
+        if self.has_no_points():
+            return np.zeros(3)
+
+        nppc = self.n_points_per_curve
+        p0 = self.points[0::nppc]
+        p1 = self.points[nppc - 1::nppc]
+
+        # Each term goes through all edges [(x1, y1, z1), (x2, y2, z2)]
+        return 0.5 * np.array([
+            sum((p0[:, 1] + p1[:, 1]) * (p1[:, 2] - p0[:, 2])),  # Add up (y1 + y2)*(z2 - z1)
+            sum((p0[:, 2] + p1[:, 2]) * (p1[:, 0] - p0[:, 0])),  # Add up (z1 + z2)*(x2 - x1)
+            sum((p0[:, 0] + p1[:, 0]) * (p1[:, 1] - p0[:, 1])),  # Add up (x1 + x2)*(y2 - y1)
+        ])
+
+    def get_unit_normal_vector(self):
+        if len(self.points) < 3:
+            return OUT
+        area_vect = self.get_area_vector()
+        area = get_norm(area_vect)
+        if area > 0:
+            return area_vect / area
+        else:
+            return get_unit_normal(
+                self.points[1] - self.points[0],
+                self.points[2] - self.points[1],
+            )
+
     # Alignment
     def align_points(self, vmobject):
         self.align_rgbas(vmobject)
@@ -910,12 +953,16 @@ class VMobject(Mobject):
         if len(stroke_width) > 1:
             stroke_width = self.stretched_style_array_matching_points(stroke_width)
 
-        data = self.get_blank_shader_data_array(len(self.points), "stroke_data")
+        points = self.get_points_without_null_curves()
-        data["point"] = self.points
+        nppc = self.n_points_per_curve
-        data["prev_point"][:3] = self.points[-3:]
+
-        data["prev_point"][3:] = self.points[:-3]
+        data = self.get_blank_shader_data_array(len(points), "stroke_data")
-        data["next_point"][:-3] = self.points[3:]
+        data["point"] = points
-        data["next_point"][-3:] = self.points[:3]
+        data["prev_point"][:nppc] = points[-nppc:]
+        data["prev_point"][nppc:] = points[:-nppc]
+        data["next_point"][:-nppc] = points[nppc:]
+        data["next_point"][-nppc:] = points[:nppc]
+        data["unit_normal"] = self.get_unit_normal_vector()
         data["stroke_width"][:, 0] = stroke_width
         data["color"] = rgbas
         data["joint_type"] = joint_type_to_code[self.joint_type]
@@ -939,27 +986,6 @@ class VMobject(Mobject):
             if sm.triangulation_locked:
                 sm.lock_triangulation(family=False)
 
-    def get_area_vector(self):
-        # Returns a vector whose length is the area bound by
-        # the polygon formed by the anchor points, pointing
-        # in a direction perpendicular to the polygon according
-        # to the right hand rule.
-        nppc = self.n_points_per_curve
-        p0 = self.points[0::nppc]
-        p1 = self.points[nppc - 1::nppc]
-
-        # Each term goes through all edges [(x1, y1, z1), (x2, y2, z2)]
-        return 0.5 * np.array([
-            sum((p0[:, 1] + p1[:, 1]) * (p1[:, 2] - p0[:, 2])),  # Add up (y1 + y2)*(z2 - z1)
-            sum((p0[:, 2] + p1[:, 2]) * (p1[:, 0] - p0[:, 0])),  # Add up (z1 + z2)*(x2 - x1)
-            sum((p0[:, 0] + p1[:, 0]) * (p1[:, 1] - p0[:, 1])),  # Add up (x1 + x2)*(y2 - y1)
-        ])
-
-    def get_unit_normal_vector(self):
-        if self.has_no_points():
-            return ORIGIN
-        return normalize(self.get_area_vector())
-
     def get_triangulation(self, normal_vector=None):
         # Figure out how to triangulate the interior to know
         # how to send the points as to the vertex shader.

manimlib/shaders/get_unit_normal.glsl

@@ -13,7 +13,7 @@ vec3 get_unit_normal(in vec3[3] points){
         if(new_cp_norm < tol){
             // We only come here if all three points line up
             // on the z-axis.
-            return vec3(1.0, 0.0, 0.0);
+            return vec3(0.0, 1.0, 0.0);
             // return k_hat;
         }
         return new_cp / new_cp_norm;

manimlib/shaders/quadratic_bezier_fill_frag.glsl

@@ -8,7 +8,7 @@ in float fill_all;  // Either 0 or 1e
 in float uv_anti_alias_width;
 
 in vec3 xyz_coords;
-in vec3 local_unit_normal;
+in vec3 global_unit_normal;
 in float orientation;
 in vec2 uv_coords;
 in vec2 uv_b2;
@@ -67,7 +67,7 @@ float sdf(){
 
 void main() {
     if (color.a == 0) discard;
-    frag_color = add_light(color, xyz_coords, local_unit_normal, light_source_position, gloss);
+    frag_color = add_light(color, xyz_coords, global_unit_normal, light_source_position, gloss);
     if (fill_all == 1.0) return;
     frag_color.a *= smoothstep(1, 0, sdf() / uv_anti_alias_width);
 }

manimlib/shaders/quadratic_bezier_fill_geom.glsl

@@ -20,7 +20,7 @@ out float fill_all;
 out float uv_anti_alias_width;
 
 out vec3 xyz_coords;
-out vec3 local_unit_normal;
+out vec3 global_unit_normal;
 out float orientation;
 // uv space is where b0 = (0, 0), b1 = (1, 0), and transform is orthogonal
 out vec2 uv_coords;
@@ -38,6 +38,7 @@ out float bezier_degree;
 void emit_vertex_wrapper(vec3 point, int index){
     color = v_color[index];
     gloss = v_gloss[index];
+    global_unit_normal = v_global_unit_normal[index];
     xyz_coords = point;
     gl_Position = get_gl_Position(xyz_coords);
     EmitVertex();
@@ -59,29 +60,29 @@ void emit_pentagon(vec3[3] points, vec3 normal){
     // Tangent vectors
     vec3 t01 = normalize(p1 - p0);
     vec3 t12 = normalize(p2 - p1);
-    // Vectors normal to the curve in the plane of the curve pointing outside the curve
+    // Vectors perpendicular to the curve in the plane of the curve pointing outside the curve
-    vec3 n01 = cross(t01, normal);
+    vec3 p0_perp = cross(t01, normal);
-    vec3 n12 = cross(t12, normal);
+    vec3 p2_perp = cross(t12, normal);
 
     bool fill_in = orientation > 0;
     float aaw = anti_alias_width;
     vec3 corners[5];
     if(fill_in){
-        // Note, straight lines will also fall into this case, and since n01 and n12
+        // Note, straight lines will also fall into this case, and since p0_perp and p2_perp
         // will point to the right of the curve, it's just what we want
         corners = vec3[5](
-            p0 + aaw * n01,
+            p0 + aaw * p0_perp,
             p0,
-            p1 + 0.5 * aaw * (n01 + n12),
+            p1 + 0.5 * aaw * (p0_perp + p2_perp),
             p2,
-            p2 + aaw * n12
+            p2 + aaw * p2_perp
         );
     }else{
         corners = vec3[5](
             p0,
-            p0 - aaw * n01,
+            p0 - aaw * p0_perp,
             p1,
-            p2 - aaw * n12,
+            p2 - aaw * p2_perp,
             p2
         );
     }
@@ -102,7 +103,7 @@ void emit_pentagon(vec3[3] points, vec3 normal){
 
 void main(){
     fill_all = v_fill_all[0];
-    local_unit_normal = get_unit_normal(vec3[3](bp[0], bp[1], bp[2]));
+    vec3 local_unit_normal = get_unit_normal(vec3[3](bp[0], bp[1], bp[2]));
     orientation = sign(dot(v_global_unit_normal[0], local_unit_normal));
 
     if(fill_all == 1){

manimlib/shaders/quadratic_bezier_fill_vert.glsl

@@ -21,7 +21,7 @@ out float v_gloss;
 
 void main(){
     bp = position_point_into_frame(point);
-    v_global_unit_normal = position_point_into_frame(unit_normal);
+    v_global_unit_normal = normalize(position_point_into_frame(unit_normal));
     v_color = color;
     v_fill_all = fill_all;
     v_gloss = gloss;

manimlib/shaders/quadratic_bezier_geometry_functions.glsl

@@ -47,7 +47,8 @@ float get_reduced_control_points(in vec3 points[3], out vec3 new_points[3]){
     vec3 v01 = (p1 - p0);
     vec3 v12 = (p2 - p1);
 
-    bool aligned = acos(dot(normalize(v01), normalize(v12))) < angle_threshold;
+    float dot_prod = clamp(dot(normalize(v01), normalize(v12)), -1, 1);
+    bool aligned = acos(dot_prod) < angle_threshold;
     bool distinct_01 = length(v01) > length_threshold;  // v01 is considered nonzero
     bool distinct_12 = length(v12) > length_threshold;  // v12 is considered nonzero
     int n_uniques = int(distinct_01) + int(distinct_12);

manimlib/shaders/quadratic_bezier_stroke_frag.glsl

@@ -4,7 +4,7 @@ uniform mat4 to_screen_space;
 uniform vec3 light_source_position;
 
 in vec3 xyz_coords;
-in vec3 unit_normal;
+in vec3 global_unit_normal;
 in vec2 uv_coords;
 in vec2 uv_b2;
 
@@ -91,10 +91,12 @@ void main() {
     if (uv_stroke_width == 0) discard;
 
     // Add lighting if needed
-    frag_color = add_light(color, xyz_coords, unit_normal, light_source_position, gloss);
+    frag_color = add_light(color, xyz_coords, global_unit_normal, light_source_position, gloss);
 
     float dist_to_curve = min_dist_to_curve(uv_coords, uv_b2, bezier_degree);
     // An sdf for the region around the curve we wish to color.
     float signed_dist = abs(dist_to_curve) - 0.5 * uv_stroke_width;
     frag_color.a *= smoothstep(0.5, -0.5, signed_dist / uv_anti_alias_width);
+
+    // frag_color.a += 0.3;
 }

manimlib/shaders/quadratic_bezier_stroke_geom.glsl

@@ -11,6 +11,7 @@ uniform float anti_alias_width;
 in vec3 bp[3];
 in vec3 prev_bp[3];
 in vec3 next_bp[3];
+in vec3 v_global_unit_normal[3];
 
 in vec4 v_color[3];
 in float v_stroke_width[3];
@@ -32,7 +33,7 @@ out float angle_to_next;
 out float bezier_degree;
 
 out vec3 xyz_coords;
-out vec3 unit_normal;
+out vec3 global_unit_normal;
 out vec2 uv_coords;
 out vec2 uv_b2;
 
@@ -51,9 +52,19 @@ const float MITER_JOINT = 3;
 #INSERT get_unit_normal.glsl
 
 
+float get_aaw_scalar(vec3 normal){
+    return min(abs(normal.z), 5);
+}
+
+
 float angle_between_vectors(vec3 v1, vec3 v2, vec3 normal){
-    vec3 nv1 = normalize(v1);
+    float v1_norm = length(v1);
-    vec3 nv2 = normalize(v2);
+    float v2_norm = length(v2);
+    if(v1_norm == 0 || v2_norm == 0) return 0;
+    vec3 nv1 = v1 / v1_norm;
+    vec3 nv2 = v2 / v2_norm;
+    // float signed_area = clamp(dot(cross(nv1, nv2), normal), -1, 1);
+    // return asin(signed_area);
     float unsigned_angle = acos(clamp(dot(nv1, nv2), -1, 1));
     float sn = sign(dot(cross(nv1, nv2), normal));
     return sn * unsigned_angle;
@@ -129,23 +140,20 @@ int get_corners(vec3 controls[3], vec3 normal, int degree, out vec3 corners[5]){
     vec3 p1 = controls[1];
     vec3 p2 = controls[2];
 
-    // Unit vectors for directions between
+    // Unit vectors for directions between control points
-    // Various control points
-    vec3 v02 = normalize(p2 - p0);
     vec3 v10 = normalize(p0 - p1);
     vec3 v12 = normalize(p2 - p1);
-    vec3 v20 = -v02;
     vec3 v01 = -v10;
     vec3 v21 = -v12;
 
-    // Find bounding points around ends
+    // 
-    vec3 p0_perp = cross(normal, v01);
+    vec3 p0_perp = cross(normal, v01);  // Pointing to the left of the curve from p0
-    vec3 p2_perp = cross(normal, v21);
+    vec3 p2_perp = cross(normal, v12);  // Pointing to the left of the curve from p2
 
     // aaw is the added width given around the polygon for antialiasing.
     // In case the normal is faced away from (0, 0, 1), the vector to the
     // camera, this is scaled up.
-    float aaw = anti_alias_width / normal.z;
+    float aaw = anti_alias_width / get_aaw_scalar(normal);
     float buff0 = 0.5 * v_stroke_width[0] + aaw;
     float buff2 = 0.5 * v_stroke_width[2] + aaw;
     float aaw0 = (1 - has_prev) * aaw;
@@ -153,16 +161,12 @@ int get_corners(vec3 controls[3], vec3 normal, int degree, out vec3 corners[5]){
 
     vec3 c0 = p0 - buff0 * p0_perp + aaw0 * v10;
     vec3 c1 = p0 + buff0 * p0_perp + aaw0 * v10;
-    vec3 c2 = p2 - buff2 * p2_perp + aaw2 * v12;
+    vec3 c2 = p2 + buff2 * p2_perp + aaw2 * v12;
-    vec3 c3 = p2 + buff2 * p2_perp + aaw2 * v12;
+    vec3 c3 = p2 - buff2 * p2_perp + aaw2 * v12;
 
     // Account for previous and next control points
-    if(has_prev == 1){
+    if(has_prev > 0) create_joint(angle_from_prev, v01, buff0, bevel_start, c0, c0, c1, c1);
-        create_joint(angle_from_prev, v01, buff0, bevel_start, c0, c0, c1, c1);
+    if(has_next > 0) create_joint(angle_to_next, v21, buff2, bevel_end, c3, c3, c2, c2);
-    }
-    if(has_next == 1){
-        create_joint(-angle_to_next, v21, buff2, bevel_end, c2, c2, c3, c3);
-    }
 
     // Linear case is the simplest
     if(degree == 1){
@@ -171,46 +175,12 @@ int get_corners(vec3 controls[3], vec3 normal, int degree, out vec3 corners[5]){
         corners = vec3[5](c0, c1, c3, c2, vec3(0.0));
         return 4;
     }
-
+    // Otherwise, form a pentagon around the curve
-    // Some admitedly complicated logic to (hopefully efficiently)
+    float orientation = sign(dot(cross(v01, v12), normal));  // Positive for ccw curves
-    // make sure corners forms a convex hull around the curve.
+    if(orientation > 0) corners = vec3[5](c0, c1, p1, c2, c3);
-    if(dot(cross(v10, v12), normal) > 0){
+    else                corners = vec3[5](c1, c0, p1, c3, c2);
-        bool change_c0 = (
+    // Replace corner[2] with convex hull point accounting for stroke width
-            // has_prev == 0 &&
+    find_intersection(corners[0], v01, corners[4], v21, normal, corners[2]);
-            dot(v21, v20) > 0 &&
-            should_motify_corner(c0, v01, c2, c3, v21, normal, buff0)
-        );
-        if(change_c0) c0 = p0 + p2_perp * buff0;
-
-        bool change_c3 = (
-            // has_next == 0 &&
-            dot(v01, v02) > 0 &&
-            should_motify_corner(c3, v21, c1, c0, v01, normal, buff2)
-        );
-        if(change_c3) c3 = p2 - p0_perp * buff2;
-
-        vec3 i12;
-        find_intersection(c1, v01, c2, v21, normal, i12);
-        corners = vec3[5](c1, c0, i12, c3, c2);
-    }else{
-        bool change_c1 = (
-            // has_prev == 0 &&
-            dot(v21, v20) > 0 &&
-            should_motify_corner(c1, v01, c3, c2, v21, normal, buff0)
-        );
-        if(change_c1) c1 = p0 - p2_perp * buff0;
-
-        bool change_c2 = (
-            // has_next == 0 &&
-            dot(v01, v02) > 0 &&
-            should_motify_corner(c2, v21, c0, c1, v01, normal, buff2)
-        );
-        if(change_c2) c2 = p2 + p0_perp * buff2;
-
-        vec3 i03;
-        find_intersection(c0, v01, c3, v21, normal, i03);
-        corners = vec3[5](c0, c1, i03, c2, c3);
-    }
     return 5;
 }
 
@@ -219,23 +189,14 @@ void set_adjascent_info(vec3 c0, vec3 tangent,
                         int degree,
                         vec3 normal,
                         vec3 adj[3],
-                        out float has,
                         out float bevel,
                         out float angle
                         ){
     float joint_type = v_joint_type[0];
-
-    has = 0;
-    bevel = 0;
-    angle = 0;
-
     vec3 new_adj[3];
     float adj_degree = get_reduced_control_points(adj, new_adj);
-    has = float(adj_degree > 0);
+    // Check if adj_degree is zero?
-    if(has == 1){
+    angle = angle_between_vectors(c0 - new_adj[1], tangent, normal);
-        vec3 adj = new_adj[1];
-        angle = angle_between_vectors(c0 - adj, tangent, normal);
-    }
     // Decide on joint type
     bool one_linear = (degree == 1 || adj_degree == 1.0);
     bool should_bevel = (
@@ -247,29 +208,35 @@ void set_adjascent_info(vec3 c0, vec3 tangent,
 
 
 void set_previous_and_next(vec3 controls[3], int degree, vec3 normal){
-    float a_tol = 1e-10;
+    float a_tol = 1e-8;
 
-    if(distance(prev_bp[2], bp[0]) < a_tol){
+    // Made as floats not bools so they can be passed to the frag shader
+    has_prev = float(distance(prev_bp[2], bp[0]) < a_tol);
+    has_next = float(distance(next_bp[0], bp[2]) < a_tol);
+
+    if(has_prev > 0){
         vec3 tangent = controls[1] - controls[0];
         set_adjascent_info(
             controls[0], tangent, degree, normal,
             vec3[3](prev_bp[0], prev_bp[1], prev_bp[2]),
-            has_prev, bevel_start, angle_from_prev
+            bevel_start, angle_from_prev
         );
     }
-    if(distance(next_bp[0], bp[2]) < a_tol){
+    if(has_next > 0){
         vec3 tangent = controls[1] - controls[2];
         set_adjascent_info(
             controls[2], tangent, degree, normal,
             vec3[3](next_bp[0], next_bp[1], next_bp[2]),
-            has_next, bevel_end, angle_to_next
+            bevel_end, angle_to_next
         );
+        angle_to_next *= -1;
     }
 }
 
 
 void main() {
-    unit_normal = get_unit_normal(vec3[3](bp[0], bp[1], bp[2]));
+    vec3 unit_normal = v_global_unit_normal[0];
+    // anti_alias_width /= cos(0.5 * acos(abs(unit_normal.z)));
 
     vec3 controls[3];
     bezier_degree = get_reduced_control_points(vec3[3](bp[0], bp[1], bp[2]), controls);
@@ -283,16 +250,13 @@ void main() {
     // Find uv conversion matrix
     mat4 xyz_to_uv = get_xyz_to_uv(controls[0], controls[1], unit_normal);
     float scale_factor = length(controls[1] - controls[0]);
-    uv_anti_alias_width = anti_alias_width / scale_factor / unit_normal.z;
+    uv_anti_alias_width = anti_alias_width / scale_factor / get_aaw_scalar(unit_normal);
     uv_b2 = (xyz_to_uv * vec4(controls[2], 1.0)).xy;
 
     // Corners of a bounding region around curve
     vec3 corners[5];
     int n_corners = get_corners(controls, unit_normal, degree, corners);
 
-    // Get style info aligned to the corners
-    float stroke_widths[5];
-    vec4 stroke_colors[5];
     int index_map[5] = int[5](0, 0, 1, 2, 2);
     if(n_corners == 4) index_map[2] = 2;
 
@@ -303,6 +267,7 @@ void main() {
         uv_stroke_width = v_stroke_width[index_map[i]] / scale_factor;
         color = v_color[index_map[i]];
         gloss = v_gloss[index_map[i]];
+        global_unit_normal = v_global_unit_normal[index_map[i]];
         gl_Position = get_gl_Position(xyz_coords);
         EmitVertex();
     }

manimlib/shaders/quadratic_bezier_stroke_vert.glsl

@@ -6,6 +6,7 @@ uniform float focal_distance;
 in vec3 point;
 in vec3 prev_point;
 in vec3 next_point;
+in vec3 unit_normal;
 
 in float stroke_width;
 in vec4 color;
@@ -16,6 +17,7 @@ in float gloss;
 out vec3 bp;
 out vec3 prev_bp;
 out vec3 next_bp;
+out vec3 v_global_unit_normal;
 
 out float v_stroke_width;
 out vec4 v_color;
@@ -33,9 +35,9 @@ void main(){
     bp = position_point_into_frame(point);
     prev_bp = position_point_into_frame(prev_point);
     next_bp = position_point_into_frame(next_point);
+    v_global_unit_normal = normalize(position_point_into_frame(unit_normal));
 
     v_stroke_width = STROKE_WIDTH_CONVERSION * stroke_width;
-    // v_stroke_width /= (1 - bp.z / focal_distance);  // Change stroke width by perspective
     v_color = color;
     v_joint_type = joint_type;
     v_gloss = gloss;

manimlib/utils/space_ops.py

@@ -3,9 +3,10 @@ import math
 import itertools as it
 from mapbox_earcut import triangulate_float32 as earcut
 
+from manimlib.constants import RIGHT
+from manimlib.constants import UP
 from manimlib.constants import OUT
 from manimlib.constants import PI
-from manimlib.constants import RIGHT
 from manimlib.constants import TAU
 from manimlib.utils.iterables import adjacent_pairs
 
@@ -192,7 +193,7 @@ def get_unit_normal(v1, v2, tol=1e-6):
         new_cp = cross(cross(v1, OUT), v1)
         new_cp_norm = get_norm(new_cp)
         if new_cp_norm < tol:
-            return RIGHT
+            return UP
         return new_cp / new_cp_norm
     return cp / cp_norm