From 714fde9591f1175e8c001f264c958b367db3d3ca Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Thu, 12 Oct 2023 11:51:22 -0700 Subject: [PATCH] Improve numerical robustness of conservative line rasterization Apply additional logic to make sure that first and last tiles of conservative line rasterization land on the right tile. Note: this fixes the assert and should avoid artifacts, but still doesn't render the poland svg correctly (there are black splotches). Fixes #376 and #378 --- shader/fine.wgsl | 17 +++++++++++++---- shader/path_count.wgsl | 15 ++++++++++++--- shader/path_tiling.wgsl | 15 ++++++++++++--- src/cpu_shader/path_count.rs | 13 +++++++++---- src/cpu_shader/path_tiling.rs | 16 ++++++++++++---- src/cpu_shader/util.rs | 15 +++++++++++++++ 6 files changed, 73 insertions(+), 18 deletions(-) diff --git a/shader/fine.wgsl b/shader/fine.wgsl index f41747d82..5d500ad04 100644 --- a/shader/fine.wgsl +++ b/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, local_id: vec2) -> array { let n_segs = fill.size_and_rule >> 1u; @@ -175,16 +179,21 @@ fn fill_path_ms(fill: CmdFill, wg_id: vec2, local_id: vec2) -> array= 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 diff --git a/shader/path_count.wgsl b/shader/path_count.wgsl index e34ff7a99..dabb6d193 100644 --- a/shader/path_count.wgsl +++ b/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]; diff --git a/shader/path_tiling.wgsl b/shader/path_tiling.wgsl index ace01ab97..fad5d725e 100644 --- a/shader/path_tiling.wgsl +++ b/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); diff --git a/src/cpu_shader/path_count.rs b/src/cpu_shader/path_count.rs index 2cee5b815..3cb8855ee 100644 --- a/src/cpu_shader/path_count.rs +++ b/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]; diff --git a/src/cpu_shader/path_tiling.rs b/src/cpu_shader/path_tiling.rs index 41549bb54..56bc2b47a 100644 --- a/src/cpu_shader/path_tiling.rs +++ b/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; diff --git a/src/cpu_shader/util.rs b/src/cpu_shader/util.rs index 2bb3279aa..eae6f8bd4 100644 --- a/src/cpu_shader/util.rs +++ b/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;