james-willis opened a new issue, #3118:
URL: https://github.com/apache/sedona/issues/3118

   ## Expected behavior
   
   Rasterizing a `LINESTRING`, and rasterizing a polygon with `allTouched = 
true`, should burn every pixel the geometry touches — this is what the 
[`RS_AsRaster` 
docs](https://sedona.apache.org/latest/api/sql/Raster-Operators/RS_AsRaster/) 
promise ("selects all pixels touched by the geometry") and what GDAL 
(`gdal_rasterize -at`, `rasterio.features.rasterize(all_touched=True)`) 
produces.
   
   ## Actual behavior
   
   Pixels that the geometry edge crosses over a short distance are silently not 
burned.
   
   A `LINESTRING` on a 6x6 unit raster:
   
   ```
   RS_AsRaster of LINESTRING (3.97 1.57, 0.31 3.24):
   
   Sedona                 GDAL / rasterio
   0 0 0 0 0 0            0 0 0 0 0 0
   0 0 0 0 0 0            0 0 0 0 0 0
   1 0 0 0 0 0            1 0 0 0 0 0
   0 1 1 0 0 0            1 1 1 1 0 0
   0 0 0 1 0 0            0 0 0 1 0 0
   0 0 0 0 0 0            0 0 0 0 0 0
   ```
   
   Sedona drops pixels (row 3, col 0) and (row 3, col 3), which the line 
clearly passes through.
   
   The same happens to polygon boundaries under `allTouched = true`. For 
`POLYGON ((3.7 1.28, 0.92 5.23, 4.26 4.15, 3.7 1.28))` on the same grid, Sedona 
omits pixel (row 3, col 4) — the boundary crosses it with a chord of 0.19 
pixel-widths — where GDAL burns it.
   
   The raster here is square with unit pixels, so this is unrelated to the 
non-square-pixel selection bug (a separate report).
   
   ## Steps to reproduce the problem
   
   ```sql
   -- LINESTRING: pixels (3,0) and (3,3) are missing
   SELECT RS_BandAsArray(RS_AsRaster(
     ST_GeomFromWKT('LINESTRING (3.97 1.57, 0.31 3.24)'),
     RS_MakeEmptyRaster(1, 'B', 6, 6, 0.0, 6.0, 1.0, -1.0, 0.0, 0.0, 0),
     'B', false, 1.0, 0.0, false), 1);
   
   -- Polygon with allTouched = true: pixel (3,4) is missing
   SELECT RS_BandAsArray(RS_AsRaster(
     ST_GeomFromWKT('POLYGON ((3.7 1.28, 0.92 5.23, 4.26 4.15, 3.7 1.28))'),
     RS_MakeEmptyRaster(1, 'B', 6, 6, 0.0, 6.0, 1.0, -1.0, 0.0, 0.0, 0),
     'B', true, 1.0, 0.0, false), 1);
   ```
   
   Comparison reference (Python):
   
   ```python
   import numpy as np, rasterio.features
   from affine import Affine
   from shapely.wkt import loads
   
   wkt = "LINESTRING (3.97 1.57, 0.31 3.24)"
   print(rasterio.features.rasterize(
       [(loads(wkt), 1)], out_shape=(6, 6), fill=0,
       transform=Affine(1, 0, 0, 0, -1, 6), all_touched=True, dtype="uint8"))
   ```
   
   ## Root cause
   
   `Rasterization.drawLineBresenham` is not a Bresenham line walk despite the 
name: it samples the segment at fixed 0.2-pixel steps and burns the pixel under 
each sample ([Rasterization.java#L204-L256 @ 
0299fa8](https://github.com/apache/sedona/blob/0299fa8a380a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java#L204-L256)):
   
   ```java
   drawLineBresenham(params, x0, y0, x1, y1, value, 0.2);
   ...
   int steps = (int) Math.ceil(distance / stepSize);   // stepSize = 0.2
   double stepX = dx / steps, stepY = dy / steps;
   for (int i = 0; i <= steps; i++) {
     params.writableRaster.setSample((int) Math.floor(x), (int) Math.floor(y), 
0, value);
     x += stepX; y += stepY;
   }
   ```
   
   A pixel is burned only if a sample lands in it, so any pixel the edge 
crosses with a chord shorter than the 0.2-pixel step can fall between two 
consecutive samples and be missed. This path is used for `LINESTRING` 
geometries (independent of `allTouched`) and for polygon-boundary pixels when 
`allTouched = true`, so it affects `RS_AsRaster`, `RS_Clip`, 
`RS_ZonalStats`/`RS_ZonalStatsAll`, and `RS_SetValues`. Introduced with the 
`allTouched` parameter in SEDONA-707 (#1788) and present on current master.
   
   The fix is exact grid traversal (integer DDA / Amanatides–Woo), visiting 
each cell the segment enters. Reducing the step size only narrows the miss 
window; it cannot close it.
   
   ## Settings
   
   Sedona version = 1.9.0 (also reproduced on current master)
   
   Apache Spark version = 4.0.0
   
   API type = SQL / Python
   
   Scala version = 2.13
   
   JRE version = 17
   
   Python version = 3.10
   


-- 
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.

To unsubscribe, e-mail: [email protected]

For queries about this service, please contact Infrastructure at:
[email protected]

Reply via email to