1use geop_core_geometry::{
2 contains::curve::curve_could_contain,
3 intersection::curve_curve_intersect,
4 nurb_curve::{NurbCurve, NurbCurve2D},
5 nurb_surface::NurbSurface3D,
6};
7use geop_core_math::{
8 geop_error::{GeopError, GeopResult},
9 scalars::Scalar,
10 vector::{Vector2, Vector3},
11};
12
13use crate::{CoedgeId, FaceId, Model, boundary::BoundaryType, contains::rng::Rng};
14
15#[derive(Clone, Copy, Debug, PartialEq, Eq)]
17pub enum PointClassification {
18 OnVertex,
20 OnCoedge,
22 Inside,
24 Outside,
27}
28
29const MAX_RAY_ATTEMPTS: usize = 64;
30
31pub fn face_contains<S: Scalar>(
50 model: &Model<S>,
51 face_id: FaceId,
52 u: S,
53 v: S,
54 max_nodes: usize,
55 epsilon: S,
56 seed: u64,
57) -> GeopResult<PointClassification> {
58 let face = &model.faces[&face_id];
59 let coedges: Vec<CoedgeId> = model.iterate_face_coedges(face_id).collect();
60 loops_contain(
61 model,
62 &face.surface,
63 &coedges,
64 Vector2::from_array([u, v]),
65 max_nodes,
66 epsilon,
67 seed,
68 )
69}
70
71pub fn loops_contain<S: Scalar>(
81 model: &Model<S>,
82 surface: &NurbSurface3D<S>,
83 coedges: &[CoedgeId],
84 query: Vector2<S>,
85 max_nodes: usize,
86 epsilon: S,
87 seed: u64,
88) -> GeopResult<PointClassification> {
89 for &coedge_id in coedges {
94 let pcurve = &model.coedges[&coedge_id].pcurve;
95 let vertex_pt = pcurve.evaluate(pcurve.domain().0)?;
96 if vertex_pt.could_be_equal(&query) {
97 return Ok(PointClassification::OnVertex);
98 }
99 }
100 for &coedge_id in coedges {
101 let pcurve = &model.coedges[&coedge_id].pcurve;
102 if curve_could_contain(pcurve, &query, max_nodes, epsilon)?.is_some() {
103 return Ok(PointClassification::OnCoedge);
104 }
105 }
106
107 let (u_lo, u_hi) = surface.domain_u();
108 let (v_lo, v_hi) = surface.domain_v();
109 let du = u_hi.sub(u_lo);
110 let dv = v_hi.sub(v_lo);
111 let diag = du.mul(du).add(dv.mul(dv)).sqrt()?;
112 let ray_length = diag.mul(S::from_f64(3.0)).add(S::ONE);
113 let t_epsilon = epsilon.div(ray_length)?;
114
115 let mut rng = Rng::new(seed);
116 let mut last_rejection = String::new();
119 'attempt: for _ in 0..MAX_RAY_ATTEMPTS {
120 let dir = rng.next_direction2::<S>();
121 let far = query.add(&dir.prod_scalar(ray_length));
122 let ray: NurbCurve2D<S> = NurbCurve::try_new(
123 1,
124 vec![
125 Vector3::from_array([query[0], query[1], S::ONE]),
126 Vector3::from_array([far[0], far[1], S::ONE]),
127 ],
128 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
129 )?;
130
131 for &coedge_id in coedges {
132 let pcurve = &model.coedges[&coedge_id].pcurve;
133 let vertex_pt = pcurve.evaluate(pcurve.domain().0)?;
134 if curve_could_contain(&ray, &vertex_pt, max_nodes, epsilon)?.is_some() {
135 last_rejection = format!(
136 "ray {ray:?} could pass through coedge {coedge_id}'s start {vertex_pt:?}"
137 );
138 continue 'attempt;
139 }
140 }
141
142 let mut count = 0usize;
143 for &coedge_id in coedges {
144 let pcurve = &model.coedges[&coedge_id].pcurve;
145 let (d0, d1) = pcurve.domain();
146 let hits = match curve_curve_intersect(&ray, pcurve, max_nodes, max_nodes, epsilon) {
154 Ok(hits) => hits.into_vec(),
155 Err(e) => {
156 last_rejection =
157 format!("ray {ray:?} x coedge {coedge_id} pcurve {pcurve:?}: {e:?}");
158 continue 'attempt;
159 }
160 };
161 for (t, mid) in hits {
162 if !t.definitely_greater(t_epsilon) {
163 continue;
164 }
165 let mid = mid.midpoint();
170 if !mid.sub(d0).abs().definitely_greater(epsilon)
171 || !mid.sub(d1).abs().definitely_greater(epsilon)
172 {
173 last_rejection = format!(
177 "ray {ray:?} grazes coedge {coedge_id}'s end at t={mid:?} (domain {d0:?}..{d1:?})"
178 );
179 continue 'attempt;
180 }
181 count += 1;
182 }
183 }
184 return Ok(if count % 2 == 1 {
185 PointClassification::Inside
186 } else {
187 PointClassification::Outside
188 });
189 }
190 Err(GeopError::new(format!(
191 "loops_contain: could not find a ray direction clear of every vertex after many attempts; \
192 the last one was rejected because {last_rejection}"
193 )))
194}
195
196const POINTS_PER_BASE: usize = 2;
200
201const MAX_HALVINGS: usize = 40;
205
206pub fn face_interior_point<S: Scalar>(
222 model: &Model<S>,
223 face_id: FaceId,
224 max_nodes: usize,
225 epsilon: S,
226 seed: u64,
227) -> GeopResult<(S, S)> {
228 let found =
229 face_interior_point_where(model, face_id, max_nodes, epsilon, seed, |_, _| Ok(true))?;
230 Ok(found.expect("the first interior point found is always accepted"))
231}
232
233pub fn face_interior_point_where<S: Scalar>(
239 model: &Model<S>,
240 face_id: FaceId,
241 max_nodes: usize,
242 epsilon: S,
243 seed: u64,
244 mut accept: impl FnMut(S, S) -> GeopResult<bool>,
245) -> GeopResult<Option<(S, S)>> {
246 let face = &model.faces[&face_id];
247 let BoundaryType::Loop(anchor) = face.outer else {
248 return Err(GeopError::new(format!(
249 "face_interior_point: face {face_id} is bounded by a bare vertex, so it has no interior to sample"
250 )));
251 };
252
253 let (u_lo, u_hi) = face.surface.domain_u();
254 let (v_lo, v_hi) = face.surface.domain_v();
255 let du = u_hi.sub(u_lo);
256 let dv = v_hi.sub(v_lo);
257 let diagonal = if dv.definitely_greater(du) { dv } else { du };
258
259 let coedges: Vec<CoedgeId> = model
267 .iterate_loop_coedges(anchor)
268 .take(model.coedges.len() + 1)
269 .collect();
270
271 let mut found_any = false;
272 'base: for &coedge_id in &coedges {
273 let pcurve = &model.get_coedge(coedge_id)?.pcurve;
274 let (t0, t1) = pcurve.domain();
275 let t = t0.add(t1).div(S::TWO)?.sharpen();
276 let Ok(base) = pcurve.evaluate(t) else {
277 continue;
278 };
279 let Ok(tangent) = pcurve.tangent(t).and_then(|d| d.normalize()) else {
280 continue;
281 };
282
283 let normals = [
288 Vector2::from_array([tangent[1].neg(), tangent[0]]),
289 Vector2::from_array([tangent[1], tangent[0].neg()]),
290 ];
291
292 let mut step = diagonal.div(S::TWO)?;
293 let mut offered = 0;
294 for _ in 0..MAX_HALVINGS {
295 for inward in normals {
296 let u = base[0].add(inward[0].mul(step)).sharpen();
301 let v = base[1].add(inward[1].mul(step)).sharpen();
302 if u.definitely_less(u_lo)
303 || u.definitely_greater(u_hi)
304 || v.definitely_less(v_lo)
305 || v.definitely_greater(v_hi)
306 {
307 continue;
308 }
309 let neighbourhood = |t: S| t.sub(epsilon).union(t.add(epsilon));
317 if matches!(
318 face_contains(
319 model,
320 face_id,
321 neighbourhood(u),
322 neighbourhood(v),
323 max_nodes,
324 epsilon,
325 seed
326 )?,
327 PointClassification::Inside
328 ) {
329 found_any = true;
330 if accept(u, v)? {
331 return Ok(Some((u, v)));
332 }
333 offered += 1;
340 if offered == POINTS_PER_BASE {
341 continue 'base;
342 }
343 break;
345 }
346 }
347 if !step.definitely_greater(epsilon) {
353 break;
354 }
355 step = step.div(S::TWO)?;
356 }
357 }
358
359 if found_any {
360 return Ok(None);
361 }
362
363 let mut u_extent = None;
370 let mut v_extent = None;
371 for &coedge_id in &coedges {
372 let Ok(coedge) = model.get_coedge(coedge_id) else {
373 continue;
374 };
375 let (t0, t1) = coedge.pcurve.domain();
376 for i in 0..=4 {
377 let Ok(frac) = S::from_ratio(i, 4) else {
378 continue;
379 };
380 let Ok(uv) = coedge.pcurve.evaluate(t0.add(t1.sub(t0).mul(frac))) else {
381 continue;
382 };
383 u_extent = Some(match u_extent {
384 None => uv[0],
385 Some(e) => S::union(e, uv[0]),
386 });
387 v_extent = Some(match v_extent {
388 None => uv[1],
389 Some(e) => S::union(e, uv[1]),
390 });
391 }
392 }
393 let loop_description: Vec<String> = coedges
398 .iter()
399 .map(|&coedge_id| {
400 let start = model.coedge_start_vertex(coedge_id).map(|v| v.point);
401 let coedge = model.get_coedge(coedge_id);
402 let geometry = coedge.as_ref().ok().map(|c| c.geometry);
403 let start_uv = coedge.and_then(|c| c.pcurve.evaluate(c.pcurve.domain().0));
404 format!("{coedge_id} ({geometry:?}): starts at {start:?}, (u, v) = {start_uv:?}")
405 })
406 .collect();
407 Err(GeopError::new(format!(
408 "face_interior_point: no point strictly inside face {face_id} was found, stepping inward from the midpoint of each of its {} outer coedges; that loop spans u={u_extent:?}, v={v_extent:?} (a loop spanning nothing encloses no area, so the face is degenerate) within the surface's domain u={:?}, v={:?}; the loop: [{}]",
409 coedges.len(),
410 (u_lo, u_hi),
411 (v_lo, v_hi),
412 loop_description.join("; ")
413 )))
414}
415
416#[cfg(test)]
417mod interior_point_tests {
418 use super::face_interior_point;
419 use crate::{Model, test_fixtures::test_cube_solid};
420 use geop_core_math::scalars::{ScalInF64, Scalar};
421
422 const MAX: usize = 20000;
423 const SEED: u64 = 99;
424
425 fn eps() -> ScalInF64 {
426 <ScalInF64 as Scalar>::from_f64(1e-4)
427 }
428
429 #[test]
431 fn cube_faces_all_have_interior_points() {
432 let mut model = Model::<ScalInF64>::new();
433 let solid = test_cube_solid(&mut model);
434 for face_id in model.solid_faces(solid).unwrap() {
435 face_interior_point(&model, face_id, MAX, eps(), SEED)
436 .unwrap_or_else(|e| panic!("face {face_id}: {e}"));
437 }
438 }
439
440 }
449
450#[cfg(test)]
451mod tests {
452 use super::{PointClassification, face_contains};
453 use crate::{
454 Coedge, CoedgeGeometry, CoedgeId, Edge, Face, FaceId, Model, Sense, ShellId, Vertex,
455 VertexId, boundary::BoundaryType, model::Curve3,
456 };
457 use geop_core_geometry::{
458 nurb_curve::{NurbCurve, NurbCurve2D},
459 nurb_surface::NurbSurface3D,
460 };
461 use geop_core_math::{
462 for_all_scalars,
463 scalars::Scalar,
464 vector::{Vector3, Vector4},
465 };
466
467 const MAX: usize = 200;
468 const EPS: f64 = 1e-3;
469 const SEED: u64 = 12345;
470
471 fn p2<S: Scalar>(x: f64, y: f64) -> Vector3<S> {
472 Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ONE])
473 }
474
475 fn line2<S: Scalar>(a: (f64, f64), b: (f64, f64)) -> NurbCurve2D<S> {
476 NurbCurve::try_new(
477 1,
478 vec![p2(a.0, a.1), p2(b.0, b.1)],
479 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
480 )
481 .unwrap()
482 }
483
484 fn polygon_face<S: Scalar>(model: &mut Model<S>, points: &[(f64, f64)]) -> FaceId {
487 let p =
488 |x: f64, y: f64| Vector4::from_array([S::from_f64(x), S::from_f64(y), S::ZERO, S::ONE]);
489 let surface = NurbSurface3D::try_new(
490 1,
491 1,
492 vec![p(0.0, 0.0), p(0.0, 1.0), p(1.0, 0.0), p(1.0, 1.0)],
493 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
494 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
495 )
496 .unwrap();
497
498 let face_id = model.insert_face(Face {
499 surface,
500 outer: BoundaryType::Vertex(VertexId(0)),
501 holes: Vec::new(),
502 shell: ShellId(999),
503 });
504
505 let n = points.len();
506 let verts: Vec<VertexId> = points
507 .iter()
508 .map(|&(x, y)| {
509 model.insert_vertex(Vertex {
510 point: Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ZERO]),
511 })
512 })
513 .collect();
514 let edges = (0..n)
515 .map(|i| {
516 model.insert_edge(Edge {
517 curve: Curve3::try_new(
518 1,
519 vec![
520 Vector4::from_array([
521 S::from_f64(points[i].0),
522 S::from_f64(points[i].1),
523 S::ZERO,
524 S::ONE,
525 ]),
526 Vector4::from_array([
527 S::from_f64(points[(i + 1) % n].0),
528 S::from_f64(points[(i + 1) % n].1),
529 S::ZERO,
530 S::ONE,
531 ]),
532 ],
533 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
534 )
535 .unwrap(),
536 start_vertex: verts[i],
537 end_vertex: verts[(i + 1) % n],
538 })
539 })
540 .collect::<Vec<_>>();
541 let coedges: Vec<CoedgeId> = (0..n)
542 .map(|i| {
543 model.insert_coedge(Coedge {
544 geometry: CoedgeGeometry::Edge(edges[i]),
545 sense: Sense::Forward,
546 pcurve: line2(points[i], points[(i + 1) % n]),
547 next: CoedgeId(0),
548 prev: CoedgeId(0),
549 face: face_id,
550 })
551 })
552 .collect();
553 for i in 0..n {
554 model.coedges.get_mut(&coedges[i]).unwrap().next = coedges[(i + 1) % n];
555 model.coedges.get_mut(&coedges[i]).unwrap().prev = coedges[(i + n - 1) % n];
556 }
557 model.faces.get_mut(&face_id).unwrap().outer = BoundaryType::Loop(coedges[0]);
558
559 face_id
560 }
561
562 fn diamond_face<S: Scalar>(model: &mut Model<S>) -> FaceId {
565 polygon_face(model, &[(1., 0.5), (0.5, 1.), (0., 0.5), (0.5, 0.)])
566 }
567
568 fn check_diamond_interior_point_is_contained<S: Scalar>() {
569 let mut model = Model::<S>::new();
570 let face_id = diamond_face(&mut model);
571 assert_eq!(
572 face_contains(
573 &model,
574 face_id,
575 S::from_f64(0.5),
576 S::from_f64(0.3),
577 MAX,
578 S::from_f64(EPS),
579 SEED
580 )
581 .unwrap(),
582 PointClassification::Inside
583 );
584 }
585 #[test]
586 fn diamond_interior_point_is_contained() {
587 for_all_scalars!(check_diamond_interior_point_is_contained);
588 }
589
590 fn check_diamond_exterior_point_is_not_contained<S: Scalar>() {
591 let mut model = Model::<S>::new();
592 let face_id = diamond_face(&mut model);
593 assert_eq!(
594 face_contains(
595 &model,
596 face_id,
597 S::from_f64(0.1),
598 S::from_f64(0.3),
599 MAX,
600 S::from_f64(EPS),
601 SEED
602 )
603 .unwrap(),
604 PointClassification::Outside
605 );
606 }
607 #[test]
608 fn diamond_exterior_point_is_not_contained() {
609 for_all_scalars!(check_diamond_exterior_point_is_not_contained);
610 }
611
612 fn check_diamond_center_hits_convex_vertex_from_inside<S: Scalar>() {
613 let mut model = Model::<S>::new();
614 let face_id = diamond_face(&mut model);
615 assert_eq!(
616 face_contains(
617 &model,
618 face_id,
619 S::from_f64(0.5),
620 S::from_f64(0.5),
621 MAX,
622 S::from_f64(EPS),
623 SEED
624 )
625 .unwrap(),
626 PointClassification::Inside
627 );
628 }
629 #[test]
630 fn diamond_center_hits_convex_vertex_from_inside() {
631 for_all_scalars!(check_diamond_center_hits_convex_vertex_from_inside);
632 }
633
634 fn check_diamond_vertex_query_is_on_vertex<S: Scalar>() {
635 let mut model = Model::<S>::new();
636 let face_id = diamond_face(&mut model);
637 assert_eq!(
638 face_contains(
639 &model,
640 face_id,
641 S::ONE,
642 S::from_f64(0.5),
643 MAX,
644 S::from_f64(EPS),
645 SEED
646 )
647 .unwrap(),
648 PointClassification::OnVertex
649 );
650 }
651 #[test]
652 fn diamond_vertex_query_is_on_vertex() {
653 for_all_scalars!(check_diamond_vertex_query_is_on_vertex);
654 }
655
656 fn check_diamond_edge_query_is_on_coedge<S: Scalar>() {
657 let mut model = Model::<S>::new();
658 let face_id = diamond_face(&mut model);
659 assert_eq!(
660 face_contains(
661 &model,
662 face_id,
663 S::from_f64(0.75),
664 S::from_f64(0.75),
665 MAX,
666 S::from_f64(EPS),
667 SEED
668 )
669 .unwrap(),
670 PointClassification::OnCoedge
671 );
672 }
673 #[test]
674 fn diamond_edge_query_is_on_coedge() {
675 for_all_scalars!(check_diamond_edge_query_is_on_coedge);
676 }
677}