Add angle-weighted smooth normals implementation (#18383) (#18552)

# Objective

Closes #18383 

## Solution

Given the 2 votes (me and @komadori ) for making angle-weighted normals
the default, I went ahead and did so, moving the face-weighted
implementation to the new `Mesh::compute_face_weighted_normals` method.
I factored out the common code between both into
`Mesh::compute_custom_smooth_normals`, which I went ahead and made
public to make it easier for users to add any other weighting methods
they might come up with.

If any of these decisions are undesirable for any reason, please let me
know and I will gladly make modifications.

## Testing & Showcase

I made a demo that exaggerates the difference at
[Waridley/bevy_smooth_normals_comparison](https://github.com/Waridley/bevy_smooth_normals_comparison).
Screenshots included in the readme.

Another way it could be demonstrated is via scaling a mesh along its
normals, like for generating outline meshes with inverted faces. I'd be
willing to make a demo for that as well.

I also edited and renamed the `compute_smooth_normals` tests to use
face-weighted normals, and added a new test for angle-weighted ones
which validates that all normals of a unit cube have each component
equal to `(±1 / √3) ± f32::EPSILON`.

## Migration Guide

#16050 already did not mention a migration guide, and it is not even in
a stable release yet. If this lands in a 0.16 RC, updating from RC1
would probably not require any changes in the vast majority of cases,
anyway. If someone really needs face-weighted normals, they can switch
to `.compute_face_weighted_normals()` or
`.with_computed_face_weighted_normals()`. And if for some reason they
really liked the old count-weighted implementation from 0.15, there is
an example in the docs for `compute_custom_smooth_normals`.

Author: Waridley

benches/benches/bevy_render/compute_normals.rs

@@ -54,7 +54,7 @@ fn compute_normals(c: &mut Criterion) {
         });
     });
 
-    c.bench_function("face_weighted_normals", |b| {
+    c.bench_function("angle_weighted_normals", |b| {
         b.iter_custom(|iters| {
             let mut total = Duration::default();
             for _ in 0..iters {
@@ -70,11 +70,23 @@ fn compute_normals(c: &mut Criterion) {
         });
     });
 
-    let new_mesh = || {
-        new_mesh()
-            .with_duplicated_vertices()
-            .with_computed_flat_normals()
-    };
+    c.bench_function("face_weighted_normals", |b| {
+        b.iter_custom(|iters| {
+            let mut total = Duration::default();
+            for _ in 0..iters {
+                let mut mesh = new_mesh();
+                black_box(mesh.attribute(Mesh::ATTRIBUTE_NORMAL));
+                let start = Instant::now();
+                mesh.compute_area_weighted_normals();
+                let end = Instant::now();
+                black_box(mesh.attribute(Mesh::ATTRIBUTE_NORMAL));
+                total += end.duration_since(start);
+            }
+            total
+        });
+    });
+
+    let new_mesh = || new_mesh().with_duplicated_vertices();
 
     c.bench_function("flat_normals", |b| {
         b.iter_custom(|iters| {

crates/bevy_mesh/src/mesh.rs

@@ -2,7 +2,7 @@ use bevy_transform::components::Transform;
 pub use wgpu_types::PrimitiveTopology;
 
 use super::{
-    face_area_normal, face_normal, generate_tangents_for_mesh, scale_normal, FourIterators,
+    generate_tangents_for_mesh, scale_normal, triangle_area_normal, triangle_normal, FourIterators,
     GenerateTangentsError, Indices, MeshAttributeData, MeshTrianglesError, MeshVertexAttribute,
     MeshVertexAttributeId, MeshVertexBufferLayout, MeshVertexBufferLayoutRef,
     MeshVertexBufferLayouts, MeshWindingInvertError, VertexAttributeValues, VertexBufferLayout,
@@ -649,11 +649,7 @@ impl Mesh {
     ///
     /// # Panics
     /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
-    /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
-    ///
-    /// FIXME: This should handle more cases since this is called as a part of gltf
-    /// mesh loading where we can't really blame users for loading meshes that might
-    /// not conform to the limitations here!
+    /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].=
     pub fn compute_normals(&mut self) {
         assert!(
             matches!(self.primitive_topology, PrimitiveTopology::TriangleList),
@@ -695,7 +691,7 @@ impl Mesh {
 
         let normals: Vec<_> = positions
             .chunks_exact(3)
-            .map(|p| face_normal(p[0], p[1], p[2]))
+            .map(|p| triangle_normal(p[0], p[1], p[2]))
             .flat_map(|normal| [normal; 3])
             .collect();
 
@@ -705,22 +701,141 @@ impl Mesh {
     /// Calculates the [`Mesh::ATTRIBUTE_NORMAL`] of an indexed mesh, smoothing normals for shared
     /// vertices.
     ///
+    /// This method weights normals by the angles of the corners of connected triangles, thus
+    /// eliminating triangle area and count as factors in the final normal. This does make it
+    /// somewhat slower than [`Mesh::compute_area_weighted_normals`] which does not need to
+    /// greedily normalize each triangle's normal or calculate corner angles.
+    ///
+    /// If you would rather have the computed normals be weighted by triangle area, see
+    /// [`Mesh::compute_area_weighted_normals`] instead. If you need to weight them in some other
+    /// way, see [`Mesh::compute_custom_smooth_normals`].
+    ///
     /// # Panics
     /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
     /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
     /// Panics if the mesh does not have indices defined.
-    ///
-    /// FIXME: This should handle more cases since this is called as a part of gltf
-    /// mesh loading where we can't really blame users for loading meshes that might
-    /// not conform to the limitations here!
     pub fn compute_smooth_normals(&mut self) {
+        self.compute_custom_smooth_normals(|[a, b, c], positions, normals| {
+            let pa = Vec3::from(positions[a]);
+            let pb = Vec3::from(positions[b]);
+            let pc = Vec3::from(positions[c]);
+
+            let ab = pb - pa;
+            let ba = pa - pb;
+            let bc = pc - pb;
+            let cb = pb - pc;
+            let ca = pa - pc;
+            let ac = pc - pa;
+
+            const EPS: f32 = f32::EPSILON;
+            let weight_a = if ab.length_squared() * ac.length_squared() > EPS {
+                ab.angle_between(ac)
+            } else {
+                0.0
+            };
+            let weight_b = if ba.length_squared() * bc.length_squared() > EPS {
+                ba.angle_between(bc)
+            } else {
+                0.0
+            };
+            let weight_c = if ca.length_squared() * cb.length_squared() > EPS {
+                ca.angle_between(cb)
+            } else {
+                0.0
+            };
+
+            let normal = Vec3::from(triangle_normal(positions[a], positions[b], positions[c]));
+
+            normals[a] += normal * weight_a;
+            normals[b] += normal * weight_b;
+            normals[c] += normal * weight_c;
+        });
+    }
+
+    /// Calculates the [`Mesh::ATTRIBUTE_NORMAL`] of an indexed mesh, smoothing normals for shared
+    /// vertices.
+    ///
+    /// This method weights normals by the area of each triangle containing the vertex. Thus,
+    /// larger triangles will skew the normals of their vertices towards their own normal more
+    /// than smaller triangles will.
+    ///
+    /// This method is actually somewhat faster than [`Mesh::compute_smooth_normals`] because an
+    /// intermediate result of triangle normal calculation is already scaled by the triangle's area.
+    ///
+    /// If you would rather have the computed normals be influenced only by the angles of connected
+    /// edges, see [`Mesh::compute_smooth_normals`] instead. If you need to weight them in some
+    /// other way, see [`Mesh::compute_custom_smooth_normals`].
+    ///
+    /// # Panics
+    /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
+    /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
+    /// Panics if the mesh does not have indices defined.
+    pub fn compute_area_weighted_normals(&mut self) {
+        self.compute_custom_smooth_normals(|[a, b, c], positions, normals| {
+            let normal = Vec3::from(triangle_area_normal(
+                positions[a],
+                positions[b],
+                positions[c],
+            ));
+            [a, b, c].into_iter().for_each(|pos| {
+                normals[pos] += normal;
+            });
+        });
+    }
+
+    /// Calculates the [`Mesh::ATTRIBUTE_NORMAL`] of an indexed mesh, smoothing normals for shared
+    /// vertices.
+    ///
+    /// This method allows you to customize how normals are weighted via the `per_triangle` parameter,
+    /// which must be a function or closure that accepts 3 parameters:
+    /// - The indices of the three vertices of the triangle as a `[usize; 3]`.
+    /// - A reference to the values of the [`Mesh::ATTRIBUTE_POSITION`] of the mesh (`&[[f32; 3]]`).
+    /// - A mutable reference to the sums of all normals so far.
+    ///
+    /// See also the standard methods included in Bevy for calculating smooth normals:
+    /// - [`Mesh::compute_smooth_normals`]
+    /// - [`Mesh::compute_area_weighted_normals`]
+    ///
+    /// An example that would weight each connected triangle's normal equally, thus skewing normals
+    /// towards the planes divided into the most triangles:
+    /// ```
+    /// # use bevy_asset::RenderAssetUsages;
+    /// # use bevy_mesh::{Mesh, PrimitiveTopology, Meshable, MeshBuilder};
+    /// # use bevy_math::{Vec3, primitives::Cuboid};
+    /// # let mut mesh = Cuboid::default().mesh().build();
+    /// mesh.compute_custom_smooth_normals(|[a, b, c], positions, normals| {
+    ///     let normal = Vec3::from(bevy_mesh::triangle_normal(positions[a], positions[b], positions[c]));
+    ///     for idx in [a, b, c] {
+    ///         normals[idx] += normal;
+    ///     }
+    /// });
+    /// ```
+    ///
+    /// # Panics
+    /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
+    /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
+    /// Panics if the mesh does not have indices defined.
+    //
+    // FIXME: This should handle more cases since this is called as a part of gltf
+    // mesh loading where we can't really blame users for loading meshes that might
+    // not conform to the limitations here!
+    //
+    // When fixed, also update "Panics" sections of
+    // - [Mesh::compute_smooth_normals]
+    // - [Mesh::with_computed_smooth_normals]
+    // - [Mesh::compute_area_weighted_normals]
+    // - [Mesh::with_computed_area_weighted_normals]
+    pub fn compute_custom_smooth_normals(
+        &mut self,
+        mut per_triangle: impl FnMut([usize; 3], &[[f32; 3]], &mut [Vec3]),
+    ) {
         assert!(
             matches!(self.primitive_topology, PrimitiveTopology::TriangleList),
-            "`compute_smooth_normals` can only work on `TriangleList`s"
+            "smooth normals can only be computed on `TriangleList`s"
         );
         assert!(
             self.indices().is_some(),
-            "`compute_smooth_normals` can only work on indexed meshes"
+            "smooth normals can only be computed on indexed meshes"
         );
 
         let positions = self
@@ -736,16 +851,8 @@ impl Mesh {
             .iter()
             .collect::<Vec<usize>>()
             .chunks_exact(3)
-            .for_each(|face| {
-                let [a, b, c] = [face[0], face[1], face[2]];
-                let normal = Vec3::from(face_area_normal(positions[a], positions[b], positions[c]));
-                [a, b, c].iter().for_each(|pos| {
-                    normals[*pos] += normal;
-                });
-            });
+            .for_each(|face| per_triangle([face[0], face[1], face[2]], positions, &mut normals));
 
-        // average (smooth) normals for shared vertices...
-        // TODO: support different methods of weighting the average
         for normal in &mut normals {
             *normal = normal.try_normalize().unwrap_or(Vec3::ZERO);
         }
@@ -786,6 +893,10 @@ impl Mesh {
     ///
     /// (Alternatively, you can use [`Mesh::compute_smooth_normals`] to mutate an existing mesh in-place)
     ///
+    /// This method weights normals by the angles of triangle corners connected to each vertex. If
+    /// you would rather have the computed normals be weighted by triangle area, see
+    /// [`Mesh::with_computed_area_weighted_normals`] instead.
+    ///
     /// # Panics
     /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
     /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
@@ -796,6 +907,25 @@ impl Mesh {
         self
     }
 
+    /// Consumes the mesh and returns a mesh with calculated [`Mesh::ATTRIBUTE_NORMAL`].
+    ///
+    /// (Alternatively, you can use [`Mesh::compute_area_weighted_normals`] to mutate an existing mesh in-place)
+    ///
+    /// This method weights normals by the area of each triangle containing the vertex. Thus,
+    /// larger triangles will skew the normals of their vertices towards their own normal more
+    /// than smaller triangles will. If you would rather have the computed normals be influenced
+    /// only by the angles of connected edges, see [`Mesh::with_computed_smooth_normals`] instead.
+    ///
+    /// # Panics
+    /// Panics if [`Mesh::ATTRIBUTE_POSITION`] is not of type `float3`.
+    /// Panics if the mesh has any other topology than [`PrimitiveTopology::TriangleList`].
+    /// Panics if the mesh does not have indices defined.
+    #[must_use]
+    pub fn with_computed_area_weighted_normals(mut self) -> Self {
+        self.compute_area_weighted_normals();
+        self
+    }
+
     /// Generate tangents for the mesh using the `mikktspace` algorithm.
     ///
     /// Sets the [`Mesh::ATTRIBUTE_TANGENT`] attribute if successful.
@@ -1587,7 +1717,7 @@ mod tests {
     }
 
     #[test]
-    fn compute_smooth_normals() {
+    fn compute_area_weighted_normals() {
         let mut mesh = Mesh::new(
             PrimitiveTopology::TriangleList,
             RenderAssetUsages::default(),
@@ -1604,7 +1734,7 @@ mod tests {
             vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
         );
         mesh.insert_indices(Indices::U16(vec![0, 1, 2, 0, 2, 3]));
-        mesh.compute_smooth_normals();
+        mesh.compute_area_weighted_normals();
         let normals = mesh
             .attribute(Mesh::ATTRIBUTE_NORMAL)
             .unwrap()
@@ -1622,7 +1752,7 @@ mod tests {
     }
 
     #[test]
-    fn compute_smooth_normals_proportionate() {
+    fn compute_area_weighted_normals_proportionate() {
         let mut mesh = Mesh::new(
             PrimitiveTopology::TriangleList,
             RenderAssetUsages::default(),
@@ -1639,7 +1769,7 @@ mod tests {
             vec![[0., 0., 0.], [2., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
         );
         mesh.insert_indices(Indices::U16(vec![0, 1, 2, 0, 2, 3]));
-        mesh.compute_smooth_normals();
+        mesh.compute_area_weighted_normals();
         let normals = mesh
             .attribute(Mesh::ATTRIBUTE_NORMAL)
             .unwrap()
@@ -1656,6 +1786,59 @@ mod tests {
         assert_eq!([1., 0., 0.], normals[3]);
     }
 
+    #[test]
+    fn compute_angle_weighted_normals() {
+        // CuboidMeshBuilder duplicates vertices (even though it is indexed)
+
+        //   5---------4
+        //  /|        /|
+        // 1-+-------0 |
+        // | 6-------|-7
+        // |/        |/
+        // 2---------3
+        let verts = vec![
+            [1.0, 1.0, 1.0],
+            [-1.0, 1.0, 1.0],
+            [-1.0, -1.0, 1.0],
+            [1.0, -1.0, 1.0],
+            [1.0, 1.0, -1.0],
+            [-1.0, 1.0, -1.0],
+            [-1.0, -1.0, -1.0],
+            [1.0, -1.0, -1.0],
+        ];
+
+        let indices = Indices::U16(vec![
+            0, 1, 2, 2, 3, 0, // front
+            5, 4, 7, 7, 6, 5, // back
+            1, 5, 6, 6, 2, 1, // left
+            4, 0, 3, 3, 7, 4, // right
+            4, 5, 1, 1, 0, 4, // top
+            3, 2, 6, 6, 7, 3, // bottom
+        ]);
+        let mut mesh = Mesh::new(
+            PrimitiveTopology::TriangleList,
+            RenderAssetUsages::default(),
+        );
+        mesh.insert_attribute(Mesh::ATTRIBUTE_POSITION, verts);
+        mesh.insert_indices(indices);
+        mesh.compute_smooth_normals();
+
+        let normals = mesh
+            .attribute(Mesh::ATTRIBUTE_NORMAL)
+            .unwrap()
+            .as_float3()
+            .unwrap();
+
+        for new in normals.iter().copied().flatten() {
+            // std impl is unstable
+            const FRAC_1_SQRT_3: f32 = 0.57735026;
+            const MIN: f32 = FRAC_1_SQRT_3 - f32::EPSILON;
+            const MAX: f32 = FRAC_1_SQRT_3 + f32::EPSILON;
+            assert!(new.abs() >= MIN, "{new} < {MIN}");
+            assert!(new.abs() <= MAX, "{new} > {MAX}");
+        }
+    }
+
     #[test]
     fn triangles_from_triangle_list() {
         let mut mesh = Mesh::new(

crates/bevy_mesh/src/vertex.rs

@@ -219,16 +219,16 @@ impl SerializedMeshAttributeData {
 /// the sum of these vectors which are then normalized, a constant multiple has
 /// no effect.
 #[inline]
-pub fn face_area_normal(a: [f32; 3], b: [f32; 3], c: [f32; 3]) -> [f32; 3] {
+pub fn triangle_area_normal(a: [f32; 3], b: [f32; 3], c: [f32; 3]) -> [f32; 3] {
     let (a, b, c) = (Vec3::from(a), Vec3::from(b), Vec3::from(c));
     (b - a).cross(c - a).into()
 }
 
 /// Compute the normal of a face made of three points: a, b, and c.
 #[inline]
-pub fn face_normal(a: [f32; 3], b: [f32; 3], c: [f32; 3]) -> [f32; 3] {
+pub fn triangle_normal(a: [f32; 3], b: [f32; 3], c: [f32; 3]) -> [f32; 3] {
     let (a, b, c) = (Vec3::from(a), Vec3::from(b), Vec3::from(c));
-    (b - a).cross(c - a).normalize().into()
+    (b - a).cross(c - a).normalize_or_zero().into()
 }
 
 /// Contains an array where each entry describes a property of a single vertex.