render: add perspective UV interpolation and derivatives (#1461)

Co-Authored-By: selune <[email protected]>
Co-Authored-By: Claude Sonnet 5.5 <[email protected]>

crates/render/src/lib.rs

@@ -14,6 +14,7 @@ mod raster;
 pub mod studio;
 pub mod surface;
 pub mod texture_input;
+pub mod texture_interpolation;
 
 pub use backend::{
     Availability, BackendDescriptor, BackendId, ExecutorDescriptor, ExecutorId, MAX_RENDER_ID_BYTES, RenderFeature, RenderFormat, RenderLimits,

crates/render/src/texture_interpolation.rs

@@ -0,0 +1,114 @@
+//! Perspective-correct unwrapped UV and analytic gradients from caller-owned barycentrics.
+//! No projection, clipping, raster coverage, instance transform or texture sampling.
+
+pub const UV_INTERPOLATION_POLICY: &str = "solidcraft.perspective-uv-differential/1";
+/// Matches the existing raster's near-edge admission; does not replace its top-left rule.
+pub const BARYCENTRIC_EDGE_TOLERANCE: f64 = 1e-9;
+const MAX_UV: f64 = 1_000_000.0;
+const MAX_DERIVATIVE: f64 = 1_000_000.0;
+const SUM_TOLERANCE: f64 = 1e-12;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum InterpolationError {
+    VertexUv,
+    ReciprocalW,
+    Barycentrics,
+    Derivatives,
+    Output,
+}
+impl std::fmt::Display for InterpolationError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.write_str(match self {
+            Self::VertexUv => "vertex UV must be finite and within one million repeats",
+            Self::ReciprocalW => "reciprocal clip-w must be finite and in [1e-12, 1e12]",
+            Self::Barycentrics => "barycentric weights exceed near-edge bounds or do not sum to one",
+            Self::Derivatives => "barycentric derivatives must be bounded, finite and sum to zero",
+            Self::Output => "interpolated UV or gradients exceed the numerical reference bounds",
+        })
+    }
+}
+impl std::error::Error for InterpolationError {}
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct UvDifferential {
+    pub uv: [f64; 2],
+    /// Unwrapped UV per actual shading-grid pixel step, not necessarily final-image pixels.
+    pub duv_dx: [f64; 2],
+    pub duv_dy: [f64; 2],
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct PerspectiveUv {
+    anchor: [f64; 2],
+    deltas: [[f64; 2]; 3],
+    reciprocal_w: [f64; 3],
+}
+impl PerspectiveUv {
+    /// Vertices must already have valid positive clip-w after caller-owned clipping.
+    /// Only a common reciprocal-w scale is removed; UVs are never wrapped or flipped.
+    pub fn new(uv: [[f64; 2]; 3], reciprocal_w: [f64; 3]) -> Result<Self, InterpolationError> {
+        if uv.iter().flatten().any(|v| !v.is_finite() || v.abs() > MAX_UV) {
+            return Err(InterpolationError::VertexUv);
+        }
+        if reciprocal_w.iter().any(|v| !v.is_finite() || !(1e-12..=1e12).contains(v)) {
+            return Err(InterpolationError::ReciprocalW);
+        }
+        let [anchor, _, _] = uv;
+        let [anchor_u, anchor_v] = anchor;
+        let scale = reciprocal_w.into_iter().fold(0.0, f64::max);
+        let deltas = uv.map(|[u, v]| [u - anchor_u, v - anchor_v]);
+        Ok(Self { anchor, deltas, reciprocal_w: reciprocal_w.map(|v| v / scale) })
+    }
+
+    /// Barycentrics and derivatives come from the same projected triangle/shading grid.
+    /// No coverage rule is applied here: weights beyond the explicit near-edge bound, malformed partitions and
+    /// unsupported gradients refuse instead of falling back to affine interpolation.
+    pub fn evaluate(&self, weights: [f64; 3], dweights_dx: [f64; 3], dweights_dy: [f64; 3]) -> Result<UvDifferential, InterpolationError> {
+        let admitted_weights = -BARYCENTRIC_EDGE_TOLERANCE..=1.0 + 2.0 * BARYCENTRIC_EDGE_TOLERANCE;
+        if weights.iter().any(|v| !v.is_finite() || !admitted_weights.contains(v)) || (weights.into_iter().sum::<f64>() - 1.0).abs() > SUM_TOLERANCE {
+            return Err(InterpolationError::Barycentrics);
+        }
+        for derivatives in [dweights_dx, dweights_dy] {
+            let scale = derivatives.into_iter().map(f64::abs).fold(1.0, f64::max);
+            if derivatives.iter().any(|v| !v.is_finite() || v.abs() > MAX_DERIVATIVE)
+                || derivatives.into_iter().sum::<f64>().abs() > SUM_TOLERANCE * scale
+            {
+                return Err(InterpolationError::Derivatives);
+            }
+        }
+        let denominator = weighted_sum(weights, self.reciprocal_w);
+        if !(denominator > 0.0) || !denominator.is_finite() {
+            return Err(InterpolationError::Output);
+        }
+        let dx_denominator = weighted_sum(dweights_dx, self.reciprocal_w);
+        let dy_denominator = weighted_sum(dweights_dy, self.reciprocal_w);
+        let evaluate_channel = |anchor, [a, b, c]: [f64; 3]| {
+            let [qa, qb, qc] = self.reciprocal_w;
+            let corrected = [a * qa, b * qb, c * qc];
+            let delta = weighted_sum(weights, corrected) / denominator;
+            // Avoid denominator squared and subtract the anchor before differentiating.
+            [
+                anchor + delta,
+                (weighted_sum(dweights_dx, corrected) - delta * dx_denominator) / denominator,
+                (weighted_sum(dweights_dy, corrected) - delta * dy_denominator) / denominator,
+            ]
+        };
+        let [[u0, v0], [u1, v1], [u2, v2]] = self.deltas;
+        let [anchor_u, anchor_v] = self.anchor;
+        let [u, du_dx, du_dy] = evaluate_channel(anchor_u, [u0, u1, u2]);
+        let [v, dv_dx, dv_dy] = evaluate_channel(anchor_v, [v0, v1, v2]);
+        let result = UvDifferential { uv: [u, v], duv_dx: [du_dx, dv_dx], duv_dy: [du_dy, dv_dy] };
+        if result.uv.iter().chain(&result.duv_dx).chain(&result.duv_dy).any(|v| !v.is_finite() || v.abs() > MAX_UV) {
+            return Err(InterpolationError::Output);
+        }
+        Ok(result)
+    }
+}
+
+fn weighted_sum(a: [f64; 3], b: [f64; 3]) -> f64 {
+    a.into_iter().zip(b).map(|(a, b)| a * b).sum()
+}
+
+#[cfg(test)]
+#[path = "texture_interpolation_tests.rs"]
+mod tests;

crates/render/src/texture_interpolation_tests.rs

@@ -0,0 +1,104 @@
+use super::*;
+const UV: [[f64; 2]; 3] = [[0., 0.], [1., 0.], [0., 1.]];
+const WEIGHTS: [f64; 3] = [0.5, 0.25, 0.25];
+const DX: [f64; 3] = [-1., 1., 0.];
+const DY: [f64; 3] = [-1., 0., 1.];
+
+fn close(actual: [f64; 2], expected: [f64; 2]) {
+    for (a, b) in actual.into_iter().zip(expected) {
+        assert!((a - b).abs() < 1e-10, "{actual:?} != {expected:?}");
+    }
+}
+
+#[test]
+fn orthographic_uv_and_derivatives_match_the_affine_plane() {
+    let result = PerspectiveUv::new(UV, [1.; 3]).unwrap().evaluate(WEIGHTS, DX, DY).unwrap();
+    close(result.uv, [0.25, 0.25]);
+    close(result.duv_dx, [1., 0.]);
+    close(result.duv_dy, [0., 1.]);
+}
+
+#[test]
+fn perspective_quotient_and_derivatives_match_closed_form() {
+    let result = PerspectiveUv::new(UV, [1., 2., 1.]).unwrap().evaluate(WEIGHTS, DX, DY).unwrap();
+    close(result.uv, [0.4, 0.2]);
+    close(result.duv_dx, [1.28, -0.16]);
+    close(result.duv_dy, [0., 0.8]);
+}
+
+#[test]
+fn common_clip_scale_and_vertex_order_preserve_the_result() {
+    for scale in [1e-6, 1., 1e6] {
+        for order in [[0, 1, 2], [2, 1, 0], [1, 2, 0]] {
+            let q = [scale, 2. * scale, scale];
+            let interpolation = PerspectiveUv::new(order.map(|i| UV[i]), order.map(|i| q[i])).unwrap();
+            let result = interpolation.evaluate(order.map(|i| WEIGHTS[i]), order.map(|i| DX[i]), order.map(|i| DY[i])).unwrap();
+            close(result.uv, [0.4, 0.2]);
+            close(result.duv_dx, [1.28, -0.16]);
+            close(result.duv_dy, [0., 0.8]);
+        }
+    }
+}
+
+#[test]
+fn constant_large_offset_has_exactly_zero_gradients() {
+    let interpolation = PerspectiveUv::new([[999_999., -999_999.]; 3], [1e-12, 1., 1e12]).unwrap();
+    let result = interpolation.evaluate(WEIGHTS, DX, DY).unwrap();
+    assert_eq!(result.uv, [999_999., -999_999.]);
+    assert_eq!(result.duv_dx, [0.; 2]);
+    assert_eq!(result.duv_dy, [0.; 2]);
+}
+
+#[test]
+fn vertices_edges_and_shading_grid_units_remain_explicit() {
+    let interpolation = PerspectiveUv::new(UV, [1., 2., 1.]).unwrap();
+    close(interpolation.evaluate([0., 1., 0.], DX, DY).unwrap().uv, [1., 0.]);
+    close(interpolation.evaluate([0.5, 0.5, 0.], DX, DY).unwrap().uv, [2. / 3., 0.]);
+    let fine = interpolation.evaluate(WEIGHTS, DX.map(|v| v / 2.), DY.map(|v| v / 4.)).unwrap();
+    close(fine.uv, [0.4, 0.2]);
+    close(fine.duv_dx, [0.64, -0.08]);
+    close(fine.duv_dy, [0., 0.2]);
+}
+
+#[test]
+fn invalid_vertices_and_reciprocal_w_are_refused() {
+    for value in [f64::NAN, f64::INFINITY, MAX_UV + 1., -MAX_UV - 1.] {
+        assert!(matches!(PerspectiveUv::new([[value, 0.]; 3], [1.; 3]), Err(InterpolationError::VertexUv)));
+    }
+    for value in [f64::NAN, f64::INFINITY, 0., -1., 1e-13, 1e13] {
+        assert!(matches!(PerspectiveUv::new(UV, [1., value, 1.]), Err(InterpolationError::ReciprocalW)));
+    }
+}
+
+#[test]
+fn malformed_partitions_and_gradients_are_refused() {
+    let interpolation = PerspectiveUv::new(UV, [1.; 3]).unwrap();
+    for weights in [[0.; 3], [1.; 3], [-0.1, 0.6, 0.5], [f64::NAN, 0., 1.], [f64::INFINITY, 0., 0.]] {
+        assert_eq!(interpolation.evaluate(weights, DX, DY), Err(InterpolationError::Barycentrics));
+    }
+    for derivatives in [[1.; 3], [f64::NAN, 0., 0.], [f64::INFINITY, 0., 0.], [-MAX_DERIVATIVE - 1., MAX_DERIVATIVE + 1., 0.]] {
+        assert_eq!(interpolation.evaluate(WEIGHTS, derivatives, DY), Err(InterpolationError::Derivatives));
+        assert_eq!(interpolation.evaluate(WEIGHTS, DX, derivatives), Err(InterpolationError::Derivatives));
+    }
+}
+
+#[test]
+fn unsupported_output_gradient_refuses_without_affine_fallback() {
+    let interpolation = PerspectiveUv::new(UV, [1e-12, 1e12, 1.]).unwrap();
+    assert_eq!(interpolation.evaluate([1., 0., 0.], DX, DY), Err(InterpolationError::Output));
+    assert!(PerspectiveUv::new(UV, [1.; 3]).unwrap().evaluate(WEIGHTS, DX, DY).is_ok());
+}
+
+#[test]
+fn near_edge_admission_keeps_exact_weights_and_refuses_beyond_bound_or_nonpositive_denominator() {
+    let e = BARYCENTRIC_EDGE_TOLERANCE;
+    let interpolation = PerspectiveUv::new(UV, [1.; 3]).unwrap();
+    let result = interpolation.evaluate([-e, 1. + 2. * e, -e], DX, DY).unwrap();
+    close(result.uv, [1. + 2. * e, -e]);
+    close(result.duv_dx, [1., 0.]);
+    close(result.duv_dy, [0., 1.]);
+    let outside = -e * 1.00001;
+    assert_eq!(interpolation.evaluate([outside, 0.5, 0.5 - outside], DX, DY), Err(InterpolationError::Barycentrics));
+    let perspective = PerspectiveUv::new(UV, [1e12, 1e-12, 1.]).unwrap();
+    assert_eq!(perspective.evaluate([-e, 1. + e, 0.], DX, DY), Err(InterpolationError::Output));
+}