cluster-cone 0.1.0

Generate cones for clusters of triangles using the algorithm from https://github.com/zeux/meshoptimizer
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144

pub use nalgebra_glm as glm;

#[derive(Debug)]
pub struct Triangle {
    pub p0: glm::Vec3,
    pub p1: glm::Vec3,
    pub p2: glm::Vec3,
    pub normal: glm::Vec3,
}

impl Triangle {
    pub fn from_points(p0: &glm::Vec3, p1: &glm::Vec3, p2: &glm::Vec3) -> Triangle {
        Triangle {
            p0: *p0,
            p1: *p1,
            p2: *p2,
            normal: (p1 - p0).cross(&(p2 - p0)).normalize(),
        }
    }
}

#[derive(Debug)]
pub struct Bounds {
    pub center: glm::Vec3,
    pub radius: f32,
    pub cone_apex: glm::Vec3,
    pub cone_axis: glm::Vec3,
    pub cone_cutoff: f32,
}

pub fn compute_bounding_sphere(points: &[glm::Vec3]) -> (glm::Vec3, f32) {
    let mut pmin = [0; 3];
    let mut pmax = [0; 3];

    for (idx, p) in points.iter().enumerate() {
        for axis in 0..3 {
            if p[axis] < points[pmin[axis]][axis] {
                pmin[axis] = idx
            }

            if p[axis] > points[pmax[axis]][axis] {
                pmax[axis] = idx
            }
        }
    }

    let mut paxisd2 = 0.0;
    let mut paxis = 0;

    for axis in 0..3 {
        let p1 = points[pmin[axis]];
        let p2 = points[pmax[axis]];

        let d2 = glm::length2(&(p2 - p1));

        if d2 > paxisd2 {
            paxisd2 = d2;
            paxis = axis;
        }
    }

    let p1 = points[pmin[paxis]];
    let p2 = points[pmax[paxis]];

    let mut center = (p1 + p2) * 0.5;
    let mut radius = paxisd2.sqrt() * 0.5;

    for p in points.iter() {
        let d2 = glm::length2(&(p - center));

        if d2 > radius * radius {
            let d = d2.sqrt();

            let k = 0.5 + (radius / d) * 0.5;

            center = center * k + p * (1.0 - k);
            radius = (radius + d) * 0.5;
        }
    }

    (center, radius)
}

pub fn compute_cone(triangles: &[Triangle]) -> Bounds {
    let mut corners = vec![];
    let mut normals = vec![];

    for tri in triangles.iter() {
        corners.push(tri.p0);
        corners.push(tri.p1);
        corners.push(tri.p2);
        normals.push(tri.normal);
    }

    let point_sphere = compute_bounding_sphere(&corners);
    let normal_sphere = compute_bounding_sphere(&normals);

    let center = point_sphere.0;
    let axis = glm::normalize(&normal_sphere.0);

    let mut min_dot = 1.0;

    for n in normals.iter() {
        let dot = glm::dot(&n, &axis);
        if dot < min_dot {
            min_dot = dot;
        }
    }

    // degenerate cluster, normal cone is larger than a hemisphere => trivial accept
    // note that if mindp is positive but close to 0, the triangle intersection code below gets less stable
    // we arbitrarily decide that if a normal cone is ~168 degrees wide or more, the cone isn't useful
    if min_dot <= 0.1 {
        return Bounds {
            center,
            radius: point_sphere.1,
            cone_apex: glm::vec3(0.0, 0.0, 0.0),
            cone_axis: glm::vec3(0.0, 0.0, 0.0),
            cone_cutoff: 0.0,
        };
    }

    let mut maxt = 0.0;
    for tri in triangles.iter() {
        for p in [tri.p0, tri.p1, tri.p2].iter() {
            let c = center - p;
            let dc = glm::dot(&c, &tri.normal);
            let dn = glm::dot(&axis, &tri.normal);
            let t = dc / dn;

            if t > maxt {
                maxt = t;
            }
        }
    }

    Bounds {
        center,
        radius: point_sphere.1,
        cone_apex: center - axis * maxt,
        cone_axis: axis,
        cone_cutoff: (1.0 - min_dot * min_dot).sqrt(),
    }
}