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 Err(GeopError::new(format!(
394 "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)",
395 coedges.len()
396 )))
397}
398
399#[cfg(test)]
400mod interior_point_tests {
401 use super::face_interior_point;
402 use crate::{Model, test_fixtures::test_cube_solid};
403 use geop_core_math::scalars::{ScalInF64, Scalar};
404
405 const MAX: usize = 20000;
406 const SEED: u64 = 99;
407
408 fn eps() -> ScalInF64 {
409 <ScalInF64 as Scalar>::from_f64(1e-4)
410 }
411
412 #[test]
414 fn cube_faces_all_have_interior_points() {
415 let mut model = Model::<ScalInF64>::new();
416 let solid = test_cube_solid(&mut model);
417 for face_id in model.solid_faces(solid).unwrap() {
418 face_interior_point(&model, face_id, MAX, eps(), SEED)
419 .unwrap_or_else(|e| panic!("face {face_id}: {e}"));
420 }
421 }
422
423 }
432
433#[cfg(test)]
434mod tests {
435 use super::{PointClassification, face_contains};
436 use crate::{
437 Coedge, CoedgeGeometry, CoedgeId, Edge, Face, FaceId, Model, Sense, ShellId, Vertex,
438 VertexId, boundary::BoundaryType, model::Curve3,
439 };
440 use geop_core_geometry::{
441 nurb_curve::{NurbCurve, NurbCurve2D},
442 nurb_surface::NurbSurface3D,
443 };
444 use geop_core_math::{
445 for_all_scalars,
446 scalars::Scalar,
447 vector::{Vector3, Vector4},
448 };
449
450 const MAX: usize = 200;
451 const EPS: f64 = 1e-3;
452 const SEED: u64 = 12345;
453
454 fn p2<S: Scalar>(x: f64, y: f64) -> Vector3<S> {
455 Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ONE])
456 }
457
458 fn line2<S: Scalar>(a: (f64, f64), b: (f64, f64)) -> NurbCurve2D<S> {
459 NurbCurve::try_new(
460 1,
461 vec![p2(a.0, a.1), p2(b.0, b.1)],
462 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
463 )
464 .unwrap()
465 }
466
467 fn polygon_face<S: Scalar>(model: &mut Model<S>, points: &[(f64, f64)]) -> FaceId {
470 let p =
471 |x: f64, y: f64| Vector4::from_array([S::from_f64(x), S::from_f64(y), S::ZERO, S::ONE]);
472 let surface = NurbSurface3D::try_new(
473 1,
474 1,
475 vec![p(0.0, 0.0), p(0.0, 1.0), p(1.0, 0.0), p(1.0, 1.0)],
476 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
477 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
478 )
479 .unwrap();
480
481 let face_id = model.insert_face(Face {
482 surface,
483 outer: BoundaryType::Vertex(VertexId(0)),
484 holes: Vec::new(),
485 shell: ShellId(999),
486 });
487
488 let n = points.len();
489 let verts: Vec<VertexId> = points
490 .iter()
491 .map(|&(x, y)| {
492 model.insert_vertex(Vertex {
493 point: Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ZERO]),
494 })
495 })
496 .collect();
497 let edges = (0..n)
498 .map(|i| {
499 model.insert_edge(Edge {
500 curve: Curve3::try_new(
501 1,
502 vec![
503 Vector4::from_array([
504 S::from_f64(points[i].0),
505 S::from_f64(points[i].1),
506 S::ZERO,
507 S::ONE,
508 ]),
509 Vector4::from_array([
510 S::from_f64(points[(i + 1) % n].0),
511 S::from_f64(points[(i + 1) % n].1),
512 S::ZERO,
513 S::ONE,
514 ]),
515 ],
516 vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
517 )
518 .unwrap(),
519 start_vertex: verts[i],
520 end_vertex: verts[(i + 1) % n],
521 })
522 })
523 .collect::<Vec<_>>();
524 let coedges: Vec<CoedgeId> = (0..n)
525 .map(|i| {
526 model.insert_coedge(Coedge {
527 geometry: CoedgeGeometry::Edge(edges[i]),
528 sense: Sense::Forward,
529 pcurve: line2(points[i], points[(i + 1) % n]),
530 next: CoedgeId(0),
531 prev: CoedgeId(0),
532 face: face_id,
533 })
534 })
535 .collect();
536 for i in 0..n {
537 model.coedges.get_mut(&coedges[i]).unwrap().next = coedges[(i + 1) % n];
538 model.coedges.get_mut(&coedges[i]).unwrap().prev = coedges[(i + n - 1) % n];
539 }
540 model.faces.get_mut(&face_id).unwrap().outer = BoundaryType::Loop(coedges[0]);
541
542 face_id
543 }
544
545 fn diamond_face<S: Scalar>(model: &mut Model<S>) -> FaceId {
548 polygon_face(model, &[(1., 0.5), (0.5, 1.), (0., 0.5), (0.5, 0.)])
549 }
550
551 fn check_diamond_interior_point_is_contained<S: Scalar>() {
552 let mut model = Model::<S>::new();
553 let face_id = diamond_face(&mut model);
554 assert_eq!(
555 face_contains(
556 &model,
557 face_id,
558 S::from_f64(0.5),
559 S::from_f64(0.3),
560 MAX,
561 S::from_f64(EPS),
562 SEED
563 )
564 .unwrap(),
565 PointClassification::Inside
566 );
567 }
568 #[test]
569 fn diamond_interior_point_is_contained() {
570 for_all_scalars!(check_diamond_interior_point_is_contained);
571 }
572
573 fn check_diamond_exterior_point_is_not_contained<S: Scalar>() {
574 let mut model = Model::<S>::new();
575 let face_id = diamond_face(&mut model);
576 assert_eq!(
577 face_contains(
578 &model,
579 face_id,
580 S::from_f64(0.1),
581 S::from_f64(0.3),
582 MAX,
583 S::from_f64(EPS),
584 SEED
585 )
586 .unwrap(),
587 PointClassification::Outside
588 );
589 }
590 #[test]
591 fn diamond_exterior_point_is_not_contained() {
592 for_all_scalars!(check_diamond_exterior_point_is_not_contained);
593 }
594
595 fn check_diamond_center_hits_convex_vertex_from_inside<S: Scalar>() {
596 let mut model = Model::<S>::new();
597 let face_id = diamond_face(&mut model);
598 assert_eq!(
599 face_contains(
600 &model,
601 face_id,
602 S::from_f64(0.5),
603 S::from_f64(0.5),
604 MAX,
605 S::from_f64(EPS),
606 SEED
607 )
608 .unwrap(),
609 PointClassification::Inside
610 );
611 }
612 #[test]
613 fn diamond_center_hits_convex_vertex_from_inside() {
614 for_all_scalars!(check_diamond_center_hits_convex_vertex_from_inside);
615 }
616
617 fn check_diamond_vertex_query_is_on_vertex<S: Scalar>() {
618 let mut model = Model::<S>::new();
619 let face_id = diamond_face(&mut model);
620 assert_eq!(
621 face_contains(
622 &model,
623 face_id,
624 S::ONE,
625 S::from_f64(0.5),
626 MAX,
627 S::from_f64(EPS),
628 SEED
629 )
630 .unwrap(),
631 PointClassification::OnVertex
632 );
633 }
634 #[test]
635 fn diamond_vertex_query_is_on_vertex() {
636 for_all_scalars!(check_diamond_vertex_query_is_on_vertex);
637 }
638
639 fn check_diamond_edge_query_is_on_coedge<S: Scalar>() {
640 let mut model = Model::<S>::new();
641 let face_id = diamond_face(&mut model);
642 assert_eq!(
643 face_contains(
644 &model,
645 face_id,
646 S::from_f64(0.75),
647 S::from_f64(0.75),
648 MAX,
649 S::from_f64(EPS),
650 SEED
651 )
652 .unwrap(),
653 PointClassification::OnCoedge
654 );
655 }
656 #[test]
657 fn diamond_edge_query_is_on_coedge() {
658 for_all_scalars!(check_diamond_edge_query_is_on_coedge);
659 }
660}