conspire 0.7.6

The Rust interface to conspire.
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

#[cfg(test)]
mod test;

use crate::{
    geometry::{Coordinate, Coordinates, Direction},
    math::{CrossProduct, FxHashSet, Quantity, SquareMatrix, Tensor, TensorRank1Vec, Vector},
    units::{Area, ReciprocalLength},
};

const D: usize = 3;
const N: usize = 3;

pub(crate) struct Jet {
    // Off until anisotropic remesh: the fitted surface normal, part of the
    // (u, v, n) curvature frame the per-vertex metric will need.
    // pub normal: Direction<D>,
    pub principal_curvatures: [Quantity<ReciprocalLength>; 2],
}

impl Jet {
    pub(crate) fn max_abs_curvature(&self) -> Quantity<ReciprocalLength> {
        self.principal_curvatures[0]
            .abs()
            .max(self.principal_curvatures[1].abs())
    }
}

pub(crate) fn fit_jet(
    center: &Coordinate<D>,
    neighbors: &Coordinates<D>,
    normal_guess: &Direction<D>,
) -> Option<Jet> {
    if neighbors.len() < 5 {
        return None;
    }
    let w = normal_guess.clone();
    let reference = if w[0].abs() < 0.9 {
        Direction::const_from([1.0, 0.0, 0.0])
    } else {
        Direction::const_from([0.0, 1.0, 0.0])
    };
    let u = w.cross(reference).normalized();
    let v = w.cross(&u);
    let mut normal_equations = SquareMatrix::zero(5);
    let mut rhs = Vector::zero(5);
    for neighbor in neighbors {
        let delta = neighbor - center;
        let (a, b, height) = (
            (&delta * &u).value(),
            (&delta * &v).value(),
            (&delta * &w).value(),
        );
        let basis = [a, b, a * a, a * b, b * b];
        for i in 0..5 {
            for j in 0..5 {
                normal_equations[i][j] += basis[i] * basis[j];
            }
            rhs[i] += basis[i] * height;
        }
    }
    let c = normal_equations.solve_lu(&rhs).ok()?;
    let (h_u, h_v) = (c[0], c[1]);
    let (h_uu, h_uv, h_vv) = (2.0 * c[2], c[3], 2.0 * c[4]);
    let area = 1.0 + h_u * h_u + h_v * h_v;
    let root = area.sqrt();
    let (big_e, big_f, big_g) = (1.0 + h_u * h_u, h_u * h_v, 1.0 + h_v * h_v);
    let (big_l, big_m, big_n) = (h_uu / root, h_uv / root, h_vv / root);
    let mean = (big_e * big_n - 2.0 * big_f * big_m + big_g * big_l) / (2.0 * area);
    let gauss = (big_l * big_n - big_m * big_m) / area;
    let spread = (mean * mean - gauss).max(0.0).sqrt();
    // Off until anisotropic remesh (see Jet::normal):
    // let normal = (&(&w - &(&u * h_u)) - &(&v * h_v)).normalized();
    Some(Jet {
        // normal,
        principal_curvatures: [Quantity::new(mean + spread), Quantity::new(mean - spread)],
    })
}

pub(crate) fn vertex_jets(
    connectivity: &[[usize; N]],
    coordinates: &Coordinates<D>,
    discarded: &FxHashSet<usize>,
) -> Vec<Option<Jet>> {
    let count = coordinates.len();
    let mut neighbors = vec![FxHashSet::default(); count];
    let mut normals = TensorRank1Vec::<D, _, Area>::zero(count);
    for &[a, b, c] in connectivity {
        for (i, j) in [(a, b), (b, c), (c, a)] {
            neighbors[i].insert(j);
            neighbors[j].insert(i);
        }
        let face = (&coordinates[b] - &coordinates[a]).cross(&coordinates[c] - &coordinates[a]);
        for vertex in [a, b, c] {
            normals[vertex] += &face;
        }
    }
    (0..count)
        .map(|vertex| {
            if discarded.contains(&vertex) {
                return None;
            }
            let mut ring = neighbors[vertex].clone();
            neighbors[vertex]
                .iter()
                .filter(|w| !discarded.contains(w))
                .for_each(|&w| ring.extend(&neighbors[w]));
            ring.remove(&vertex);
            ring.retain(|w| !discarded.contains(w));
            let points = ring.iter().map(|&w| coordinates[w].clone()).collect();
            fit_jet(
                &coordinates[vertex],
                &points,
                &normals[vertex].clone().normalized(),
            )
        })
        .collect()
}