1pub fn uniform_knot_vector(n: usize, degree: usize) -> Vec<f64> {
7 let m = n + degree + 1;
8 (0..m).map(|i| i as f64 / (m - 1) as f64).collect()
9}
10pub fn clamped_knot_vector(n: usize, degree: usize) -> Vec<f64> {
12 let m = n + degree + 1;
13 let mut knots = vec![0.0f64; m];
14 for (i, knot) in knots.iter_mut().enumerate().take(m) {
15 if i <= degree {
16 *knot = 0.0;
17 } else if i >= n {
18 *knot = (n - degree) as f64;
19 } else {
20 *knot = (i - degree) as f64;
21 }
22 }
23 let max_k = knots[m - 1];
24 if max_k > 0.0 {
25 for k in &mut knots {
26 *k /= max_k;
27 }
28 }
29 knots
30}
31#[cfg(test)]
32mod tests {
33 use super::*;
34 use crate::fractal_geometry::IteratedFunctionSystem;
35 use crate::fractal_geometry::LSystem;
36 use crate::procedural_geometry::CurveParametric;
37 use crate::procedural_geometry::DelaunayTri;
38 use crate::procedural_geometry::DelaunayTriangulation2d;
39 use crate::procedural_geometry::FractalGeometry;
40 use crate::procedural_geometry::NoiseGeometry;
41 use crate::procedural_geometry::PerlinNoise2d;
42 use crate::procedural_geometry::PerlinNoise3d;
43 use crate::procedural_geometry::ReactionDiffusion;
44 use crate::procedural_geometry::SpaceColonization;
45 use crate::procedural_geometry::SurfaceParametric;
46 use crate::procedural_geometry::TurtleState;
47 use crate::procedural_geometry::VoronoiTessellation2d;
48 use crate::procedural_geometry::WorleyNoise;
49 pub(super) const EPS: f64 = 1e-6;
50 #[test]
51 fn test_koch_fractal_dimension_approx() {
52 let koch = FractalGeometry::koch_snowflake(5, 4);
53 let dim = koch.box_counting_dimension(10);
54 assert!(
55 dim > 1.1 && dim < 1.5,
56 "Koch dim = {:.6}, expected ~1.26",
57 dim
58 );
59 }
60 #[test]
61 fn test_sierpinski_fractal_dimension_approx() {
62 let sier = FractalGeometry::sierpinski_triangle(7);
63 let dim = sier.hausdorff_dimension();
64 assert!(
65 dim > 0.5 && dim < 2.0,
66 "Sierpinski dim = {:.6} out of expected range",
67 dim
68 );
69 }
70 #[test]
71 fn test_fractal_dim_line_is_one() {
72 let pts: Vec<[f64; 2]> = (0..500).map(|i| [i as f64 / 499.0, 0.5]).collect();
73 let fg = FractalGeometry::new(pts);
74 let dim = fg.box_counting_dimension(8);
75 assert!(dim > 0.8 && dim < 1.2, "Line dim = {:.6}", dim);
76 }
77 #[test]
78 fn test_fractal_dim_plane_fill_near_two() {
79 let pts: Vec<[f64; 2]> = (0..50)
80 .flat_map(|i| (0..50).map(move |j| [i as f64 / 49.0, j as f64 / 49.0]))
81 .collect();
82 let fg = FractalGeometry::new(pts);
83 let dim = fg.box_counting_dimension(8);
84 assert!(dim > 1.1, "Dense grid dim = {:.6} should be > 1.1", dim);
85 }
86 #[test]
87 fn test_koch_snowflake_point_count() {
88 let k = FractalGeometry::koch_snowflake(3, 2);
89 let expected = 3 * 4usize.pow(3) * 2;
90 assert_eq!(k.points.len(), expected);
91 }
92 #[test]
93 fn test_sierpinski_triangle_non_empty() {
94 let s = FractalGeometry::sierpinski_triangle(4);
95 assert!(!s.points.is_empty());
96 assert_eq!(s.points.len(), 3usize.pow(4) * 3);
97 }
98 #[test]
99 fn test_lsystem_string_grows_geometrically() {
100 let ls = LSystem::koch_curve();
101 let s1 = ls.evolve(1).string.len();
102 let s2 = ls.evolve(2).string.len();
103 let s3 = ls.evolve(3).string.len();
104 assert!(s2 > s1, "s2={} should be > s1={}", s2, s1);
105 assert!(s3 > s2, "s3={} should be > s2={}", s3, s2);
106 }
107 #[test]
108 fn test_lsystem_axiom_preserved_gen0() {
109 let ls = LSystem::sierpinski();
110 let state = ls.evolve(0);
111 assert_eq!(state.string, ls.axiom);
112 assert_eq!(state.generation, 0);
113 }
114 #[test]
115 fn test_lsystem_dragon_curve_generates_segments() {
116 let ls = LSystem::dragon_curve();
117 let state = ls.evolve(6);
118 let segs = ls.interpret(&state);
119 assert!(!segs.is_empty(), "dragon curve should have segments");
120 }
121 #[test]
122 fn test_lsystem_plant_branching_has_brackets() {
123 let ls = LSystem::plant_branching();
124 let state = ls.evolve(3);
125 let open = state.string.chars().filter(|&c| c == '[').count();
126 let close = state.string.chars().filter(|&c| c == ']').count();
127 assert_eq!(
128 open, close,
129 "unbalanced brackets: open={}, close={}",
130 open, close
131 );
132 }
133 #[test]
134 fn test_lsystem_generation_count() {
135 let ls = LSystem::koch_curve();
136 let state = ls.evolve(5);
137 assert_eq!(state.generation, 5);
138 }
139 #[test]
140 fn test_lsystem_turtle_state_default() {
141 let t = TurtleState::default();
142 assert_eq!(t.x, 0.0);
143 assert_eq!(t.y, 0.0);
144 assert_eq!(t.angle, 0.0);
145 }
146 #[test]
147 fn test_perlin2d_range() {
148 let noise = PerlinNoise2d::new(42);
149 for i in 0..100 {
150 let x = i as f64 * 0.137;
151 let y = i as f64 * 0.251;
152 let v = noise.noise(x, y);
153 assert!((-1.5..=1.5).contains(&v), "Perlin2D out of range: {:.6}", v);
154 }
155 }
156 #[test]
157 fn test_perlin3d_range() {
158 let noise = PerlinNoise3d::new(123);
159 for i in 0..100 {
160 let x = i as f64 * 0.1;
161 let y = i as f64 * 0.17;
162 let z = i as f64 * 0.23;
163 let v = noise.noise(x, y, z);
164 assert!((-1.5..=1.5).contains(&v), "Perlin3D out of range: {:.6}", v);
165 }
166 }
167 #[test]
168 fn test_perlin2d_deterministic() {
169 let n1 = PerlinNoise2d::new(99);
170 let n2 = PerlinNoise2d::new(99);
171 for i in 0..10 {
172 let x = i as f64 * 0.3;
173 let y = i as f64 * 0.4;
174 assert!((n1.noise(x, y) - n2.noise(x, y)).abs() < EPS);
175 }
176 }
177 #[test]
178 fn test_fbm_different_seeds_differ() {
179 let n1 = PerlinNoise2d::new(1);
180 let n2 = PerlinNoise2d::new(999);
181 let v1 = n1.fbm(1.5, 2.3, 4, 2.0, 0.5);
182 let v2 = n2.fbm(1.5, 2.3, 4, 2.0, 0.5);
183 assert!(
184 (v1 - v2).abs() > 1e-8,
185 "Different seeds should produce different noise"
186 );
187 }
188 #[test]
189 fn test_fbm_nonnegative_max_val() {
190 let noise = PerlinNoise2d::new(7);
191 let v1 = noise.fbm(1.0, 1.0, 1, 2.0, 0.5);
192 let v4 = noise.fbm(1.0, 1.0, 4, 2.0, 0.5);
193 assert!(v1.is_finite(), "fBm 1 octave should be finite");
194 assert!(v4.is_finite(), "fBm 4 octaves should be finite");
195 }
196 #[test]
197 fn test_worley_f1_nonneg() {
198 let w = WorleyNoise::new(20, 42);
199 for i in 0..50 {
200 let x = i as f64 / 49.0;
201 let y = (50 - i) as f64 / 49.0;
202 assert!(w.f1(x, y) >= 0.0);
203 }
204 }
205 #[test]
206 fn test_worley_f2_minus_f1_nonneg() {
207 let w = WorleyNoise::new(20, 42);
208 for i in 0..50 {
209 let x = i as f64 / 49.0;
210 let y = (50 - i) as f64 / 49.0;
211 assert!(w.f2_minus_f1(x, y) >= -EPS);
212 }
213 }
214 #[test]
215 fn test_heightmap_fbm_dimensions() {
216 let hm = NoiseGeometry::heightmap_fbm(16, 16, 4.0, 4, 77);
217 assert_eq!(hm.len(), 16);
218 assert_eq!(hm[0].len(), 16);
219 }
220 #[test]
221 fn test_space_colonization_grows() {
222 let attractors: Vec<[f64; 3]> = (0..50)
223 .map(|i| [i as f64 * 0.05, 1.0 + i as f64 * 0.03, 0.0])
224 .collect();
225 let mut sc = SpaceColonization::new(attractors, [0.0, 0.0, 0.0], 0.1, 2.0, 0.15);
226 let initial = sc.nodes.len();
227 sc.grow();
228 assert!(sc.nodes.len() >= initial);
229 }
230 #[test]
231 fn test_space_colonization_root_has_no_parent() {
232 let attractors = vec![[0.0, 1.0, 0.0]];
233 let sc = SpaceColonization::new(attractors, [0.0, 0.0, 0.0], 0.1, 5.0, 0.15);
234 assert!(sc.nodes[0].parent.is_none());
235 }
236 #[test]
237 fn test_space_colonization_grow_until_done() {
238 let attractors: Vec<[f64; 3]> = (0..10).map(|i| [i as f64 * 0.1, 0.5, 0.0]).collect();
239 let mut sc = SpaceColonization::new(attractors, [0.5, 0.0, 0.0], 0.15, 1.0, 0.2);
240 sc.grow_until_done(100);
241 assert!(sc.attractors.is_empty() || sc.nodes.len() > 1);
242 }
243 #[test]
244 fn test_reaction_diffusion_step_runs() {
245 let mut rd = ReactionDiffusion::new(20, 20, 0.2, 0.1, 0.055, 0.062);
246 rd.seed_center();
247 rd.step(5, 1.0);
248 for &u in &rd.u {
249 assert!((0.0..=1.0).contains(&u), "U out of range: {:.6}", u);
250 }
251 }
252 #[test]
253 fn test_reaction_diffusion_v_in_range() {
254 let mut rd = ReactionDiffusion::new(20, 20, 0.2, 0.1, 0.055, 0.062);
255 rd.seed_center();
256 rd.step(10, 1.0);
257 for &v in &rd.v {
258 assert!((0.0..=1.0).contains(&v), "V out of range: {:.6}", v);
259 }
260 }
261 #[test]
262 fn test_reaction_diffusion_grid_size() {
263 let rd = ReactionDiffusion::new(10, 15, 0.2, 0.1, 0.055, 0.062);
264 assert_eq!(rd.u.len(), 10 * 15);
265 assert_eq!(rd.v.len(), 10 * 15);
266 }
267 #[test]
268 fn test_ifs_barnsley_fern_generates_points() {
269 let ifs = IteratedFunctionSystem::barnsley_fern();
270 let pts = ifs.generate(1000, 42);
271 assert_eq!(pts.len(), 1000);
272 }
273 #[test]
274 fn test_ifs_barnsley_fern_y_range() {
275 let ifs = IteratedFunctionSystem::barnsley_fern();
276 let pts = ifs.generate(500, 1);
277 for &[_, y] in &pts {
278 assert!(
279 (-0.5..=11.0).contains(&y),
280 "y={:.6} out of expected range",
281 y
282 );
283 }
284 }
285 #[test]
286 fn test_ifs_empty_transforms() {
287 let ifs = IteratedFunctionSystem::new(vec![], vec![]);
288 let pts = ifs.generate(10, 42);
289 assert!(pts.is_empty());
290 }
291 #[test]
292 fn test_voronoi_cell_areas_sum() {
293 let voronoi = VoronoiTessellation2d::random_seeds(10, [0.0, 1.0, 0.0, 1.0], 42);
294 let areas = voronoi.cell_areas(50);
295 let total: f64 = areas.iter().sum();
296 assert!((total - 1.0).abs() < 0.05, "Total area = {:.6}", total);
297 }
298 #[test]
299 fn test_voronoi_nearest_seed_is_correct() {
300 let seeds = vec![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]];
301 let voronoi = VoronoiTessellation2d::new(seeds, [-0.5, 1.5, -0.5, 1.5]);
302 assert_eq!(voronoi.nearest_seed(0.9, 0.1), 1);
303 }
304 #[test]
305 fn test_voronoi_no_negative_areas() {
306 let voronoi = VoronoiTessellation2d::random_seeds(8, [0.0, 1.0, 0.0, 1.0], 77);
307 let areas = voronoi.cell_areas(30);
308 for &a in &areas {
309 assert!(a >= 0.0);
310 }
311 }
312 #[test]
313 fn test_voronoi_neighbors_symmetric() {
314 let voronoi = VoronoiTessellation2d::random_seeds(6, [0.0, 1.0, 0.0, 1.0], 55);
315 let neighbors = voronoi.neighbors(30);
316 for (i, nbrs) in neighbors.iter().enumerate() {
317 for &j in nbrs {
318 assert!(
319 neighbors[j].contains(&i),
320 "neighbor relation not symmetric: {} neighbors {} but not vice versa",
321 i,
322 j
323 );
324 }
325 }
326 }
327 #[test]
328 fn test_delaunay_circumcircle_property() {
329 let points = vec![[0.0, 0.0], [1.0, 0.0], [0.5, 0.866], [0.5, 0.3], [0.2, 0.5]];
330 let dt = DelaunayTriangulation2d::new(points);
331 assert!(dt.verify_delaunay(), "Delaunay property violated");
332 }
333 #[test]
334 fn test_delaunay_triangle_count() {
335 let points = vec![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
336 let dt = DelaunayTriangulation2d::new(points);
337 assert!(!dt.triangles.is_empty(), "Should have at least 1 triangle");
338 }
339 #[test]
340 fn test_delaunay_triangle_indices_valid() {
341 let points: Vec<[f64; 2]> = (0..8)
342 .map(|i| {
343 let a = i as f64 * std::f64::consts::PI / 4.0;
344 [a.cos(), a.sin()]
345 })
346 .collect();
347 let dt = DelaunayTriangulation2d::new(points.clone());
348 for tri in &dt.triangles {
349 for &idx in &tri.indices {
350 assert!(idx < points.len(), "Triangle index {} out of bounds", idx);
351 }
352 }
353 }
354 #[test]
355 fn test_delaunay_circumcircle_center() {
356 let points = vec![[0.0, 0.0], [1.0, 0.0], [0.5, 0.866_025_4]];
357 let tri = DelaunayTri::new(0, 1, 2);
358 let ([cx, cy], _r2) = tri.circumcircle(&points);
359 assert!((cx - 0.5).abs() < 0.01, "circumcenter x = {:.6}", cx);
360 assert!(cy > 0.0 && cy < 1.0, "circumcenter y = {:.6}", cy);
361 }
362 #[test]
363 fn test_bezier_endpoint_interpolation() {
364 let pts = vec![
365 [0.0, 0.0, 0.0],
366 [1.0, 2.0, 0.0],
367 [3.0, 1.0, 0.0],
368 [4.0, 0.0, 0.0],
369 ];
370 let curve = CurveParametric::new(pts.clone());
371 let p0 = curve.bezier_de_casteljau(0.0);
372 let p1 = curve.bezier_de_casteljau(1.0);
373 assert!((p0[0] - pts[0][0]).abs() < EPS, "start x mismatch");
374 assert!((p1[0] - pts[3][0]).abs() < EPS, "end x mismatch");
375 }
376 #[test]
377 fn test_bezier_midpoint_on_curve() {
378 let pts = vec![
379 [0.0, 0.0, 0.0],
380 [0.0, 1.0, 0.0],
381 [1.0, 1.0, 0.0],
382 [1.0, 0.0, 0.0],
383 ];
384 let curve = CurveParametric::new(pts);
385 let mid = curve.bezier_de_casteljau(0.5);
386 assert!(mid[0] >= -0.01 && mid[0] <= 1.01);
387 assert!(mid[1] >= -0.01 && mid[1] <= 1.01);
388 }
389 #[test]
390 fn test_bezier_bernstein_matches_de_casteljau() {
391 let pts = vec![
392 [0.0, 0.0, 0.0],
393 [1.0, 2.0, 0.0],
394 [2.0, 2.0, 0.0],
395 [3.0, 0.0, 0.0],
396 ];
397 let curve = CurveParametric::new(pts);
398 for i in 0..=10 {
399 let t = i as f64 / 10.0;
400 let dc = curve.bezier_de_casteljau(t);
401 let bern = curve.bezier_bernstein(t);
402 assert!((dc[0] - bern[0]).abs() < EPS, "x mismatch at t={}", t);
403 assert!((dc[1] - bern[1]).abs() < EPS, "y mismatch at t={}", t);
404 }
405 }
406 #[test]
407 fn test_bspline_endpoint_at_clamped_knots() {
408 let pts = vec![
409 [0.0, 0.0, 0.0],
410 [1.0, 1.0, 0.0],
411 [2.0, 0.0, 0.0],
412 [3.0, 1.0, 0.0],
413 ];
414 let n = pts.len();
415 let degree = 3;
416 let knots = clamped_knot_vector(n, degree);
417 let curve = CurveParametric::new(pts.clone());
418 let start = curve.bspline_de_boor(0.0, degree, &knots);
419 let end = curve.bspline_de_boor(1.0, degree, &knots);
420 assert!((start[0] - pts[0][0]).abs() < EPS, "start x mismatch");
421 assert!((end[0] - pts[n - 1][0]).abs() < EPS, "end x mismatch");
422 }
423 #[test]
424 fn test_bezier_arc_length_positive() {
425 let pts = vec![[0.0, 0.0, 0.0], [1.0, 1.0, 0.0], [2.0, 0.0, 0.0]];
426 let curve = CurveParametric::new(pts);
427 let len = curve.bezier_arc_length(50);
428 assert!(len > 0.0, "arc length should be positive");
429 }
430 #[test]
431 fn test_bezier_sample_count() {
432 let pts = vec![[0.0, 0.0, 0.0], [1.0, 1.0, 0.0], [2.0, 0.0, 0.0]];
433 let curve = CurveParametric::new(pts);
434 let samples = curve.sample_bezier(20);
435 assert_eq!(samples.len(), 20);
436 }
437 #[test]
438 fn test_tensor_bezier_corner_interpolation() {
439 let net = vec![
440 vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]],
441 vec![[0.0, 1.0, 0.0], [1.0, 1.0, 1.0], [2.0, 1.0, 0.0]],
442 vec![[0.0, 2.0, 0.0], [1.0, 2.0, 0.0], [2.0, 2.0, 0.0]],
443 ];
444 let surf = SurfaceParametric::new(net);
445 let p00 = surf.tensor_bezier(0.0, 0.0);
446 let p11 = surf.tensor_bezier(1.0, 1.0);
447 assert!((p00[0] - 0.0).abs() < EPS);
448 assert!((p11[0] - 2.0).abs() < EPS);
449 }
450 #[test]
451 fn test_coons_patch_corners() {
452 let net = vec![
453 vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]],
454 vec![[0.0, 1.0, 0.0], [1.0, 1.0, 0.0]],
455 ];
456 let surf = SurfaceParametric::new(net);
457 let p00 = surf.coons_patch(0.0, 0.0);
458 let p11 = surf.coons_patch(1.0, 1.0);
459 assert!((p00[0] - 0.0).abs() < EPS && (p00[1] - 0.0).abs() < EPS);
460 assert!((p11[0] - 1.0).abs() < EPS && (p11[1] - 1.0).abs() < EPS);
461 }
462 #[test]
463 fn test_surface_sample_dimensions() {
464 let net = vec![
465 vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]],
466 vec![[0.0, 1.0, 0.0], [1.0, 1.0, 0.0]],
467 ];
468 let surf = SurfaceParametric::new(net);
469 let samples = surf.sample(5, 7);
470 assert_eq!(samples.len(), 5);
471 assert_eq!(samples[0].len(), 7);
472 }
473 #[test]
474 fn test_surface_normal_unit_length() {
475 let net = vec![
476 vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]],
477 vec![[0.0, 1.0, 0.0], [1.0, 1.0, 0.0], [2.0, 1.0, 0.0]],
478 vec![[0.0, 2.0, 0.0], [1.0, 2.0, 0.0], [2.0, 2.0, 0.0]],
479 ];
480 let surf = SurfaceParametric::new(net);
481 let n = surf.normal(0.5, 0.5);
482 let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
483 assert!((len - 1.0).abs() < 1e-4, "normal length = {:.6}", len);
484 }
485 #[test]
486 fn test_uniform_knot_vector_length() {
487 let n = 5;
488 let degree = 3;
489 let knots = uniform_knot_vector(n, degree);
490 assert_eq!(knots.len(), n + degree + 1);
491 }
492 #[test]
493 fn test_clamped_knot_vector_endpoints() {
494 let n = 6;
495 let degree = 3;
496 let knots = clamped_knot_vector(n, degree);
497 assert_eq!(knots.len(), n + degree + 1);
498 assert!((knots[0] - 0.0).abs() < EPS);
499 assert!((knots[knots.len() - 1] - 1.0).abs() < EPS);
500 }
501 #[test]
502 fn test_clamped_knot_vector_monotone() {
503 let knots = clamped_knot_vector(5, 2);
504 for i in 1..knots.len() {
505 assert!(knots[i] >= knots[i - 1], "knot not monotone at {}", i);
506 }
507 }
508}