1use axiolid_contracts::{GeomError, GeomResult};
169use axiolid_core::{Point3, Scalar, Vec3};
170use axiolid_mesh::TriMesh;
171
172use crate::loft::{self, Frame, Station};
173use crate::profile::Rings;
174use crate::revolve::rotate;
175use crate::sweep;
176
177mod band;
178pub use band::MitreBand;
179
180#[derive(Debug, Clone, Copy, PartialEq)]
182#[non_exhaustive]
183pub enum PathPiece {
184 Segment {
186 start: Point3,
188 end: Point3,
190 },
191 Arc {
196 centre: Point3,
198 axis: Vec3,
200 start: Point3,
202 angle: Scalar,
204 },
205}
206
207impl PathPiece {
208 #[must_use]
210 pub fn start_point(&self) -> Point3 {
211 match *self {
212 Self::Segment { start, .. } | Self::Arc { start, .. } => start,
213 }
214 }
215
216 #[must_use]
218 pub fn end_point(&self) -> Point3 {
219 self.point(1.0)
220 }
221
222 #[must_use]
224 pub fn point(&self, s: Scalar) -> Point3 {
225 match *self {
226 Self::Segment { start, end } => start + (end - start) * s,
227 Self::Arc {
228 centre,
229 axis,
230 start,
231 angle,
232 } => rotate(start, centre, axis, angle * s),
233 }
234 }
235
236 fn tangent_at(&self, point: Point3) -> Vec3 {
239 match *self {
240 Self::Segment { start, end } => (end - start).normalize_or_zero(),
241 Self::Arc {
242 centre,
243 axis,
244 angle,
245 ..
246 } => (axis.cross(point - centre) * angle.signum()).normalize_or_zero(),
247 }
248 }
249
250 #[must_use]
252 pub fn start_tangent(&self) -> Vec3 {
253 self.tangent_at(self.start_point())
254 }
255
256 #[must_use]
258 pub fn end_tangent(&self) -> Vec3 {
259 self.tangent_at(self.end_point())
260 }
261
262 #[must_use]
264 pub fn length(&self) -> Scalar {
265 match *self {
266 Self::Segment { start, end } => (end - start).length(),
267 Self::Arc {
268 centre,
269 start,
270 angle,
271 ..
272 } => (start - centre).length() * angle.abs(),
273 }
274 }
275
276 #[must_use]
278 pub fn reversed(&self) -> Self {
279 match *self {
280 Self::Segment { start, end } => Self::Segment {
281 start: end,
282 end: start,
283 },
284 Self::Arc {
285 centre,
286 axis,
287 angle,
288 ..
289 } => Self::Arc {
290 centre,
291 axis,
292 start: self.end_point(),
293 angle: -angle,
294 },
295 }
296 }
297
298 #[must_use]
301 pub fn between(&self, from: Scalar, to: Scalar) -> Self {
302 match *self {
303 Self::Segment { .. } => Self::Segment {
304 start: self.point(from),
305 end: self.point(to),
306 },
307 Self::Arc {
308 centre,
309 axis,
310 angle,
311 ..
312 } => Self::Arc {
313 centre,
314 axis,
315 start: self.point(from),
316 angle: angle * (to - from),
317 },
318 }
319 }
320
321 fn check(&self, radius: Scalar) -> GeomResult<()> {
323 let finite = |p: Point3| p.is_finite();
324 match *self {
325 Self::Segment { start, end } => {
326 if !(finite(start) && finite(end)) {
327 return Err(GeomError::InvalidInput(
328 "swept disk directrix segment must be finite".to_owned(),
329 ));
330 }
331 if start == end {
332 return Err(GeomError::Degenerate(
333 "swept disk directrix segment has no length".to_owned(),
334 ));
335 }
336 }
337 Self::Arc {
338 centre,
339 axis,
340 start,
341 angle,
342 } => {
343 if !(finite(centre) && finite(axis) && finite(start) && angle.is_finite()) {
344 return Err(GeomError::InvalidInput(
345 "swept disk directrix arc must be finite".to_owned(),
346 ));
347 }
348 let bend = (start - centre).length();
349 if (axis.length() - 1.0).abs() > 1e-9
350 || bend == 0.0
351 || (start - centre).dot(axis).abs() > 1e-9 * bend
352 {
353 return Err(GeomError::InvalidInput(
354 "swept disk directrix arc needs a unit axis normal to a non-zero \
355 radius"
356 .to_owned(),
357 ));
358 }
359 if angle == 0.0 || angle.abs() > core::f64::consts::TAU * (1.0 + 1e-12) {
360 return Err(GeomError::InvalidInput(format!(
361 "swept disk directrix arc must turn by a non-zero angle of at most \
362 a full turn, got {angle}"
363 )));
364 }
365 if bend < radius {
366 return Err(GeomError::InvalidInput(format!(
367 "swept disk radius {radius} is not below the bend radius {bend}: \
368 the tube would fold through the bend's axis"
369 )));
370 }
371 if bend == radius {
372 return Err(horn_torus("bend", radius));
373 }
374 }
375 }
376 Ok(())
377 }
378
379 fn steps(&self, budget: Scalar) -> GeomResult<usize> {
382 match *self {
383 Self::Segment { .. } => Ok(1),
384 Self::Arc {
385 centre,
386 start,
387 angle,
388 ..
389 } => crate::revolve::steps((start - centre).length(), angle, budget).map_err(|_| {
390 GeomError::BudgetExceeded {
391 resource: "swept disk bend steps",
392 }
393 }),
394 }
395 }
396}
397
398#[must_use]
402pub fn joint_tolerance(chord: Scalar, radius: Scalar) -> Scalar {
403 2.0 * (chord / (8.0 * radius)).min(1.0).asin()
404}
405
406fn turn(a: Vec3, b: Vec3) -> Scalar {
408 a.cross(b).length().atan2(a.dot(b))
409}
410
411#[derive(Debug, Clone, Copy, PartialEq)]
413#[non_exhaustive]
414pub enum MitreSide {
415 Segment {
418 corner: Point3,
420 direction: Vec3,
422 },
423 Arc(MitreBand),
425}
426
427#[derive(Debug, Clone, Copy, PartialEq)]
430#[non_exhaustive]
431pub struct ArcMitre {
432 pub joint: usize,
435 pub normal: Vec3,
438 pub sides: [MitreSide; 2],
440}
441
442#[derive(Debug, Clone, PartialEq)]
446#[non_exhaustive]
447pub struct PipePlan {
448 pub pieces: Vec<PathPiece>,
450 pub closed: bool,
452 pub arc_mitres: Vec<ArcMitre>,
454 pub twist: Scalar,
460}
461
462pub fn swept_disk_along_pieces(
472 pieces: &[PathPiece],
473 radius: Scalar,
474 inner_radius: Option<Scalar>,
475 fillet_radius: Option<Scalar>,
476 chord: Scalar,
477) -> GeomResult<TriMesh> {
478 Pipe::new(pieces, false, radius, inner_radius, fillet_radius, chord)?.mesh()
479}
480
481pub fn swept_disk_along_loop(
487 pieces: &[PathPiece],
488 radius: Scalar,
489 inner_radius: Option<Scalar>,
490 fillet_radius: Option<Scalar>,
491 chord: Scalar,
492) -> GeomResult<TriMesh> {
493 Pipe::new(pieces, true, radius, inner_radius, fillet_radius, chord)?.mesh()
494}
495
496pub fn pipe_plan(
500 pieces: &[PathPiece],
501 closed: bool,
502 radius: Scalar,
503 inner_radius: Option<Scalar>,
504 fillet_radius: Option<Scalar>,
505 chord: Scalar,
506) -> GeomResult<PipePlan> {
507 Ok(Pipe::new(pieces, closed, radius, inner_radius, fillet_radius, chord)?.plan())
508}
509
510const MAX_REFINEMENTS: usize = 16;
512
513const MAX_TWIST_STEPS: usize = 4096;
515
516const ROUNDING_TWIST: Scalar = 1e-3;
520
521const STRIP_OFFSET: Scalar = 1.0 / 16.0;
524
525struct Pipe {
527 path: Vec<PathPiece>,
528 closed: bool,
529 radius: Scalar,
530 inner: Option<Scalar>,
531 chord: Scalar,
532 mitres: Vec<Option<Vec3>>,
536 sharpest: Scalar,
538 bands: Vec<[Option<MitreBand>; 2]>,
540 twist: Scalar,
543 rate: Scalar,
544}
545
546struct Built {
549 stations: Vec<Station>,
550 end: Frame,
551 turned: Scalar,
552}
553
554#[derive(Clone, Copy)]
556enum Close<'a> {
557 Open,
559 To(&'a Station),
561 Own,
563}
564
565struct Sampled {
568 points: Vec<Point3>,
569 tangents: Vec<Vec3>,
570 twists: Vec<Scalar>,
571}
572
573impl Pipe {
574 fn new(
575 pieces: &[PathPiece],
576 closed: bool,
577 radius: Scalar,
578 inner: Option<Scalar>,
579 fillet_radius: Option<Scalar>,
580 chord: Scalar,
581 ) -> GeomResult<Self> {
582 sweep::check_radii(radius, inner)?;
583 if !(chord.is_finite() && chord > 0.0) {
584 return Err(GeomError::InvalidInput(format!(
585 "swept disk chord budget must be positive and finite, got {chord}"
586 )));
587 }
588 if let Some(f) = fillet_radius {
589 if !(f.is_finite() && f > 0.0) {
590 return Err(GeomError::InvalidInput(format!(
591 "swept disk fillet radius must be positive and finite, got {f}"
592 )));
593 }
594 }
595 if pieces.is_empty() {
596 return Err(GeomError::InvalidInput(
597 "a swept disk directrix needs at least one piece".to_owned(),
598 ));
599 }
600 let tolerance = joint_tolerance(chord, radius);
601 let path = match fillet_radius {
602 Some(f) => fillet(pieces, f, radius, tolerance, closed)?,
603 None => pieces.to_vec(),
604 };
605 for piece in &path {
606 piece.check(radius)?;
607 }
608 let m = path.len();
609 let joints = if closed { m } else { m - 1 };
610 let mut mitres: Vec<Option<Vec3>> = Vec::with_capacity(joints);
611 let mut sharpest: Scalar = 0.0;
612 for k in 0..joints {
613 let (a, b) = (path[k], path[(k + 1) % m]);
614 let gap = (b.start_point() - a.end_point()).length();
615 if gap > 0.25 * chord {
616 return Err(GeomError::InvalidInput(format!(
617 "swept disk directrix pieces leave a {gap} gap at joint {k}, more than a \
618 quarter of the {chord} chord budget"
619 )));
620 }
621 let (u, v) = (a.end_tangent(), b.start_tangent());
622 let theta = turn(u, v);
623 if theta <= tolerance {
624 mitres.push(None);
625 continue;
626 }
627 let straight = matches!(
628 (a, b),
629 (PathPiece::Segment { .. }, PathPiece::Segment { .. })
630 );
631 if fillet_radius.is_some() && !straight {
632 return Err(GeomError::InvalidInput(format!(
633 "swept disk directrix turns a corner of {theta} rad beside an arc at joint \
634 {k}, above the {tolerance} rad a tangent-continuous joint may turn at this \
635 chord budget: only corners between two straight segments are filleted"
636 )));
637 }
638 let bisector = u + v;
639 if theta >= core::f64::consts::PI - 1e-9 || bisector.length() <= 1e-9 {
640 return Err(GeomError::InvalidInput(format!(
641 "swept disk directrix reverses at joint {k}: a mitre at half angle would \
642 hold both legs"
643 )));
644 }
645 sharpest = sharpest.max(theta);
646 mitres.push(Some(bisector.normalize()));
647 }
648 let mut bands: Vec<[Option<MitreBand>; 2]> = vec![[None, None]; m];
649 for (k, normal) in mitres.iter().enumerate() {
650 let Some(n) = *normal else { continue };
651 for (i, end) in [(k, 1), ((k + 1) % m, 0)] {
652 let PathPiece::Arc {
653 centre,
654 axis,
655 angle,
656 ..
657 } = path[i]
658 else {
659 continue;
660 };
661 let (toward, corner, normal) = if end == 1 {
664 (axis * angle.signum(), path[i].end_point(), n)
665 } else {
666 (-axis * angle.signum(), path[i].start_point(), -n)
667 };
668 let (back, ahead) = band::reach(centre, toward, corner, normal, radius)
669 .ok_or_else(|| {
670 GeomError::InvalidInput(format!(
671 "swept disk directrix mitres a corner beside an arc at joint {k}, but \
672 the mitre plane does not clear the bend's tube within a quarter turn \
673 of the corner: the bend is too tight for a corner this sharp"
674 ))
675 })?;
676 bands[i][end] = Some(MitreBand {
677 centre,
678 axis: toward,
679 corner,
680 normal,
681 back,
682 ahead,
683 });
684 }
685 }
686 for (i, piece) in path.iter().enumerate() {
687 let PathPiece::Arc { angle, .. } = *piece else {
688 continue;
689 };
690 let need: Scalar = bands[i].iter().flatten().map(|b| b.back).sum();
691 if need > angle.abs() {
692 return Err(GeomError::InvalidInput(format!(
693 "swept disk directrix arc {i} turns {} rad, but the mitre bands beside it \
694 need {need} rad of it: the mitres reach past the arc",
695 angle.abs()
696 )));
697 }
698 }
699 let mut pipe = Self {
700 path,
701 closed,
702 radius,
703 inner,
704 chord,
705 mitres,
706 sharpest,
707 bands,
708 twist: 0.0,
709 rate: 0.0,
710 };
711 pipe.check_segment_reach()?;
712 if closed {
713 pipe.twist = pipe.closing_twist();
714 if pipe.twist.abs() * radius > ROUNDING_TWIST * chord {
715 let clear: Scalar = (0..m).map(|i| pipe.clear(i).1).sum();
716 if clear <= 0.0 {
717 return Err(GeomError::InvalidInput(format!(
718 "swept disk along a closed directrix: a frame carried round the loop \
719 comes back turned by {} rad about the tangent, and no part of the \
720 loop is clear of its mitres to take that twist back",
721 pipe.twist
722 )));
723 }
724 pipe.rate = pipe.twist / clear;
725 }
726 }
727 Ok(pipe)
728 }
729
730 fn before(&self, i: usize) -> Option<usize> {
732 match i.checked_sub(1) {
733 Some(k) => Some(k),
734 None => self.closed.then(|| self.path.len() - 1),
735 }
736 }
737
738 fn after(&self, i: usize) -> Option<usize> {
740 (self.closed || i + 1 < self.path.len()).then_some(i)
741 }
742
743 fn mitre_ends(&self, i: usize) -> [Option<Vec3>; 2] {
745 [
746 self.before(i).and_then(|k| self.mitres[k]),
747 self.after(i).and_then(|k| self.mitres[k]),
748 ]
749 }
750
751 fn check_segment_reach(&self) -> GeomResult<()> {
753 let radius = self.radius;
754 for (i, piece) in self.path.iter().enumerate() {
755 let PathPiece::Segment { start, end } = *piece else {
756 continue;
757 };
758 let ends = self.mitre_ends(i);
759 if ends == [None, None] {
760 continue;
761 }
762 let u = (end - start).normalize();
763 let lean = |n: Option<Vec3>| n.map_or(u, |n| n / u.dot(n));
764 let g = lean(ends[1]) - lean(ends[0]);
765 let reach = radius * (g - u * u.dot(g)).length();
766 let length = (end - start).length();
767 if reach >= length {
768 return Err(GeomError::InvalidInput(format!(
769 "swept disk directrix segment {i} is {length} long, but the mitres at its \
770 ends reach {reach} along it at the disk radius {radius}: the tube would cut \
771 through itself"
772 )));
773 }
774 }
775 Ok(())
776 }
777
778 fn hold(&self, u: Vec3, n: Vec3) -> Scalar {
782 let c = u.dot(n).abs();
783 self.radius * (1.0 - c * c).max(0.0).sqrt() / c
784 }
785
786 fn clear(&self, i: usize) -> (Scalar, Scalar) {
791 match self.path[i] {
792 PathPiece::Segment { start, end } => {
793 let u = (end - start).normalize();
794 let [z0, z1] = self
795 .mitre_ends(i)
796 .map(|n| n.map_or(0.0, |n| self.hold(u, n)));
797 let clear = (end - start).length() - z0 - z1;
798 if clear > 0.0 {
799 (z0, clear)
800 } else {
801 (0.0, 0.0)
802 }
803 }
804 PathPiece::Arc {
805 centre,
806 start,
807 angle,
808 ..
809 } => {
810 let banded: Scalar = self.bands[i].iter().flatten().map(|b| b.back).sum();
811 (
812 0.0,
813 (start - centre).length() * (angle.abs() - banded).max(0.0),
814 )
815 }
816 }
817 }
818
819 fn closing_twist(&self) -> Scalar {
824 let m = self.path.len();
825 let t0 = self.path[0].start_tangent();
826 let x0 = seed_reference(t0).normalize();
827 let mut x = x0;
828 for (i, piece) in self.path.iter().enumerate() {
829 if let PathPiece::Arc { axis, angle, .. } = *piece {
830 x = rotate(x, Point3::ZERO, axis, angle);
831 }
832 x = least_rotation(
833 x,
834 piece.end_tangent(),
835 self.path[(i + 1) % m].start_tangent(),
836 );
837 }
838 let x = (x - t0 * t0.dot(x)).normalize();
839 x.cross(x0).dot(t0).atan2(x.dot(x0))
840 }
841
842 fn plan(&self) -> PipePlan {
843 let m = self.path.len();
844 let side = |i: usize, end: usize| -> Option<MitreSide> {
845 let piece = &self.path[i];
846 match piece {
847 PathPiece::Segment { .. } => Some(MitreSide::Segment {
848 corner: if end == 1 {
849 piece.end_point()
850 } else {
851 piece.start_point()
852 },
853 direction: piece.start_tangent(),
854 }),
855 PathPiece::Arc { .. } => self.bands[i][end].map(MitreSide::Arc),
856 }
857 };
858 let arc_mitres = self
859 .mitres
860 .iter()
861 .enumerate()
862 .filter_map(|(k, normal)| {
863 let normal = (*normal)?;
864 let sides = [side(k, 1)?, side((k + 1) % m, 0)?];
865 matches!(sides, [MitreSide::Arc(_), _] | [_, MitreSide::Arc(_)]).then_some(
866 ArcMitre {
867 joint: k,
868 normal,
869 sides,
870 },
871 )
872 })
873 .collect();
874 PipePlan {
875 pieces: self.path.clone(),
876 closed: self.closed,
877 arc_mitres,
878 twist: self.twist,
879 }
880 }
881
882 fn mesh(&self) -> GeomResult<TriMesh> {
883 let rings = sweep::disk_rings(
884 self.radius,
885 self.inner,
886 0.5 * self.chord * (0.5 * self.sharpest).cos(),
887 )?;
888 let m = self.path.len();
889 let mut stations: Vec<Station> = Vec::new();
890 let mut carried: Option<Frame> = None;
891 let mut turned: Scalar = 0.0;
892 for i in 0..m {
893 let start = self.start_frame(i, carried.as_ref())?;
894 let closing = self.closed && i + 1 == m;
895 let close = match (closing, stations.first()) {
896 (false, _) => Close::Open,
897 (true, Some(first)) => Close::To(first),
898 (true, None) => Close::Own,
899 };
900 let built = self.piece_stations(i, &rings, &start, stations.last(), close, turned)?;
901 let skip = usize::from(!stations.is_empty());
904 let keep = built.stations.len() - usize::from(closing);
905 stations.extend(built.stations.into_iter().take(keep).skip(skip));
906 carried = Some(built.end);
907 turned = built.turned;
908 }
909 let mut mesh = loft::loft(&rings, &stations, self.closed)?;
910 if matches!(crate::extrude::outward_orientation(&mesh), Some(false)) {
913 for triangle in mesh.indices.as_chunks_mut::<3>().0 {
914 triangle.swap(1, 2);
915 }
916 }
917 Ok(mesh)
918 }
919
920 fn start_frame(&self, i: usize, carried: Option<&Frame>) -> GeomResult<Frame> {
923 let piece = &self.path[i];
924 let tangent = piece.start_tangent();
925 let seed = carried.map_or_else(
926 || seed_reference(tangent),
927 |f| least_rotation(f.x, f.x.cross(f.y), tangent),
928 );
929 Frame::from_reference(piece.start_point(), tangent, seed)
930 }
931
932 fn piece_stations(
933 &self,
934 i: usize,
935 rings: &Rings,
936 start: &Frame,
937 joint: Option<&Station>,
938 close: Close<'_>,
939 turned: Scalar,
940 ) -> GeomResult<Built> {
941 match self.path[i] {
942 PathPiece::Segment { .. } => {
943 self.segment_stations(i, rings, start, joint, close, turned)
944 }
945 PathPiece::Arc { .. } => self.arc_stations(i, rings, start, joint, close, turned),
946 }
947 }
948
949 fn segment_stations(
953 &self,
954 i: usize,
955 rings: &Rings,
956 start: &Frame,
957 joint: Option<&Station>,
958 close: Close<'_>,
959 turned: Scalar,
960 ) -> GeomResult<Built> {
961 let piece = self.path[i];
962 let length = piece.length();
963 let tangent = piece.start_tangent();
964 let (held, clear) = self.clear(i);
965 let tail = if clear > 0.0 {
966 length - held - clear
967 } else {
968 0.0
969 };
970 let twist = self.rate * clear;
971 let sample = |budget: Scalar| -> GeomResult<Sampled> {
972 let mut at = vec![0.0];
973 if twist != 0.0 {
974 let n = twist_steps(twist, budget, self.radius)?;
975 if held > 0.0 {
976 at.push(held);
977 }
978 at.extend((1..n).map(|k| held + clear * k as Scalar / n as Scalar));
979 if tail > 0.0 {
980 at.push(length - tail);
981 }
982 }
983 at.push(length);
984 Ok(Sampled {
985 points: at.iter().map(|s| piece.point(s / length)).collect(),
986 tangents: vec![tangent; at.len()],
987 twists: at
988 .iter()
989 .map(|s| turned + self.rate * (s - held).clamp(0.0, clear))
990 .collect(),
991 })
992 };
993 let (frames, placed) = bounded_piece(
994 sample,
995 rings,
996 start.x,
997 self.radius,
998 0.5 * self.chord,
999 joint,
1000 self.mitre_ends(i),
1001 close,
1002 )?;
1003 Ok(Built {
1004 stations: placed,
1005 end: last_frame(frames)?,
1006 turned: turned + twist,
1007 })
1008 }
1009
1010 fn arc_stations(
1013 &self,
1014 i: usize,
1015 rings: &Rings,
1016 start: &Frame,
1017 joint: Option<&Station>,
1018 close: Close<'_>,
1019 turned: Scalar,
1020 ) -> GeomResult<Built> {
1021 let piece = self.path[i];
1022 let PathPiece::Arc {
1023 centre,
1024 start: from,
1025 angle,
1026 axis,
1027 } = piece
1028 else {
1029 return Err(not_an_arc());
1030 };
1031 let total = angle.abs();
1032 let along = axis * angle.signum();
1033 let reach_radius = (from - centre).length() + self.radius;
1034 let share = 0.5 * self.chord;
1035 let band_tol = 0.5 * share;
1036 let offset = STRIP_OFFSET * self.chord;
1037 let [first, last] = self.bands[i];
1038 let (head, tail) = (first.map_or(0.0, |b| b.back), last.map_or(0.0, |b| b.back));
1039 let place = |frame: &Frame| loft::place(rings, |p| loft::at(&twisted(frame, turned), p));
1040 let carried = |frame: &Frame, by: Scalar| Frame {
1041 origin: rotate(frame.origin, centre, along, by),
1042 x: rotate(frame.x, Point3::ZERO, along, by),
1043 y: rotate(frame.y, Point3::ZERO, along, by),
1044 };
1045 let mut out: Vec<Station> = Vec::new();
1046 let mut frame = start.clone();
1047 if let Some(band) = first {
1048 let k = self.before(i).ok_or_else(not_an_arc)?;
1049 let shared = self.shared_ring(k, &twisted(start, turned), rings)?;
1050 frame = carried(start, head);
1051 let far = place(&frame);
1052 let strip = matches!(self.path[k], PathPiece::Segment { .. });
1055 let mut toward = band.stations_toward(
1056 &far,
1057 &shared,
1058 strip.then_some(offset),
1059 reach_radius,
1060 band_tol,
1061 )?;
1062 toward.reverse();
1063 out.extend(toward);
1064 out.push(far);
1065 }
1066 let middle = total - head - tail;
1067 let mut twist = 0.0;
1068 if middle > 0.0 {
1069 let part = piece.between(head / total, 1.0 - tail / total);
1070 twist = self.rate * part.length();
1071 let sample = |budget: Scalar| -> GeomResult<Sampled> {
1072 let n = part
1073 .steps(budget)?
1074 .max(twist_steps(twist, budget, self.radius)?);
1075 let points: Vec<Point3> = (0..=n)
1076 .map(|k| part.point(k as Scalar / n as Scalar))
1077 .collect();
1078 let tangents = points.iter().map(|p| part.tangent_at(*p)).collect();
1079 let twists = (0..=n)
1080 .map(|k| turned + twist * k as Scalar / n as Scalar)
1081 .collect();
1082 Ok(Sampled {
1083 points,
1084 tangents,
1085 twists,
1086 })
1087 };
1088 let (joint, close) = (
1089 if first.is_some() { out.last() } else { joint },
1090 if last.is_some() { Close::Open } else { close },
1091 );
1092 let (frames, placed) = bounded_piece(
1093 sample,
1094 rings,
1095 frame.x,
1096 self.radius,
1097 share,
1098 joint,
1099 [None, None],
1100 close,
1101 )?;
1102 let skip = usize::from(first.is_some());
1103 out.extend(placed.into_iter().skip(skip));
1104 frame = last_frame(frames)?;
1105 } else if out.is_empty() {
1106 out.push(place(start));
1107 }
1108 let turned_end = turned + twist;
1109 if let Some(band) = last {
1110 let k = self.after(i).ok_or_else(not_an_arc)?;
1111 let end = carried(&frame, tail);
1112 let target = match close {
1113 Close::To(first) => first.clone(),
1114 Close::Own => out[0].clone(),
1115 Close::Open => {
1116 let next = self.start_frame((k + 1) % self.path.len(), Some(&end))?;
1117 self.shared_ring(k, &twisted(&next, turned_end), rings)?
1118 }
1119 };
1120 let far = out.last().cloned().ok_or_else(not_an_arc)?;
1121 out.extend(band.stations_toward(
1122 &far,
1123 &target,
1124 Some(offset),
1125 reach_radius,
1126 band_tol,
1127 )?);
1128 frame = end;
1129 }
1130 Ok(Built {
1131 stations: out,
1132 end: frame,
1133 turned: turned_end,
1134 })
1135 }
1136
1137 fn shared_ring(&self, k: usize, frame: &Frame, rings: &Rings) -> GeomResult<Station> {
1143 let next = (k + 1) % self.path.len();
1144 let normal = self.mitres[k].ok_or_else(not_an_arc)?;
1145 let square = loft::place(rings, |p| loft::at(frame, p));
1146 match (self.path[k], self.path[next]) {
1147 (PathPiece::Arc { .. }, PathPiece::Arc { .. }) => {
1148 self.bands[next][0].ok_or_else(not_an_arc)?.cut(&square)
1149 }
1150 _ => Ok(mitred(&square, frame, normal)),
1151 }
1152 }
1153}
1154
1155fn not_an_arc() -> GeomError {
1157 GeomError::Degenerate("swept disk: a mitre band without its arc or corner".to_owned())
1158}
1159
1160fn last_frame(frames: Vec<Frame>) -> GeomResult<Frame> {
1161 frames
1162 .into_iter()
1163 .last()
1164 .ok_or_else(|| GeomError::Degenerate("swept disk: a piece without a station".to_owned()))
1165}
1166
1167fn twisted(frame: &Frame, angle: Scalar) -> Frame {
1169 if angle == 0.0 {
1170 return frame.clone();
1171 }
1172 let (s, c) = angle.sin_cos();
1173 Frame {
1174 origin: frame.origin,
1175 x: frame.x * c + frame.y * s,
1176 y: frame.y * c - frame.x * s,
1177 }
1178}
1179
1180fn twist_steps(twist: Scalar, budget: Scalar, radius: Scalar) -> GeomResult<usize> {
1183 if twist == 0.0 {
1184 return Ok(1);
1185 }
1186 let per = 2.0 * (1.0 - (budget / radius).min(1.0)).acos();
1187 let n = (twist.abs() / per).ceil();
1188 if n.is_nan() || n > MAX_TWIST_STEPS as Scalar {
1189 return Err(GeomError::BudgetExceeded {
1190 resource: "swept disk twist steps",
1191 });
1192 }
1193 Ok((n as usize).max(1))
1194}
1195
1196#[allow(clippy::too_many_arguments)]
1204fn bounded_piece(
1205 sample: impl Fn(Scalar) -> GeomResult<Sampled>,
1206 rings: &Rings,
1207 seed: Vec3,
1208 reach: Scalar,
1209 share: Scalar,
1210 joint: Option<&Station>,
1211 ends: [Option<Vec3>; 2],
1212 close: Close<'_>,
1213) -> GeomResult<(Vec<Frame>, Vec<Station>)> {
1214 let mut b = share;
1215 for _ in 0..MAX_REFINEMENTS {
1216 let Sampled {
1217 points,
1218 tangents,
1219 twists,
1220 } = sample(b)?;
1221 let frames = sweep::frames_carried(&points, &tangents, seed)?;
1222 let placing: Vec<Frame> = frames
1223 .iter()
1224 .zip(&twists)
1225 .map(|(f, a)| twisted(f, *a))
1226 .collect();
1227 let mut placed: Vec<Station> = placing
1228 .iter()
1229 .map(|f| loft::place(rings, |p| loft::at(f, p)))
1230 .collect();
1231 let last = placed.len() - 1;
1232 for (k, normal) in [(0, ends[0]), (last, ends[1])] {
1233 if let Some(n) = normal {
1234 placed[k] = mitred(&placed[k], &placing[k], n);
1235 }
1236 }
1237 let shift = joint.map_or(0.0, |j| displacement(j, &placed[0]));
1238 let closing = match close {
1239 Close::Open => 0.0,
1240 Close::To(station) => displacement(station, &placed[last]),
1241 Close::Own => displacement(&placed[0], &placed[last]),
1242 };
1243 let worst = placing
1244 .windows(2)
1245 .zip(placed.windows(2))
1246 .enumerate()
1247 .map(|(k, (f, s))| {
1248 let moved = if k == 0 { shift } else { 0.0 };
1249 let closed = if k + 1 == last { closing } else { 0.0 };
1250 sweep::span_bound(&f[0], &f[1], b, reach)
1251 + loft::span_deviation(&s[0], &s[1])
1252 + moved
1253 + closed
1254 })
1255 .fold(0.0, Scalar::max);
1256 if worst <= share {
1257 return Ok((frames, placed));
1258 }
1259 b *= 0.5;
1260 }
1261 Err(GeomError::BudgetExceeded {
1262 resource: "swept disk directrix refinement",
1263 })
1264}
1265
1266fn mitred(station: &Station, frame: &Frame, normal: Vec3) -> Station {
1270 let tangent = frame.x.cross(frame.y);
1271 let along = tangent.dot(normal);
1272 Station {
1273 loops: station
1274 .loops
1275 .iter()
1276 .map(|ring| {
1277 ring.iter()
1278 .map(|q| *q - tangent * ((*q - frame.origin).dot(normal) / along))
1279 .collect()
1280 })
1281 .collect(),
1282 }
1283}
1284
1285fn horn_torus(what: &str, radius: Scalar) -> GeomError {
1288 GeomError::InvalidInput(format!(
1289 "swept disk radius {radius} equals the {what} radius: the format rule permits a \
1290 fillet radius equal to the disk radius, but the bend is then a horn torus whose \
1291 inner wall degenerates to one point on the bend's axis, a boundary that is not a \
1292 two-manifold and that no closed mesh bounds without meeting itself there"
1293 ))
1294}
1295
1296fn displacement(a: &Station, b: &Station) -> Scalar {
1298 a.loops
1299 .iter()
1300 .zip(&b.loops)
1301 .flat_map(|(p, q)| p.iter().zip(q).map(|(u, v)| (*u - *v).length()))
1302 .fold(0.0, Scalar::max)
1303}
1304
1305fn least_rotation(v: Vec3, from: Vec3, to: Vec3) -> Vec3 {
1307 let axis = from.cross(to);
1308 let sin = axis.length();
1309 if sin == 0.0 {
1310 return v;
1311 }
1312 rotate(v, Point3::ZERO, axis / sin, turn(from, to))
1313}
1314
1315fn seed_reference(t: Vec3) -> Vec3 {
1317 let candidate = if t.x.abs() < 0.9 { Vec3::X } else { Vec3::Y };
1318 candidate - t * t.dot(candidate)
1319}
1320
1321fn fillet(
1325 pieces: &[PathPiece],
1326 fillet: Scalar,
1327 radius: Scalar,
1328 tolerance: Scalar,
1329 closed: bool,
1330) -> GeomResult<Vec<PathPiece>> {
1331 struct Leg {
1334 piece: PathPiece,
1335 dir: Vec3,
1336 trim_start: Scalar,
1337 trim_end: Scalar,
1338 }
1339 let mut legs: Vec<Leg> = Vec::with_capacity(pieces.len());
1340 for piece in pieces {
1341 piece.check(radius)?;
1342 legs.push(Leg {
1343 piece: *piece,
1344 dir: piece.start_tangent(),
1345 trim_start: 0.0,
1346 trim_end: 0.0,
1347 });
1348 }
1349 let m = legs.len();
1350 let mut arcs: Vec<Option<PathPiece>> = vec![None; m];
1352 let first = if closed { 0 } else { 1 };
1353 for k in first..m {
1354 let prev = (k + m - 1) % m;
1355 if prev == k {
1356 continue;
1357 }
1358 let (PathPiece::Segment { .. }, PathPiece::Segment { start: corner, .. }) =
1359 (legs[prev].piece, legs[k].piece)
1360 else {
1361 continue;
1362 };
1363 let (u, v) = (legs[prev].dir, legs[k].dir);
1364 let theta = turn(u, v);
1365 if theta <= tolerance {
1366 continue;
1367 }
1368 if radius > fillet {
1369 return Err(GeomError::InvalidInput(format!(
1370 "swept disk radius {radius} is above the fillet radius {fillet}: the \
1371 tube would fold through the bend's axis"
1372 )));
1373 }
1374 if radius == fillet {
1375 return Err(horn_torus("fillet", radius));
1376 }
1377 let normal = u.cross(v);
1378 if normal.length() <= 1e-12 {
1379 return Err(GeomError::InvalidInput(format!(
1380 "swept disk directrix reverses at corner {k}; no fillet rounds it"
1381 )));
1382 }
1383 let cut = fillet * (0.5 * theta).tan();
1384 legs[prev].trim_end = cut;
1385 legs[k].trim_start = cut;
1386 let tangent_point = corner - u * cut;
1387 let inward = (v - u * u.dot(v)).normalize();
1388 arcs[k] = Some(PathPiece::Arc {
1389 centre: tangent_point + inward * fillet,
1390 axis: normal.normalize(),
1391 start: tangent_point,
1392 angle: theta,
1393 });
1394 }
1395 let mut out = Vec::with_capacity(2 * m);
1396 for (k, leg) in legs.iter().enumerate() {
1397 if k > 0 {
1398 if let Some(arc) = arcs[k] {
1399 out.push(arc);
1400 }
1401 }
1402 let length = leg.piece.length();
1403 let needed = leg.trim_start + leg.trim_end;
1404 if needed == 0.0 {
1405 out.push(leg.piece);
1406 continue;
1407 }
1408 let slack = 1e-9 * length;
1411 if needed > length + slack {
1412 return Err(GeomError::InvalidInput(format!(
1413 "swept disk fillet radius {fillet} does not fit: segment {k} is {length} \
1414 long but its fillets need {needed} of it"
1415 )));
1416 }
1417 if needed < length - slack {
1418 out.push(
1419 leg.piece
1420 .between(leg.trim_start / length, 1.0 - leg.trim_end / length),
1421 );
1422 }
1423 }
1424 if closed {
1425 if let Some(arc) = arcs[0] {
1426 out.push(arc);
1427 }
1428 }
1429 Ok(out)
1430}
1431
1432#[cfg(test)]
1433mod tests {
1434 use super::*;
1435
1436 #[test]
1437 fn a_fillet_is_tangent_to_both_segments() {
1438 let corner = Point3::new(1.0, 0.0, 0.0);
1439 let pieces = [
1440 PathPiece::Segment {
1441 start: Point3::ZERO,
1442 end: corner,
1443 },
1444 PathPiece::Segment {
1445 start: corner,
1446 end: Point3::new(1.0, 1.0, 1.0),
1447 },
1448 ];
1449 let path = fillet(&pieces, 0.2, 0.05, 1e-6, false).unwrap();
1450 assert_eq!(path.len(), 3);
1451 for pair in path.windows(2) {
1452 assert!((pair[0].end_point() - pair[1].start_point()).length() < 1e-12);
1453 assert!(turn(pair[0].end_tangent(), pair[1].start_tangent()) < 1e-12);
1454 }
1455 }
1456
1457 #[test]
1458 fn a_mitre_ring_lies_on_both_cylinders_and_on_the_bisector_plane() {
1459 let corner = Point3::new(1.0, 0.2, -0.3);
1460 let (u, v) = (
1461 Vec3::new(1.0, 0.0, 0.0),
1462 Vec3::new(0.2, 0.7, -0.4).normalize(),
1463 );
1464 let pieces = [
1465 PathPiece::Segment {
1466 start: corner - u,
1467 end: corner,
1468 },
1469 PathPiece::Segment {
1470 start: corner,
1471 end: corner + v,
1472 },
1473 ];
1474 let (r, inner) = (0.05, 0.03);
1475 let mesh = swept_disk_along_pieces(&pieces, r, Some(inner), None, 1e-3).unwrap();
1476 let per = mesh.positions.len() / 3;
1478 let n = (u + v).normalize();
1479 let off = |p: Point3, axis: Vec3| {
1480 let d = p - corner;
1481 (d - axis * d.dot(axis)).length()
1482 };
1483 for p in &mesh.positions[per..2 * per] {
1484 assert!((*p - corner).dot(n).abs() < 1e-12);
1485 let (a, b) = (off(*p, u), off(*p, v));
1486 assert!((a - b).abs() < 1e-12, "{a} {b}");
1487 assert!((a - r).abs() < 1e-12 || (a - inner).abs() < 1e-12, "{a}");
1488 }
1489 }
1490
1491 #[test]
1492 fn an_arc_walked_backwards_retraces_itself() {
1493 let arc = PathPiece::Arc {
1494 centre: Point3::new(0.0, 1.0, 0.0),
1495 axis: Vec3::Z,
1496 start: Point3::ZERO,
1497 angle: 1.0,
1498 };
1499 let back = arc.reversed();
1500 assert!((back.end_point() - arc.start_point()).length() < 1e-12);
1501 assert!((back.start_tangent() + arc.end_tangent()).length() < 1e-12);
1502 let part = arc.between(0.25, 0.75);
1503 assert!((part.length() - 0.5).abs() < 1e-12);
1504 assert!((part.start_point() - arc.point(0.25)).length() < 1e-12);
1505 }
1506}