ph_faces/
lib.rs

1#![doc(html_root_url = "https://docs.rs/polyhedron-faces/0.4.3")]
2//! polyhedron faces for Rust
3//!
4
5pub mod polyhedron;
6pub use polyhedron::*;
7
8use num::Float;
9// use qm::v::TVector;
10use qm::m::{TMatrix, m3::Matrix3};
11pub use qm::{prec_eq, prec_eq_f};
12
13/// sum of vec [F; 2] without trait Sum
14/// when use += need trait Float + std::ops::AddAssign &lt; [F; 2] &gt;
15pub fn sum_f2<F: Float>(uvs: &Vec<[F; 2]>) -> Vec<F> {
16  uvs.iter().fold(vec![<F>::from(0).unwrap(); 2], |s, p|
17    s.iter().zip(p.iter()).map(|(&q, &p)| q + p).collect())
18}
19
20/// avg of vec [F; 2]
21pub fn avg_f2<F: Float>(uvs: &Vec<[F; 2]>) -> Vec<F> {
22  let n = <F>::from(uvs.len()).unwrap();
23  sum_f2(uvs).iter().map(|&p| p / n).collect()
24}
25
26/// center indexed uv [F; 2]
27pub fn center_indexed_uv<F: Float>(idx: &[u16], uvs: &Vec<[F; 2]>) -> [F; 2] {
28  let p = avg_f2(&idx.iter().map(|&i| uvs[i as usize]).collect());
29  p.as_slice().try_into().unwrap()
30}
31
32/// sum of vec [F; 3] without trait Sum
33/// when use += need trait Float + std::ops::AddAssign &lt; [F; 3] &gt;
34pub fn sum_f3<F: Float>(vs: &Vec<[F; 3]>) -> Vec<F> {
35  vs.iter().fold(vec![<F>::from(0).unwrap(); 3], |s, p|
36    s.iter().zip(p.iter()).map(|(&q, &p)| q + p).collect())
37}
38
39/// avg of vec [F; 3]
40pub fn avg_f3<F: Float>(vs: &Vec<[F; 3]>) -> Vec<F> {
41  let n = <F>::from(vs.len()).unwrap();
42  sum_f3(vs).iter().map(|&v| v / n).collect()
43}
44
45/// center indexed [F; 3]
46pub fn center_indexed<F: Float>(idx: &[u16], vtx: &Vec<[F; 3]>) -> [F; 3] {
47  let p = avg_f3(&idx.iter().map(|&i| vtx[i as usize]).collect());
48  p.as_slice().try_into().unwrap()
49}
50
51/// sum of vec [F; 4] without trait Sum
52/// when use += need trait Float + std::ops::AddAssign &lt; [F; 4] &gt;
53pub fn sum_f4<F: Float>(vs: &Vec<[F; 4]>) -> Vec<F> {
54  vs.iter().fold(vec![<F>::from(0).unwrap(); 4], |s, p|
55    s.iter().zip(p.iter()).map(|(&q, &p)| q + p).collect())
56}
57
58/// avg of vec [F; 4]
59pub fn avg_f4<F: Float>(vs: &Vec<[F; 4]>) -> Vec<F> {
60  let n = <F>::from(vs.len()).unwrap();
61  sum_f4(vs).iter().map(|&v| v / n).collect()
62}
63
64/// center indexed [F; 4]
65pub fn center_indexed_f4<F: Float>(idx: &[u16], vtx: &Vec<[F; 4]>) -> [F; 4] {
66  let p = avg_f4(&idx.iter().map(|&i| vtx[i as usize]).collect());
67  p.as_slice().try_into().unwrap()
68}
69
70/// divide internally
71pub fn divide_int<F: Float>(p: &[F; 3], q: &[F; 3], m: i32, n: i32) -> [F; 3] {
72  let mf = <F>::from(m).unwrap();
73  let nf = <F>::from(n).unwrap();
74  let sf = mf + nf;
75  p.iter().zip(q.iter()).map(|(&a, &b)|
76    (nf * a + mf * b) / sf).collect::<Vec<_>>().as_slice().try_into().unwrap()
77}
78
79/// divide externally
80pub fn divide_ext<F: Float>(p: &[F; 3], q: &[F; 3], m: i32, n: i32) -> [F; 3] {
81  divide_int(p, q, m, -n)
82}
83
84/// solve quadratic equation
85pub fn solve<F: Float>(a: F, b: F, c: F) -> [F; 2] {
86  let z = <F>::from(2).unwrap();
87  let d = b * b - <F>::from(4).unwrap() * a * c;
88  let r = d.sqrt();
89  [(-b - r) / (z * a), (-b + r) / (z * a)]
90}
91
92/// solve quadratic equation (b2 = b / 2)
93pub fn solve2<F: Float>(a: F, b2: F, c: F) -> [F; 2] {
94  let d = b2 * b2 - a * c;
95  let r = d.sqrt();
96  [(-b2 - r) / a, (-b2 + r) / a]
97}
98
99/// sol
100pub fn sol<F: Float>(r: &[F], s: F, e: F) -> F {
101  for i in 0..r.len() {
102    if r[i] >= s && r[i] <= e { return r[i]; }
103  }
104  <F>::from(0).unwrap()
105}
106
107/// calc cg f2 x axis
108pub fn calc_cg_f2_x<F: Float>(vs: &Vec<[F; 2]>) -> Vec<F> {
109  let o = <F>::from(0).unwrap();
110  let z = <F>::from(2).unwrap();
111  let mut stk = (1..vs.len()).into_iter().map(|vn| {
112    let (a, b) = (vs[vn - 1], vs[vn]);
113    [(b[0] - a[0]) * (b[1] + a[1]) / z, o]
114  }).collect::<Vec<_>>();
115  stk.push([o, o]);
116  for vn in 0..vs.len() {
117    stk[vn][1] = if vn == 0 { o } else { stk[vn - 1][0] + stk[vn - 1][1] };
118  }
119  let middle = stk[stk.len() - 1][1] / z;
120  let mut n = 0;
121  for vn in 0..vs.len() { if stk[vn][1] >= middle { n = vn; break; } }
122  let m = (middle - stk[n - 1][1]) / (stk[n][1] - stk[n - 1][1]);
123  let (h0, h1) = (vs[n - 1][0], vs[n][0]);
124  let (w0, w1) = (vs[n - 1][1], vs[n][1]);
125  let h = h1 - h0;
126  // vec![m * h + h0, o] // fast approximation (not accurate)
127  let r = solve2(w1 - w0, h * w0, -m * h * h * (w1 + w0));
128  vec![sol(&r, o, h) + h0, o] // TODO: y = o
129}
130
131/// calc cg f3 CAUTION not accurate
132/// (depends on density of vertices because no care of mass volume)
133/// - vs: length &ge; 3
134/// - p: precision for equality
135pub fn calc_cg_f3<F: Float>(vs: &Vec<[F; 3]>, p: F) -> Vec<F> {
136  let mut vtmp = Vec::<[F; 3]>::new();
137  for i in 0..vs.len() {
138    let mut found = false;
139    for v in vtmp.iter() { if prec_eq(v, p, &vs[i]) { found = true; break; } }
140    if !found { vtmp.push(vs[i]); }
141  }
142  round_prec(&avg_f3(&vtmp), p, <F>::from(0).unwrap()) // not accurate
143}
144
145/// calc cg with volume
146/// - idx: index of triangles on each faces
147/// - vtx: length &ge; 3
148/// - p: precision for equality
149/// when use += need trait Float + std::ops::AddAssign &lt; [F; 3] &gt;
150pub fn calc_cg_with_volume<F: Float + std::fmt::Debug>(
151  idx: &Vec<Vec<[u16; 3]>>, vtx: &Vec<[F; 3]>, p: F) -> (Vec<F>, F)
152  where F: std::iter::Sum {
153  let o = <F>::from(0).unwrap();
154  let mut m_total = o; // 6 * volume
155  let mut moi = [o, o, o];
156  for f in idx.iter() {
157    for t in f.iter() {
158      let vs = (0..t.len()).into_iter().map(|i|
159        vtx[t[i] as usize]).collect::<Vec<_>>();
160      let c = calc_cg_o(&vs);
161      // let m = vs[2].dot(&vs[0].cross(&vs[1])); // iter::Sum
162      let m = Matrix3::rowmajor3(vs).det(); // iter::Sum
163      m_total = m_total + m; // +=
164      for j in 0..moi.len() { moi[j] = moi[j] + m * c[j]; } // +=
165    }
166  }
167  let cg = moi.into_iter().map(|p| p / m_total).collect::<Vec<_>>();
168  // println!("CG: {:?}", cg);
169  (round_prec(&cg, p, <F>::from(0).unwrap()), m_total / <F>::from(6).unwrap())
170}
171
172/// calc cg
173/// - idx: index of triangles on each faces
174/// - vtx: length &ge; 3
175/// - p: precision for equality
176/// when use += need trait Float + std::ops::AddAssign &lt; [F; 3] &gt;
177pub fn calc_cg<F: Float + std::fmt::Debug>(
178  idx: &Vec<Vec<[u16; 3]>>, vtx: &Vec<[F; 3]>, p: F) -> Vec<F>
179  where F: std::iter::Sum {
180  let (cg, _vol) = calc_cg_with_volume(idx, vtx, p);
181  cg
182}
183
184/// calc cg skip o
185/// - vs: length = 3 (It means: (o + vs[0] + vs[1] + vs[2]) / 4)
186pub fn calc_cg_o<F: Float>(vs: &Vec<[F; 3]>) -> Vec<F> {
187  let n = <F>::from(4).unwrap(); // always 4
188  sum_f3(vs).iter().map(|&v| v / n).collect()
189}
190
191/// adjust cg with volume
192pub fn adjust_cg_with_volume<F: Float + std::fmt::Debug>(
193  idx: &Vec<Vec<[u16; 3]>>, vtx: &mut Vec<[F; 3]>, p: F) -> (Vec<F>, F)
194  where F: std::iter::Sum {
195  let (cg, vol) = calc_cg_with_volume(idx, vtx, p);
196  // println!("cg: {:?}", cg);
197  translate(vtx, &[-cg[0], -cg[1], -cg[2]]);
198/*
199  // double check after adjust cg
200  // TODO: calc cg and volume are not accurate (prec 1e-3) for some polyhedron
201  let (cg0, vol0) = calc_cg_with_volume(idx, vtx, p); // re calc for test
202  let q = <F>::from(1e3).unwrap();
203  assert!(prec_eq_f(vol0 / q, p, vol / q)); // expect
204  assert_eq!(f_to_f32(&cg0), &[0.0, 0.0, 0.0]); // expect
205*/
206  (cg, vol)
207}
208
209/// adjust cg
210pub fn adjust_cg<F: Float + std::fmt::Debug>(
211  idx: &Vec<Vec<[u16; 3]>>, vtx: &mut Vec<[F; 3]>, p: F) -> Vec<F>
212  where F: std::iter::Sum {
213  let (cg, _vol) = adjust_cg_with_volume(idx, vtx, p);
214  cg
215}
216
217/// translate
218pub fn translate<F: Float + std::fmt::Debug>(vtx: &mut Vec<[F; 3]>, o: &[F]) {
219  for v in vtx.iter_mut() {
220    *v = [v[0] + o[0], v[1] + o[1], v[2] + o[2]];
221  }
222}
223
224/// round precision
225pub fn round_prec<F: Float>(v: &[F], e: F, q: F) -> Vec<F> {
226  let o = <F>::from(0).unwrap();
227  v.iter().map(|&p| {
228    if (p - q).abs() >= e { p - q } else { o }
229  }).collect()
230}
231
232/// f_to_f32
233pub fn f_to_f32<F: Float>(v: &[F]) -> Vec<f32> {
234  v.iter().map(|i| i.to_f32().unwrap()).collect()
235}
236
237/// f_to_f64
238pub fn f_to_f64<F: Float>(v: &[F]) -> Vec<f64> {
239  v.iter().map(|i| i.to_f64().unwrap()).collect()
240}
241
242/// tests
243#[cfg(test)]
244mod tests {
245  use super::{prec_eq, f_to_f32};
246  use super::polyhedron::TUV;
247  // use super::polyhedron::tetra::*;
248  use super::tetra::*; // short cut
249  use super::cube::*;
250  use super::octa::*;
251  use super::pipe::*;
252  use super::polyhedron; // polyhedron::pin::Pin
253
254  /// [-- --nocapture] [-- --show-output]
255  #[test]
256  fn test_tetra() {
257    let tetra32_e = Tetra::new(1.0f32);
258    let tetra32 = tetra32_e.ph;
259    assert_eq!(tetra32.tri.len(), 4);
260    assert_eq!(tetra32.vtx.len(), 4);
261    let s32 = format!("{:?}", tetra32.vtx[2]);
262    println!("{}", s32);
263    assert_eq!(s32, "[0.0, 0.0, 0.61237246]"); // TODO: prec_eq
264
265    let tetra64_e = Tetra::new(1.0f64);
266    let tetra64 = tetra64_e.ph;
267    assert_eq!(tetra64.tri.len(), 4);
268    assert_eq!(tetra64.vtx.len(), 4);
269    let s64 = format!("{:?}", tetra64.vtx[2]);
270    println!("{}", s64);
271    assert_eq!(s64, "[0.0, 0.0, 0.6123724356957946]"); // TODO: prec_eq
272  }
273
274  #[test]
275  fn test_cube() {
276    let cube32_e = Cube::new(1.0f32);
277    let cube32 = cube32_e.ph;
278    assert_eq!(cube32.tri.len(), 6);
279    assert_eq!(cube32.tri[0].len(), 2);
280    assert_eq!(cube32.vtx.len(), 6 * 4);
281    let s32 = format!("{:?}", cube32.vtx[23]);
282    println!("{}", s32);
283    assert_eq!(s32, "[-1.0, -1.0, -1.0]");
284
285    let uv32 = format!("{:?}", cube32.with_uv(false)[0][0][0].puv());
286    println!("{}", uv32);
287    assert_eq!(uv32, "([1.0, -1.0, 1.0], [0.5, 0.0])");
288
289    let uv32 = format!("{:?}", cube32.with_uv(false)[5][1][2].puv());
290    println!("{}", uv32);
291    assert_eq!(uv32, "([-1.0, -1.0, -1.0], [1.0, 0.5])");
292
293    let uv32 = format!("{:?}", cube32.with_uv(true)[0][0][0].puv());
294    println!("{}", uv32);
295    assert_eq!(uv32, "([1.0, -1.0, 1.0], [0.25, 0.5])"); // not use 0.50
296
297    let uv32 = format!("{:?}", cube32.with_uv(true)[5][1][2].puv());
298    println!("{}", uv32);
299    assert_eq!(uv32, "([-1.0, -1.0, -1.0], [0.5, 0.75])"); // not use 0.50
300  }
301
302  #[test]
303  fn test_cube_center() {
304    let cube32_e = CubeCenter::new(1.0f32);
305    let cube32 = cube32_e.ph;
306    assert_eq!(cube32.tri.len(), 6);
307    assert_eq!(cube32.tri[0].len(), 4);
308    assert_eq!(cube32.vtx.len(), 6 * 4 + 6);
309    let s32 = format!("{:?}", cube32.vtx[23]);
310    println!("{}", s32);
311    assert_eq!(s32, "[-1.0, -1.0, -1.0]");
312
313    let s32 = format!("{:?}", cube32.vtx[29]);
314    println!("{}", s32);
315    assert_eq!(s32, "[0.0, 0.0, -1.0]");
316
317    let uv32 = format!("{:?}", cube32.with_uv(false)[0][0][0].puv());
318    println!("{}", uv32);
319    assert_eq!(uv32, "([1.0, 0.0, 0.0], [0.5, 0.5])");
320
321    let uv32 = format!("{:?}", cube32.with_uv(false)[5][3][2].puv());
322    println!("{}", uv32);
323    assert_eq!(uv32, "([-1.0, 1.0, -1.0], [0.5, 0.0])");
324
325    let uv32 = format!("{:?}", cube32.with_uv(true)[0][0][0].puv());
326    println!("{}", uv32);
327    assert_eq!(uv32, "([1.0, 0.0, 0.0], [0.375, 0.375])");
328
329    let uv32 = format!("{:?}", cube32.with_uv(true)[5][3][2].puv());
330    println!("{}", uv32);
331    assert_eq!(uv32, "([-1.0, 1.0, -1.0], [0.75, 0.75])");
332  }
333
334  #[test]
335  fn test_octa() {
336    let octa32_e = Octa::new(1.0f32);
337    let octa32 = octa32_e.ph;
338    assert_eq!(octa32.tri.len(), 8);
339    assert_eq!(octa32.vtx.len(), 6);
340    let s32 = format!("{:?}", octa32.vtx[0]);
341    println!("{}", s32);
342    assert_eq!(s32, "[1.0, 0.0, 0.0]");
343  }
344
345  #[test]
346  fn test_tube() {
347    let tube32_e = Tube::new(0.5, 0.4, 1.0, 6);
348    let tube32 = tube32_e.ph;
349    assert_eq!(tube32.tri.len(), 4 * (6*4));
350    assert_eq!(tube32.vtx.len(), (2*2+1) * (6*4));
351  }
352
353  #[test]
354  fn test_halfpipe() {
355//    let cmp = [0.0, 0.5, 0.2 - 0.05825041967286819]; // [0, l/2, idm/2 - cg]
356    let cmp = [0.0, 0.5, 0.2 - 0.06758362864954912]; // [0, l/2, idm/2 - cg]
357    let halfpipe32_e = HalfPipe::new(4.712388980, 0.5, 0.4, 1.0, 6); // 3pi/2
358/*
359    let cmp = [0.0, 0.5, 0.2]; // [0, l/2, idm/2]
360    let halfpipe32_e = HalfPipe::new(6.283185307, 0.5, 0.4, 1.0, 6); // 2pi
361*/
362    let halfpipe32 = halfpipe32_e.ph;
363    assert_eq!(halfpipe32.tri.len(), 4 * (6*4) + 2);
364    assert_eq!(halfpipe32.vtx.len(), 4 * (6*4+1));
365    let middle32 = halfpipe32.vtx[4 * (6*4/2) + 3]; // middle idm top
366    println!("{:?}", middle32);
367    assert!(prec_eq(&f_to_f32(&middle32), 1e-6, &cmp));
368  }
369
370  #[test]
371  fn test_pin() {
372    let pin32_e = polyhedron::pin::Pin::new(1.0f32, 8, 6);
373    let pin32 = pin32_e.ph;
374    assert_eq!(pin32.tri.len(), 432);
375    assert_eq!(pin32.vtx.len(), 410); // (8*2+1) * (6*4) + 2 (bottom, top)
376    let btm32 = pin32.vtx[(8*2+1) * (6*4)]; // bottom
377    println!("{:?}", btm32);
378//    assert!(prec_eq(&f_to_f32(&btm32), 1e-6, &[0.0, -5.8672757, 0.0]));
379    assert!(prec_eq(&f_to_f32(&btm32), 1e-6, &[0.0, -5.779917, 0.0]));
380  }
381}
ph_faces/lib.rs

ph_faces/
lib.rs
