fix(vector): bound hostile geometry and raster work (#8769)

* fix(vector): bound hostile geometry raster work

Guard geometry and tracing operations against malformed dimensions and degenerate inputs with bounded work and safe failures.

🤖 Generated with Codebuff
Co-Authored-By: Codebuff <noreply@codebuff.com>

* fix(vector): address review (knot cap on split, relative identity tol, non-vacuous test)

* revert(geom): splitting past the quilt limit is intentional (smart_map test)

* fix(geom): cap identity tolerance at 1e-6 px

---------

Co-authored-by: Codebuff <noreply@codebuff.com>

crates/geom/src/warp.rs

@@ -474,38 +474,50 @@ impl BezierMesh {
 
     /// Control points per row.
     pub fn nx(&self) -> usize {
-        3 * (self.us.len() - 1) + 1
+        3usize.saturating_mul(self.us.len().saturating_sub(1)).saturating_add(1)
     }
     /// Control point rows.
     pub fn ny(&self) -> usize {
-        3 * (self.vs.len() - 1) + 1
+        3usize.saturating_mul(self.vs.len().saturating_sub(1)).saturating_add(1)
     }
     pub fn point(&self, i: usize, j: usize) -> [f64; 2] {
-        self.points[j * self.nx() + i]
+        self.points.get(j.saturating_mul(self.nx()).saturating_add(i)).copied().unwrap_or([0.0; 2])
     }
 
     /// Checks the invariants (knots and point count); malformed data from files is rejected.
     pub fn is_valid(&self) -> bool {
-        let knots_ok = |k: &[f64]| k.len() >= 2 && k[0] == 0.0 && k[k.len() - 1] == 1.0 && k.windows(2).all(|w| w[1] > w[0]);
-        knots_ok(&self.us)
-            && knots_ok(&self.vs)
-            && self.points.len() == self.nx() * self.ny()
-            && self.points.iter().all(|p| p[0].is_finite() && p[1].is_finite())
+        let knots_ok = |k: &[f64]| (2..=65).contains(&k.len()) && k[0] == 0.0 && k[k.len() - 1] == 1.0 && k.windows(2).all(|w| w[1] > w[0]);
+        if !knots_ok(&self.us) || !knots_ok(&self.vs) {
+            return false;
+        }
+        let Some(nx) = (self.us.len() - 1).checked_mul(3).and_then(|v| v.checked_add(1)) else {
+            return false;
+        };
+        let Some(ny) = (self.vs.len() - 1).checked_mul(3).and_then(|v| v.checked_add(1)) else {
+            return false;
+        };
+        let Some(expected) = nx.checked_mul(ny) else {
+            return false;
+        };
+        self.points.len() == expected && self.points.iter().all(|p| p[0].is_finite() && p[1].is_finite())
     }
 
     fn locate(knots: &[f64], s: f64) -> (usize, f64) {
-        let n = knots.len() - 1;
-        let s = s.clamp(0.0, 1.0);
-        let mut i = 0;
-        while i + 1 < n && s > knots[i + 1] {
-            i += 1;
+        if knots.len() < 2 {
+            return (0, 0.0);
         }
+        let n = knots.len() - 1;
+        let s = if s.is_finite() { s.clamp(0.0, 1.0) } else { 0.0 };
+        let i = knots[1..n].partition_point(|&k| k < s);
         let w = knots[i + 1] - knots[i];
         (i, if w > 0.0 { (s - knots[i]) / w } else { 0.0 })
     }
 
     /// Evaluates the surface at normalised `(s, t)`.
     pub fn eval(&self, s: f64, t: f64) -> [f64; 2] {
+        if self.us.len() < 2 || self.vs.len() < 2 {
+            return [0.0; 2];
+        }
         let (pi, a) = Self::locate(&self.us, s);
         let (pj, b) = Self::locate(&self.vs, t);
         let (ba, bb) = (bernstein(a), bernstein(b));
@@ -513,7 +525,10 @@ impl BezierMesh {
         let mut out = [0.0; 2];
         for (j, wb) in bb.iter().enumerate() {
             for (i, wa) in ba.iter().enumerate() {
-                let p = self.points[(3 * pj + j) * nx + 3 * pi + i];
+                let idx = (3 * pj + j).saturating_mul(nx).saturating_add(3 * pi + i);
+                let Some(&p) = self.points.get(idx) else {
+                    return [0.0; 2];
+                };
                 let k = wa * wb;
                 out[0] += p[0] * k;
                 out[1] += p[1] * k;
@@ -949,8 +964,29 @@ impl Warp {
             WarpStyle::Custom => match &self.mesh {
                 None => true,
                 Some(m) => {
-                    let id = BezierMesh::fit(&|s, t| self.box_point(s, t), m.us.clone(), m.vs.clone());
-                    m.points.iter().zip(&id.points).all(|(a, b)| (a[0] - b[0]).abs() < 1e-9 && (a[1] - b[1]).abs() < 1e-9)
+                    if !m.is_valid() {
+                        return false;
+                    }
+                    let nx = m.nx();
+                    let (pc, pr) = (m.us.len() - 1, m.vs.len() - 1);
+                    for pj in 0..pr {
+                        for j in 0..4 {
+                            let v = m.vs[pj] + (m.vs[pj + 1] - m.vs[pj]) * j as f64 / 3.0;
+                            for pi in 0..pc {
+                                for i in 0..4 {
+                                    let u = m.us[pi] + (m.us[pi + 1] - m.us[pi]) * i as f64 / 3.0;
+                                    let expected = self.box_point(u, v);
+                                    let actual = m.points[(3 * pj + j) * nx + 3 * pi + i];
+                                    // Relative to the magnitude so the fit's rounding at large origins isn't rejected.
+                                    let tol = |e: f64| (1e-9 * e.abs().max(1.0)).min(1e-6);
+                                    if (actual[0] - expected[0]).abs() >= tol(expected[0]) || (actual[1] - expected[1]).abs() >= tol(expected[1]) {
+                                        return false;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                    true
                 }
             },
             _ => self.style_warp().is_none(),
@@ -1384,4 +1420,38 @@ mod tests {
         let ob = Warp::none(B).output_bounds();
         assert!((ob[0] - 10.0).abs() < 1e-9 && (ob[3] - 70.0).abs() < 1e-9);
     }
+
+    #[test]
+    fn is_identity_rejects_malformed_custom_mesh_without_allocating() {
+        let us: Vec<f64> = (0..=1000).map(|i| i as f64 / 1000.0).collect();
+        let vs: Vec<f64> = (0..=1000).map(|i| i as f64 / 1000.0).collect();
+        let malformed = BezierMesh { us, vs, points: Vec::new() };
+        assert!(!malformed.is_valid());
+        let w = Warp::custom(malformed.clone(), B);
+        assert!(!w.is_identity());
+        assert_eq!(malformed.eval(0.5, 0.5), [0.0, 0.0]);
+    }
+
+    #[test]
+    fn locate_binary_search_matches_linear_scan() {
+        let knots = [0.0, 0.1, 0.25, 0.5, 0.85, 1.0];
+        let linear = |s: f64| {
+            let n = knots.len() - 1;
+            let s = if s.is_finite() { s.clamp(0.0, 1.0) } else { 0.0 };
+            let mut i = 0;
+            while i + 1 < n && s > knots[i + 1] {
+                i += 1;
+            }
+            let w = knots[i + 1] - knots[i];
+            (i, if w > 0.0 { (s - knots[i]) / w } else { 0.0 })
+        };
+        for k in -10..=110 {
+            let s = k as f64 / 100.0;
+            assert_eq!(BezierMesh::locate(&knots, s), linear(s), "s={s}");
+        }
+        for &s in &knots {
+            assert_eq!(BezierMesh::locate(&knots, s), linear(s), "knot={s}");
+        }
+        assert_eq!(BezierMesh::locate(&knots, f64::NAN), linear(f64::NAN));
+    }
 }

crates/vector/src/lib.rs

@@ -156,11 +156,14 @@ fn tile_bands(r: Rect) -> Vec<Rect> {
     if r.is_empty() {
         return out;
     }
-    let mut y = r.y0.div_euclid(TILE_SIZE) * TILE_SIZE;
-    while y < r.y1 {
-        let band = Rect::new(r.x0, y.max(r.y0), r.x1, (y + TILE_SIZE).min(r.y1));
+    let step = i64::from(TILE_SIZE);
+    let (y0, y1) = (i64::from(r.y0), i64::from(r.y1));
+    let mut y = y0.div_euclid(step) * step;
+    while y < y1 {
+        let next = (y + step).min(y1);
+        let band = Rect::new(r.x0, y.max(y0) as i32, r.x1, next as i32);
         out.push(band);
-        y += TILE_SIZE;
+        y = next;
     }
     out
 }
@@ -170,9 +173,15 @@ fn tile_bands(r: Rect) -> Vec<Rect> {
 fn band_tiles(band: Rect, data: &[f32], channels: usize, keep: impl Fn(&[f32]) -> bool) -> Vec<(Rect, Vec<f32>)> {
     let w = band.width() as usize;
     let mut out = Vec::new();
-    let mut x = band.x0.div_euclid(TILE_SIZE) * TILE_SIZE;
-    while x < band.x1 {
-        let tr = Rect::new(x.max(band.x0), band.y0, (x + TILE_SIZE).min(band.x1), band.y1);
+    if band.is_empty() {
+        return out;
+    }
+    let step = i64::from(TILE_SIZE);
+    let (x0, x1) = (i64::from(band.x0), i64::from(band.x1));
+    let mut x = x0.div_euclid(step) * step;
+    while x < x1 {
+        let next = (x + step).min(x1);
+        let tr = Rect::new(x.max(x0) as i32, band.y0, next as i32, band.y1);
         let tw = tr.width() as usize;
         let mut buf = Vec::with_capacity(tw * tr.height() as usize * channels);
         let mut any = false;
@@ -187,7 +196,7 @@ fn band_tiles(band: Rect, data: &[f32], channels: usize, keep: impl Fn(&[f32]) -
         if any {
             out.push((tr, buf));
         }
-        x += TILE_SIZE;
+        x = next;
     }
     out
 }

crates/vector/src/raster.rs

@@ -19,6 +19,11 @@ use photocraft_geom::Rect;
 /// a 16-bit step (7.6e-6), so 8- and 16-bit output is unchanged by dropping it.
 const RESIDUE: f32 = 1e-6;
 
+/// Maximum sub-band cuts within a single pixel row before adjacent cuts are coalesced, and
+/// cumulative edge-evaluation budget per row to prevent quadratic work on degenerate inputs.
+const MAX_ROW_CUTS: usize = 256;
+const MAX_ROW_EDGE_WORK: usize = 262_144;
+
 #[derive(Clone, Copy, Debug)]
 struct Edge {
     x0: f64,
@@ -235,9 +240,36 @@ impl Rasterizer {
             }
             st.cuts.sort_by(f64::total_cmp);
             st.cuts.dedup();
+            if st.cuts.len() > MAX_ROW_CUTS {
+                let step = (st.cuts.len() - 1).div_ceil(MAX_ROW_CUTS - 1);
+                let last = st.cuts[st.cuts.len() - 1];
+                let mut thinned = Vec::with_capacity(MAX_ROW_CUTS);
+                for (idx, &c) in st.cuts.iter().enumerate() {
+                    if idx % step == 0 {
+                        thinned.push(c);
+                    }
+                }
+                if thinned.last().copied() != Some(last) {
+                    thinned.push(last);
+                }
+                st.cuts = thinned;
+            }
+            st.active.sort_unstable_by(|&a, &b| self.edges[a].y0.total_cmp(&self.edges[b].y0));
+            st.sub_active.clear();
+            st.work = 0;
+            let mut next_sub = 0usize;
             let cuts = std::mem::take(&mut st.cuts);
             for c in cuts.windows(2) {
-                self.band(c[0], c[1], rect.x0, &mut st, 0);
+                let (a, b) = (c[0], c[1]);
+                while next_sub < st.active.len() && self.edges[st.active[next_sub]].y0 < b {
+                    let ei = st.active[next_sub];
+                    if self.edges[ei].y1 > a {
+                        st.sub_active.push(ei);
+                    }
+                    next_sub += 1;
+                }
+                st.sub_active.retain(|&ei| self.edges[ei].y1 > a);
+                self.band(a, b, rect.x0, &mut st, 0);
             }
             st.cuts = cuts;
             let row_out = &mut out[row * w..(row + 1) * w];
@@ -259,13 +291,14 @@ impl Rasterizer {
         }
         st.band.clear();
         let mid = 0.5 * (a + b);
-        for &i in &st.active {
+        for &i in &st.sub_active {
             let e = &self.edges[i];
             if e.y0 <= mid && e.y1 >= mid {
                 let (xa, xb) = (e.x_at(a), e.x_at(b));
                 st.band.push(BandEdge { xa, xb, xm: 0.5 * (xa + xb), idx: i as u32 });
             }
         }
+        st.work = st.work.saturating_add(st.band.len());
         if st.band.is_empty() {
             if self.inverted {
                 add_full(&mut st.acc, b - a);
@@ -275,7 +308,7 @@ impl Rasterizer {
         st.band.sort_by(|p, q| p.xm.total_cmp(&q.xm));
         // An inversion of the order at the band's top or bottom means two edges cross inside
         // it: split there so the order is constant in every piece.
-        if depth < 24 {
+        if depth < 24 && st.work < MAX_ROW_EDGE_WORK {
             for k in 1..st.band.len() {
                 let (p, q) = (st.band[k - 1], st.band[k]);
                 let top = p.xa - q.xa;
@@ -330,16 +363,26 @@ struct BandEdge {
 
 struct RowState {
     active: Vec<usize>,
+    sub_active: Vec<usize>,
     cuts: Vec<f64>,
     band: Vec<BandEdge>,
     wind: Vec<i32>,
     /// Cell accumulator: `width + 1` cells (the extra one absorbs spill-over).
     acc: Vec<f32>,
+    work: usize,
 }
 
 impl RowState {
     fn new(comps: usize, width: usize) -> Self {
-        RowState { active: Vec::new(), cuts: Vec::new(), band: Vec::new(), wind: vec![0; comps.max(1)], acc: vec![0.0; width + 1] }
+        RowState {
+            active: Vec::new(),
+            sub_active: Vec::new(),
+            cuts: Vec::new(),
+            band: Vec::new(),
+            wind: vec![0; comps.max(1)],
+            acc: vec![0.0; width + 1],
+            work: 0,
+        }
     }
 }
 
@@ -519,4 +562,23 @@ mod tests {
 
         assert_eq!(r.pixel_bounds(), Some(Rect::new(1, 2, 5, 6)));
     }
+
+    #[test]
+    fn many_short_segments_in_single_row_complete_quickly() {
+        let mut r = Rasterizer::new(false);
+        let n = 4000;
+        let mut poly = Vec::with_capacity(n * 2);
+        for i in 0..n {
+            let y = 0.01 + 0.98 * (i as f64 / n as f64);
+            let x = if i % 2 == 0 { 1.0 } else { 3.0 };
+            poly.push((x, y));
+        }
+        poly.push((5.0, 0.99));
+        poly.push((5.0, 0.01));
+        r.add_component(&[poly], PathOp::Combine, FillRule::NonZero);
+
+        let cov = r.render(Rect::new(0, 0, 8, 2));
+        assert_eq!(cov.len(), 16);
+        assert!(cov.iter().all(|v| v.is_finite() && (0.0..=1.0).contains(v)));
+    }
 }

crates/vector/src/shapes.rs

@@ -142,7 +142,7 @@ pub fn fillet(prev: [f64; 2], v: [f64; 2], next: [f64; 2], r: f64) -> Option<Fil
     let (u2, l2) = unit(v, next)?;
     let cos = (u1[0] * u2[0] + u1[1] * u2[1]).clamp(-1.0, 1.0);
     let half = cos.acos() / 2.0;
-    if half < 1e-6 || half > std::f64::consts::FRAC_PI_2 - 1e-6 {
+    if !(1e-6..=std::f64::consts::FRAC_PI_2 - 1e-6).contains(&half) {
         return None;
     }
     let r = r.min(l1.min(l2) / 2.0 * half.tan());

crates/vector/src/tests.rs

@@ -570,3 +570,16 @@ fn thin_antialiased_edges_still_count_toward_bounds() {
     let c = r.render(Rect::new(0, 0, 102, 10));
     assert!(c[100] > 5e-5 && c[101] == 0.0, "{} {}", c[100], c[101]);
 }
+
+#[test]
+fn tile_bands_and_band_tiles_near_i32_max_do_not_overflow() {
+    let r = Rect::new(0, i32::MAX - 10, 10, i32::MAX);
+    let bands = tile_bands(r);
+    assert_eq!(bands, vec![r]);
+
+    let band = Rect::new(i32::MAX - 10, 0, i32::MAX, 2);
+    let data = vec![1.0f32; (band.width() * band.height()) as usize];
+    let tiles = band_tiles(band, &data, 1, |p| p[0] > 0.0);
+    assert_eq!(tiles.len(), 1);
+    assert_eq!(tiles[0].0, band);
+}

crates/vector/src/trace.rs

@@ -37,14 +37,26 @@ pub fn trace_mask(values: &[f32], rect: Rect, outside: f32, tolerance: f64) -> P
 
 /// Closed iso-contours at 0.5 through pixel centres.
 pub fn contours(values: &[f32], rect: Rect, outside: f32) -> Vec<Vec<P>> {
+    if rect.is_empty() {
+        return Vec::new();
+    }
     let w = rect.width() as i64;
     let h = rect.height() as i64;
-    let get = |x: i64, y: i64| -> f32 { if x < 0 || y < 0 || x >= w || y >= h { outside } else { values[(y * w + x) as usize] } };
+    let Some(need) = (w as usize).checked_mul(h as usize) else {
+        return Vec::new();
+    };
+    if values.len() < need {
+        return Vec::new();
+    }
+    let clean = |v: f32| if v.is_finite() { v } else { 0.0 };
+    let outside = clean(outside);
+    let get = |x: i64, y: i64| -> f32 { if x < 0 || y < 0 || x >= w || y >= h { outside } else { clean(values[(y * w + x) as usize]) } };
     // Edge ids: horizontal edge between (x,y)-(x+1,y) = (x, y, 0); vertical (x,y)-(x,y+1) = (x, y, 1).
     let lvl = 0.5f32;
     let point = |a: (i64, i64), b: (i64, i64)| -> P {
         let (va, vb) = (get(a.0, a.1), get(b.0, b.1));
-        let t = if (vb - va).abs() > 1e-9 { f64::from((lvl - va) / (vb - va)).clamp(0.0, 1.0) } else { 0.5 };
+        let raw_t = if (vb - va).abs() > 1e-9 { f64::from((lvl - va) / (vb - va)) } else { 0.5 };
+        let t = if raw_t.is_finite() { raw_t.clamp(0.0, 1.0) } else { 0.5 };
         let (ax, ay) = (a.0 as f64 + 0.5, a.1 as f64 + 0.5);
         let (bx, by) = (b.0 as f64 + 0.5, b.1 as f64 + 0.5);
         (f64::from(rect.x0) + ax + (bx - ax) * t, f64::from(rect.y0) + ay + (by - ay) * t)
@@ -156,14 +168,18 @@ fn norm(a: P) -> P {
 /// an open polyline, where the walk stops at the ends.
 fn walk(lp: &[P], i: usize, window: f64, fwd: bool, closed: bool) -> P {
     let n = lp.len();
+    let max_steps = n.min(256);
     let mut acc = 0.0;
     let mut cur = i;
-    for _ in 0..n {
+    for _ in 0..max_steps {
         if !closed && (if fwd { cur + 1 >= n } else { cur == 0 }) {
             break;
         }
         let nx = if fwd { (cur + 1) % n } else { (cur + n - 1) % n };
         let d = len(sub(lp[nx], lp[cur]));
+        if !d.is_finite() || !window.is_finite() {
+            break;
+        }
         if acc + d >= window {
             let t = (window - acc) / d.max(1e-12);
             return add(lp[cur], mul(sub(lp[nx], lp[cur]), t));
@@ -178,6 +194,9 @@ fn walk(lp: &[P], i: usize, window: f64, fwd: bool, closed: bool) -> P {
 /// ends of an open polyline are never corners.
 fn corners(lp: &[P], window: f64, closed: bool) -> Vec<usize> {
     let n = lp.len();
+    if n < 3 || !window.is_finite() || lp.iter().any(|p| !p.0.is_finite() || !p.1.is_finite()) {
+        return Vec::new();
+    }
     let cos_thresh = (60f64).to_radians().cos();
     let turn: Vec<f64> = (0..n)
         .map(|i| {
@@ -190,21 +209,33 @@ fn corners(lp: &[P], window: f64, closed: bool) -> Vec<usize> {
         })
         .collect();
     let mut out = Vec::new();
+    let max_neighbor_steps = n.min(256);
+    let mut work_budget = n.saturating_mul(64).max(1024);
     for i in 0..n {
         if turn[i] >= cos_thresh {
             continue;
         }
+        if work_budget == 0 {
+            break;
+        }
         // Local minimum of the cosine within ±window (ties: first index).
         let mut is_min = true;
         let mut acc = 0.0;
         let mut j = i;
-        while acc < window {
+        let mut steps = 0usize;
+        while acc < window && steps < max_neighbor_steps && work_budget > 0 {
             if !closed && j + 1 >= n {
                 break;
             }
             let nx = (j + 1) % n;
-            acc += len(sub(lp[nx], lp[j]));
+            let d = len(sub(lp[nx], lp[j]));
+            if !d.is_finite() {
+                break;
+            }
+            acc += d;
             j = nx;
+            steps += 1;
+            work_budget -= 1;
             if j == i {
                 break;
             }
@@ -215,13 +246,20 @@ fn corners(lp: &[P], window: f64, closed: bool) -> Vec<usize> {
         }
         acc = 0.0;
         j = i;
-        while is_min && acc < window {
+        steps = 0;
+        while is_min && acc < window && steps < max_neighbor_steps && work_budget > 0 {
             if !closed && j == 0 {
                 break;
             }
             let pv = (j + n - 1) % n;
-            acc += len(sub(lp[j], lp[pv]));
+            let d = len(sub(lp[j], lp[pv]));
+            if !d.is_finite() {
+                break;
+            }
+            acc += d;
             j = pv;
+            steps += 1;
+            work_budget -= 1;
             if j == i {
                 break;
             }
@@ -239,6 +277,9 @@ fn corners(lp: &[P], window: f64, closed: bool) -> Vec<usize> {
 /// Fits a closed loop with cubic Béziers; `window` is the corner detection window.
 fn fit_closed(lp: &[P], tol: f64, window: f64) -> Subpath {
     let n = lp.len();
+    if n < 3 || lp.iter().any(|p| !p.0.is_finite() || !p.1.is_finite()) {
+        return Subpath { closed: true, knots: Vec::new(), op: PathOp::Combine };
+    }
     let cs = corners(lp, window, true);
     let mut segs: Vec<[P; 4]> = Vec::new();
     let mut corner_at_start: Vec<bool> = Vec::new();
@@ -617,4 +658,39 @@ mod tests {
         let sp = fit_trace(&long, -5.0, false).unwrap();
         assert!(sp.knots.len() >= 2 && sp.knots.iter().all(|k| k.anchor.x.is_finite() && k.out_ctrl.y.is_finite()));
     }
+
+    #[test]
+    fn corners_have_bounded_work_for_many_tiny_segments() {
+        let points: Vec<P> = (0..20_000).map(|i| (i as f64 * 0.01, if i % 2 == 0 { 0.0 } else { 1e-5 })).collect();
+        let selected = corners(&points, 100.0, false);
+        assert!(selected.len() <= points.len());
+    }
+
+    #[test]
+    fn contours_reject_empty_rect_and_clean_non_finite_values() {
+        assert!(contours(&[], Rect::new(0, 0, 0, 1_000_000_000), 0.0).is_empty());
+        assert!(contours(&[], Rect::new(0, 0, 1_000_000_000, 0), 0.0).is_empty());
+        assert!(contours(&[1.0], Rect::new(0, 0, 10, 10), 0.0).is_empty());
+
+        let (w, h) = (6, 6);
+        let mut v = vec![0.0f32; w * h];
+        for y in 1..5 {
+            for x in 1..5 {
+                v[y * w + x] = f32::INFINITY;
+            }
+        }
+        v[2 * w + 2] = f32::NAN;
+        v[w + 1] = 1.0;
+        v[w + 2] = 1.0;
+        v[2 * w + 1] = 1.0;
+        let path = trace_mask(&v, Rect::new(0, 0, w as i32, h as i32), f32::NAN, 1.0);
+        assert!(path.subpaths.iter().any(|sp| !sp.knots.is_empty()));
+        for sp in &path.subpaths {
+            for k in &sp.knots {
+                assert!(k.anchor.x.is_finite() && k.anchor.y.is_finite());
+                assert!(k.in_ctrl.x.is_finite() && k.in_ctrl.y.is_finite());
+                assert!(k.out_ctrl.x.is_finite() && k.out_ctrl.y.is_finite());
+            }
+        }
+    }
 }