efd/
tests.rs

1#![doc(hidden)]
2use crate::*;
3#[cfg(test)]
4use alloc::{vec, vec::Vec};
5#[cfg(test)]
6use approx::assert_abs_diff_eq;
7#[cfg(test)]
8use core::iter::zip;
9pub use util::dist_err as curve_diff;
10
11/// Epsilon for curve difference.
12pub const EPS: f64 = 2.2e-14;
13pub const RES: usize = 1000;
14
15#[test]
16fn error() {
17    let coeff = vec![
18        na::matrix![10., 20.; 20., 10.],
19        na::matrix![ 3., 4.; 4., 3.],
20    ];
21    let a = Efd2::from_coeffs_unchecked(coeff);
22    let coeff = vec![na::matrix![10., 20.; 20., 10.]];
23    let b = Efd2::from_coeffs_unchecked(coeff);
24    assert_eq!(a.square_err(&b), 50.);
25    assert_eq!(b.square_err(&a), 50.);
26}
27
28#[test]
29fn efd2d() {
30    let efd = Efd2::from_curve(CURVE2D, false);
31    assert!(!efd.is_open());
32    // Test starting point
33    let curve = CURVE2D
34        .iter()
35        .cycle()
36        .skip(CURVE2D.len() / 2)
37        .take(CURVE2D.len())
38        .copied()
39        .collect::<Vec<_>>();
40    let efd_half = Efd2::from_curve_harmonic(curve, false, efd.harmonic());
41    assert_abs_diff_eq!(efd.err(&efd_half), 0., epsilon = 1e-12);
42    assert_eq!(efd.harmonic(), 8);
43    // Test rotation
44    for ang in 0..6 {
45        let ang = core::f64::consts::TAU * ang as f64 / 6.;
46        let curve = GeoVar::from_rot(na::UnitComplex::new(ang)).transform(CURVE2D);
47        let efd_rot = Efd2::from_curve_harmonic(curve, false, efd.harmonic());
48        assert_abs_diff_eq!(efd_rot.err(&efd), 0., epsilon = 1e-12);
49    }
50    // Test transformation
51    let geo = efd.as_geo();
52    assert_abs_diff_eq!(geo.trans()[0], -1.248409055632358);
53    assert_abs_diff_eq!(geo.trans()[1], 55.26080122817753);
54    assert_abs_diff_eq!(geo.rot().angle(), -2.49925101855502);
55    assert_abs_diff_eq!(geo.scale(), 48.16765830752243);
56    // Test reconstruction
57    let target = efd.recon_norm_by(&get_norm_t(CURVE2D, false));
58    let curve = efd.as_geo().inverse().transform(CURVE2D);
59    assert_abs_diff_eq!(curve_diff(target, curve), 0., epsilon = 0.0029);
60}
61
62#[test]
63fn efd2d_open() {
64    let efd = Efd2::from_curve(CURVE2D_OPEN, true);
65    assert!(efd.is_open());
66    assert_eq!(efd.harmonic(), 14);
67    // Test rotation
68    for ang in 0..6 {
69        let ang = core::f64::consts::TAU * ang as f64 / 6.;
70        let curve = GeoVar::from_rot(na::UnitComplex::new(ang)).transform(CURVE2D_OPEN);
71        let efd_rot = Efd2::from_curve_harmonic(curve, true, efd.harmonic());
72        assert_abs_diff_eq!(efd_rot.err(&efd), 0., epsilon = 1e-12);
73    }
74    // Test transformation
75    let geo = efd.as_geo();
76    assert_abs_diff_eq!(geo.trans()[0], 43.03456791427352);
77    assert_abs_diff_eq!(geo.trans()[1], 48.107208358019015);
78    assert_abs_diff_eq!(geo.rot().angle(), 2.7330524299596815);
79    assert_abs_diff_eq!(geo.scale(), 33.930916934329495);
80    // Test reconstruction
81    let target = efd.recon_norm_by(&get_norm_t(CURVE2D_OPEN, true));
82    let curve = efd.as_geo().inverse().transform(CURVE2D_OPEN);
83    assert_abs_diff_eq!(curve_diff(target, curve), 0., epsilon = 0.0411);
84}
85
86#[test]
87fn efd3d() {
88    let efd = Efd3::from_curve(CURVE3D, false);
89    assert!(!efd.is_open());
90    // Test starting point
91    let curve = CURVE3D
92        .iter()
93        .cycle()
94        .skip(CURVE3D.len() / 2)
95        .take(CURVE3D.len())
96        .copied()
97        .collect::<Vec<_>>();
98    let efd_half = Efd3::from_curve_harmonic(curve, false, efd.harmonic());
99    assert_abs_diff_eq!(efd.err(&efd_half), 0., epsilon = 1e-12);
100    assert_eq!(efd.harmonic(), 5);
101    // Test rotation
102    for ang in 0..6 {
103        let ang = core::f64::consts::TAU * ang as f64 / 6.;
104        let curve = GeoVar::from_rot(na::UnitQuaternion::new(na::matrix![1.; 1.; 0.] * ang))
105            .transform(CURVE3D);
106        let efd_rot = Efd3::from_curve_harmonic(curve, false, efd.harmonic());
107        assert_abs_diff_eq!(efd.err(&efd_rot), 0., epsilon = 1e-12);
108    }
109    // Test transformation
110    let geo = efd.as_geo();
111    assert_abs_diff_eq!(geo.trans()[0], 0.7239345388499508);
112    assert_abs_diff_eq!(geo.trans()[1], 0.09100107896533066);
113    assert_abs_diff_eq!(geo.trans()[2], 0.49979194975846675);
114    assert_abs_diff_eq!(geo.rot().angle(), 2.9160714030359416);
115    assert_abs_diff_eq!(geo.scale(), 0.5629099155595344);
116    // Test reconstruction
117    let target = efd.recon_norm_by(&get_norm_t(CURVE3D, false));
118    let curve = efd.as_geo().inverse().transform(CURVE3D);
119    assert_abs_diff_eq!(curve_diff(target, curve), 0., epsilon = 0.0013);
120}
121
122#[test]
123fn posed_efd2d_open() {
124    let efd = PosedEfd2::from_angles(CURVE2D_POSE, ANGLE2D_POSE, true);
125    // Test rotation
126    for ang in 0..6 {
127        let ang = core::f64::consts::TAU * ang as f64 / 6.;
128        let curve = GeoVar::from_rot(na::UnitComplex::new(ang)).transform(CURVE2D_POSE);
129        let angles = ANGLE2D_POSE.iter().map(|a| a + ang).collect::<Vec<_>>();
130        let efd_rot = PosedEfd2::from_angles(curve, &angles, true);
131        assert_abs_diff_eq!(efd.err(&efd_rot), 0., epsilon = 1e-12);
132    }
133}
134
135#[test]
136#[cfg(feature = "std")]
137fn plot2d_harmonic() -> Result<(), Box<dyn std::error::Error>> {
138    const N: usize = 8;
139    let mut efd = Efd2::from_curve_harmonic(ORIGIN_CURVE, false, N);
140    *efd.as_geo_mut() *= GeoVar2::from_scale(20.);
141    plot2d(&efd, format!("img/2dh{N}.svg"))?;
142    for n in (1..N).rev() {
143        efd.set_harmonic(n);
144        plot2d(&efd, format!("img/2dh{n}.svg"))?;
145    }
146    Ok(())
147}
148
149#[test]
150#[cfg(feature = "std")]
151fn plot2d_closed() -> Result<(), Box<dyn std::error::Error>> {
152    let efd = Efd::from_coeffs_unchecked(vec![
153        na::matrix![12., 35.; 35., 13.],
154        na::matrix![5., 21.; 21., 5.],
155        na::matrix![1., 12.; 12., 1.],
156    ]);
157    plot2d(&efd, "img/2d.svg")
158}
159
160#[test]
161#[cfg(feature = "std")]
162fn plot2d_open() -> Result<(), Box<dyn std::error::Error>> {
163    let efd = Efd::from_coeffs_unchecked(vec![
164        na::matrix![35., 0.; 8., 0.],
165        na::matrix![10., 0.; 24., 0.],
166        na::matrix![5., 0.; -8., 0.],
167    ]);
168    plot2d(&efd, "img/2d_open.svg")
169}
170
171#[test]
172#[cfg(feature = "std")]
173fn plot3d_closed() -> Result<(), Box<dyn std::error::Error>> {
174    let efd = Efd::from_coeffs_unchecked(vec![
175        na::matrix![12., 22.; 35., 5.; 20., 21.],
176        na::matrix![21., 12.; 5., 12.; 1., 1.],
177        na::matrix![3., 3.; 7., 5.; 12., 21.],
178    ]);
179    plot3d(&efd, "img/3d.svg")
180}
181
182#[test]
183#[cfg(feature = "std")]
184fn plot3d_open() -> Result<(), Box<dyn std::error::Error>> {
185    let efd = Efd::from_coeffs_unchecked(vec![
186        na::matrix![16., 0.; 35., 0.; 27., 0.],
187        na::matrix![21., 0.; 8., 0.; 16., 0.],
188        na::matrix![3., 0.; 7., 0.; 12., 0.],
189    ]);
190    plot3d(&efd, "img/3d_open.svg")
191}
192
193#[cfg(all(test, feature = "std"))]
194fn get_area<const D: usize>(pts: &[[f64; D]]) -> [[f64; 2]; D] {
195    pts.iter()
196        .fold([[f64::INFINITY, f64::NEG_INFINITY]; D], |mut b, c| {
197            zip(&mut b, c).for_each(|(b, c)| *b = [b[0].min(*c), b[1].max(*c)]);
198            b
199        })
200}
201
202#[cfg(all(test, feature = "std"))]
203fn plot2d(efd: &Efd2, path: impl AsRef<std::path::Path>) -> Result<(), Box<dyn std::error::Error>> {
204    use plotters::prelude::*;
205
206    fn bounding_box(pts: &[[f64; 2]]) -> [f64; 4] {
207        let [[x_min, x_max], [y_min, y_max]] = get_area(pts);
208        let dx = (x_max - x_min).abs();
209        let dy = (y_max - y_min).abs();
210        if dx > dy {
211            let cen = (y_min + y_max) * 0.5;
212            let r = dx * 0.5;
213            let mg = dx * 0.1;
214            [x_min - mg, x_max + mg, cen - r - mg, cen + r + mg]
215        } else {
216            let cen = (x_min + x_max) * 0.5;
217            let r = dy * 0.5;
218            let mg = dy * 0.1;
219            [cen - r - mg, cen + r + mg, y_min - mg, y_max + mg]
220        }
221    }
222
223    let curve = efd.recon(360);
224    let [x_min, x_max, y_min, y_max] = bounding_box(&curve);
225    let b = SVGBackend::new(&path, (1200, 1200));
226    let root = b.into_drawing_area();
227    root.fill(&WHITE)?;
228    let mut chart = ChartBuilder::on(&root).build_cartesian_2d(x_min..x_max, y_min..y_max)?;
229    let p0 = curve[0];
230    chart.draw_series([Circle::new((p0[0], p0[1]), 3, BLACK.filled())])?;
231    for (p, color) in [((10., 0.), RED), ((0., 10.), BLUE)] {
232        chart.draw_series(LineSeries::new([(0., 0.), p], color.stroke_width(10)))?;
233    }
234    let geo0 = efd.as_geo();
235    let mut c0 = [0.; 2];
236    for m in efd.coeffs_iter() {
237        let geo = geo0 * GeoVar2::new(c0, na::UnitComplex::new(0.), 1.);
238        const N: usize = 100;
239        const N_F: f64 = N as f64;
240        let ellipse = (0..N)
241            .map(|i| {
242                let t = i as f64 * std::f64::consts::TAU / N_F;
243                m * na::matrix![t.cos(); t.sin()]
244            })
245            .map(|c| geo.transform_pt(c.data.0[0]).into());
246        let p1 = c0;
247        zip(&mut c0, &m.column(0)).for_each(|(c, u)| *c += u);
248        let p1 = geo0.transform_pt(p1).into();
249        let p2 = geo0.transform_pt(c0).into();
250        chart.draw_series([Circle::new(p2, 5, RED.filled())])?;
251        chart.draw_series(LineSeries::new([p1, p2], RED.stroke_width(5)))?;
252        chart.draw_series(LineSeries::new(ellipse, RED.stroke_width(7)))?;
253    }
254    chart.draw_series(LineSeries::new(
255        curve.into_iter().map(|c| c.into()),
256        BLACK.stroke_width(10),
257    ))?;
258    Ok(())
259}
260
261#[cfg(all(test, feature = "std"))]
262fn plot3d(efd: &Efd3, path: impl AsRef<std::path::Path>) -> Result<(), Box<dyn std::error::Error>> {
263    use plotters::prelude::*;
264
265    fn bounding_box(pts: &[[f64; 3]]) -> [f64; 6] {
266        let mut b = get_area(pts);
267        let center = b.map(|[min, max]| (min + max) * 0.5);
268        let width = b
269            .iter()
270            .map(|[min, max]| (max - min).abs())
271            .fold(0., f64::max);
272        for ([min, max], c) in zip(&mut b, center) {
273            *min = c - width * 0.5;
274            *max = c + width * 0.5;
275        }
276        let [[x_min, x_max], [y_min, y_max], [z_min, z_max]] = b;
277        [x_min, x_max, y_min, y_max, z_min, z_max]
278    }
279
280    let curve = efd.recon(360);
281    let [x_min, x_max, y_min, y_max, z_min, z_max] = bounding_box(&curve);
282    let b = SVGBackend::new(&path, (1200, 1200));
283    let root = b.into_drawing_area();
284    root.fill(&WHITE)?;
285    let mut chart =
286        ChartBuilder::on(&root).build_cartesian_3d(x_min..x_max, y_min..y_max, z_min..z_max)?;
287    chart.with_projection(|mut pb| {
288        pb.yaw = 0.3;
289        pb.scale = 1.4;
290        pb.into_matrix()
291    });
292    let p0 = curve[0];
293    chart.draw_series([Circle::new((p0[0], p0[1], p0[2]), 3, BLACK.filled())])?;
294    for (p, color) in [
295        ((10., 0., 0.), RED),
296        ((0., 10., 0.), BLUE),
297        ((0., 0., 10.), GREEN),
298    ] {
299        chart.draw_series(LineSeries::new([(0., 0., 0.), p], color.stroke_width(10)))?;
300    }
301    let geo0 = efd.as_geo();
302    let mut c0 = [0.; 3];
303    for m in efd.coeffs_iter() {
304        let geo = geo0 * GeoVar3::new(c0, na::UnitQuaternion::identity(), 1.);
305        const N: usize = 100;
306        const N_F: f64 = N as f64;
307        let ellipse = (0..N)
308            .map(|i| {
309                let t = i as f64 * std::f64::consts::TAU / N_F;
310                m * na::matrix![t.cos(); t.sin()]
311            })
312            .map(|c| geo.transform_pt(c.data.0[0]).into());
313        let p1 = c0;
314        zip(&mut c0, &m.column(0)).for_each(|(c, u)| *c += u);
315        let p1 = geo0.transform_pt(p1).into();
316        let p2 = geo0.transform_pt(c0).into();
317        chart.draw_series([Circle::new(p2, 5, RED.filled())])?;
318        chart.draw_series(LineSeries::new([p1, p2], RED.stroke_width(5)))?;
319        chart.draw_series(LineSeries::new(ellipse, RED.stroke_width(7)))?;
320    }
321    chart.draw_series(LineSeries::new(
322        curve.into_iter().map(|c| c.into()),
323        BLACK.stroke_width(10),
324    ))?;
325    Ok(())
326}
327
328pub const CURVE2D: &[[f64; 2]] = &[
329    [14.928108089437242, 90.01002059789568],
330    [-3.25371009238094, 85.46456605244113],
331    [-16.763462931659024, 76.52439024390245],
332    [-39.6173464560173, 57.055475143350215],
333    [-49.46583130450215, 35.085778173653246],
334    [-27.739072687756586, 14.939024390243903],
335    [-2.117346456017304, 19.17668726456234],
336    [17.958411119740273, 37.7372933251684],
337    [26.291744453073605, 57.81305090092597],
338    [43.71598687731603, 68.41911150698658],
339    [47.12507778640693, 80.5403236281987],
340    [38.41295657428572, 90.38880847668355],
341    [27.80689596822512, 91.1463842342593],
342];
343pub const CURVE2D_OPEN: &[[f64; 2]] = &[
344    [0.028607755880487345, 47.07692307692308],
345    [6.182453909726641, 52.76923076923077],
346    [14.797838525111256, 57.07692307692308],
347    [24.643992371265103, 58.61538461538461],
348    [41.10553083280357, 59.07692307692308],
349    [50.18245390972664, 56.76923076923077],
350    [60.6439923712651, 51.53846153846154],
351    [65.41322314049587, 46.0],
352    [68.79783852511126, 36.92307692307692],
353    [67.41322314049587, 25.384615384615383],
354    [60.6439923712651, 18.153846153846153],
355];
356#[rustfmt::skip]
357pub const CURVE3D: &[[f64; 3]] = &[
358    [0.6999100262927096,0.43028119112952473,0.5700737247541725],
359    [0.6296590344074013,0.4512760872295425,0.6323601770225047],
360    [0.5638138739696974,0.46183051089148464,0.6847090584540213],
361    [0.5058369206546075,0.4615280764519245,0.7287803814245077],
362    [0.45702863707701546,0.4510622388540364,0.7665948614304103],
363    [0.41723197726966244,0.4317766784472379,0.799678921250722],
364    [0.3856541747672426,0.40523683682222406,0.8288871838597307],
365    [0.3614034251971859,0.3729768425814701,0.8545617819407205],
366    [0.3437381466626122,0.33639334903869056,0.8767460300745514],
367    [0.3321447389761846,0.29671906154191235,0.8953422088163436],
368    [0.3263372629749841,0.25502826620566654,0.9102002934684913],
369    [0.32622858809088523,0.21225076122841705,0.9211539082423659],
370    [0.3318929803861661,0.16918542545203516,0.928021196624841],
371    [0.34352480290022414,0.1265118028240236,0.9305834049340109],
372    [0.36139084538662236,0.08480105307198857,0.9285501808026447],
373    [0.38577056486989364,0.04452862638492877,0.9215195454857344],
374    [0.41687739559371884,0.0060909757379047635,0.9089423177834546],
375    [0.45475516434546037,-0.030172412132808968,0.8901052556002697],
376    [0.49914754370187137,-0.06395289690096129,0.8641537806399759],
377    [0.5493470560271839,-0.09494849409160033,0.8301822664355119],
378    [0.6040449685807827,-0.1228552958162816,0.7874238072487618],
379    [0.6612249903953387,-0.14739768525858138,0.7355647044666407],
380    [0.7181712990997067,-0.1684251622294386,0.6751762361616609],
381    [0.7716940290779204,-0.18608109658863914,0.6081629312757029],
382    [0.8186957920221034,-0.20095221972437277,0.5379176568148215],
383    [0.8570911228971989,-0.21393997152327254,0.46864111603170266],
384    [0.88661962597743,-0.2256744734813517,0.403703691896119],
385    [0.9087058112647335,-0.2359493070484034,0.3443568978212587],
386    [0.9253987983439792,-0.24384980367577586,0.29012813939844717],
387    [0.9384640530551117,-0.24830445184260166,0.24006274246229356],
388    [0.9491216566898215,-0.2484529388271337,0.19349216519158882],
389    [0.9581400770939057,-0.2437369857356437,0.1501994489037393],
390    [0.9659871776934699,-0.2338686278496317,0.11033692690815355],
391    [0.9729414021915203,-0.21877677851516128,0.07430847249221559],
392    [0.979158591144934,-0.1985626427303512,0.04267704651223003],
393    [0.984708076435491,-0.17346680001041803,0.0161019718314504],
394    [0.9895891445905367,-0.1438444290548901,-0.004701610196814504],
395    [0.9937341248474196,-0.11014503640716575,-0.01898841932855333],
396    [0.997000753888979,-0.07289400668205587,-0.02599924103976281],
397    [0.9991540450179144,-0.03267407497300495,-0.024971967263275538],
398    [0.9998362404365794,0.009894685600921771,-0.015152145277328999],
399    [0.9985225360959517,0.05417737861443478,0.004190054292864898],
400    [0.9944609178696164,0.09954962956012486,0.03372468064136638],
401    [0.9865989401537969,0.14542428454053347,0.07398857177482451],
402    [0.9735136458469561,0.19127348741773503,0.12524230260109778],
403    [0.9533914378420906,0.2366305939786505,0.18721839718014333],
404    [0.9241589851285424,0.2810406966408692,0.25873982499375203],
405    [0.8839235581077262,0.32391579444519214,0.33727985639690555],
406    [0.8318395864591331,0.36427712952733915,0.4187422540212782],
407    [0.7691918202384619,0.4004903064827955,0.49794724428553494],
408];
409pub const CURVE2D_POSE: &[[f64; 2]] = &[
410    [18.8, 12.1],
411    [13.3, 18.1],
412    [6.3, 19.8],
413    [-0.4, 17.1],
414    [-2.7, 10.3],
415    [-1.1, 6.0],
416    [0.2, 1.7],
417    [3.4, -2.2],
418    [7.8, -4.9],
419];
420pub const ANGLE2D_POSE: &[f64] = &[-0.9, 0., 0.7, 1.5, 2.8, -2.3, -2., -1.9, -2.1];
421pub const ORIGIN_CURVE: &[[f64; 2]] = &[
422    [0., 0.],
423    [1., 1.],
424    [1., 2.],
425    [0., 3.],
426    [1., 3.],
427    [2., 3.],
428    [3., 3.],
429    [2., 2.],
430    [2., 1.],
431    [3., 0.],
432    [3., -1.],
433    [3., -2.],
434    [2., -2.],
435    [1., -2.],
436    [0., -2.],
437    [0., -1.],
438    [0., 0.],
439];
efd/tests.rs

efd/
tests.rs
