Compute better smooth normals for cheaper, maybe (#16050)

# Objective

Avoid a premature normalize operation and get better smooth normals for
it.

## Inspiration

@IceSentry suggested `face_normal()` could have its normalize removed
based on [this article](https://iquilezles.org/articles/normals/) in PR
#16039.

## Solution

I did not want to change `face_normal()` to return a vector that's not
normalized. The name "normal" implies it'll be normalized. Instead I
added the `face_area_normal()` function, whose result is not normalized.
Its magnitude is equal two times the triangle's area. I've noted why
this is the case in its doc comment.

I changed `compute_smooth_normals()` from computing normals from
adjacent faces with equal weight to use the area of the faces as a
weight. This has the benefit of being cheaper computationally and
hopefully produces better normals.

The `compute_flat_normals()` is unchanged and still uses
`face_normal()`.

## Testing

One test was added which shows the bigger triangle having an effect on
the normal, but the previous test that uses the same size triangles is
unchanged.

**WARNING:** No visual test has been done yet. No example exists that
demonstrates the compute_smooth_normals(). Perhaps there's a good model
to demonstrate what the differences are. I would love to have some input
on this.

I'd suggest @IceSentry and @stepancheg to review this PR.

## Further Considerations

It's possible weighting normals by their area is not definitely better
than unweighted. It's possible there may be aesthetic reasons to prefer
one over the other. In such a case, we could offer two variants:
weighted or unweighted. Or we could offer another function perhaps like
this: `compute_smooth_normals_with_weights(|normal, area| 1.0)` which
would restore the original unweighted sum of normals.

---

## Showcase

Smooth normal calculation now weights adjacent face normals by their
area.

## Migration Guide

crates/bevy_mesh/src/mesh.rs

@@ -2,10 +2,11 @@ use bevy_transform::components::Transform;
 pub use wgpu::PrimitiveTopology;
 
 use super::{
-    face_normal, generate_tangents_for_mesh, scale_normal, FourIterators, GenerateTangentsError,
-    Indices, MeshAttributeData, MeshTrianglesError, MeshVertexAttribute, MeshVertexAttributeId,
-    MeshVertexBufferLayout, MeshVertexBufferLayoutRef, MeshVertexBufferLayouts,
-    MeshWindingInvertError, VertexAttributeValues, VertexBufferLayout, VertexFormatSize,
+    face_area_normal, face_normal, generate_tangents_for_mesh, scale_normal, FourIterators,
+    GenerateTangentsError, Indices, MeshAttributeData, MeshTrianglesError, MeshVertexAttribute,
+    MeshVertexAttributeId, MeshVertexBufferLayout, MeshVertexBufferLayoutRef,
+    MeshVertexBufferLayouts, MeshWindingInvertError, VertexAttributeValues, VertexBufferLayout,
+    VertexFormatSize,
 };
 use alloc::collections::BTreeMap;
 use bevy_asset::{Asset, Handle, RenderAssetUsages};
@@ -698,7 +699,7 @@ impl Mesh {
             .chunks_exact(3)
             .for_each(|face| {
                 let [a, b, c] = [face[0], face[1], face[2]];
-                let normal = Vec3::from(face_normal(positions[a], positions[b], positions[c]));
+                let normal = Vec3::from(face_area_normal(positions[a], positions[b], positions[c]));
                 [a, b, c].iter().for_each(|pos| {
                     normals[*pos] += normal;
                 });
@@ -1402,6 +1403,41 @@ mod tests {
         assert_eq!([1., 0., 0.], normals[3]);
     }
 
+    #[test]
+    fn compute_smooth_normals_proportionate() {
+        let mut mesh = Mesh::new(
+            PrimitiveTopology::TriangleList,
+            RenderAssetUsages::default(),
+        );
+
+        //  z      y
+        //  |    /
+        //  3---2..
+        //  | /    \
+        //  0-------1---x
+
+        mesh.insert_attribute(
+            Mesh::ATTRIBUTE_POSITION,
+            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();
+        let normals = mesh
+            .attribute(Mesh::ATTRIBUTE_NORMAL)
+            .unwrap()
+            .as_float3()
+            .unwrap();
+        assert_eq!(4, normals.len());
+        // 0
+        assert_eq!(Vec3::new(1., 0., 2.).normalize().to_array(), normals[0]);
+        // 1
+        assert_eq!([0., 0., 1.], normals[1]);
+        // 2
+        assert_eq!(Vec3::new(1., 0., 2.).normalize().to_array(), normals[2]);
+        // 3
+        assert_eq!([1., 0., 0.], normals[3]);
+    }
+
     #[test]
     fn triangles_from_triangle_list() {
         let mut mesh = Mesh::new(

crates/bevy_mesh/src/vertex.rs

@@ -138,7 +138,29 @@ pub(crate) struct MeshAttributeData {
     pub(crate) values: VertexAttributeValues,
 }
 
-pub(crate) fn face_normal(a: [f32; 3], b: [f32; 3], c: [f32; 3]) -> [f32; 3] {
+/// Compute a vector whose direction is the normal of the triangle formed by
+/// points a, b, c, and whose magnitude is double the area of the triangle. This
+/// is useful for computing smooth normals where the contributing normals are
+/// proportionate to the areas of the triangles as [discussed
+/// here](https://iquilezles.org/articles/normals/).
+///
+/// Question: Why double the area? Because the area of a triangle _A_ is
+/// determined by this equation:
+///
+/// _A = |(b - a) x (c - a)| / 2_
+///
+/// By computing _2 A_ we avoid a division operation, and when calculating the
+/// 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] {
+    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] {
     let (a, b, c) = (Vec3::from(a), Vec3::from(b), Vec3::from(c));
     (b - a).cross(c - a).normalize().into()
 }