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 mode = carry_verdict_mode();
163    let carried = mode != crate::env_read::Mode::Off
164        && carries_clean_verdict(solid, transform, reverse_orientation);
165    if carried && mode == crate::env_read::Mode::On {
166        CARRIED_VERDICTS.with(|count| count.set(count.get() + 1));
167        if crate::profile::profile_enabled() {
168            profile_line!("transform.verdict carried=1");
169        }
170        return Ok(result);
171    }
172    let issues = result.validate();
173    if mode == crate::env_read::Mode::Check && crate::profile::profile_enabled() {
174        // What carrying would have answered, against the copy's own verdict.
175        profile_line!(
176            "transform.verdict check carried={} issues={} disagree={}",
177            u8::from(carried),
178            issues.len(),
179            u8::from(carried && !issues.is_empty())
180        );
181    }
182    if issues.is_empty() {
183        Ok(result)
184    } else {
185        Err(format!("transformed BREP is invalid: {issues:?}"))
186    }
187}
188
189/// Whether a rigid placement may keep its source's clean validation verdict
190/// instead of validating the moved copy. A proper rigid map (the shared
191/// `RIGID_TOL` band components are placed under) changes no distance, angle,
192/// parameter or topology. It does change one input of `validate`: the
193/// identity and pcurve bands are sized from the AXIS-ALIGNED bounding box of
194/// the vertices, whose diagonal moves under a rotation, so a copy read with
195/// its own (narrower) bands could report an issue its source does not. The
196/// carried verdict is the source's, read in the source's frame. The source's
197/// verdict comes from the content-keyed validate memo: the first placement of
198/// a part validates the part once, each further placement of the same
199/// content reads the memo. A source with issues is never carried; its copy is
200/// validated as before. `BREP_TRANSFORM_CARRY_VERDICT`: off by default, `1`
201/// carries, `check` validates every copy anyway and prints (`BREP_PROFILE`)
202/// whether carrying would have disagreed.
203fn carries_clean_verdict(
204    solid: &BrepSolid,
205    transform: AffineTransform,
206    reverse_orientation: bool,
207) -> bool {
208    !reverse_orientation
209        && transform.is_rigid(crate::RIGID_TOL)
210        && solid.validate_memoised().is_empty()
211}
212
213thread_local! {
214    /// Placements on this thread that kept their source's verdict.
215    static CARRIED_VERDICTS: std::cell::Cell<u64> = const { std::cell::Cell::new(0) };
216}
217
218/// How many placements on this thread have kept their source's verdict.
219pub(crate) fn carried_verdicts() -> u64 {
220    CARRIED_VERDICTS.with(|count| count.get())
221}
222
223fn carry_verdict_mode() -> crate::env_read::Mode {
224    crate::mode_switch!("BREP_TRANSFORM_CARRY_VERDICT", false)
225}
226
227/// Reflect an exact BREP across the plane through `plane_point` with unit
228/// normal `n = plane_normal.normalized()`. The reflection is the affine map
229/// `p' = R p + t` with linear part `R = I - 2 n nᵀ` and translation
230/// `t = 2 (plane_point·n) n`, so points on the plane map to themselves. A
231/// reflection has negative determinant (it flips handedness), so faces are
232/// re-oriented via the existing transform with `reverse_orientation = true`,
233/// keeping outward normals.
234pub fn mirror_brep(
235    solid: &BrepSolid,
236    plane_point: Vec3,
237    plane_normal: Vec3,
238) -> Result<BrepSolid, String> {
239    let n = plane_normal.normalized()?;
240    let (nx, ny, nz) = (n.x, n.y, n.z);
241    let d = plane_point.dot(n);
242    let (tx, ty, tz) = (2.0 * d * nx, 2.0 * d * ny, 2.0 * d * nz);
243    let matrix = [
244        1.0 - 2.0 * nx * nx,
245        -2.0 * nx * ny,
246        -2.0 * nx * nz,
247        tx,
248        -2.0 * ny * nx,
249        1.0 - 2.0 * ny * ny,
250        -2.0 * ny * nz,
251        ty,
252        -2.0 * nz * nx,
253        -2.0 * nz * ny,
254        1.0 - 2.0 * nz * nz,
255        tz,
256        0.0,
257        0.0,
258        0.0,
259        1.0,
260    ];
261    let transform = AffineTransform::new(matrix)?;
262    transform_brep(solid, transform, true)
263}
264
265
266pub(crate) fn multiply_elements(a: &[f64; 16], b: &[f64; 16]) -> [f64; 16] {
267    multiply_elements_seeded(a, b, 0.0)
268}
269
270fn multiply_elements_seeded(a: &[f64; 16], b: &[f64; 16], seed: f64) -> [f64; 16] {
271    let mut out = [0.0f64; 16];
272    for row in 0..4 {
273        for col in 0..4 {
274            let mut sum = seed;
275            for k in 0..4 {
276                sum += a[row * 4 + k] * b[k * 4 + col];
277            }
278            out[row * 4 + col] = sum;
279        }
280    }
281    out
282}
283