bertt / mapbox-vector-tile-cs

A .NET library for decoding a Mapbox vector tile
MIT License
74 stars 14 forks source link

Signed area calculation is incorrect ... #12

Closed Tom-Cuthill closed 4 years ago

Tom-Cuthill commented 7 years ago

It appears that the code to implement the Surveyors / Shoelace formula in the SignedArea method is incomplete. It is missing the boundary conditions when you consider the last point's X needs to match with the first point's Y and vice-versa.

Consider the simple triangle with vertices: [1,1], [2,2], [3,1]. The area of a triangle is a half of the base times the height. That gives an area of 1. Your algorithm gives an area of 2.

A proper implementation would look like this:

private double SignedArea(List points) { double sum = 0.0; for (int i = 0; i < points.Count; ++i) { sum += points[i].X (((i + 1) == points.Count) ? points[0].Y : points[i + 1].Y); sum -= points[i].Y (((i + 1) == points.Count) ? points[0].X : points[i + 1].X); } return 0.5 * sum; }

Note as well, that vector tiles have an origin at the top left of the tile (ie. the Y axis is inverted). As a result your tests for clockwise / counterclockwise are the opposite of what they should be.

According to the 2.1 spec (https://github.com/mapbox/vector-tile-spec/tree/master/2.1) :

4.3.4.4. Polygon Geometry Type ... An exterior ring is DEFINED as a linear ring having a positive area as calculated by applying the surveyor's formula to the vertices of the polygon in tile coordinates. In the tile coordinate system (with the Y axis positive down and X axis positive to the right) this makes the exterior ring's winding order appear clockwise.

An interior ring is DEFINED as a linear ring having a negative area as calculated by applying the surveyor's formula to the vertices of the polygon in tile coordinates. In the tile coordinate system (with the Y axis positive down and X axis positive to the right) this makes the interior ring's winding order appear counterclockwise.

Hope this helps, Tom

bertt commented 7 years ago

Hi, I've added a testcase for the signed area https://github.com/bertt/mapbox-vector-tile-cs/commit/100197fcc094090ce9b1f319e67b1fd4563e23ba

The signed area function assumes the polygon is closed (first point is the same as last point), maybe I should add a check on that.

In your simple triangle case it results in area == -1 (when polygon closed), its negative because polygon is defined clockwise. So I think the area calculation is correct?

I'll take an extra check on the CW/CCW issue you mention later this week.

Tom-Cuthill commented 7 years ago

Hi. In the definition seen in the 4.3.4.4 above, it is a positive area that is clockwise, not a negative area. This is because the Y axis is inverted. Your code claims a negative area is clockwise, but it is the opposite that is true.

    public bool IsCW()
    {
        return SignedArea() < 0;
    }

Here is your calculation for the signed area: for (var i = 0; i < points.Count-1; i++) { sum = sum + (points[i].X points[i + 1].Y - (points[i].Y points[i + 1].X)); }

It is missing out the term for point[Count-1].X. This X value has to loop back to the point[0].Y. That is why it's called the shoelace formula, the X and Y values cross over, and the last X and Y values go back to the beginning.

shoelace

Your code misses out the XnY1 and -X1Yn terms, where n is point.Count - 1. This corrects the sign of the area and also the area itself.

bertt commented 7 years ago

can you create a pr for the code changes involved?