Fix up fill shaders to work when being viewed from different orientations, along with several other little bugs

manimlib/mobject/types/vectorized_mobject.py

@@ -19,10 +19,12 @@ from manimlib.utils.iterables import make_even
 from manimlib.utils.iterables import stretch_array_to_length
 from manimlib.utils.iterables import stretch_array_to_length_with_interpolation
 from manimlib.utils.iterables import listify
-from manimlib.utils.space_ops import cross2d
-from manimlib.utils.space_ops import get_norm
 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 z_to_vector
 from manimlib.utils.shaders import get_shader_info
 
 
@@ -64,6 +66,7 @@ class VMobject(Mobject):
         "triangulation_locked": False,
         "fill_dtype": [
             ('point', np.float32, (3,)),
+            ('unit_normal', np.float32, (3,)),
             ('color', np.float32, (4,)),
             ('fill_all', np.float32, (1,)),
             ('gloss', np.float32, (1,)),
@@ -924,7 +927,6 @@ class VMobject(Mobject):
         for mob in mobs:
             mob.triangulation_locked = False
             mob.saved_triangulation = mob.get_triangulation()
-            mob.saved_orientation = mob.get_orientation()
             mob.triangulation_locked = True
         return self
 
@@ -937,26 +939,33 @@ class VMobject(Mobject):
             if sm.triangulation_locked:
                 sm.lock_triangulation(family=False)
 
-    def get_signed_polygonal_area(self):
+    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]
-        # Add up (x1 + x2)*(y2 - y1) for all edges (x1, y1), (x2, y2)
-        return sum((p0[:, 0] + p1[:, 0]) * (p1[:, 1] - p0[:, 1]))
 
-    def get_orientation(self):
-        if self.triangulation_locked:
-            return self.saved_orientation
+        # 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 0
-        return np.sign(self.get_signed_polygonal_area())
+            return ORIGIN
+        return normalize(self.get_area_vector())
 
-    def get_triangulation(self, orientation=None):
+    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.
         # First triangles come directly from the points
-        if orientation is None:
-            orientation = self.get_orientation()
+        if normal_vector is None:
+            normal_vector = self.get_unit_normal_vector()
 
         if self.triangulation_locked:
             return self.saved_triangulation
@@ -964,8 +973,9 @@ class VMobject(Mobject):
         if len(self.points) <= 1:
             return []
 
-        # Otherwise, compute from scratch
-        points = self.points
+        # Rotate points such that unit normal vector is OUT
+        # TODO, 99% of the time this does nothing.  Do a check for that?
+        points = np.dot(self.points, z_to_vector(normal_vector))
         indices = np.arange(len(points), dtype=int)
 
         b0s = points[0::3]
@@ -974,9 +984,8 @@ class VMobject(Mobject):
         v01s = b1s - b0s
         v12s = b2s - b1s
 
-        # TODO, account for 3d
         crosses = cross2d(v01s, v12s)
-        convexities = orientation * np.sign(crosses)
+        convexities = np.sign(crosses)
 
         atol = self.tolerance_for_point_equality
         end_of_loop = np.zeros(len(b0s), dtype=bool)
@@ -1003,23 +1012,21 @@ class VMobject(Mobject):
 
     def get_fill_shader_data(self):
         points = self.points
-        orientation = self.get_orientation()
-        tri_indices = self.get_triangulation(orientation)
+        unit_normal = self.get_unit_normal_vector()
+        tri_indices = self.get_triangulation(unit_normal)
 
         # TODO, best way to enable multiple colors?
         rgbas = self.get_fill_rgbas()[:1]
 
         data = self.get_blank_shader_data_array(len(tri_indices), "fill_data")
         data["point"] = points[tri_indices]
+        data["unit_normal"] = unit_normal
         data["color"] = rgbas
         # Assume the triangulation is such that the first n_points points
         # are on the boundary, and the rest are in the interior
         data["fill_all"][:len(points)] = 0
         data["fill_all"][len(points):] = 1
         data["gloss"] = self.gloss
-        # Always send points in a positively oriented way
-        if orientation < 0:
-            data["point"][:len(points)] = points[::-1]
         return data
 
 

manimlib/shaders/add_light.glsl

@@ -1,13 +1,19 @@
 vec4 add_light(vec4 raw_color, vec3 point, vec3 unit_normal, vec3 light_coords, float gloss){
     if(gloss == 0.0) return raw_color;
 
+    // TODO, do we actually want this?  For VMobjects its nice to just choose whichever unit normal
+    // is pointing towards the camera.
+    if(unit_normal.z < 0){
+        unit_normal *= -1;
+    }
+
     float camera_distance = 6;
     // Assume everything has already been rotated such that camera is in the z-direction
     vec3 to_camera = vec3(0, 0, camera_distance) - point;
     vec3 to_light = light_coords - point;
     vec3 light_reflection = -to_light + 2 * unit_normal * dot(to_light, unit_normal);
     float dot_prod = dot(normalize(light_reflection), normalize(to_camera));
-    float shine = gloss * exp(-2 * pow(1 - dot_prod, 2));
+    float shine = gloss * exp(-3 * pow(1 - dot_prod, 2));
     return vec4(
         mix(raw_color.rgb, vec3(1.0), shine),
         raw_color.a

manimlib/shaders/get_gl_Position.glsl

@@ -10,5 +10,5 @@ vec4 get_gl_Position(vec3 point){
     // the z-coordiante of gl_Position only matter for z-indexing.  The reason
     // for thie line is to avoid agressive clipping of distant points.
     if(point.z < 0) point.z *= 0.1;
-    return vec4(point, 1);
+    return vec4(point.xy, -point.z, 1);
 }

manimlib/shaders/get_unit_normal.glsl

@@ -1,15 +1,22 @@
-
-vec3 get_unit_normal(in vec3 point0, in vec3 point1, in vec3 point2){
-    vec3 cp = cross(point1 - point0, point2 - point1);
-    if(length(cp) == 0){
-        return vec3(0.0, 0.0, 1.0);
-    }else{
-        if(cp.z < 0){
-            // After re-orienting, camera will always sit in the positive
-            // z-direction.  We always want normal vectors pointing towards
-            // the camera.
-            cp *= -1;
+vec3 get_unit_normal(in vec3[3] points){
+    float tol = 1e-6;
+    vec3 v1 = normalize(points[1] - points[0]);
+    vec3 v2 = normalize(points[2] - points[0]);
+    vec3 cp = cross(v1, v2);
+    float cp_norm = length(cp);
+    if(cp_norm < tol){
+        // Three points form a line, so find a normal vector
+        // to that line in the plane shared with the z-axis
+        vec3 k_hat = vec3(0.0, 0.0, 1.0);
+        vec3 new_cp = cross(cross(v2, k_hat), v2);
+        float new_cp_norm = length(new_cp);
+        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 k_hat;
         }
-        return normalize(cp);
+        return new_cp / new_cp_norm;
     }
+    return cp / cp_norm;
 }

manimlib/shaders/quadratic_bezier_fill_frag.glsl

@@ -8,7 +8,8 @@ in float fill_all;  // Either 0 or 1e
 in float uv_anti_alias_width;
 
 in vec3 xyz_coords;
-in vec3 unit_normal;
+in vec3 local_unit_normal;
+in float orientation;
 in vec2 uv_coords;
 in vec2 uv_b2;
 in float bezier_degree;
@@ -29,24 +30,33 @@ float modify_distance_for_endpoints(vec2 p, float dist, float t){
 
 
 float sdf(){
-    // For really flat curves, just take the distance to the curve
-    if(bezier_degree < 2 || abs(uv_b2.y / uv_b2.x) < uv_anti_alias_width){
+    if(bezier_degree < 2){
+        return abs(uv_coords[1]);
+    }
+    float u2 = uv_b2.x;
+    float v2 = uv_b2.y;
+    // For really flat curves, just take the distance to x-axis
+    if(abs(v2 / u2) < 0.1 * uv_anti_alias_width){
+        return abs(uv_coords[1]);
+    }
+    // For flat-ish curves, take the curve
+    else if(abs(v2 / u2) < 0.5 * uv_anti_alias_width){
         return min_dist_to_curve(uv_coords, uv_b2, bezier_degree);
     }
+    // I know, I don't love this amount of arbitrary-seeming branching either,
+    // but a number of strange dimples and bugs pop up otherwise.
+
     // This converts uv_coords to yet another space where the bezier points sit on
     // (0, 0), (1/2, 0) and (1, 1), so that the curve can be expressed implicityly
     // as y = x^2.
-    float u2 = uv_b2.x;
-    float v2 = uv_b2.y;
     mat2 to_simple_space = mat2(
         v2, 0,
         2 - u2, 4 * v2
     );
     vec2 p = to_simple_space * uv_coords;
-
     // Sign takes care of whether we should be filling the inside or outside of curve.
-    float Fp = sign(v2) * (p.x * p.x - p.y);
-
+    float sn = orientation * sign(v2);
+    float Fp = sn * (p.x * p.x - p.y);
     vec2 grad = vec2(
         -2 * p.x * v2,  // del C / del u
         4 * v2 - 4 * p.x * (2 - u2)  // del C / del v
@@ -57,7 +67,7 @@ float sdf(){
 
 void main() {
     if (color.a == 0) discard;
-    frag_color = add_light(color, xyz_coords, unit_normal, light_source_position, gloss);
+    frag_color = add_light(color, xyz_coords, local_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

@@ -9,21 +9,23 @@ uniform float aspect_ratio;
 uniform float focal_distance;
 
 in vec3 bp[3];
+in vec3 v_global_unit_normal[3];
 in vec4 v_color[3];
 in float v_fill_all[3];
 in float v_gloss[3];
 
 out vec4 color;
+out float gloss;
 out float fill_all;
 out float uv_anti_alias_width;
 
 out vec3 xyz_coords;
-out vec3 unit_normal;
+out vec3 local_unit_normal;
+out float orientation;
 // uv space is where b0 = (0, 0), b1 = (1, 0), and transform is orthogonal
 out vec2 uv_coords;
 out vec2 uv_b2;
 out float bezier_degree;
-out float gloss;
 
 // To my knowledge, there is no notion of #include for shaders,
 // so to share functionality between this and others, the caller
@@ -32,13 +34,19 @@ out float gloss;
 #INSERT get_gl_Position.glsl
 #INSERT get_unit_normal.glsl
 
+
+void emit_vertex_wrapper(vec3 point, int index){
+    color = v_color[index];
+    gloss = v_gloss[index];
+    xyz_coords = point;
+    gl_Position = get_gl_Position(xyz_coords);
+    EmitVertex();
+}
+
+
 void emit_simple_triangle(){
     for(int i = 0; i < 3; i++){
-        color = v_color[i];
-        gloss = v_gloss[i];
-        xyz_coords = bp[i];
-        gl_Position = get_gl_Position(bp[i]);
-        EmitVertex();
+        emit_vertex_wrapper(bp[i], i);
     }
     EndPrimitive();
 }
@@ -51,42 +59,42 @@ 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
+    // Vectors normal to the curve in the plane of the curve pointing outside the curve
     vec3 n01 = cross(t01, normal);
     vec3 n12 = cross(t12, normal);
 
-    // Assume you always fill in to the left of the curve
-    float orient = sign(dot(cross(t01, t12), normal));
-    bool fill_in = (orient > 0);
-
-    float aaw = anti_alias_width / normal.z;
-    vec3 nudge1 = fill_in ? 0.5 * aaw * (n01 + n12) : vec3(0);
-    vec3 corners[5] = vec3[5](
+    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
+        // will point to the right of the curve, it's just what we want
+        corners = vec3[5](
             p0 + aaw * n01,
             p0,
-        p1 + nudge1,
+            p1 + 0.5 * aaw * (n01 + n12),
             p2,
             p2 + aaw * n12
         );
-
-    int coords_index_map[5] = int[5](0, 1, 2, 3, 4);
-    if(!fill_in) coords_index_map = int[5](1, 0, 2, 4, 3);
+    }else{
+        corners = vec3[5](
+            p0,
+            p0 - aaw * n01,
+            p1,
+            p2 - aaw * n12,
+            p2
+        );
+    }
 
     mat4 xyz_to_uv = get_xyz_to_uv(p0, p1, normal);
     uv_b2 = (xyz_to_uv * vec4(p2, 1)).xy;
     uv_anti_alias_width = anti_alias_width / length(p1 - p0);
 
     for(int i = 0; i < 5; i++){
-        vec3 corner = corners[coords_index_map[i]];
-        xyz_coords = corner;
+        vec3 corner = corners[i];
         uv_coords = (xyz_to_uv * vec4(corner, 1)).xy;
-        // I haven't a clue why an index map doesn't work just
-        // as well here, but for some reason it doesn't.
-        int j = int(sign(i - 1) + 1);  // Maps 0, 1, 2, 3, 4 onto 0, 0, 1, 2, 2
-        color = v_color[j];
-        gloss = v_gloss[j];
-        gl_Position = get_gl_Position(corner);
-        EmitVertex();
+        int j = int(sign(i - 1) + 1);  // Maps i = [0, 1, 2, 3, 4] onto j = [0, 0, 1, 2, 2]
+        emit_vertex_wrapper(corner, j);
     }
     EndPrimitive();
 }
@@ -94,7 +102,8 @@ void emit_pentagon(vec3[3] points, vec3 normal){
 
 void main(){
     fill_all = v_fill_all[0];
-    unit_normal = get_unit_normal(bp[0], bp[1], bp[2]);
+    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){
         emit_simple_triangle();
@@ -103,7 +112,9 @@ void main(){
 
     vec3 new_bp[3];
     bezier_degree = get_reduced_control_points(vec3[3](bp[0], bp[1], bp[2]), new_bp);
-    if(bezier_degree == 0) return;  // Don't emit any vertices
-    emit_pentagon(new_bp, unit_normal);
+    if(bezier_degree >= 1){
+        emit_pentagon(new_bp, local_unit_normal);
+    }
+    // Don't emit any vertices for bezier_degree 0
 }
 

manimlib/shaders/quadratic_bezier_fill_vert.glsl

@@ -3,11 +3,13 @@
 uniform mat4 to_screen_space;
 
 in vec3 point;
+in vec3 unit_normal;
 in vec4 color;
 in float fill_all;  // Either 0 or 1
 in float gloss;
 
 out vec3 bp;  // Bezier control point
+out vec3 v_global_unit_normal;
 out vec4 v_color;
 out float v_fill_all;
 out float v_gloss;
@@ -19,6 +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_color = color;
     v_fill_all = fill_all;
     v_gloss = gloss;

manimlib/shaders/quadratic_bezier_stroke_geom.glsl

@@ -269,8 +269,7 @@ void set_previous_and_next(vec3 controls[3], int degree, vec3 normal){
 
 
 void main() {
-    unit_normal = get_unit_normal(bp[0], bp[1], bp[2]);
-    // unit_normal = vec3(0, 0, 1);
+    unit_normal = get_unit_normal(vec3[3](bp[0], bp[1], bp[2]));
 
     vec3 controls[3];
     bezier_degree = get_reduced_control_points(vec3[3](bp[0], bp[1], bp[2]), controls);

manimlib/utils/space_ops.py

@@ -136,25 +136,11 @@ def z_to_vector(vector):
     Returns some matrix in SO(3) which takes the z-axis to the
     (normalized) vector provided as an argument
     """
-    norm = get_norm(vector)
-    if norm == 0:
+    cp = cross(OUT, vector)
+    if get_norm(cp) == 0:
         return np.identity(3)
-    v = np.array(vector) / norm
-    phi = np.arccos(v[2])
-    if any(v[:2]):
-        # projection of vector to unit circle
-        axis_proj = v[:2] / get_norm(v[:2])
-        theta = np.arccos(axis_proj[0])
-        if axis_proj[1] < 0:
-            theta = -theta
-    else:
-        theta = 0
-    phi_down = np.array([
-        [math.cos(phi), 0, math.sin(phi)],
-        [0, 1, 0],
-        [-math.sin(phi), 0, math.cos(phi)]
-    ])
-    return np.dot(rotation_about_z(theta), phi_down)
+    angle = np.arccos(np.dot(OUT, normalize(vector)))
+    return rotation_matrix(angle, axis=cp)
 
 
 def angle_of_vector(vector):
@@ -196,8 +182,19 @@ def cross(v1, v2):
     ])
 
 
-def get_unit_normal(v1, v2):
-    return normalize(cross(v1, v2))
+def get_unit_normal(v1, v2, tol=1e-6):
+    v1 = normalize(v1)
+    v2 = normalize(v2)
+    cp = cross(v1, v2)
+    cp_norm = get_norm(cp)
+    if cp_norm < tol:
+        # Vectors align, so find a normal to them in the plane shared with the z-axis
+        new_cp = cross(cross(v1, OUT), v1)
+        new_cp_norm = get_norm(new_cp)
+        if new_cp_norm < tol:
+            return RIGHT
+        return new_cp / new_cp_norm
+    return cp / cp_norm
 
 
 ###