-
-
Notifications
You must be signed in to change notification settings - Fork 110
/
surface_surface.rs
178 lines (148 loc) · 5.47 KB
/
surface_surface.rs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
use fj_math::{Line, Point, Scalar, Vector};
use crate::objects::{CurveKind, Surface};
/// The intersection between two surfaces
#[derive(Clone, Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
pub struct SurfaceSurfaceIntersection {
/// The intersection curves, in the coordinates of the input surfaces
pub local_intersection_curves: [CurveKind<2>; 2],
/// The intersection curve, in global coordinates
pub global_intersection_curve: CurveKind<3>,
}
impl SurfaceSurfaceIntersection {
/// Compute the intersection between two surfaces
pub fn compute(a: &Surface, b: &Surface) -> Option<Self> {
// Algorithm from Real-Time Collision Detection by Christer Ericson. See
// section 5.4.4, Intersection of Two Planes.
//
// Adaptations were made to get the intersection curves in local
// coordinates for each surface.
let a_parametric = PlaneParametric::extract_from_surface(a);
let b_parametric = PlaneParametric::extract_from_surface(b);
let a = PlaneConstantNormal::from_parametric_plane(&a_parametric);
let b = PlaneConstantNormal::from_parametric_plane(&b_parametric);
let direction = a.normal.cross(&b.normal);
let denom = direction.dot(&direction);
if denom == Scalar::ZERO {
// Comparing `denom` against zero looks fishy. It's probably better
// to compare it against an epsilon value, but I don't know how
// large that epsilon should be.
//
// I'll just leave it like that, until we had the opportunity to
// collect some experience with this code.
// - @hannobraun
return None;
}
let origin = (b.normal * a.distance - a.normal * b.distance)
.cross(&direction)
/ denom;
let origin = Point { coords: origin };
let line = Line { origin, direction };
let curve_a = project_line_into_plane(&line, &a_parametric);
let curve_b = project_line_into_plane(&line, &b_parametric);
let curve_global = CurveKind::Line(Line { origin, direction });
Some(Self {
local_intersection_curves: [curve_a, curve_b],
global_intersection_curve: curve_global,
})
}
}
/// A plane in parametric form
struct PlaneParametric {
pub origin: Point<3>,
pub u: Vector<3>,
pub v: Vector<3>,
}
impl PlaneParametric {
pub fn extract_from_surface(surface: &Surface) -> Self {
let Surface::SweptCurve(surface) = surface;
let line = match surface.curve {
CurveKind::Line(line) => line,
_ => todo!("Only plane-plane intersection is currently supported."),
};
Self {
origin: line.origin,
u: line.direction,
v: surface.path,
}
}
}
/// A plane in constant-normal form
struct PlaneConstantNormal {
pub distance: Scalar,
pub normal: Vector<3>,
}
impl PlaneConstantNormal {
/// Extract a plane in constant-normal form from a `Surface`
///
/// Panics, if the given `Surface` is not a plane.
pub fn from_parametric_plane(plane: &PlaneParametric) -> Self {
// Convert plane from parametric form to three-point form.
let a = plane.origin;
let b = plane.origin + plane.u;
let c = plane.origin + plane.v;
// Convert plane from three-point form to constant-normal form. See
// Real-Time Collision Detection by Christer Ericson, section 3.6, Planes
// and Halfspaces.
let normal = (b - a).cross(&(c - a)).normalize();
let distance = normal.dot(&a.coords);
PlaneConstantNormal { distance, normal }
}
}
fn project_line_into_plane(
line: &Line<3>,
plane: &PlaneParametric,
) -> CurveKind<2> {
let line_origin_relative_to_plane = line.origin - plane.origin;
let line_origin_in_plane = Vector::from([
plane
.u
.scalar_projection_onto(&line_origin_relative_to_plane),
plane
.v
.scalar_projection_onto(&line_origin_relative_to_plane),
]);
let line_direction_in_plane = Vector::from([
plane.u.scalar_projection_onto(&line.direction),
plane.v.scalar_projection_onto(&line.direction),
]);
let line = Line {
origin: Point {
coords: line_origin_in_plane,
},
direction: line_direction_in_plane,
};
CurveKind::Line(line)
}
#[cfg(test)]
mod tests {
use fj_math::Transform;
use crate::{
algorithms::TransformObject,
objects::{CurveKind, Surface},
};
use super::SurfaceSurfaceIntersection;
#[test]
fn plane_plane() {
let xy = Surface::xy_plane();
let xz = Surface::xz_plane();
// Coincident and parallel planes don't have an intersection curve.
assert_eq!(SurfaceSurfaceIntersection::compute(&xy, &xy), None);
assert_eq!(
SurfaceSurfaceIntersection::compute(
&xy,
&xy.transform(&Transform::translation([0., 0., 1.]))
),
None,
);
let expected_xy = CurveKind::u_axis();
let expected_xz = CurveKind::u_axis();
let expected_global = CurveKind::x_axis();
assert_eq!(
SurfaceSurfaceIntersection::compute(&xy, &xz),
Some(SurfaceSurfaceIntersection {
local_intersection_curves: [expected_xy, expected_xz],
global_intersection_curve: expected_global,
})
);
}
}