chull/
wrapper.rs

1use super::convex::{ConvexHull, ErrorKind};
2use super::util::det;
3use num_bigint::{BigInt, ToBigInt};
4use num_traits::Float;
5
6#[derive(Clone, Debug)]
7/// A wrapper that holds a given points as a bigint.
8/// This can compute a convex hull robustly, although the computation is slower.
9pub struct ConvexHullWrapper<T: Float + ToBigInt> {
10    pub inner: ConvexHull<BigInt>,
11    pub conversion_factor: T,
12}
13
14impl<T: Float + ToBigInt> ConvexHullWrapper<T> {
15    pub fn try_new(points: &[Vec<T>], max_iter: Option<usize>) -> Result<Self, ErrorKind> {
16        let conversion_factor =T::epsilon();
17        let int_points:Vec<_> = to_bigint_points(points, conversion_factor);
18        let convex_hull = ConvexHull::try_new(&int_points, 0, max_iter)?;
19        Ok(
20            Self{
21                inner: convex_hull,
22                conversion_factor,
23            }
24        )
25    }
26    pub fn add_points(
27        &mut self,
28        points: &[Vec<T>],
29    ) -> Result<(), ErrorKind>{
30        let int_points = to_bigint_points(points, self.conversion_factor);
31        self.inner.add_points(&int_points, 0)
32    }
33    pub fn vertices_indices(&self) -> (Vec<Vec<T>>, Vec<usize>) {
34        let (int_v, i) = self.inner.vertices_indices();
35        let v = to_float_points(&int_v, self.conversion_factor);
36        (v,i)
37    }
38    pub fn volume(&self) -> T {
39        let dim = self.inner.points[0].len();
40        let (c_hull_vertices, c_hull_indices) = self.vertices_indices();
41        let mut reference_point = c_hull_vertices[c_hull_indices[0]].to_vec();
42        reference_point.push(T::one());
43        let mut volume = T::zero();
44        for i in (dim..c_hull_indices.len()).step_by(dim) {
45            let mut mat = Vec::new();
46            for j in 0..dim {
47                let mut row = c_hull_vertices[c_hull_indices[i + j]].to_vec();
48                row.push(T::one());
49                mat.push(row);
50            }
51            mat.push(reference_point.to_vec());
52            volume = volume + det(&mat);
53        }
54        let factorial = {
55            let mut result = T::one();
56            let mut m = T::one() + T::one();
57            let mut n = dim;
58            while n > 1 {
59                result = result * m.clone();
60                n = n - 1;
61                m = m + T::one();
62            }
63            result
64        };
65        volume / factorial
66    }
67    pub fn support_point(&self, direction: &[T]) -> Result<Vec<T>, ErrorKind> {
68        let d = to_bigint_points(&[direction.to_vec()], self.conversion_factor);
69        let p = to_float_points(&[self.inner.support_point(&d[0])?], self.conversion_factor);
70        Ok(p[0].to_vec())
71    }
72}
73
74fn to_bigint_points<T: Float+ToBigInt>(points: &[Vec<T>], conversion_factor: T) ->Vec<Vec<BigInt>>{
75    let mut int_points = Vec::new();
76    for point in points {
77        int_points.push(
78            point
79                .iter()
80                .map(|&p| p / conversion_factor)
81                .map(|p| p.to_bigint().expect("cannot convert to bigint"))
82                .collect::<Vec<_>>(),
83        );
84    }
85    int_points
86}
87
88fn to_float_points<T: Float+ToBigInt>(points: &[Vec<BigInt>], conversion_factor: T) ->Vec<Vec<T>>{
89    let mut float_points = Vec::new();
90    for point in points {
91        float_points.push(
92            point
93                .iter()
94                .map(|p| T::from(p.clone()).expect("cannot convert to bigint") * conversion_factor)
95                .collect::<Vec<_>>(),
96        );
97    }
98    float_points
99}
100
101#[test]
102fn wrapper_cube_test() {
103    let p1 = vec![1.0, 1.0, 1.0];
104    let p2 = vec![1.0, 1.0, -1.0];
105    let p3 = vec![1.0, -1.0, 1.0];
106    let p4 = vec![1.0, -1.0, -1.0];
107    let p5 = vec![-1.0, 1.0, 1.0];
108    let p6 = vec![-1.0, 1.0, -1.0];
109    let p7 = vec![-1.0, -1.0, 1.0];
110    let p8 = vec![-1.0, -1.0, -1.0];
111    let (_v, i) = ConvexHullWrapper::try_new(&[p1, p2, p3, p4, p5, p6, p7, p8], None)
112        .unwrap()
113        .vertices_indices();
114    assert_eq!(i.len(), 6 * 2 * 3);
115}
116
117#[test]
118fn wrapper_cube_volume_test() {
119    let p1 = vec![2.0, 2.0, 2.0];
120    let p2 = vec![2.0, 2.0, 0.0];
121    let p3 = vec![2.0, 0.0, 2.0];
122    let p4 = vec![2.0, 0.0, 0.0];
123    let p5 = vec![0.0, 2.0, 2.0];
124    let p6 = vec![0.0, 2.0, 0.0];
125    let p7 = vec![0.0, 0.0, 2.0];
126    let p8 = vec![0.0, 0.0, 0.0];
127    let cube = ConvexHullWrapper::try_new(&[p1, p2, p3, p4, p5, p6, p7, p8], None).unwrap();
128    assert_eq!(cube.volume(), 8.0);
129}
chull/wrapper.rs

chull/
wrapper.rs
