Skip to main content

brep_kernel/brep/
transform_topology.rs

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    /// Whether this affine map is a finite, proper rigid placement.
11    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    /// Compose row-major transforms, applying `b` first.
32    pub fn compose(&self, b: &AffineTransform) -> Result<AffineTransform, String> {
33        // Iterator Sum starts at -0.0; the raw STEP loop starts at +0.0.
34        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    /// Invert a rigid map as `[Rᵀ | −Rᵀ·t]`. The caller must verify rigidity.
55    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]; // Rᵀ
63            }
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
126/// Transform exact BREP geometry without changing entity IDs or parameter
127/// ranges. `reverse_orientation` is required for negative-determinant maps.
128pub 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
170/// Reflect an exact BREP across the plane through `plane_point` with unit
171/// normal `n = plane_normal.normalized()`. The reflection is the affine map
172/// `p' = R p + t` with linear part `R = I - 2 n nᵀ` and translation
173/// `t = 2 (plane_point·n) n`, so points on the plane map to themselves. A
174/// reflection has negative determinant (it flips handedness), so faces are
175/// re-oriented via the existing transform with `reverse_orientation = true`,
176/// keeping outward normals.
177pub 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