Skip to content

Commit

Permalink
Reduce number of transform invocations in area calc
Browse files Browse the repository at this point in the history
When calculating an area with a transformation, we only really have to
transform each point exactly once. However, the previous implementation
was naively calculating each transformation twice. For simple
transformations, this may not matter much. However, the transformation
is user-supplied, so could be arbitrarily expensive to compute.

The bounds of the for loop in the shoelace formula implementation have
been reduced by 1 from from `[0, n)` to `[0, n-1)`. This is because the
shoelace formula applies to rings, and the start and end points of rings
are always the same. So the final iteration in the previous
implementation was redundant.
  • Loading branch information
peterstace committed Dec 1, 2021
1 parent 7c7cb33 commit 3fb4846
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 5 deletions.
12 changes: 12 additions & 0 deletions geom/attr_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -642,6 +642,18 @@ func TestTransformedArea(t *testing.T) {
}
}

func TestTransformedAreaInvocationCount(t *testing.T) {
g := geomFromWKT(t, "POLYGON((0 0,0 1,1 0,0 0))")
var count int
g.Area(WithTransform(func(xy XY) XY {
count++
return xy
}))

// Each of the 4 points making up the polygon get transformed once each.
expectIntEq(t, count, 4)
}

func TestCentroid(t *testing.T) {
for i, tt := range []struct {
input string
Expand Down
19 changes: 14 additions & 5 deletions geom/type_polygon.go
Original file line number Diff line number Diff line change
Expand Up @@ -398,13 +398,22 @@ func signedAreaOfLinearRing(lr LineString, transform func(XY) XY) float64 {
var sum float64
seq := lr.Coordinates()
n := seq.Length()
for i := 0; i < n; i++ {
pt0 := seq.GetXY(i)
pt1 := seq.GetXY((i + 1) % n)
if n == 0 {
return 0
}

nthPt := func(i int) XY {
pt := seq.GetXY(i)
if transform != nil {
pt0 = transform(pt0)
pt1 = transform(pt1)
pt = transform(pt)
}
return pt
}

pt1 := nthPt(0)
for i := 0; i < n-1; i++ {
pt0 := pt1
pt1 = nthPt(i + 1)
sum += (pt1.X + pt0.X) * (pt1.Y - pt0.Y)
}
return sum / 2
Expand Down

0 comments on commit 3fb4846

Please sign in to comment.