Click the display to select spheres (#25)

On the incoming branch, you can select a sphere by clicking it in the display. Holding *shift* while clicking enables multiple selection. These controls match the ones already implemented in the outline view.

Since the selection routine is now used in multiple places, the incoming branch factors it out into the `AppState::select` method.

Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: glen/dyna3#25
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>

app-proto/src/assembly.rs

@@ -1,4 +1,4 @@
-use nalgebra::{DMatrix, DVector};
+use nalgebra::{DMatrix, DVector, Vector3};
 use rustc_hash::FxHashMap;
 use slab::Slab;
 use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
@@ -65,6 +65,49 @@ impl Element {
             column_index: 0
         }
     }
+    
+    // the smallest positive depth, represented as a multiple of `dir`, where
+    // the line generated by `dir` hits the element (which is assumed to be a
+    // sphere). returns `None` if the line misses the sphere. this function
+    // should be kept synchronized with `sphere_cast` in `inversive.frag`, which
+    // does essentially the same thing on the GPU side
+    pub fn cast(&self, dir: Vector3<f64>, assembly_to_world: &DMatrix<f64>) -> Option<f64> {
+        // if `a/b` is less than this threshold, we approximate
+        // `a*u^2 + b*u + c` by the linear function `b*u + c`
+        const DEG_THRESHOLD: f64 = 1e-9;
+        
+        let rep = self.representation.with_untracked(|rep| assembly_to_world * rep);
+        let a = -rep[3] * dir.norm_squared();
+        let b = rep.rows_range(..3).dot(&dir);
+        let c = -rep[4];
+        
+        let adjust = 4.0*a*c/(b*b);
+        if adjust < 1.0 {
+            // as long as `b` is non-zero, the linear approximation of
+            //
+            //   a*u^2 + b*u + c
+            //
+            // at `u = 0` will reach zero at a finite depth `u_lin`. the root of
+            // the quadratic adjacent to `u_lin` is stored in `lin_root`. if
+            // both roots have the same sign, `lin_root` will be the one closer
+            // to `u = 0`
+            let square_rect_ratio = 1.0 + (1.0 - adjust).sqrt();
+            let lin_root = -(2.0*c)/b / square_rect_ratio;
+            if a.abs() > DEG_THRESHOLD * b.abs() {
+                if lin_root > 0.0 {
+                    Some(lin_root)
+                } else {
+                    let other_root = -b/(2.*a) * square_rect_ratio;
+                    (other_root > 0.0).then_some(other_root)
+                }
+            } else {
+                (lin_root > 0.0).then_some(lin_root)
+            }
+        } else {
+            // the line through `dir` misses the sphere completely
+            None
+        }
+    }
 }