cluster_cone/
lib.rs

1pub use nalgebra_glm as glm;
2
3#[derive(Debug)]
4pub struct Triangle {
5    pub p0: glm::Vec3,
6    pub p1: glm::Vec3,
7    pub p2: glm::Vec3,
8    pub normal: glm::Vec3,
9}
10
11impl Triangle {
12    pub fn from_points(p0: &glm::Vec3, p1: &glm::Vec3, p2: &glm::Vec3) -> Triangle {
13        Triangle {
14            p0: *p0,
15            p1: *p1,
16            p2: *p2,
17            normal: (p1 - p0).cross(&(p2 - p0)).normalize(),
18        }
19    }
20}
21
22#[derive(Debug)]
23pub struct Bounds {
24    pub center: glm::Vec3,
25    pub radius: f32,
26    pub cone_apex: glm::Vec3,
27    pub cone_axis: glm::Vec3,
28    pub cone_cutoff: f32,
29}
30
31pub fn compute_bounding_sphere(points: &[glm::Vec3]) -> (glm::Vec3, f32) {
32    let mut pmin = [0; 3];
33    let mut pmax = [0; 3];
34
35    for (idx, p) in points.iter().enumerate() {
36        for axis in 0..3 {
37            if p[axis] < points[pmin[axis]][axis] {
38                pmin[axis] = idx
39            }
40
41            if p[axis] > points[pmax[axis]][axis] {
42                pmax[axis] = idx
43            }
44        }
45    }
46
47    let mut paxisd2 = 0.0;
48    let mut paxis = 0;
49
50    for axis in 0..3 {
51        let p1 = points[pmin[axis]];
52        let p2 = points[pmax[axis]];
53
54        let d2 = glm::length2(&(p2 - p1));
55
56        if d2 > paxisd2 {
57            paxisd2 = d2;
58            paxis = axis;
59        }
60    }
61
62    let p1 = points[pmin[paxis]];
63    let p2 = points[pmax[paxis]];
64
65    let mut center = (p1 + p2) * 0.5;
66    let mut radius = paxisd2.sqrt() * 0.5;
67
68    for p in points.iter() {
69        let d2 = glm::length2(&(p - center));
70
71        if d2 > radius * radius {
72            let d = d2.sqrt();
73
74            let k = 0.5 + (radius / d) * 0.5;
75
76            center = center * k + p * (1.0 - k);
77            radius = (radius + d) * 0.5;
78        }
79    }
80
81    (center, radius)
82}
83
84pub fn compute_cone(triangles: &[Triangle]) -> Bounds {
85    let mut corners = vec![];
86    let mut normals = vec![];
87
88    for tri in triangles.iter() {
89        corners.push(tri.p0);
90        corners.push(tri.p1);
91        corners.push(tri.p2);
92        normals.push(tri.normal);
93    }
94
95    let point_sphere = compute_bounding_sphere(&corners);
96    let normal_sphere = compute_bounding_sphere(&normals);
97
98    let center = point_sphere.0;
99    let axis = glm::normalize(&normal_sphere.0);
100
101    let mut min_dot = 1.0;
102
103    for n in normals.iter() {
104        let dot = glm::dot(&n, &axis);
105        if dot < min_dot {
106            min_dot = dot;
107        }
108    }
109
110    // degenerate cluster, normal cone is larger than a hemisphere => trivial accept
111    // note that if mindp is positive but close to 0, the triangle intersection code below gets less stable
112    // we arbitrarily decide that if a normal cone is ~168 degrees wide or more, the cone isn't useful
113    if min_dot <= 0.1 {
114        return Bounds {
115            center,
116            radius: point_sphere.1,
117            cone_apex: glm::vec3(0.0, 0.0, 0.0),
118            cone_axis: glm::vec3(0.0, 0.0, 0.0),
119            cone_cutoff: 0.0,
120        };
121    }
122
123    let mut maxt = 0.0;
124    for tri in triangles.iter() {
125        for p in [tri.p0, tri.p1, tri.p2].iter() {
126            let c = center - p;
127            let dc = glm::dot(&c, &tri.normal);
128            let dn = glm::dot(&axis, &tri.normal);
129            let t = dc / dn;
130
131            if t > maxt {
132                maxt = t;
133            }
134        }
135    }
136
137    Bounds {
138        center,
139        radius: point_sphere.1,
140        cone_apex: center - axis * maxt,
141        cone_axis: axis,
142        cone_cutoff: (1.0 - min_dot * min_dot).sqrt(),
143    }
144}
cluster_cone/lib.rs

cluster_cone/
lib.rs
