Improve numerical robustness of conservative line rasterization

shader/fine.wgsl

@@ -77,6 +77,10 @@ fn span(a: f32, b: f32) -> u32 {
 
 let SEG_SIZE = 5u;
 
+// See cpu_shaders/util.rs for explanation of these.
+let ONE_MINUS_ULP: f32 = 0.99999994;
+let ROBUST_EPSILON: f32 = 2e-7;
+
 // New multisampled algorithm.
 fn fill_path_ms(fill: CmdFill, wg_id: vec2<u32>, local_id: vec2<u32>) -> array<f32, PIXELS_PER_THREAD> {
     let n_segs = fill.size_and_rule >> 1u;
@@ -175,16 +179,21 @@ fn fill_path_ms(fill: CmdFill, wg_id: vec2<u32>, local_id: vec2<u32>) -> array<f
             // One alternative is to compute it in a separate dispatch.
             let dx = abs(xy1.x - xy0.x);
             let dy = xy1.y - xy0.y;
-            // TODO: apply numerical robustness and optimization
-            let dy_dxdy = dy / (dx + dy);
-            let a = dx / (dx + dy);
+            let idxdy = 1.0 / (dx + dy);
+            var a = dx * idxdy;
             let is_positive_slope = xy1.x >= xy0.x;
             let sign = select(-1.0, 1.0, is_positive_slope);
             let xt0 = floor(xy0.x * sign);
             let c = xy0.x * sign - xt0;
             let y0i = floor(xy0.y);
             let ytop = y0i + 1.0;
-            let b = dy_dxdy * c + a * (ytop - xy0.y);
+            let b = min((dy * c + dx * (ytop - xy0.y)) * idxdy, ONE_MINUS_ULP);
+            let count_x = span(xy0.x, xy1.x) - 1u;
+            let count = count_x + span(xy0.y, xy1.y);
+            let robust_err = floor(a * (f32(count) - 1.0) + b) - f32(count_x);
+            if robust_err != 0.0 {
+                a -= ROBUST_EPSILON * sign(robust_err);
+            }
             let x0i = i32(xt0 * sign + 0.5 * (sign - 1.0));
             // Use line equation to plot pixel coordinates
 

shader/path_count.wgsl

@@ -36,6 +36,10 @@ fn span(a: f32, b: f32) -> u32 {
     return u32(max(ceil(max(a, b)) - floor(min(a, b)), 1.0));
 }
 
+// See cpu_shaders/util.rs for explanation of these.
+let ONE_MINUS_ULP: f32 = 0.99999994;
+let ROBUST_EPSILON: f32 = 2e-7;
+
 // Note regarding clipping to bounding box:
 //
 // We have to do the backdrop bumps for all tiles to the left of the bbox.
@@ -57,7 +61,8 @@ fn main(
         let xy1 = select(line.p0, line.p1, is_down);
         let s0 = xy0 * TILE_SCALE;
         let s1 = xy1 * TILE_SCALE;
-        count = span(s0.x, s1.x) + span(s0.y, s1.y) - 1u;
+        let count_x = span(s0.x, s1.x) - 1u;
+        count = count_x + span(s0.y, s1.y);
         let line_ix = global_id.x;
 
         let dx = abs(s1.x - s0.x);
@@ -72,14 +77,18 @@ fn main(
             return;
         }
         let idxdy = 1.0 / (dx + dy);
-        let a = dx * idxdy;
+        var a = dx * idxdy;
         let is_positive_slope = s1.x >= s0.x;
         let sign = select(-1.0, 1.0, is_positive_slope);
         let xt0 = floor(s0.x * sign);
         let c = s0.x * sign - xt0;
         let y0 = floor(s0.y);
         let ytop = select(y0 + 1.0, ceil(s0.y), s0.y == s1.y);
-        let b = (dy * c + dx * (ytop - s0.y)) * idxdy;
+        let b = min((dy * c + dx * (ytop - s0.y)) * idxdy, ONE_MINUS_ULP);
+        let robust_err = floor(a * (f32(count) - 1.0) + b) - f32(count_x);
+        if robust_err != 0.0 {
+            a -= ROBUST_EPSILON * sign(robust_err);
+        }
         let x0 = xt0 * sign + select(-1.0, 0.0, is_positive_slope);
 
         let path = paths[line.path_ix];

shader/path_tiling.wgsl

@@ -29,6 +29,10 @@ fn span(a: f32, b: f32) -> u32 {
     return u32(max(ceil(max(a, b)) - floor(min(a, b)), 1.0));
 }
 
+// See cpu_shaders/util.rs for explanation of these.
+let ONE_MINUS_ULP: f32 = 0.99999994;
+let ROBUST_EPSILON: f32 = 2e-7;
+
 // One invocation for each tile that is to be written.
 // Total number of invocations = bump.seg_counts
 @compute @workgroup_size(256)
@@ -49,20 +53,25 @@ fn main(
         var xy1 = select(line.p0, line.p1, is_down);
         let s0 = xy0 * TILE_SCALE;
         let s1 = xy1 * TILE_SCALE;
-        let count = span(s0.x, s1.x) + span(s0.y, s1.y) - 1u;
+        let count_x = span(s0.x, s1.x) - 1u;
+        let count = count_x + span(s0.y, s1.y);
         let dx = abs(s1.x - s0.x);
         let dy = s1.y - s0.y;
         // Division by zero can't happen because zero-length lines
         // have already been discarded in the path_count stage.
         let idxdy = 1.0 / (dx + dy);
-        let a = dx * idxdy;
+        var a = dx * idxdy;
         let is_positive_slope = s1.x >= s0.x;
         let sign = select(-1.0, 1.0, is_positive_slope);
         let xt0 = floor(s0.x * sign);
         let c = s0.x * sign - xt0;
         let y0i = floor(s0.y);
         let ytop = select(y0i + 1.0, ceil(s0.y), s0.y == s1.y);
-        let b = (dy * c + dx * (ytop - s0.y)) * idxdy;
+        let b = min((dy * c + dx * (ytop - s0.y)) * idxdy, ONE_MINUS_ULP);
+        let robust_err = floor(a * (f32(count) - 1.0) + b) - f32(count_x);
+        if robust_err != 0.0 {
+            a -= ROBUST_EPSILON * sign(robust_err);
+        }
         let x0i = i32(xt0 * sign + 0.5 * (sign - 1.0));
         let z = floor(a * f32(seg_within_line) + b);
         let x = x0i + i32(sign * z);

src/cpu_shader/path_count.rs

@@ -5,7 +5,7 @@ use vello_encoding::{BumpAllocators, LineSoup, Path, SegmentCount, Tile};
 
 use crate::cpu_dispatch::CpuBinding;
 
-use super::util::{span, Vec2};
+use super::util::{span, Vec2, ONE_MINUS_ULP, ROBUST_EPSILON};
 
 const TILE_SCALE: f32 = 1.0 / 16.0;
 
@@ -24,7 +24,8 @@ fn path_count_main(
         let (xy0, xy1) = if is_down { (p0, p1) } else { (p1, p0) };
         let s0 = xy0 * TILE_SCALE;
         let s1 = xy1 * TILE_SCALE;
-        let count = span(s0.x, s1.x) + span(s0.y, s1.y) - 1;
+        let count_x = span(s0.x, s1.x) - 1;
+        let count = count_x + span(s0.y, s1.y);
 
         let dx = (s1.x - s0.x).abs();
         let dy = s1.y - s0.y;
@@ -35,14 +36,18 @@ fn path_count_main(
             continue;
         }
         let idxdy = 1.0 / (dx + dy);
-        let a = dx * idxdy;
+        let mut a = dx * idxdy;
         let is_positive_slope = s1.x >= s0.x;
         let sign = if is_positive_slope { 1.0 } else { -1.0 };
         let xt0 = (s0.x * sign).floor();
         let c = s0.x * sign - xt0;
         let y0 = s0.y.floor();
         let ytop = if s0.y == s1.y { s0.y.ceil() } else { y0 + 1.0 };
-        let b = (dy * c + dx * (ytop - s0.y)) * idxdy;
+        let b = ((dy * c + dx * (ytop - s0.y)) * idxdy).min(ONE_MINUS_ULP);
+        let robust_err = (a * (count as f32 - 1.0) + b).floor() - count_x as f32;
+        if robust_err != 0.0 {
+            a -= ROBUST_EPSILON.copysign(robust_err);
+        }
         let x0 = xt0 * sign + if is_positive_slope { 0.0 } else { -1.0 };
 
         let path = paths[line.path_ix as usize];

src/cpu_shader/path_tiling.rs

@@ -3,7 +3,10 @@
 
 use vello_encoding::{BumpAllocators, LineSoup, Path, PathSegment, SegmentCount, Tile};
 
-use crate::cpu_dispatch::CpuBinding;
+use crate::{
+    cpu_dispatch::CpuBinding,
+    cpu_shader::util::{ONE_MINUS_ULP, ROBUST_EPSILON},
+};
 
 use super::util::{span, Vec2};
 
@@ -33,19 +36,24 @@ fn path_tiling_main(
         let (mut xy0, mut xy1) = if is_down { (p0, p1) } else { (p1, p0) };
         let s0 = xy0 * TILE_SCALE;
         let s1 = xy1 * TILE_SCALE;
-        let count = span(s0.x, s1.x) + span(s0.y, s1.y) - 1;
+        let count_x = span(s0.x, s1.x) - 1;
+        let count = count_x + span(s0.y, s1.y);
 
         let dx = (s1.x - s0.x).abs();
         let dy = s1.y - s0.y;
         let idxdy = 1.0 / (dx + dy);
-        let a = dx * idxdy;
+        let mut a = dx * idxdy;
         let is_positive_slope = s1.x >= s0.x;
         let sign = if is_positive_slope { 1.0 } else { -1.0 };
         let xt0 = (s0.x * sign).floor();
         let c = s0.x * sign - xt0;
         let y0 = s0.y.floor();
         let ytop = if s0.y == s1.y { s0.y.ceil() } else { y0 + 1.0 };
-        let b = (dy * c + dx * (ytop - s0.y)) * idxdy;
+        let b = ((dy * c + dx * (ytop - s0.y)) * idxdy).min(ONE_MINUS_ULP);
+        let robust_err = (a * (count as f32 - 1.0) + b).floor() - count_x as f32;
+        if robust_err != 0.0 {
+            a -= ROBUST_EPSILON.copysign(robust_err);
+        }
         let x0 = xt0 * sign + if is_positive_slope { 0.0 } else { -1.0 };
         let z = (a * seg_within_line as f32 + b).floor();
         let x = x0 as i32 + (sign * z) as i32;

src/cpu_shader/util.rs

@@ -111,3 +111,18 @@ pub fn read_draw_tag_from_scene(config: &ConfigUniform, scene: &[u32], ix: u32)
         DRAWTAG_NOP
     }
 }
+
+/// The largest floating point value strictly less than 1.
+///
+/// This value is used to limit the value of b so that its floor is strictly less
+/// than 1. That guarantees that floor(a * i + b) == 0 for i == 0, which lands on
+/// the correct first tile.
+pub const ONE_MINUS_ULP: f32 = 0.99999994;
+
+/// An epsilon to be applied in path numerical robustness.
+///
+/// When floor(a * (n - 1) + b) does not match the expected value (the width in
+/// grid cells minus one), this delta is applied to a to push it in the correct
+/// direction. The theory is that a is not off by more than a few ulp, and it's
+/// always in the range of 0..1.
+pub const ROBUST_EPSILON: f32 = 2e-7;