vtk-pure-rs 0.2.0

Pure Rust visualization toolkit — data structures, filters, I/O, rendering
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

use crate::data::{AnyDataArray, DataArray, PolyData};

/// Compute circumradius for each triangle.
///
/// R = (a*b*c)/(4*area). Adds "Circumradius" cell data.
/// Also computes the ratio R/r (circumradius/inradius) as a quality metric.
pub fn circumradius(input: &PolyData) -> PolyData {
    let mut radii = Vec::new();
    let mut ratios = Vec::new();

    for cell in input.polys.iter() {
        if cell.len()<3 { radii.push(0.0); ratios.push(0.0); continue; }
        let v0=input.points.get(cell[0] as usize);
        let v1=input.points.get(cell[1] as usize);
        let v2=input.points.get(cell[2] as usize);

        let a=dist(v0,v1); let b=dist(v1,v2); let c=dist(v2,v0);
        let s=(a+b+c)*0.5;
        let area=(s*(s-a)*(s-b)*(s-c)).max(0.0).sqrt();

        let cr = if area>1e-15{a*b*c/(4.0*area)}else{0.0};
        let ir = if s>1e-15{area/s}else{0.0};
        let ratio = if ir>1e-15{cr/ir}else{0.0};

        radii.push(cr);
        ratios.push(ratio);
    }

    let mut pd=input.clone();
    pd.cell_data_mut().add_array(AnyDataArray::F64(DataArray::from_vec("Circumradius", radii, 1)));
    pd.cell_data_mut().add_array(AnyDataArray::F64(DataArray::from_vec("CircumInRatio", ratios, 1)));
    pd
}

fn dist(a:[f64;3],b:[f64;3])->f64{((a[0]-b[0]).powi(2)+(a[1]-b[1]).powi(2)+(a[2]-b[2]).powi(2)).sqrt()}

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

    #[test]
    fn equilateral_circumradius() {
        let h=(3.0f64).sqrt()/2.0;
        let mut pd = PolyData::new();
        pd.points.push([0.0,0.0,0.0]); pd.points.push([1.0,0.0,0.0]); pd.points.push([0.5,h,0.0]);
        pd.polys.push_cell(&[0,1,2]);

        let result=circumradius(&pd);
        let arr=result.cell_data().get_array("Circumradius").unwrap();
        let mut buf=[0.0f64];
        arr.tuple_as_f64(0,&mut buf);
        // Equilateral triangle circumradius = a/sqrt(3)
        assert!((buf[0]-1.0/3.0f64.sqrt()).abs()<0.01);
    }

    #[test]
    fn ratio_equilateral() {
        let h=(3.0f64).sqrt()/2.0;
        let mut pd = PolyData::new();
        pd.points.push([0.0,0.0,0.0]); pd.points.push([1.0,0.0,0.0]); pd.points.push([0.5,h,0.0]);
        pd.polys.push_cell(&[0,1,2]);

        let result=circumradius(&pd);
        let arr=result.cell_data().get_array("CircumInRatio").unwrap();
        let mut buf=[0.0f64];
        arr.tuple_as_f64(0,&mut buf);
        assert!((buf[0]-2.0).abs()<0.01); // equilateral R/r = 2
    }

    #[test]
    fn empty_input() {
        let pd = PolyData::new();
        let result=circumradius(&pd);
        assert_eq!(result.polys.num_cells(), 0);
    }
}