Experiment: drop polygons using adjacent normals.

If any polygon's normal is different from its neighbors, drop it.

src/output.rs

@@ -478,6 +478,50 @@ impl Mesh {
             pl.report_status(0.9 * percent_complete);
         }
 
+        let mut point_normals = vec![vec![]; triangulated_surface.num_vertices()];
+
+        triangulated_surface.inner_faces().for_each(|f| {
+            let vertices = f.vertices();
+            let v0 = vertices[0].data().track_i;
+            let v1 = vertices[1].data().track_i;
+            let v2 = vertices[2].data().track_i;
+            // TODO: use 3D points directly
+            let polygon = Polygon::new(camera_i, [v0, v1, v2]);
+            let polygon_normal = self.polygon_normal(&polygon).normalize();
+            f.vertices()
+                .iter()
+                .for_each(|vertex| point_normals[vertex.index()].push(polygon_normal));
+        });
+        let point_dominant_normal = point_normals
+            .par_iter()
+            .map(|normals| {
+                if normals.len() < 2 {
+                    return None;
+                }
+                let (min_angle_cos, dominant_direction) = normals
+                    .iter()
+                    .enumerate()
+                    .take(normals.len() - 1)
+                    .flat_map(|(i, normal_i)| {
+                        normals
+                            .iter()
+                            .skip(i + 1)
+                            .map(|normal_j| {
+                                let angle_cos = normal_i.dot(&normal_j);
+                                let dominant_direction = (normal_i + normal_j).normalize();
+                                (angle_cos, dominant_direction)
+                            })
+                            .reduce(|a, b| if a.0 > b.0 { a } else { b })
+                    })
+                    .reduce(|a, b| if a.0 > b.0 { a } else { b })?;
+
+                if min_angle_cos > self.min_angle_cos {
+                    Some(dominant_direction)
+                } else {
+                    None
+                }
+            })
+            .collect::<Vec<_>>();
         let mut new_polygons = triangulated_surface
             .inner_faces()
             .par_bridge()
@@ -491,11 +535,17 @@ impl Mesh {
                     return None;
                 }
 
+                let polygon_normal = self.polygon_normal(&polygon).normalize();
                 // Discard polygons that are too steep.
-                if self
-                    .max_angle_cos(&polygon)
-                    .map_or(true, |angle_cos| angle_cos < self.min_angle_cos)
-                {
+                // TODO: precalculate normals once, then read a yes/no value
+                let angle_too_steep = vertices.iter().any(|vertex| {
+                    if let Some(point_normal) = point_dominant_normal[vertex.index()] {
+                        point_normal.dot(&polygon_normal) < self.min_angle_cos
+                    } else {
+                        true
+                    }
+                });
+                if angle_too_steep {
                     None
                 } else {
                     Some(polygon)