Kontinuation commented on code in PR #1788: URL: https://github.com/apache/sedona/pull/1788#discussion_r1960925114
########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { Review Comment: Why not use a traditional Bresenham's line algorithm using integer arithmetics? Using a `stepSize = 0.2` here introduces lots of repeated setSample operations. ########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { + + double dx = x1 - x0; + double dy = y1 - y0; + + // Compute the number of steps based on the larger of dx or dy + double distance = Math.sqrt(dx * dx + dy * dy); + int steps = (int) Math.ceil(distance / stepSize); + + // Calculate the step increment for each axis + double stepX = dx / steps; + double stepY = dy / steps; + + // Start stepping through the line + double x = x0; + double y = y0; + + for (int i = 0; i <= steps; i++) { + int rasterX = (int) (Math.floor(x)); + int rasterY = (int) (Math.floor(y)); + + // Only write if within raster bounds + if (rasterX >= 0 + && rasterX < params.writableRaster.getWidth() + && rasterY >= 0 + && rasterY < params.writableRaster.getHeight()) { + params.writableRaster.setSample(rasterX, rasterY, 0, value); + } + + // Increment by fractional steps + x += stepX; + y += stepY; + } + } + + private static LineSegment clipSegmentToRasterBounds( + Coordinate p1, Coordinate p2, Envelope2D geomExtent) { + double minX = geomExtent.getMinX(); + double maxX = geomExtent.getMaxX(); + double minY = geomExtent.getMinY(); + double maxY = geomExtent.getMaxY(); + + double x1 = p1.x, y1 = p1.y; + double x2 = p2.x, y2 = p2.y; + + boolean p1Inside = isInsideBounds(x1, y1, minX, maxX, minY, maxY); + boolean p2Inside = isInsideBounds(x2, y2, minX, maxX, minY, maxY); + + if (p1Inside && p2Inside) { + // Both points inside: no clipping needed + return new LineSegment(p1, p2); + } + + if (!p1Inside && !p2Inside) { + // Both points outside: no clipping needed + return null; + } + + double dx = x2 - x1; + double dy = y2 - y1; + + // Clip using parametric line equation + double[] tValues = {0, 1}; // Stores valid segment proportions + + // Clip against minX and maxX + if (dx != 0) { + double tMin = (minX - x1) / dx; + double tMax = (maxX - x1) / dx; + if (dx < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // Clip against minY and maxY + if (dy != 0) { + double tMin = (minY - y1) / dy; + double tMax = (maxY - y1) / dy; + if (dy < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // If tValues are invalid (segment is outside), return null + if (tValues[0] > tValues[1]) { + return null; // No valid clipped segment + } + + // Compute new clipped endpoints + Coordinate newP1 = new Coordinate(x1 + tValues[0] * dx, y1 + tValues[0] * dy); + Coordinate newP2 = new Coordinate(x1 + tValues[1] * dx, y1 + tValues[1] * dy); + + return new LineSegment(newP1, newP2); + } + + // Helper function to check if a point is inside the bounding box + private static boolean isInsideBounds( + double x, double y, double minX, double maxX, double minY, double maxY) { + return x >= minX && x <= maxX && y >= minY && y <= maxY; + } + + private static Envelope2D rasterizeGeomExtent( + Geometry geom, GridCoverage2D raster, double[] metadata, boolean allTouched) { + + if (Objects.equals(geom.getGeometryType(), "MultiLineString")) { + allTouched = true; + } + if (Objects.equals(geom.getGeometryType(), "MultiPoint")) { + allTouched = true; + } + + Envelope2D rasterExtent = + JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D()); + + // Compute the aligned min/max values + double alignedMinX = + (Math.floor((rasterExtent.getMinX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMinY = + (Math.floor((rasterExtent.getMinY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + double alignedMaxX = + (Math.ceil((rasterExtent.getMaxX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMaxY = + (Math.ceil((rasterExtent.getMaxY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + + // For points and LineStrings at intersection of 2 or more pixels, + // extend search grid by 1 pixel in each direction + if (alignedMaxX == alignedMinX) { + alignedMinX -= metadata[4]; + alignedMaxX += metadata[4]; + } + if (alignedMaxY == alignedMinY) { + alignedMinY += metadata[5]; + alignedMaxY -= metadata[5]; + } + + // Get the extent of the original raster + double originalMinX = raster.getEnvelope().getMinimum(0); + double originalMinY = raster.getEnvelope().getMinimum(1); + double originalMaxX = raster.getEnvelope().getMaximum(0); + double originalMaxY = raster.getEnvelope().getMaximum(1); + + // Extend the aligned extent by 1 pixel if allTouched is true and if any geometry extent meets + // the aligned extent + if (allTouched) { + alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0; + alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0; + alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0; + alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0; + } + + alignedMinX = Math.max(alignedMinX, originalMinX); + alignedMinY = Math.max(alignedMinY, originalMinY); + alignedMaxX = Math.min(alignedMaxX, originalMaxX); + alignedMaxY = Math.min(alignedMaxY, originalMaxY); + + // Create the aligned raster extent + Envelope2D alignedRasterExtent = + new Envelope2D( + rasterExtent.getCoordinateReferenceSystem(), + alignedMinX, + alignedMinY, + alignedMaxX - alignedMinX, + alignedMaxY - alignedMinY); + + return alignedRasterExtent; + } + + private static RasterizationParams calculateRasterizationParams( + GridCoverage2D raster, + boolean useGeometryExtent, + double[] metadata, + Envelope2D geomExtent, + String pixelType) { + + double upperLeftX = 0; + double upperLeftY = 0; + if (useGeometryExtent) { + upperLeftX = geomExtent.getMinX(); + upperLeftY = geomExtent.getMaxY(); + } else { + upperLeftX = metadata[0]; + upperLeftY = metadata[1]; + } + + WritableRaster writableRaster; + if (useGeometryExtent) { + int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4])); + int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5])); + + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), geomExtentWidth, geomExtentHeight, 1, null); + } else { + int rasterWidth = (int) raster.getGridGeometry().getGridRange2D().getWidth(); + int rasterHeight = (int) raster.getGridGeometry().getGridRange2D().getHeight(); + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null); + } + + return new RasterizationParams( + writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY); + } + + private static void validateRasterMetadata(double[] rasterMetadata) { + if (rasterMetadata[4] < 0) { + throw new IllegalArgumentException("ScaleX cannot be negative"); + } + if (rasterMetadata[5] > 0) { + throw new IllegalArgumentException("ScaleY must be negative"); + } + if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) { + throw new IllegalArgumentException("SkewX and SkewY must be zero"); + } + } + + // New condensed Rasterization parameters + private static class RasterizationParams { + WritableRaster writableRaster; + String pixelType; + double scaleX; + double scaleY; + double upperLeftX; + double upperLeftY; + + RasterizationParams( + WritableRaster writableRaster, + String pixelType, + double scaleX, + double scaleY, + double upperLeftX, + double upperLeftY) { + this.writableRaster = writableRaster; + this.pixelType = pixelType; + this.scaleX = scaleX; + this.scaleY = scaleY; + this.upperLeftX = upperLeftX; + this.upperLeftY = upperLeftY; + } + } + + public static void rasterizePolygon( + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) { + if (!(geom instanceof Polygon)) { + throw new IllegalArgumentException("Only Polygon geometry is supported"); + } + + Geometry clippedGeom = + Functions.intersection(JTS.toGeometry((BoundingBox) geomExtent), Functions.buffer(geom, 0)); + + Polygon polygon = (Polygon) clippedGeom; + + // Compute scanline X-intercepts + Map<Double, TreeSet<Double>> scanlineIntersections = + computeScanlineIntersections(polygon, params, value, geomExtent, allTouched); + + // Process intersections to get startXs and endXs for each scanline + Map<Integer, List<int[]>> scanlineFillRanges = computeFillRanges(scanlineIntersections); + + // Burn values between startX and endX pairs + fillPolygon(scanlineFillRanges, params, value); + } + + /** Computes scanline intersections by iterating over polygon edges. */ + private static Map<Double, TreeSet<Double>> computeScanlineIntersections( + Polygon polygon, + RasterizationParams params, + double value, + Envelope2D geomExtent, + boolean allTouched) { + + Map<Double, TreeSet<Double>> scanlineIntersections = new HashMap<>(); + List<LineString> allRings = new ArrayList<>(); + allRings.add(polygon.getExteriorRing()); + for (int i = 0; i < polygon.getNumInteriorRing(); i++) { + allRings.add(polygon.getInteriorRingN(i)); + } + + for (LineString ring : allRings) { + Coordinate[] coords = ring.getCoordinates(); + int numPoints = coords.length; + + if (allTouched) { + rasterizeLineString(ring, params, value, geomExtent); + } + + for (int i = 0; i < numPoints - 1; i++) { + Coordinate worldP1 = coords[i]; + Coordinate worldP2 = coords[i + 1]; + + // Ensure worldP1.y ≤ worldP2.y + if (worldP1.y > worldP2.y) { + Coordinate temp = worldP1; + worldP1 = worldP2; + worldP2 = temp; + } + + if (worldP1.y == worldP2.y) { + continue; // Skip horizontal edges + } Review Comment: Should we burn the pixels in between `worldP1.x` and `worldP2.x` instead of simply ignoring the horizontal edge? ########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { + + double dx = x1 - x0; + double dy = y1 - y0; + + // Compute the number of steps based on the larger of dx or dy + double distance = Math.sqrt(dx * dx + dy * dy); + int steps = (int) Math.ceil(distance / stepSize); + + // Calculate the step increment for each axis + double stepX = dx / steps; + double stepY = dy / steps; + + // Start stepping through the line + double x = x0; + double y = y0; + + for (int i = 0; i <= steps; i++) { + int rasterX = (int) (Math.floor(x)); + int rasterY = (int) (Math.floor(y)); + + // Only write if within raster bounds + if (rasterX >= 0 + && rasterX < params.writableRaster.getWidth() + && rasterY >= 0 + && rasterY < params.writableRaster.getHeight()) { + params.writableRaster.setSample(rasterX, rasterY, 0, value); + } + + // Increment by fractional steps + x += stepX; + y += stepY; + } + } + + private static LineSegment clipSegmentToRasterBounds( + Coordinate p1, Coordinate p2, Envelope2D geomExtent) { + double minX = geomExtent.getMinX(); + double maxX = geomExtent.getMaxX(); + double minY = geomExtent.getMinY(); + double maxY = geomExtent.getMaxY(); + + double x1 = p1.x, y1 = p1.y; + double x2 = p2.x, y2 = p2.y; + + boolean p1Inside = isInsideBounds(x1, y1, minX, maxX, minY, maxY); + boolean p2Inside = isInsideBounds(x2, y2, minX, maxX, minY, maxY); + + if (p1Inside && p2Inside) { + // Both points inside: no clipping needed + return new LineSegment(p1, p2); + } + + if (!p1Inside && !p2Inside) { + // Both points outside: no clipping needed + return null; + } + + double dx = x2 - x1; + double dy = y2 - y1; + + // Clip using parametric line equation + double[] tValues = {0, 1}; // Stores valid segment proportions + + // Clip against minX and maxX + if (dx != 0) { + double tMin = (minX - x1) / dx; + double tMax = (maxX - x1) / dx; + if (dx < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // Clip against minY and maxY + if (dy != 0) { + double tMin = (minY - y1) / dy; + double tMax = (maxY - y1) / dy; + if (dy < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // If tValues are invalid (segment is outside), return null + if (tValues[0] > tValues[1]) { + return null; // No valid clipped segment + } + + // Compute new clipped endpoints + Coordinate newP1 = new Coordinate(x1 + tValues[0] * dx, y1 + tValues[0] * dy); + Coordinate newP2 = new Coordinate(x1 + tValues[1] * dx, y1 + tValues[1] * dy); + + return new LineSegment(newP1, newP2); + } + + // Helper function to check if a point is inside the bounding box + private static boolean isInsideBounds( + double x, double y, double minX, double maxX, double minY, double maxY) { + return x >= minX && x <= maxX && y >= minY && y <= maxY; + } + + private static Envelope2D rasterizeGeomExtent( + Geometry geom, GridCoverage2D raster, double[] metadata, boolean allTouched) { + + if (Objects.equals(geom.getGeometryType(), "MultiLineString")) { + allTouched = true; + } + if (Objects.equals(geom.getGeometryType(), "MultiPoint")) { + allTouched = true; + } + + Envelope2D rasterExtent = + JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D()); + + // Compute the aligned min/max values + double alignedMinX = + (Math.floor((rasterExtent.getMinX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMinY = + (Math.floor((rasterExtent.getMinY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + double alignedMaxX = + (Math.ceil((rasterExtent.getMaxX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMaxY = + (Math.ceil((rasterExtent.getMaxY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + + // For points and LineStrings at intersection of 2 or more pixels, + // extend search grid by 1 pixel in each direction + if (alignedMaxX == alignedMinX) { + alignedMinX -= metadata[4]; + alignedMaxX += metadata[4]; + } + if (alignedMaxY == alignedMinY) { + alignedMinY += metadata[5]; + alignedMaxY -= metadata[5]; + } + + // Get the extent of the original raster + double originalMinX = raster.getEnvelope().getMinimum(0); + double originalMinY = raster.getEnvelope().getMinimum(1); + double originalMaxX = raster.getEnvelope().getMaximum(0); + double originalMaxY = raster.getEnvelope().getMaximum(1); + + // Extend the aligned extent by 1 pixel if allTouched is true and if any geometry extent meets + // the aligned extent + if (allTouched) { + alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0; + alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0; + alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0; + alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0; + } + + alignedMinX = Math.max(alignedMinX, originalMinX); + alignedMinY = Math.max(alignedMinY, originalMinY); + alignedMaxX = Math.min(alignedMaxX, originalMaxX); + alignedMaxY = Math.min(alignedMaxY, originalMaxY); + + // Create the aligned raster extent + Envelope2D alignedRasterExtent = + new Envelope2D( + rasterExtent.getCoordinateReferenceSystem(), + alignedMinX, + alignedMinY, + alignedMaxX - alignedMinX, + alignedMaxY - alignedMinY); + + return alignedRasterExtent; + } + + private static RasterizationParams calculateRasterizationParams( + GridCoverage2D raster, + boolean useGeometryExtent, + double[] metadata, + Envelope2D geomExtent, + String pixelType) { + + double upperLeftX = 0; + double upperLeftY = 0; + if (useGeometryExtent) { + upperLeftX = geomExtent.getMinX(); + upperLeftY = geomExtent.getMaxY(); + } else { + upperLeftX = metadata[0]; + upperLeftY = metadata[1]; + } + + WritableRaster writableRaster; + if (useGeometryExtent) { + int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4])); + int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5])); + + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), geomExtentWidth, geomExtentHeight, 1, null); + } else { + int rasterWidth = (int) raster.getGridGeometry().getGridRange2D().getWidth(); + int rasterHeight = (int) raster.getGridGeometry().getGridRange2D().getHeight(); + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null); + } + + return new RasterizationParams( + writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY); + } + + private static void validateRasterMetadata(double[] rasterMetadata) { + if (rasterMetadata[4] < 0) { + throw new IllegalArgumentException("ScaleX cannot be negative"); + } + if (rasterMetadata[5] > 0) { + throw new IllegalArgumentException("ScaleY must be negative"); + } + if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) { + throw new IllegalArgumentException("SkewX and SkewY must be zero"); + } + } + + // New condensed Rasterization parameters + private static class RasterizationParams { + WritableRaster writableRaster; + String pixelType; + double scaleX; + double scaleY; + double upperLeftX; + double upperLeftY; + + RasterizationParams( + WritableRaster writableRaster, + String pixelType, + double scaleX, + double scaleY, + double upperLeftX, + double upperLeftY) { + this.writableRaster = writableRaster; + this.pixelType = pixelType; + this.scaleX = scaleX; + this.scaleY = scaleY; + this.upperLeftX = upperLeftX; + this.upperLeftY = upperLeftY; + } + } + + public static void rasterizePolygon( + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) { + if (!(geom instanceof Polygon)) { + throw new IllegalArgumentException("Only Polygon geometry is supported"); + } + + Geometry clippedGeom = + Functions.intersection(JTS.toGeometry((BoundingBox) geomExtent), Functions.buffer(geom, 0)); + + Polygon polygon = (Polygon) clippedGeom; + + // Compute scanline X-intercepts + Map<Double, TreeSet<Double>> scanlineIntersections = + computeScanlineIntersections(polygon, params, value, geomExtent, allTouched); + + // Process intersections to get startXs and endXs for each scanline + Map<Integer, List<int[]>> scanlineFillRanges = computeFillRanges(scanlineIntersections); + + // Burn values between startX and endX pairs + fillPolygon(scanlineFillRanges, params, value); + } + + /** Computes scanline intersections by iterating over polygon edges. */ + private static Map<Double, TreeSet<Double>> computeScanlineIntersections( + Polygon polygon, + RasterizationParams params, + double value, + Envelope2D geomExtent, + boolean allTouched) { + + Map<Double, TreeSet<Double>> scanlineIntersections = new HashMap<>(); + List<LineString> allRings = new ArrayList<>(); + allRings.add(polygon.getExteriorRing()); + for (int i = 0; i < polygon.getNumInteriorRing(); i++) { + allRings.add(polygon.getInteriorRingN(i)); + } + + for (LineString ring : allRings) { + Coordinate[] coords = ring.getCoordinates(); + int numPoints = coords.length; + + if (allTouched) { + rasterizeLineString(ring, params, value, geomExtent); + } + + for (int i = 0; i < numPoints - 1; i++) { + Coordinate worldP1 = coords[i]; + Coordinate worldP2 = coords[i + 1]; + + // Ensure worldP1.y ≤ worldP2.y + if (worldP1.y > worldP2.y) { + Coordinate temp = worldP1; + worldP1 = worldP2; + worldP2 = temp; + } + + if (worldP1.y == worldP2.y) { + continue; // Skip horizontal edges + } + + double yStart = Math.floor((worldP1.y - params.upperLeftY) / params.scaleY); + double yEnd = Math.round((worldP2.y - params.upperLeftY) / params.scaleY); + + // Contain y range within geomExtent; Use centroid y line as scan line + yStart = Math.max(0.5, Math.abs(yStart) - 0.5); + yEnd = Math.min((params.writableRaster.getHeight() - 0.5), Math.abs(yEnd) + 0.5); + + double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x); + double p1X = (worldP1.x - params.upperLeftX) / params.scaleX; + double p1Y = (params.upperLeftY - worldP1.y) / params.scaleY; + + for (double y = yStart; y >= yEnd; y--) { + double xIntercept = p1X + ((p1Y - y) / slope); Review Comment: We can prune the bounds of this loop to `[max(yStart, worldP1.y), min(yEnd, worldP2.y)]`: if y is not within [worldP1.y, worldP2.y], then it won't intersect with the segment. This will skip inspecting most of the horizontal scan lines. ########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { + + double dx = x1 - x0; + double dy = y1 - y0; + + // Compute the number of steps based on the larger of dx or dy + double distance = Math.sqrt(dx * dx + dy * dy); + int steps = (int) Math.ceil(distance / stepSize); + + // Calculate the step increment for each axis + double stepX = dx / steps; + double stepY = dy / steps; + + // Start stepping through the line + double x = x0; + double y = y0; + + for (int i = 0; i <= steps; i++) { + int rasterX = (int) (Math.floor(x)); + int rasterY = (int) (Math.floor(y)); + + // Only write if within raster bounds + if (rasterX >= 0 + && rasterX < params.writableRaster.getWidth() + && rasterY >= 0 + && rasterY < params.writableRaster.getHeight()) { + params.writableRaster.setSample(rasterX, rasterY, 0, value); + } + + // Increment by fractional steps + x += stepX; + y += stepY; + } + } + + private static LineSegment clipSegmentToRasterBounds( + Coordinate p1, Coordinate p2, Envelope2D geomExtent) { + double minX = geomExtent.getMinX(); + double maxX = geomExtent.getMaxX(); + double minY = geomExtent.getMinY(); + double maxY = geomExtent.getMaxY(); + + double x1 = p1.x, y1 = p1.y; + double x2 = p2.x, y2 = p2.y; + + boolean p1Inside = isInsideBounds(x1, y1, minX, maxX, minY, maxY); + boolean p2Inside = isInsideBounds(x2, y2, minX, maxX, minY, maxY); + + if (p1Inside && p2Inside) { + // Both points inside: no clipping needed + return new LineSegment(p1, p2); + } + + if (!p1Inside && !p2Inside) { + // Both points outside: no clipping needed + return null; + } + + double dx = x2 - x1; + double dy = y2 - y1; + + // Clip using parametric line equation + double[] tValues = {0, 1}; // Stores valid segment proportions + + // Clip against minX and maxX + if (dx != 0) { + double tMin = (minX - x1) / dx; + double tMax = (maxX - x1) / dx; + if (dx < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // Clip against minY and maxY + if (dy != 0) { + double tMin = (minY - y1) / dy; + double tMax = (maxY - y1) / dy; + if (dy < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // If tValues are invalid (segment is outside), return null + if (tValues[0] > tValues[1]) { + return null; // No valid clipped segment + } + + // Compute new clipped endpoints + Coordinate newP1 = new Coordinate(x1 + tValues[0] * dx, y1 + tValues[0] * dy); + Coordinate newP2 = new Coordinate(x1 + tValues[1] * dx, y1 + tValues[1] * dy); + + return new LineSegment(newP1, newP2); + } + + // Helper function to check if a point is inside the bounding box + private static boolean isInsideBounds( + double x, double y, double minX, double maxX, double minY, double maxY) { + return x >= minX && x <= maxX && y >= minY && y <= maxY; + } + + private static Envelope2D rasterizeGeomExtent( + Geometry geom, GridCoverage2D raster, double[] metadata, boolean allTouched) { + + if (Objects.equals(geom.getGeometryType(), "MultiLineString")) { + allTouched = true; + } + if (Objects.equals(geom.getGeometryType(), "MultiPoint")) { + allTouched = true; + } + + Envelope2D rasterExtent = + JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D()); + + // Compute the aligned min/max values + double alignedMinX = + (Math.floor((rasterExtent.getMinX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMinY = + (Math.floor((rasterExtent.getMinY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + double alignedMaxX = + (Math.ceil((rasterExtent.getMaxX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMaxY = + (Math.ceil((rasterExtent.getMaxY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + + // For points and LineStrings at intersection of 2 or more pixels, + // extend search grid by 1 pixel in each direction + if (alignedMaxX == alignedMinX) { + alignedMinX -= metadata[4]; + alignedMaxX += metadata[4]; + } + if (alignedMaxY == alignedMinY) { + alignedMinY += metadata[5]; + alignedMaxY -= metadata[5]; + } + + // Get the extent of the original raster + double originalMinX = raster.getEnvelope().getMinimum(0); + double originalMinY = raster.getEnvelope().getMinimum(1); + double originalMaxX = raster.getEnvelope().getMaximum(0); + double originalMaxY = raster.getEnvelope().getMaximum(1); + + // Extend the aligned extent by 1 pixel if allTouched is true and if any geometry extent meets + // the aligned extent + if (allTouched) { + alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0; + alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0; + alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0; + alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0; + } + + alignedMinX = Math.max(alignedMinX, originalMinX); + alignedMinY = Math.max(alignedMinY, originalMinY); + alignedMaxX = Math.min(alignedMaxX, originalMaxX); + alignedMaxY = Math.min(alignedMaxY, originalMaxY); + + // Create the aligned raster extent + Envelope2D alignedRasterExtent = + new Envelope2D( + rasterExtent.getCoordinateReferenceSystem(), + alignedMinX, + alignedMinY, + alignedMaxX - alignedMinX, + alignedMaxY - alignedMinY); + + return alignedRasterExtent; + } + + private static RasterizationParams calculateRasterizationParams( + GridCoverage2D raster, + boolean useGeometryExtent, + double[] metadata, + Envelope2D geomExtent, + String pixelType) { + + double upperLeftX = 0; + double upperLeftY = 0; + if (useGeometryExtent) { + upperLeftX = geomExtent.getMinX(); + upperLeftY = geomExtent.getMaxY(); + } else { + upperLeftX = metadata[0]; + upperLeftY = metadata[1]; + } + + WritableRaster writableRaster; + if (useGeometryExtent) { + int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4])); + int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5])); + + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), geomExtentWidth, geomExtentHeight, 1, null); + } else { + int rasterWidth = (int) raster.getGridGeometry().getGridRange2D().getWidth(); + int rasterHeight = (int) raster.getGridGeometry().getGridRange2D().getHeight(); + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null); + } + + return new RasterizationParams( + writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY); + } + + private static void validateRasterMetadata(double[] rasterMetadata) { + if (rasterMetadata[4] < 0) { + throw new IllegalArgumentException("ScaleX cannot be negative"); + } + if (rasterMetadata[5] > 0) { + throw new IllegalArgumentException("ScaleY must be negative"); + } + if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) { + throw new IllegalArgumentException("SkewX and SkewY must be zero"); + } + } + + // New condensed Rasterization parameters + private static class RasterizationParams { + WritableRaster writableRaster; + String pixelType; + double scaleX; + double scaleY; + double upperLeftX; + double upperLeftY; + + RasterizationParams( + WritableRaster writableRaster, + String pixelType, + double scaleX, + double scaleY, + double upperLeftX, + double upperLeftY) { + this.writableRaster = writableRaster; + this.pixelType = pixelType; + this.scaleX = scaleX; + this.scaleY = scaleY; + this.upperLeftX = upperLeftX; + this.upperLeftY = upperLeftY; + } + } + + public static void rasterizePolygon( + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) { + if (!(geom instanceof Polygon)) { + throw new IllegalArgumentException("Only Polygon geometry is supported"); + } + + Geometry clippedGeom = + Functions.intersection(JTS.toGeometry((BoundingBox) geomExtent), Functions.buffer(geom, 0)); + + Polygon polygon = (Polygon) clippedGeom; + + // Compute scanline X-intercepts + Map<Double, TreeSet<Double>> scanlineIntersections = + computeScanlineIntersections(polygon, params, value, geomExtent, allTouched); + + // Process intersections to get startXs and endXs for each scanline + Map<Integer, List<int[]>> scanlineFillRanges = computeFillRanges(scanlineIntersections); + + // Burn values between startX and endX pairs + fillPolygon(scanlineFillRanges, params, value); + } + + /** Computes scanline intersections by iterating over polygon edges. */ + private static Map<Double, TreeSet<Double>> computeScanlineIntersections( + Polygon polygon, + RasterizationParams params, + double value, + Envelope2D geomExtent, + boolean allTouched) { + + Map<Double, TreeSet<Double>> scanlineIntersections = new HashMap<>(); + List<LineString> allRings = new ArrayList<>(); + allRings.add(polygon.getExteriorRing()); + for (int i = 0; i < polygon.getNumInteriorRing(); i++) { + allRings.add(polygon.getInteriorRingN(i)); + } + + for (LineString ring : allRings) { + Coordinate[] coords = ring.getCoordinates(); + int numPoints = coords.length; + + if (allTouched) { + rasterizeLineString(ring, params, value, geomExtent); + } + + for (int i = 0; i < numPoints - 1; i++) { + Coordinate worldP1 = coords[i]; + Coordinate worldP2 = coords[i + 1]; + + // Ensure worldP1.y ≤ worldP2.y + if (worldP1.y > worldP2.y) { + Coordinate temp = worldP1; + worldP1 = worldP2; + worldP2 = temp; + } + + if (worldP1.y == worldP2.y) { + continue; // Skip horizontal edges + } + + double yStart = Math.floor((worldP1.y - params.upperLeftY) / params.scaleY); + double yEnd = Math.round((worldP2.y - params.upperLeftY) / params.scaleY); + + // Contain y range within geomExtent; Use centroid y line as scan line + yStart = Math.max(0.5, Math.abs(yStart) - 0.5); + yEnd = Math.min((params.writableRaster.getHeight() - 0.5), Math.abs(yEnd) + 0.5); + + double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x); Review Comment: A better way is to use (worldP2.y - worldP1.y) as the denominator, since the case `worldP1.y == worldP2.y` is handled specially so `worldP2.y - worldP1.y` will never be 0 when we reach here. We are using the inverse of slope rather than slope here: ```java xIntercept = (y - worldP1.y) * (worldP2.x - worldP1.x) / (worldP2.y - worldP1.y) + worldP1.x ``` ########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } Review Comment: I think we can apply some arithmetic to figure out which pixels to render, instead of iterating all the pixels and test each of them by calling `intersects`. ########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { + + double dx = x1 - x0; + double dy = y1 - y0; + + // Compute the number of steps based on the larger of dx or dy + double distance = Math.sqrt(dx * dx + dy * dy); + int steps = (int) Math.ceil(distance / stepSize); + + // Calculate the step increment for each axis + double stepX = dx / steps; + double stepY = dy / steps; + + // Start stepping through the line + double x = x0; + double y = y0; + + for (int i = 0; i <= steps; i++) { + int rasterX = (int) (Math.floor(x)); + int rasterY = (int) (Math.floor(y)); + + // Only write if within raster bounds + if (rasterX >= 0 + && rasterX < params.writableRaster.getWidth() + && rasterY >= 0 + && rasterY < params.writableRaster.getHeight()) { + params.writableRaster.setSample(rasterX, rasterY, 0, value); + } + + // Increment by fractional steps + x += stepX; + y += stepY; + } + } + + private static LineSegment clipSegmentToRasterBounds( + Coordinate p1, Coordinate p2, Envelope2D geomExtent) { + double minX = geomExtent.getMinX(); + double maxX = geomExtent.getMaxX(); + double minY = geomExtent.getMinY(); + double maxY = geomExtent.getMaxY(); + + double x1 = p1.x, y1 = p1.y; + double x2 = p2.x, y2 = p2.y; + + boolean p1Inside = isInsideBounds(x1, y1, minX, maxX, minY, maxY); + boolean p2Inside = isInsideBounds(x2, y2, minX, maxX, minY, maxY); + + if (p1Inside && p2Inside) { + // Both points inside: no clipping needed + return new LineSegment(p1, p2); + } + + if (!p1Inside && !p2Inside) { + // Both points outside: no clipping needed + return null; + } + + double dx = x2 - x1; + double dy = y2 - y1; + + // Clip using parametric line equation + double[] tValues = {0, 1}; // Stores valid segment proportions + + // Clip against minX and maxX + if (dx != 0) { + double tMin = (minX - x1) / dx; + double tMax = (maxX - x1) / dx; + if (dx < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // Clip against minY and maxY + if (dy != 0) { + double tMin = (minY - y1) / dy; + double tMax = (maxY - y1) / dy; + if (dy < 0) { + double temp = tMin; + tMin = tMax; + tMax = temp; + } + tValues[0] = Math.max(tValues[0], tMin); + tValues[1] = Math.min(tValues[1], tMax); + } + + // If tValues are invalid (segment is outside), return null + if (tValues[0] > tValues[1]) { + return null; // No valid clipped segment + } + + // Compute new clipped endpoints + Coordinate newP1 = new Coordinate(x1 + tValues[0] * dx, y1 + tValues[0] * dy); + Coordinate newP2 = new Coordinate(x1 + tValues[1] * dx, y1 + tValues[1] * dy); + + return new LineSegment(newP1, newP2); + } + + // Helper function to check if a point is inside the bounding box + private static boolean isInsideBounds( + double x, double y, double minX, double maxX, double minY, double maxY) { + return x >= minX && x <= maxX && y >= minY && y <= maxY; + } + + private static Envelope2D rasterizeGeomExtent( + Geometry geom, GridCoverage2D raster, double[] metadata, boolean allTouched) { + + if (Objects.equals(geom.getGeometryType(), "MultiLineString")) { + allTouched = true; + } + if (Objects.equals(geom.getGeometryType(), "MultiPoint")) { + allTouched = true; + } + + Envelope2D rasterExtent = + JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D()); + + // Compute the aligned min/max values + double alignedMinX = + (Math.floor((rasterExtent.getMinX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMinY = + (Math.floor((rasterExtent.getMinY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + double alignedMaxX = + (Math.ceil((rasterExtent.getMaxX() + metadata[0]) / metadata[4]) * metadata[4]) + - metadata[0]; + double alignedMaxY = + (Math.ceil((rasterExtent.getMaxY() + metadata[1]) / -metadata[5]) * -metadata[5]) + - metadata[1]; + + // For points and LineStrings at intersection of 2 or more pixels, + // extend search grid by 1 pixel in each direction + if (alignedMaxX == alignedMinX) { + alignedMinX -= metadata[4]; + alignedMaxX += metadata[4]; + } + if (alignedMaxY == alignedMinY) { + alignedMinY += metadata[5]; + alignedMaxY -= metadata[5]; + } + + // Get the extent of the original raster + double originalMinX = raster.getEnvelope().getMinimum(0); + double originalMinY = raster.getEnvelope().getMinimum(1); + double originalMaxX = raster.getEnvelope().getMaximum(0); + double originalMaxY = raster.getEnvelope().getMaximum(1); + + // Extend the aligned extent by 1 pixel if allTouched is true and if any geometry extent meets + // the aligned extent + if (allTouched) { + alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0; + alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0; + alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0; + alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0; + } + + alignedMinX = Math.max(alignedMinX, originalMinX); + alignedMinY = Math.max(alignedMinY, originalMinY); + alignedMaxX = Math.min(alignedMaxX, originalMaxX); + alignedMaxY = Math.min(alignedMaxY, originalMaxY); + + // Create the aligned raster extent + Envelope2D alignedRasterExtent = + new Envelope2D( + rasterExtent.getCoordinateReferenceSystem(), + alignedMinX, + alignedMinY, + alignedMaxX - alignedMinX, + alignedMaxY - alignedMinY); + + return alignedRasterExtent; + } + + private static RasterizationParams calculateRasterizationParams( + GridCoverage2D raster, + boolean useGeometryExtent, + double[] metadata, + Envelope2D geomExtent, + String pixelType) { + + double upperLeftX = 0; + double upperLeftY = 0; + if (useGeometryExtent) { + upperLeftX = geomExtent.getMinX(); + upperLeftY = geomExtent.getMaxY(); + } else { + upperLeftX = metadata[0]; + upperLeftY = metadata[1]; + } + + WritableRaster writableRaster; + if (useGeometryExtent) { + int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4])); + int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5])); + + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), geomExtentWidth, geomExtentHeight, 1, null); + } else { + int rasterWidth = (int) raster.getGridGeometry().getGridRange2D().getWidth(); + int rasterHeight = (int) raster.getGridGeometry().getGridRange2D().getHeight(); + writableRaster = + RasterFactory.createBandedRaster( + RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null); + } + + return new RasterizationParams( + writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY); + } + + private static void validateRasterMetadata(double[] rasterMetadata) { + if (rasterMetadata[4] < 0) { + throw new IllegalArgumentException("ScaleX cannot be negative"); + } + if (rasterMetadata[5] > 0) { + throw new IllegalArgumentException("ScaleY must be negative"); + } + if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) { + throw new IllegalArgumentException("SkewX and SkewY must be zero"); + } + } + + // New condensed Rasterization parameters + private static class RasterizationParams { + WritableRaster writableRaster; + String pixelType; + double scaleX; + double scaleY; + double upperLeftX; + double upperLeftY; + + RasterizationParams( + WritableRaster writableRaster, + String pixelType, + double scaleX, + double scaleY, + double upperLeftX, + double upperLeftY) { + this.writableRaster = writableRaster; + this.pixelType = pixelType; + this.scaleX = scaleX; + this.scaleY = scaleY; + this.upperLeftX = upperLeftX; + this.upperLeftY = upperLeftY; + } + } + + public static void rasterizePolygon( + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) { + if (!(geom instanceof Polygon)) { + throw new IllegalArgumentException("Only Polygon geometry is supported"); + } + + Geometry clippedGeom = + Functions.intersection(JTS.toGeometry((BoundingBox) geomExtent), Functions.buffer(geom, 0)); + + Polygon polygon = (Polygon) clippedGeom; + + // Compute scanline X-intercepts + Map<Double, TreeSet<Double>> scanlineIntersections = + computeScanlineIntersections(polygon, params, value, geomExtent, allTouched); + + // Process intersections to get startXs and endXs for each scanline + Map<Integer, List<int[]>> scanlineFillRanges = computeFillRanges(scanlineIntersections); + + // Burn values between startX and endX pairs + fillPolygon(scanlineFillRanges, params, value); + } + + /** Computes scanline intersections by iterating over polygon edges. */ + private static Map<Double, TreeSet<Double>> computeScanlineIntersections( + Polygon polygon, + RasterizationParams params, + double value, + Envelope2D geomExtent, + boolean allTouched) { + + Map<Double, TreeSet<Double>> scanlineIntersections = new HashMap<>(); + List<LineString> allRings = new ArrayList<>(); + allRings.add(polygon.getExteriorRing()); + for (int i = 0; i < polygon.getNumInteriorRing(); i++) { + allRings.add(polygon.getInteriorRingN(i)); + } + + for (LineString ring : allRings) { + Coordinate[] coords = ring.getCoordinates(); + int numPoints = coords.length; + + if (allTouched) { + rasterizeLineString(ring, params, value, geomExtent); + } + + for (int i = 0; i < numPoints - 1; i++) { + Coordinate worldP1 = coords[i]; + Coordinate worldP2 = coords[i + 1]; + + // Ensure worldP1.y ≤ worldP2.y + if (worldP1.y > worldP2.y) { + Coordinate temp = worldP1; + worldP1 = worldP2; + worldP2 = temp; + } + + if (worldP1.y == worldP2.y) { + continue; // Skip horizontal edges + } + + double yStart = Math.floor((worldP1.y - params.upperLeftY) / params.scaleY); + double yEnd = Math.round((worldP2.y - params.upperLeftY) / params.scaleY); + + // Contain y range within geomExtent; Use centroid y line as scan line + yStart = Math.max(0.5, Math.abs(yStart) - 0.5); + yEnd = Math.min((params.writableRaster.getHeight() - 0.5), Math.abs(yEnd) + 0.5); + + double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x); Review Comment: There will be division by zero problem when the segment is an vertical edge (`worldP2.x == worldP1.x`). -- 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]
