Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fractional-Pixel Polygon Rasterizer #1873

Merged
merged 5 commits into from
Dec 7, 2016
Merged

Fractional-Pixel Polygon Rasterizer #1873

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9513863
Add New Callback Type To Support Partial Cells
Dec 2, 2016
e4065f0
Fractional-Pixel Polygon Rasterizer
Dec 5, 2016
9454ae1
Fractional-Pixel Rasterizer Unit Tests
Dec 5, 2016
75cab0d
Use Bresenham-Inspired Approach For Line Drawing
Dec 5, 2016
30af5a3
PR 1873 Review Comments
Dec 6, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
195 changes: 195 additions & 0 deletions ...er-test/src/test/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizerSpec.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,195 @@
/*
* Copyright 2016 Azavea
*
* Licensed 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 geotrellis.raster.rasterize.polygon

import geotrellis.raster._
import geotrellis.raster.rasterize._
import geotrellis.vector._

import org.scalatest.{FunSpec, Matchers}

import scala.math.round


class FractionalRasterizerSpec extends FunSpec with Matchers {

describe("Fractional-Pixel Polygon Rasterizer") {
val e = Extent(0, 0, 3, 3)
val re = RasterExtent(e, 3, 3)

it("should correctly report full pixels") {
val poly = Polygon(Point(1,1), Point(1,2), Point(2,2), Point(2,1), Point(1,1))
var actual = Double.NaN
val expected = 1.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
if (col == 1 && row == 1) actual = p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

actual should be (expected)
}

it("should correctly not report disjoint pixels") {
val poly = Polygon(Point(1,1), Point(1,2), Point(2,2), Point(2,1), Point(1,1))
var actual = 0.0
val expected = 0.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
if (col != 1 && row != 1) actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

actual should be (expected)
}

it("should correctly report partial pixels") {
val poly = Polygon(Point(1,1), Point(1,2), Point(2,2), Point(1,1))
var actual = Double.NaN
val expected = 0.5
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
if (col == 1 && row == 1) actual = p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

actual should be (expected)
}

it("should handle sub-pixel activity") {
val poly = Polygon(Point(1.2,1.2), Point(1.2,1.8), Point(1.8,1.8), Point(1.8,1.2), Point(1.2,1.2))
var actual = Double.NaN
val expected = 0.36
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
if (col == 1 && row == 1) actual = p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should handle a mix of partial and complete pixels") {
val poly = Polygon(Point(0.1,0.1), Point(0.1,2.9), Point(2.9,2.9), Point(2.9,0.1), Point(0.1,0.1))
var actual = 0.0
val expected = 2.8 * 2.8
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (1/6)") {
val re = RasterExtent(e, 30, 30)
val poly = Polygon(Point(0,0), Point(3,0), Point(3,3), Point(0,0))
var actual = 0.0
val expected = (30*30)/2.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (2/6)") {
val re = RasterExtent(e, 30, 30)
val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
var actual = 0.0
val expected = (30*30)/2.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (3/6)") {
val re = RasterExtent(Extent(0, 0, 3, 3.1), 30, 30)
val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
var actual = 0.0
val expected = 435.483871
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (4/6)") {
val re = RasterExtent(Extent(0, 0, 3.1, 3), 30, 30)
val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
var actual = 0.0
val expected = 435.483871
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (5/6)") {
val re = RasterExtent(Extent(0, 0, 3, 3.1), 30, 30)
val poly = Polygon(Point(0,0), Point(3,0), Point(0,3.1), Point(0,0))
var actual = 0.0
val expected = (30*30)/2.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

it("should efficiently handle a long diagonal line (6/6)") {
val re = RasterExtent(Extent(0, 0, 3.1, 3), 30, 30)
val poly = Polygon(Point(0,0), Point(3.1,0), Point(0,3), Point(0,0))
var actual = 0.0
val expected = (30*30)/2.0
val cb = new FractionCallback {
def callback(col: Int, row: Int, p: Double): Unit =
actual += p
}

FractionalRasterizer.foreachCellByPolygon(poly, re)(cb)

round(actual * 1000000) should be (round(expected * 1000000))
}

}
}
8 changes: 8 additions & 0 deletions raster/src/main/scala/geotrellis/raster/rasterize/package.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,4 +21,12 @@ package object rasterize {
/** Callback for given row and column (compatible with previous definition). */
type Callback = (Int, Int) => Unit

/**
* A call back which accepts a col-row-fraction triple. The
* fraction is the fraction of the pixel which is covered by the
* query object.
*/
trait FractionCallback {
def callback(col: Int, row: Int, fraction: Double): Unit
}
}
175 changes: 175 additions & 0 deletions raster/src/main/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizer.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,175 @@
package geotrellis.raster.rasterize.polygon

import geotrellis.raster._
import geotrellis.raster.rasterize._
import geotrellis.vector._

import com.vividsolutions.jts.geom.Envelope
import spire.syntax.cfor._

import scala.collection.mutable
import scala.math.{min, max, ceil, floor, abs}


object FractionalRasterizer {

private type Segment = (Double, Double, Double, Double)

private def polygonToEdges(poly: Polygon, re: RasterExtent): Seq[Segment] = {

val arrayBuffer = mutable.ArrayBuffer.empty[Segment]

/** Find the outer ring's segments */
val coords = poly.jtsGeom.getExteriorRing.getCoordinates
cfor(1)(_ < coords.length, _ + 1) { ci =>
val coord1 = coords(ci - 1)
val coord2 = coords(ci)

val col1 = re.mapXToGridDouble(coord1.x)
val row1 = re.mapYToGridDouble(coord1.y)
val col2 = re.mapXToGridDouble(coord2.x)
val row2 = re.mapYToGridDouble(coord2.y)

val segment =
if (col1 < col2) (col1, row1, col2, row2)
else (col2, row2, col1, row1)

arrayBuffer += segment
}

/** Find the segments for the holes */
cfor(0)(_ < poly.numberOfHoles, _ + 1) { i =>
val coords = poly.jtsGeom.getInteriorRingN(i).getCoordinates
cfor(1)(_ < coords.length, _ + 1) { ci =>
val coord1 = coords(ci - 1)
val coord2 = coords(ci)

val col1 = re.mapXToGridDouble(coord1.x)
val row1 = re.mapYToGridDouble(coord1.y)
val col2 = re.mapXToGridDouble(coord2.x)
val row2 = re.mapYToGridDouble(coord2.y)

val segment =
if (col1 < col2) (col1, row1, col2, row2)
else (col2, row2, col1, row1)

arrayBuffer += segment
}
}

arrayBuffer
}

private def renderEdge(
edge: Segment,
re: RasterExtent,
poly: Polygon,
set: mutable.Set[(Int, Int)],
cb: FractionCallback
): Unit = {
val (x0, y0, x1, y1) = edge
val m = (y1 - y0) / (x1 - x0)

// Grid coordinates
val colMin = floor(min(x0, x1)).toInt
val rowMin = floor(min(y0, y1)).toInt
val colMax = ceil(max(x0, x1)).toInt
val rowMax = ceil(max(y0, y1)).toInt

// Map coordinates
val xmin = re.gridColToMap(colMin) - re.cellwidth/2
val ymin = re.gridRowToMap(rowMax) + re.cellheight/2
val xmax = re.gridColToMap(colMax) - re.cellwidth/2
val ymax = re.gridRowToMap(rowMin) + re.cellheight/2

// Envelope around the edge (in map space)
val envelope = Polygon(
Point(xmin, ymin),
Point(xmin, ymax),
Point(xmax, ymax),
Point(xmax, ymin),
Point(xmin, ymin)
).jtsGeom

// Intersection of envelope and polygon (in map space)
val localPoly = poly.jtsGeom.intersection(envelope)

if (abs(m) <= 1) { // The edge is mostly horizontal
var x = colMin; while (x <= colMax) {
val _y = floor(m * (x + 0.5 - x0) + y0).toInt
var i = -1; while (i <= 1) {
val y = _y + i
val pair = (x, y)
val pixelMinX = re.gridColToMap(x+0) - re.cellwidth/2
val pixelMaxX = re.gridColToMap(x+1) - re.cellwidth/2
val pixelMinY = re.gridRowToMap(y+0) + re.cellheight/2
val pixelMaxY = re.gridRowToMap(y+1) + re.cellheight/2
val pixel = Polygon(
Point(pixelMinX, pixelMinY),
Point(pixelMinX, pixelMaxY),
Point(pixelMaxX, pixelMaxY),
Point(pixelMaxX, pixelMinY),
Point(pixelMinX, pixelMinY)
).jtsGeom
val fraction = (pixel.intersection(localPoly)).getArea / pixel.getArea

if (fraction > 0.0) {
if (!set.contains(pair)) {
set += ((x, y))
cb.callback(x, y, fraction)
}
}
i += 1
}
x += 1
}
} else { // The edge is mostly vertical
val m = (x1 - x0) / (y1 - y0)
var y = rowMin; while (y <= rowMax) {
val _x = floor(m * (y + 0.5 - y0) + x0).toInt
var i = -1; while (i <= 1) {
val x = _x + i
val pair = (x, y)
val pixelMinX = re.gridColToMap(x+0) - re.cellwidth/2
val pixelMaxX = re.gridColToMap(x+1) - re.cellwidth/2
val pixelMinY = re.gridRowToMap(y+0) + re.cellheight/2
val pixelMaxY = re.gridRowToMap(y+1) + re.cellheight/2
val pixel = Polygon(
Point(pixelMinX, pixelMinY),
Point(pixelMinX, pixelMaxY),
Point(pixelMaxX, pixelMaxY),
Point(pixelMaxX, pixelMinY),
Point(pixelMinX, pixelMinY)
).jtsGeom
val fraction = (pixel.intersection(localPoly)).getArea / pixel.getArea

if (fraction > 0.0) {
if (!set.contains(pair)) {
set += ((x, y))
cb.callback(x, y, fraction)
}
}
i += 1
}
y += 1
}
}
}

def foreachCellByPolygon(
poly: Polygon,
re: RasterExtent
)(cb: FractionCallback): Unit = {
val seen = mutable.Set.empty[(Int, Int)]
val option = Rasterizer.Options(includePartial = false, sampleType = PixelIsArea)

polygonToEdges(poly, re)
.foreach({ edge => renderEdge(edge, re, poly, seen, cb) })

PolygonRasterizer.foreachCellByPolygon(poly, re) {(col: Int, row: Int) =>
val pair = (col, row)
if (!seen.contains(pair)) cb.callback(col, row, 1.0)
}
}

}
Loading