Code shared from the Rust Playground

playground.rs

use nalgebra::{DMatrix, DVector};
use ndarray::{Array1, Array2, Axis};

struct Simplex;

impl Simplex {
    pub fn volume_lagrange(vertices: &Array2<f64>) -> f64 {
        let (num_pts, dim) = vertices.dim();
        assert_eq!(num_pts, dim + 1);
        // Explicitly annotate Vec<f64> to satisfy the compiler
        let mut data: Vec<f64> = Vec::with_capacity((dim + 1) * (dim + 1));
        for row in vertices.rows() {
            data.extend(row.iter());
            data.push(1.0);
        }
        let matrix = DMatrix::from_row_slice(dim + 1, dim + 1, &data);
        let n_fact: f64 = (1..=dim).product::<usize>() as f64;
        matrix.determinant().abs() / n_fact
    }

    pub fn volume_reduced(vertices: &Array2<f64>) -> f64 {
        let (num_pts, dim) = vertices.dim();
        let v0 = vertices.index_axis(Axis(0), 0);
        let mut reduced_data: Vec<f64> = Vec::with_capacity(dim * dim);
        for i in 1..num_pts {
            let vi = vertices.index_axis(Axis(0), i);
            let diff = &vi - &v0;
            reduced_data.extend(diff.iter());
        }
        let matrix = DMatrix::from_row_slice(dim, dim, &reduced_data);
        let n_fact: f64 = (1..=dim).product::<usize>() as f64;
        matrix.determinant().abs() / n_fact
    }

    pub fn signed_volume(vertices: &Array2<f64>) -> f64 {
        let (num_pts, dim) = vertices.dim();
        let v0 = vertices.index_axis(Axis(0), 0);
        let mut reduced_data: Vec<f64> = Vec::with_capacity(dim * dim);
        for i in 1..num_pts {
            let vi = vertices.index_axis(Axis(0), i);
            let diff = &vi - &v0;
            reduced_data.extend(diff.iter());
        }
        let matrix = DMatrix::from_row_slice(dim, dim, &reduced_data);
        let n_fact: f64 = (1..=dim).product::<usize>() as f64;
        matrix.determinant() / n_fact
    }

    pub fn centroid(vertices: &Array2<f64>) -> Array1<f64> {
        let num_vertices = vertices.nrows() as f64;
        vertices.sum_axis(Axis(0)) / num_vertices
    }

    pub fn circumcenter(vertices: &Array2<f64>) -> Option<Array1<f64>> {
        let (num_pts, dim) = vertices.dim();
        let v0 = vertices.index_axis(Axis(0), 0);

        let mut a_data: Vec<f64> = Vec::with_capacity(dim * dim);
        let mut b_data: Vec<f64> = Vec::with_capacity(dim);

        for i in 1..num_pts {
            let vi = vertices.index_axis(Axis(0), i);
            let diff = &vi - &v0;
            a_data.extend(diff.iter().map(|&x| 2.0 * x));
            let norm_sq = diff.iter().map(|&x| x * x).sum::<f64>();
            b_data.push(norm_sq);
        }

        let a_mat = DMatrix::from_row_slice(dim, dim, &a_data);
        let b_vec = DVector::from_vec(b_data);

        a_mat.lu().solve(&b_vec).map(|c_prime| {
            let mut result = v0.to_owned();
            for (j, &val) in c_prime.iter().enumerate() {
                result[j] += val;
            }
            result
        })
    }

    pub fn is_point_in_circumsphere(vertices: &Array2<f64>, p: &Array1<f64>) -> f64 {
        let (num_pts, dim) = vertices.dim();
        let matrix_dim = dim + 2;
        let mut data: Vec<f64> = Vec::with_capacity(matrix_dim * matrix_dim);

        for row in vertices.rows() {
            data.extend(row.iter());
            data.push(row.iter().map(|&x| x * x).sum());
            data.push(1.0);
        }
        data.extend(p.iter());
        data.push(p.iter().map(|&x| x * x).sum());
        data.push(1.0);

        let matrix = DMatrix::from_row_slice(matrix_dim, matrix_dim, &data);
        let det = matrix.determinant();
        let orientation = Self::signed_volume(vertices).signum();
        (det * orientation).signum()
    }
}

fn main() {
    let tetrahedron = ndarray::array![
        [0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0]
    ];

    println!("--- Simplex Stats ---");
    println!("Volume:       {:.6}", Simplex::volume_reduced(&tetrahedron));
    println!("Centroid:     {:?}", Simplex::centroid(&tetrahedron));
    
    if let Some(cc) = Simplex::circumcenter(&tetrahedron) {
        println!("Circumcenter: {:?}", cc);
    }

    let p_inside = ndarray::array![0.2, 0.2, 0.2];
    println!("In-Sphere (Inside):  {}", Simplex::is_point_in_circumsphere(&tetrahedron, &p_inside));
}

#[cfg(test)]
mod tests {
    use super::*;
    use ndarray::array;

    #[test]
    fn test_circumcenter_2d() {
        let triangle = array![[0.0, 0.0], [2.0, 0.0], [0.0, 2.0]];
        let cc = Simplex::circumcenter(&triangle).expect("Should find circumcenter");
        assert!((cc[0] - 1.0).abs() < 1e-10);
        assert!((cc[1] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_in_sphere_logic() {
        let triangle = array![[0.0, 0.0], [2.0, 0.0], [0.0, 2.0]];
        let p_inside = array![1.0, 0.5];
        let p_outside = array![3.0, 3.0];
        assert_eq!(Simplex::is_point_in_circumsphere(&triangle, &p_inside), 1.0);
        assert_eq!(Simplex::is_point_in_circumsphere(&triangle, &p_outside), -1.0);
    }
}