Replaced peak filter with improved correlation.

Instead of dealing with noise/errors, having good data out of the
correlation algorithm is a much better option:
* Increase tolerance for fundamental matrix RANSAC.
  This should prevent overfitting and give more points for P3P.
* Do not grow corridor as much as affine projection. Stick to ground
  truth found on lower scales.
* Decrease correlation threshold and instead use cross-check filter to
  to exclude accidental matches.

All these changes have significantly reduced errors with sky and other
objects, to the point that now there's too much data available.

README.md

@@ -25,7 +25,7 @@ Download a release distribution from [releases](/zlogic/cybervision/releases).
 Run cybervision:
 
 ```shell
-cybervision [--scale=<scale>] [--focal-length=<focal-length>] [--mode=<cpu|gpu>] [--interpolation=<none|delaunay>] [--min-angle-cos=<cos-value>] [--projection=<parallel|perspective>] [--mesh=<plain|vertex-colors|texture-coordinates>] [--no-bundle-adjustment] <img1> <img2> [<imgn>] <output>
+cybervision [--scale=<scale>] [--focal-length=<focal-length>] [--mode=<cpu|gpu>] [--interpolation=<none|delaunay>] [--projection=<parallel|perspective>] [--mesh=<plain|vertex-colors|texture-coordinates>] [--no-bundle-adjustment] <img1> <img2> [<imgn>] <output>
 ```
 
 `--scale=<scale>` is an optional argument to specify a depth scale, for example `--scale=-10.0`.
@@ -41,10 +41,6 @@ If not specified, EXIF metadata will be used.
 `--interpolation=<none|delaunay>` is an optional argument to specify an interpolation mode, for example `--interpolation=none` or `--interpolation=delaunay`.
 `none` means that interpolation is disabled, `delaunay` uses [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation).
 
-`--min-angle-cos=<cos-value>` is an optional argument to specify the minimum cos value between a polygon normal and the camera forward vector, for example `--min-angle-cos=0.1`.
-Increasing this ensures that steep slopes will be discarded; setting this to 0.0 will keep all polygons and produce an "extruded" look.
-If not specified, will use the default 0.1 value.
-
 `--projection=<parallel|perspective>` is an optional argument to specify a projection mode, for example `--projection=parallel` or `--projection=perspective`.
 `parallel` projection should be used for images from a scanning electron microscope, `perspective` should be used for photos from a regular camera.
 

src/correlation/gpu/mod.rs

@@ -22,7 +22,7 @@ use crate::data::Grid;
 use nalgebra::Matrix3;
 
 use crate::correlation::{
-    CorrelationDirection, CorrelationParameters, HardwareMode, ProjectionMode, CORRIDOR_MIN_RANGE,
+    CorrelationDirection, CorrelationParameters, HardwareMode, ProjectionMode,
     CROSS_CHECK_SEARCH_AREA, KERNEL_SIZE, NEIGHBOR_DISTANCE,
 };
 
@@ -115,6 +115,7 @@ pub struct GpuContext<'a> {
     corridor_segment_length: usize,
     search_area_segment_length: usize,
     corridor_size: usize,
+    corridor_min_range: f64,
     corridor_extend_range: f64,
 
     device_context: &'a mut DefaultDeviceContext,
@@ -149,6 +150,7 @@ impl GpuContext<'_> {
             min_stdev: params.min_stdev,
             correlation_threshold: params.correlation_threshold,
             corridor_size: params.corridor_size,
+            corridor_min_range: params.corridor_min_range,
             corridor_extend_range: params.corridor_extend_range,
             fundamental_matrix,
             img1_dimensions,
@@ -269,7 +271,7 @@ impl GpuContext<'_> {
             min_stdev: self.min_stdev,
             neighbor_distance: NEIGHBOR_DISTANCE as u32,
             extend_range: self.corridor_extend_range as f32,
-            min_range: CORRIDOR_MIN_RANGE as f32,
+            min_range: self.corridor_min_range as f32,
         };
 
         let device = self.device_context.device_mut()?;

src/correlation/mod.rs

@@ -18,16 +18,17 @@ const KERNEL_WIDTH: usize = KERNEL_SIZE * 2 + 1;
 const KERNEL_POINT_COUNT: usize = KERNEL_WIDTH * KERNEL_WIDTH;
 
 const THRESHOLD_AFFINE: f32 = 0.6;
-const THRESHOLD_PERSPECTIVE: f32 = 0.8;
+const THRESHOLD_PERSPECTIVE: f32 = 0.5;
 const MIN_STDEV_AFFINE: f32 = 1.0;
-const MIN_STDEV_PERSPECTIVE: f32 = 3.0;
+const MIN_STDEV_PERSPECTIVE: f32 = 1.0;
 const CORRIDOR_SIZE_AFFINE: usize = 2;
 const CORRIDOR_SIZE_PERSPECTIVE: usize = 4;
 const NEIGHBOR_DISTANCE: usize = 10;
 const CORRIDOR_EXTEND_RANGE_AFFINE: f64 = 1.0;
-const CORRIDOR_EXTEND_RANGE_PERSPECTIVE: f64 = 1.0;
-const CORRIDOR_MIN_RANGE: f64 = 2.5;
-const CROSS_CHECK_SEARCH_AREA: usize = 2;
+const CORRIDOR_EXTEND_RANGE_PERSPECTIVE: f64 = 0.25;
+const CORRIDOR_MIN_RANGE_AFFINE: f64 = 2.5;
+const CORRIDOR_MIN_RANGE_PERSPECTIVE: f64 = 0.75;
+const CROSS_CHECK_SEARCH_AREA: usize = 4;
 
 type Match = (Point2D<u32>, f32);
 
@@ -66,6 +67,7 @@ pub struct PointCorrelations<'a> {
     min_stdev: f32,
     corridor_size: usize,
     correlation_threshold: f32,
+    corridor_min_range: f64,
     corridor_extend_range: f64,
     fundamental_matrix: Matrix3<f64>,
     gpu_context: Option<GpuContext<'a>>,
@@ -102,30 +104,39 @@ struct CorrelationParameters {
     min_stdev: f32,
     corridor_size: usize,
     correlation_threshold: f32,
+    corridor_min_range: f64,
     corridor_extend_range: f64,
 }
 
 impl CorrelationParameters {
     fn for_projection(projection_mode: &ProjectionMode) -> CorrelationParameters {
-        let (min_stdev, correlation_threshold, corridor_size, corridor_extend_range) =
-            match projection_mode {
-                ProjectionMode::Affine => (
-                    MIN_STDEV_AFFINE,
-                    THRESHOLD_AFFINE,
-                    CORRIDOR_SIZE_AFFINE,
-                    CORRIDOR_EXTEND_RANGE_AFFINE,
-                ),
-                ProjectionMode::Perspective => (
-                    MIN_STDEV_PERSPECTIVE,
-                    THRESHOLD_PERSPECTIVE,
-                    CORRIDOR_SIZE_PERSPECTIVE,
-                    CORRIDOR_EXTEND_RANGE_PERSPECTIVE,
-                ),
-            };
+        let (
+            min_stdev,
+            correlation_threshold,
+            corridor_size,
+            corridor_min_range,
+            corridor_extend_range,
+        ) = match projection_mode {
+            ProjectionMode::Affine => (
+                MIN_STDEV_AFFINE,
+                THRESHOLD_AFFINE,
+                CORRIDOR_SIZE_AFFINE,
+                CORRIDOR_MIN_RANGE_AFFINE,
+                CORRIDOR_EXTEND_RANGE_AFFINE,
+            ),
+            ProjectionMode::Perspective => (
+                MIN_STDEV_PERSPECTIVE,
+                THRESHOLD_PERSPECTIVE,
+                CORRIDOR_SIZE_PERSPECTIVE,
+                CORRIDOR_MIN_RANGE_PERSPECTIVE,
+                CORRIDOR_EXTEND_RANGE_PERSPECTIVE,
+            ),
+        };
         CorrelationParameters {
             min_stdev,
             corridor_size,
             correlation_threshold,
+            corridor_min_range,
             corridor_extend_range,
         }
     }
@@ -182,6 +193,7 @@ impl PointCorrelations<'_> {
             min_stdev: params.min_stdev,
             corridor_size: params.corridor_size,
             correlation_threshold: params.correlation_threshold,
+            corridor_min_range: params.corridor_min_range,
             corridor_extend_range: params.corridor_extend_range,
             fundamental_matrix,
             gpu_context,
@@ -518,7 +530,7 @@ impl PointCorrelations<'_> {
 
         let corridor_center = mid_corridor.round() as usize;
         let corridor_length =
-            (CORRIDOR_MIN_RANGE + range_stdev * self.corridor_extend_range).round() as usize;
+            (self.corridor_min_range + range_stdev * self.corridor_extend_range).round() as usize;
         let corridor_start = corridor_center
             .saturating_sub(corridor_length)
             .clamp(corridor_start, corridor_end);

src/fundamentalmatrix.rs

@@ -20,11 +20,11 @@ const RANSAC_K_PERSPECTIVE: usize = 1_000_000;
 const RANSAC_N_AFFINE: usize = 4;
 const RANSAC_N_PERSPECTIVE: usize = 7;
 const RANSAC_T_AFFINE: f64 = 0.1;
-const RANSAC_T_PERSPECTIVE: f64 = 1.0 / 1000.0;
+const RANSAC_T_PERSPECTIVE: f64 = 10.0 / 1000.0;
 const RANSAC_D_AFFINE: usize = 10;
 const RANSAC_D_PERSPECTIVE: usize = 200;
 const RANSAC_D_EARLY_EXIT_AFFINE: usize = 1000;
-const RANSAC_D_EARLY_EXIT_PERSPECTIVE: usize = 10_000;
+const RANSAC_D_EARLY_EXIT_PERSPECTIVE: usize = 50_000;
 const RANSAC_CHECK_INTERVAL: usize = 50_000;
 const RANSAC_RANK_EPSILON_AFFINE: f64 = 0.001;
 const RANSAC_RANK_EPSILON_PERSPECTIVE: f64 = 0.001;

src/main.rs

@@ -41,7 +41,6 @@ pub struct Args {
     focal_length: Option<u32>,
     mode: HardwareMode,
     interpolation: InterpolationMode,
-    min_angle_cos: f64,
     no_bundle_adjustment: bool,
     projection: ProjectionMode,
     mesh: Mesh,
@@ -58,7 +57,6 @@ Options:\
 \n      --focal-length=<FOCAL_LENGTH>    Focal length in 35mm equivalent\
 \n      --mode=<MODE>                    Hardware mode [default: gpu] [possible values: gpu, gpu-low-power, cpu]\
 \n      --interpolation=<INTERPOLATION>  Interpolation mode [default: delaunay] [possible values: delaunay, none]\
-\n      --min-angle-cos=<COS_VALUE>      Minimum cos value for polygon normals [default: 0.2]\
 \n      --no-bundle-adjustment           Skip bundle adjustment [if unspecified, bundle adjustment will be applied]\
 \n      --projection=<PROJECTION>        Projection mode [default: perspective] [possible values: parallel, perspective]\
 \n      --mesh=<MESH>                    Mesh options [default: vertex-colors] [possible values: plain, vertex-colors, texture-coordinates]\
@@ -70,7 +68,6 @@ impl Args {
             focal_length: None,
             mode: HardwareMode::Gpu,
             interpolation: InterpolationMode::Delaunay,
-            min_angle_cos: 0.1,
             no_bundle_adjustment: false,
             projection: ProjectionMode::Perspective,
             mesh: Mesh::VertexColors,
@@ -135,11 +132,6 @@ impl Args {
                             exit(2);
                         }
                     };
-                } else if name == "--min-angle-cos" {
-                    match value.parse() {
-                        Ok(min_angle_cos) => args.min_angle_cos = min_angle_cos,
-                        Err(err) => fail_with_error(name, value, &err),
-                    };
                 } else if name == "--projection" {
                     match value {
                         "perspective" => args.projection = ProjectionMode::Perspective,

src/output.rs

@@ -335,7 +335,6 @@ impl PolygonIndex {
 
 struct Mesh {
     interpolation: InterpolationMode,
-    min_angle_cos: f64,
     points: triangulation::Surface,
     polygons: Vec<Polygon>,
     polygon_index: PolygonIndex,
@@ -351,7 +350,6 @@ impl Mesh {
         let point_normals = vec![Vector3::zeros(); surface.tracks_len()];
         let mut surface = Mesh {
             interpolation: configuration.interpolation,
-            min_angle_cos: configuration.min_angle_cos,
             points: surface,
             polygons: vec![],
             polygon_index: PolygonIndex::new(),
@@ -382,39 +380,6 @@ impl Mesh {
         Some((point0, point1, point2))
     }
 
-    #[inline]
-    fn max_angle_cos(&self, polygon: &Polygon) -> Option<f64> {
-        let (point0, point1, point2) = self.get_polygon_points(polygon)?;
-        let polygon_normal = (point1 - point0).cross(&(point2 - point0)).normalize();
-        let point0_projections = self.points.get_camera_points(polygon.vertices[0]);
-        let point1_projections = self.points.get_camera_points(polygon.vertices[1]);
-        let point2_projections = self.points.get_camera_points(polygon.vertices[2]);
-        let projections_count = point0_projections
-            .len()
-            .min(point1_projections.len())
-            .min(point2_projections.len());
-
-        let affine_projection = self.points.cameras_len() == 0;
-
-        (0..projections_count)
-            .filter_map(|camera_i| {
-                if point0_projections[camera_i].is_none()
-                    || point1_projections[camera_i].is_none()
-                    || point2_projections[camera_i].is_none()
-                {
-                    return None;
-                }
-                let look_direction = if affine_projection {
-                    Vector3::new(0.0, 0.0, 1.0)
-                } else {
-                    self.points.camera_look_direction(camera_i)
-                };
-                let cos_angle = look_direction.dot(&polygon_normal);
-                Some(cos_angle)
-            })
-            .reduce(|a, b| a.max(b))
-    }
-
     fn polygon_normal(&self, polygon: &Polygon) -> Vector3<f64> {
         let (point0, point1, point2) = if let Some(points) = self.get_polygon_points(polygon) {
             points
@@ -478,50 +443,6 @@ 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()
@@ -535,21 +456,7 @@ impl Mesh {
                     return None;
                 }
 
-                let polygon_normal = self.polygon_normal(&polygon).normalize();
-                // Discard polygons that are too steep.
-                // 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)
-                }
+                Some(polygon)
             })
             .collect::<Vec<_>>();
         drop(triangulated_surface);
@@ -689,7 +596,6 @@ impl Mesh {
 
 pub struct MeshConfiguration {
     pub interpolation: InterpolationMode,
-    pub min_angle_cos: f64,
     pub vertex_mode: VertexMode,
 }
 

src/reconstruction.rs

@@ -186,7 +186,6 @@ impl SourceImage {
 
 struct ImageReconstruction {
     interpolation_mode: output::InterpolationMode,
-    min_angle_cos: f64,
     projection_mode: fundamentalmatrix::ProjectionMode,
     vertex_mode: output::VertexMode,
     triangulation: triangulation::Triangulation,
@@ -231,7 +230,6 @@ pub fn reconstruct(args: &Args) -> Result<(), ReconstructionError> {
         crate::ProjectionMode::Perspective => out_scale,
     };
     let focal_length = args.focal_length;
-    let min_angle_cos = args.min_angle_cos;
 
     let triangulation_projection = match args.projection {
         crate::ProjectionMode::Parallel => triangulation::ProjectionMode::Affine,
@@ -247,7 +245,6 @@ pub fn reconstruct(args: &Args) -> Result<(), ReconstructionError> {
 
     let mut reconstruction_task = ImageReconstruction {
         interpolation_mode,
-        min_angle_cos,
         projection_mode,
         vertex_mode,
         triangulation,
@@ -779,7 +776,6 @@ impl ImageReconstruction {
             output_filename,
             output::MeshConfiguration {
                 interpolation: self.interpolation_mode,
-                min_angle_cos: self.min_angle_cos,
                 vertex_mode: self.vertex_mode,
             },
             Some(&pb),