Fixed bounding sphere generation

Cargo.lock

@@ -40,9 +40,9 @@ checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693"
 
 [[package]]
 name = "cfg-if"
-version = "0.1.10"
+version = "1.0.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822"
+checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
 
 [[package]]
 name = "cloudabi"
@@ -61,18 +61,18 @@ checksum = "a06f77d526c1a601b7c4cdd98f54b5eaabffc14d5f2f0296febdc7f357c6d3ba"
 
 [[package]]
 name = "generic-array"
-version = "0.12.3"
+version = "0.12.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c68f0274ae0e023facc3c97b2e00f076be70e254bc851d972503b328db79b2ec"
+checksum = "ffdf9f34f1447443d37393cc6c2b8313aebddcd96906caf34e54c68d8e57d7bd"
 dependencies = [
  "typenum",
 ]
 
 [[package]]
 name = "getrandom"
-version = "0.1.15"
+version = "0.1.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fc587bc0ec293155d5bfa6b9891ec18a1e330c234f896ea47fbada4cadbe47e6"
+checksum = "8fc3cb4d91f53b50155bdcfd23f6a4c39ae1969c2ae85982b135750cccaf5fce"
 dependencies = [
  "cfg-if",
  "libc",
@@ -81,9 +81,9 @@ dependencies = [
 
 [[package]]
 name = "libc"
-version = "0.2.80"
+version = "0.2.94"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4d58d1b70b004888f764dfbf6a26a3b0342a1632d33968e4a179d8011c760614"
+checksum = "18794a8ad5b29321f790b55d93dfba91e125cb1a9edbd4f8e3150acc771c1a5e"
 
 [[package]]
 name = "libm"
@@ -93,9 +93,9 @@ checksum = "c7d73b3f436185384286bd8098d17ec07c9a7d2388a6599f824d8502b529702a"
 
 [[package]]
 name = "matrixmultiply"
-version = "0.2.3"
+version = "0.2.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d4f7ec66360130972f34830bfad9ef05c6610a43938a467bcc9ab9369ab3478f"
+checksum = "916806ba0031cd542105d916a97c8572e1fa6dd79c9c51e7eb43a09ec2dd84c1"
 dependencies = [
  "rawpointer",
 ]
@@ -336,9 +336,9 @@ dependencies = [
 
 [[package]]
 name = "typenum"
-version = "1.12.0"
+version = "1.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "373c8a200f9e67a0c95e62a4f52fbf80c23b4381c05a17845531982fa99e6b33"
+checksum = "879f6906492a7cd215bfa4cf595b600146ccfac0c79bcbd1f3000162af5e8b06"
 
 [[package]]
 name = "wasi"

src/object/triangle.rs

@@ -160,6 +160,30 @@ impl Surface for TriangleMesh {
 
     // Uses Welzl's algorithm to solve the bounding sphere problem
     fn bound(&self) -> Bound {
+        fn smallest_sphere_plane(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f32) {
+            if points.len() == 0 || boundary.len() == 3 {
+                match boundary.len() {
+                    0 => (Point3::new(0.0, 0.0, 0.0), 0.0),
+                    1 => (*boundary[0], 0.0),
+                    2 => { let half_span = 0.5 * (boundary[1] - boundary[0]);
+                            (*boundary[0] + half_span, half_span.norm()) },
+                    3 => triangle_sphere(boundary[0], boundary[1], boundary[2]),
+                    _ => unreachable!()
+                }
+            } else {
+                let removed = points[0];
+                let points = Vec::from(&points[1..]);
+
+                let bound = smallest_sphere(points.clone(), boundary.clone());
+                if distance(&bound.0, removed) < bound.1 { return bound; }
+
+                let mut boundary = boundary.clone();
+                boundary.push(removed);
+
+                smallest_sphere_plane(points, boundary)
+            }
+        }
+
         fn triangle_sphere(point1: &Point3f, point2: &Point3f, point3: &Point3f) -> (Point3f, f32) {
             let a = point3 - point1;
             let b = point2 - point1;
@@ -181,22 +205,30 @@ impl Surface for TriangleMesh {
                                             point4.to_homogeneous().transpose()]);
 
             let a = matrix.determinant() * 2.0;
-            let mut matrix_mut = matrix.clone();
 
-            let squares = Vector4::new(point1.coords.norm_squared(), point2.coords.norm_squared(), point3.coords.norm_squared(), point4.coords.norm_squared());
+            if (a != 0.0) {
-            matrix_mut.set_column(0, &squares);
+                let mut matrix_mut = matrix.clone();
-            let center_x = matrix_mut.determinant();
 
-            matrix_mut.set_column(1, &matrix.index((.., 0)));
+                let squares = Vector4::new(point1.coords.norm_squared(), point2.coords.norm_squared(), point3.coords.norm_squared(), point4.coords.norm_squared());
-            let center_y = -matrix_mut.determinant();
+                matrix_mut.set_column(0, &squares);
+                let center_x = matrix_mut.determinant();
 
-            matrix_mut.set_column(2, &matrix.index((.., 1)));
+                matrix_mut.set_column(1, &matrix.index((.., 0)));
-            let center_z = matrix_mut.determinant();
+                let center_y = -matrix_mut.determinant();
 
-            let center = Point3::new(center_x / a, center_y / a, center_z / a);
+                matrix_mut.set_column(2, &matrix.index((.., 1)));
-            let radius = distance(point1, &center);
+                let center_z = matrix_mut.determinant();
 
-            (center, radius)
+                let center = Point3::new(center_x / a, center_y / a, center_z / a);
+                let radius = distance(point1, &center);
+
+                (center, radius)
+            } else {
+                let points = vec![point1, point2, point3, point4];
+                let boundary = Vec::new();
+
+                smallest_sphere_plane(points, boundary)
+            }
         }
 
         fn smallest_sphere(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f32) {