Add normal_vertex and normal_triangle function

examples/normals.rs

@@ -0,0 +1,20 @@
+extern crate startin;
+
+fn main() {
+    let mut pts: Vec<[f64; 3]> = Vec::new();
+    pts.push([0.0, 0.0, 1.0]);
+    pts.push([10.0, 1.0, 2.0]);
+    pts.push([10.0, 12.0, 3.0]);
+    pts.push([0.0, 10.0, 4.0]);
+    pts.push([5.0, 5.0, 4.]);
+
+    let mut dt = startin::Triangulation::new();
+    dt.insert(&pts, startin::InsertionStrategy::AsIs);
+
+    println!("*****");
+    println!("Number of points in DT: {}", dt.number_of_vertices());
+    println!("Number of triangles in DT: {}", dt.number_of_triangles());
+
+    println!("{:?}", dt.normal_vertex(2));
+    println!("{:?}", dt.normal_vertex(5));
+}

src/geom/mod.rs

@@ -15,6 +15,20 @@ pub fn area_triangle(a: &[f64], b: &[f64], c: &[f64]) -> f64 {
     det3x3t(a, b, c) / 2.0
 }
 
+pub fn normal_triangle(a: &[f64], b: &[f64], c: &[f64], normalise: bool) -> Vec<f64> {
+    let v0: Vec<f64> = vec![a[0] - b[0], a[1] - b[1], a[2] - b[2]];
+    let v1: Vec<f64> = vec![a[0] - c[0], a[1] - c[1], a[2] - c[2]];
+    let nx = (v0[1] * v1[2]) - (v0[2] * v1[1]);
+    let ny = (v0[2] * v1[0]) - (v0[0] * v1[2]);
+    let nz = (v0[0] * v1[1]) - (v0[1] * v1[0]);
+    if normalise {
+        let norm = ((nx * nx) + (ny * ny) + (nz * nz)).sqrt();
+        vec![nx / norm, ny / norm, nz / norm]
+    } else {
+        vec![nx, ny, nz]
+    }
+}
+
 pub fn circle_centre(a: &[f64], b: &[f64], c: &[f64]) -> Vec<f64> {
     //-- nicked from http://www.ambrsoft.com/trigocalc/circle3d.htm
     let val_a = det3x3t(&[a[0], a[1], 1.0], &[b[0], b[1], 1.0], &[c[0], c[1], 1.0]);

src/lib.rs

@@ -902,6 +902,51 @@ impl Triangulation {
         }
     }
 
+    /// Returns the normal of the vertex with ID `vi`.
+    pub fn normal_vertex(&self, vi: usize) -> Result<Vec<f64>, StartinError> {
+        match self.incident_triangles_to_vertex(vi) {
+            Err(why) => Err(why),
+            Ok(trs) => {
+                let mut avgns: Vec<f64> = vec![0.0, 0.0, 0.0];
+                let mut no_finite_tr = 0;
+                for tr in &trs {
+                    // println!("{:?} {:?}", tr, self.is_finite(tr));
+                    if self.is_finite(tr) == false {
+                        continue;
+                    }
+                    no_finite_tr += 1;
+                    let n = geom::normal_triangle(
+                        &self.stars[tr.v[0]].pt,
+                        &self.stars[tr.v[1]].pt,
+                        &self.stars[tr.v[2]].pt,
+                        true,
+                    );
+                    for j in 0..3 {
+                        avgns[j] += n[j];
+                    }
+                }
+                for j in 0..3 {
+                    avgns[j] /= no_finite_tr as f64;
+                }
+                let norm =
+                    ((avgns[0] * avgns[0]) + (avgns[1] * avgns[1]) + (avgns[2] * avgns[2])).sqrt();
+                Ok(vec![avgns[0] / norm, avgns[1] / norm, avgns[2] / norm])
+            }
+        }
+    }
+
+    pub fn normal_triangle(&self, tr: &Triangle) -> Result<Vec<f64>, StartinError> {
+        match self.is_triangle(tr) {
+            false => Err(StartinError::TriangleNotPresent),
+            true => Ok(geom::normal_triangle(
+                &self.stars[tr.v[0]].pt,
+                &self.stars[tr.v[1]].pt,
+                &self.stars[tr.v[2]].pt,
+                true,
+            )),
+        }
+    }
+
     /// Returns the degree of the vertex with ID `vi`.
     pub fn degree(&self, vi: usize) -> Result<usize, StartinError> {
         match self.is_vertex_removed(vi) {