eliminate inconsistent line intersections from path_coarse.comp

The finite precision of floating point computations can lead the coarse
renderer into inconsistent tile intersections, which implies impossible line
segments such as lines with gaps or double intersections. The winding number
algorithm is sensitive to these errors which show up as incorrectly filled
paths.

This change forces all intersections to be consistent.
First, the floating point top edge intersection test is removed; top edge intersections are
completely determined by left edge intersections.
Then, left edge intersections are inserted from the tile with the last top edge
intersection. The next top edge is then fixed to be the last tile with a left
edge intersection.

More details in the patch comments.

Fixes #23

Signed-off-by: Elias Naur <mail@eliasnaur.com>

piet-gpu/shader/path_coarse.comp

@@ -171,9 +171,9 @@ void main() {
                 float a = (p0.x - (p0.y - 0.5 * float(TILE_HEIGHT_PX)) * b) * SX;
 
                 int x0 = int(floor(xmin * SX));
-                int x1 = int(ceil(xmax * SX));
+                int x1 = int(floor(xmax * SX) + 1);
                 int y0 = int(floor(ymin * SY));
-                int y1 = int(ceil(ymax * SY));
+                int y1 = int(floor(ymax * SY) + 1);
 
                 x0 = clamp(x0, bbox.x, bbox.z);
                 y0 = clamp(y0, bbox.y, bbox.w);
@@ -187,19 +187,30 @@ void main() {
                 // Consider using subgroups to aggregate atomic add.
                 uint tile_offset = atomicAdd(alloc, n_tile_alloc * TileSeg_size);
                 TileSeg tile_seg;
+
+                int xray = int(floor(p0.x*SX));
+                int last_xray = int(floor(p1.x*SX));
+                if (p0.y > p1.y) {
+                    int tmp = xray;
+                    xray = last_xray;
+                    last_xray = tmp;
+                }
                 for (int y = y0; y < y1; y++) {
-                    float tile_y0 = float(y * TILE_HEIGHT_PX);
+                    int xbackdrop = max(xray + 1, bbox.x);
-                    if (tag == PathSeg_FillCubic && min(p0.y, p1.y) <= tile_y0) {
+                    if (tag == PathSeg_FillCubic && y > y0 && xbackdrop < bbox.z) {
-                        int xray = max(int(ceil(xc - 0.5 * b)), bbox.x);
-                        if (xray < bbox.z) {
                         int backdrop = p1.y < p0.y ? 1 : -1;
-                            TileRef tile_ref = Tile_index(path.tiles, uint(base + xray));
+                        TileRef tile_ref = Tile_index(path.tiles, uint(base + xbackdrop));
                         uint tile_el = tile_ref.offset >> 2;
                         atomicAdd(tile[tile_el + 1], backdrop);
                     }
-                    }
+
                     int xx0 = clamp(int(floor(xc - c)), x0, x1);
                     int xx1 = clamp(int(ceil(xc + c)), x0, x1);
+                    xx1 = max(xx1, xray + 1);
+
+                    // next_xray is the xray for the next scanline; it is derived
+                    // by left edge intersections computed below.
+                    int next_xray = xray;
                     for (int x = xx0; x < xx1; x++) {
                         float tile_x0 = float(x * TILE_WIDTH_PX);
                         TileRef tile_ref = Tile_index(path.tiles, uint(base + x));
@@ -209,8 +220,10 @@ void main() {
                         tile_seg.vector = p1 - p0;
                         float y_edge = 0.0;
                         if (tag == PathSeg_FillCubic) {
+                            float tile_y0 = float(y * TILE_HEIGHT_PX);
                             y_edge = mix(p0.y, p1.y, (tile_x0 - p0.x) / dx);
                             if (min(p0.x, p1.x) < tile_x0 && y_edge >= tile_y0 && y_edge < tile_y0 + TILE_HEIGHT_PX) {
+                                // Left edge intersection.
                                 vec2 p = vec2(tile_x0, y_edge);
                                 if (p0.x > p1.x) {
                                     tile_seg.vector = p - p0;
@@ -218,7 +231,25 @@ void main() {
                                     tile_seg.origin = p;
                                     tile_seg.vector = p1 - p;
                                 }
-                            } else {
+                                // kernel4 uses sign(vector.x) for the sign of the intersection backdrop.
+                                // Nudge zeroes towards the intended sign.
+                                if (tile_seg.vector.x == 0) {
+                                    tile_seg.vector.x += sign(p1.x - p0.x)*1e-9;
+                                }
+                                // Move next_xray consistently with previous intersections.
+                                if (x > next_xray && next_xray >= xray) {
+                                    next_xray = x;
+                                } else if (x <= next_xray && next_xray <= xray) {
+                                    next_xray = x - 1;
+                                }
+                            }
+                            // Force last xray on the last scanline for consistency with later
+                            // line segments.
+                            if (y == y1 - 1) {
+                                next_xray = last_xray;
+                            }
+                            // Drop inconsistent intersections.
+                            if (x <= min(xray, next_xray) || max(xray, next_xray) < x) {
                                 y_edge = 1e9;
                             }
                         }
@@ -229,6 +260,7 @@ void main() {
                     }
                     xc += b;
                     base += stride;
+                    xray = next_xray;
                 }
 
                 n_out += 1;