1use crate::topology::BrepSolid;
2use crate::{NurbsCurve, NurbsSurface, Vec3, Vec4};
3
4#[derive(Clone, Copy, Debug)]
5pub struct AffineTransform {
6 pub elements: [f64; 16],
7}
8
9impl AffineTransform {
10 pub fn is_rigid(&self, tolerance: f64) -> bool {
12 let m = &self.elements;
13 if !tolerance.is_finite() || tolerance < 0.0 || m.iter().any(|x| !x.is_finite()) {
14 return false;
15 }
16 if m[12].abs() > tolerance || m[13].abs() > tolerance || m[14].abs() > tolerance
17 || (m[15] - 1.0).abs() > tolerance {
18 return false;
19 }
20 let rows = [[m[0], m[1], m[2]], [m[4], m[5], m[6]], [m[8], m[9], m[10]]];
21 for i in 0..3 {
22 for j in i..3 {
23 let dot: f64 = (0..3).map(|k| rows[i][k] * rows[j][k]).sum();
24 let expected = if i == j { 1.0 } else { 0.0 };
25 if (dot - expected).abs() > tolerance { return false; }
26 }
27 }
28 (self.determinant3() - 1.0).abs() <= tolerance
29 }
30
31 pub fn compose(&self, b: &AffineTransform) -> Result<AffineTransform, String> {
33 AffineTransform::new(multiply_elements_seeded(&self.elements, &b.elements, -0.0))
35 }
36
37 pub fn new(elements: [f64; 16]) -> Result<Self, String> {
38 if elements.iter().any(|value| !value.is_finite()) {
39 return Err("AffineTransform: matrix must be finite".into());
40 }
41 if elements[12].abs() > 1e-12
42 || elements[13].abs() > 1e-12
43 || elements[14].abs() > 1e-12
44 || (elements[15] - 1.0).abs() > 1e-12
45 {
46 return Err("AffineTransform: projective matrices are unsupported".into());
47 }
48 if (Self { elements }).determinant3().abs() <= 1e-14 {
49 return Err("AffineTransform: singular matrix".into());
50 }
51 Ok(Self { elements })
52 }
53
54 pub fn rigid_inverse(&self) -> Result<Self, String> {
56 let m = &self.elements;
57 let r = [[m[0], m[1], m[2]], [m[4], m[5], m[6]], [m[8], m[9], m[10]]];
58 let t = [m[3], m[7], m[11]];
59 let mut out = [0.0f64; 16];
60 for row in 0..3 {
61 for col in 0..3 {
62 out[row * 4 + col] = r[col][row]; }
64 out[row * 4 + 3] = -(0..3).map(|k| r[k][row] * t[k]).sum::<f64>();
65 }
66 out[15] = 1.0;
67 AffineTransform::new(out)
68 }
69
70 pub fn determinant3(self) -> f64 {
71 let m = self.elements;
72 m[0] * (m[5] * m[10] - m[6] * m[9]) - m[1] * (m[4] * m[10] - m[6] * m[8])
73 + m[2] * (m[4] * m[9] - m[5] * m[8])
74 }
75
76 pub fn point(self, point: Vec3) -> Vec3 {
77 let m = self.elements;
78 Vec3::new(
79 m[0] * point.x + m[1] * point.y + m[2] * point.z + m[3],
80 m[4] * point.x + m[5] * point.y + m[6] * point.z + m[7],
81 m[8] * point.x + m[9] * point.y + m[10] * point.z + m[11],
82 )
83 }
84
85 fn homogeneous(self, point: Vec4) -> Result<Vec4, String> {
86 Ok(Vec4::from_point(self.point(point.point()?), point.w))
87 }
88}
89
90pub(crate) fn transform_curve(
91 curve: &NurbsCurve,
92 transform: AffineTransform,
93) -> Result<NurbsCurve, String> {
94 NurbsCurve::new(
95 curve.degree,
96 curve.knots.clone(),
97 curve
98 .control_points
99 .iter()
100 .map(|point| transform.homogeneous(*point))
101 .collect::<Result<_, _>>()?,
102 )
103}
104
105pub(crate) fn transform_surface(
106 surface: &NurbsSurface,
107 transform: AffineTransform,
108) -> Result<NurbsSurface, String> {
109 NurbsSurface::new(
110 surface.degree_u,
111 surface.degree_v,
112 surface.knots_u.clone(),
113 surface.knots_v.clone(),
114 surface
115 .control_points
116 .iter()
117 .map(|row| {
118 row.iter()
119 .map(|point| transform.homogeneous(*point))
120 .collect::<Result<Vec<_>, _>>()
121 })
122 .collect::<Result<_, _>>()?,
123 )
124}
125
126pub fn transform_brep(
129 solid: &BrepSolid,
130 transform: AffineTransform,
131 reverse_orientation: bool,
132) -> Result<BrepSolid, String> {
133 let determinant = transform.determinant3();
134 if determinant < 0.0 && !reverse_orientation {
135 return Err("transformSolid: reflection requires orientation reversal".into());
136 }
137 if determinant > 0.0 && reverse_orientation {
138 return Err("transformSolid: orientation reversal requires a reflection".into());
139 }
140 let mut result = solid.clone();
141 for vertex in &mut result.vertices {
142 vertex.point = transform.point(vertex.point);
143 }
144 for edge in &mut result.edges {
145 edge.curve = transform_curve(&edge.curve, transform)?;
146 }
147 for shell in &mut result.shells {
148 for face in &mut shell.faces {
149 face.surface = transform_surface(&face.surface, transform)?;
150 if reverse_orientation {
151 face.same_sense = !face.same_sense;
152 for loop_record in &mut face.loops {
153 loop_record.coedges.reverse();
154 for coedge in &mut loop_record.coedges {
155 coedge.forward = !coedge.forward;
156 coedge.pcurve = coedge.pcurve.reversed()?;
157 }
158 }
159 }
160 }
161 }
162 let issues = result.validate();
163 if issues.is_empty() {
164 Ok(result)
165 } else {
166 Err(format!("transformed BREP is invalid: {issues:?}"))
167 }
168}
169
170pub fn mirror_brep(
178 solid: &BrepSolid,
179 plane_point: Vec3,
180 plane_normal: Vec3,
181) -> Result<BrepSolid, String> {
182 let n = plane_normal.normalized()?;
183 let (nx, ny, nz) = (n.x, n.y, n.z);
184 let d = plane_point.dot(n);
185 let (tx, ty, tz) = (2.0 * d * nx, 2.0 * d * ny, 2.0 * d * nz);
186 let matrix = [
187 1.0 - 2.0 * nx * nx,
188 -2.0 * nx * ny,
189 -2.0 * nx * nz,
190 tx,
191 -2.0 * ny * nx,
192 1.0 - 2.0 * ny * ny,
193 -2.0 * ny * nz,
194 ty,
195 -2.0 * nz * nx,
196 -2.0 * nz * ny,
197 1.0 - 2.0 * nz * nz,
198 tz,
199 0.0,
200 0.0,
201 0.0,
202 1.0,
203 ];
204 let transform = AffineTransform::new(matrix)?;
205 transform_brep(solid, transform, true)
206}
207
208
209pub(crate) fn multiply_elements(a: &[f64; 16], b: &[f64; 16]) -> [f64; 16] {
210 multiply_elements_seeded(a, b, 0.0)
211}
212
213fn multiply_elements_seeded(a: &[f64; 16], b: &[f64; 16], seed: f64) -> [f64; 16] {
214 let mut out = [0.0f64; 16];
215 for row in 0..4 {
216 for col in 0..4 {
217 let mut sum = seed;
218 for k in 0..4 {
219 sum += a[row * 4 + k] * b[k * 4 + col];
220 }
221 out[row * 4 + col] = sum;
222 }
223 }
224 out
225}
226