12

I am working with geospatial shapes and looking at the centroid algorithm here,

http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon

I have implemented the code in C# like this (which is just this adapted),

Finding the centroid of a polygon?

class Program
{
    static void Main(string[] args)
    {
        List<Point> vertices = new List<Point>();

        vertices.Add(new Point() { X = 1, Y = 1 });
        vertices.Add(new Point() { X = 1, Y = 10 });
        vertices.Add(new Point() { X = 2, Y = 10 });
        vertices.Add(new Point() { X = 2, Y = 2 });
        vertices.Add(new Point() { X = 10, Y = 2 });
        vertices.Add(new Point() { X = 10, Y = 1 });
        vertices.Add(new Point() { X = 1, Y = 1 });

        Point centroid = Compute2DPolygonCentroid(vertices);
    }

    static Point Compute2DPolygonCentroid(List<Point> vertices)
    {
        Point centroid = new Point() { X = 0.0, Y = 0.0 };
        double signedArea = 0.0;
        double x0 = 0.0; // Current vertex X
        double y0 = 0.0; // Current vertex Y
        double x1 = 0.0; // Next vertex X
        double y1 = 0.0; // Next vertex Y
        double a = 0.0;  // Partial signed area

        // For all vertices except last
        int i=0;
        for (i = 0; i < vertices.Count - 1; ++i)
        {
            x0 = vertices[i].X;
            y0 = vertices[i].Y;
            x1 = vertices[i+1].X;
            y1 = vertices[i+1].Y;
            a = x0*y1 - x1*y0;
            signedArea += a;
            centroid.X += (x0 + x1)*a;
            centroid.Y += (y0 + y1)*a;
        }

        // Do last vertex
        x0 = vertices[i].X;
        y0 = vertices[i].Y;
        x1 = vertices[0].X;
        y1 = vertices[0].Y;
        a = x0*y1 - x1*y0;
        signedArea += a;
        centroid.X += (x0 + x1)*a;
        centroid.Y += (y0 + y1)*a;

        signedArea *= 0.5;
        centroid.X /= (6*signedArea);
        centroid.Y /= (6*signedArea);

        return centroid;
    }
}

public class Point
{
    public double X { get; set; }
    public double Y { get; set; }
}

The problem is that this algorithm when I have this shape (which is an L shape),

(1,1) (1,10) (2,10) (2,2) (10,2) (10,1) (1,1)

It gives me the result (3.62, 3.62). Which is OK, except that point is outside the shape. Is there another algorithm around that takes this into account?

Basically a person is going to be drawing a shape on a map. This shape might span multiple roads (so could be an L shape) and I want to work out the centre of the shape. This is so I can work out the road name at that point. It doesn't make sense to me for it to be outside the shape if they have drawn a long skinny L shape.

Community
  • 1
  • 1
peter
  • 13,009
  • 22
  • 82
  • 142
  • 4
    The centroid of a polygon doesn't have to be inside it. That is only guaranteed to apply for *convex* polygons. – Hong Ooi Mar 22 '12 at 02:51
  • Yes the algorithm is correct I agree, but is there another algorithm that will ensure that a point within the polygon is calculated? Ideally the result for the above shape would be (1.5, 1.5). – peter Mar 22 '12 at 02:57
  • If it makes sense for your problem, you can represent roads as thick lines or even union or rectangles (there by bringing notion of axis in). Your center is mid-point of this axis. – Om Deshmane Mar 22 '12 at 04:03
  • Don't use a point to find the name of the road. There would be no guarantee the point would be close enough to the road being outlined. The right thing to do is retrieve the road line segments in that area and figure out which road has the most length in the polygon. – CodeSlinger Apr 26 '12 at 18:06

5 Answers5

20

This answer is inspired by the answer by Jer2654 and this source: http://coding-experiments.blogspot.com/2009/09/xna-quest-for-centroid-of-polygon.html

  /// <summary>
  /// Method to compute the centroid of a polygon. This does NOT work for a complex polygon.
  /// </summary>
  /// <param name="poly">points that define the polygon</param>
  /// <returns>centroid point, or PointF.Empty if something wrong</returns>
  public static PointF GetCentroid(List<PointF> poly)
  {
     float accumulatedArea = 0.0f;
     float centerX = 0.0f;
     float centerY = 0.0f;

     for (int i = 0, j = poly.Count - 1; i < poly.Count; j = i++)
     {
        float temp = poly[i].X * poly[j].Y - poly[j].X * poly[i].Y;
        accumulatedArea += temp;
        centerX += (poly[i].X + poly[j].X) * temp;
        centerY += (poly[i].Y + poly[j].Y) * temp;
     }

     if (Math.Abs(accumulatedArea) < 1E-7f)
        return PointF.Empty;  // Avoid division by zero

     accumulatedArea *= 3f;
     return new PointF(centerX / accumulatedArea, centerY / accumulatedArea);
  }
RenniePet
  • 11,420
  • 7
  • 80
  • 106
  • if (Math.Abs(accumulatedArea) < 1E-7f) – Koray Feb 23 '17 at 12:37
  • 1
    @Koray Thank you! Don't quite understand the need for this but I assume you are right, and I've included your suggestion in my answer - and in my own usage of this code, also a Java version I use in an Android app. – RenniePet Feb 23 '17 at 16:38
  • 3
    Your code works just fine, however if the some coordinates are negative, this is required. If it is guaranteed that every coordinate is positive, there is no need. Thank you. – Koray Feb 24 '17 at 10:24
7
public static Point GetCentroid( Point[ ] nodes, int count )
{
    int x = 0, y = 0, area = 0, k;
    Point a, b = nodes[ count - 1 ];

    for( int i = 0; i < count; i++ )
    {
        a = nodes[ i ];

        k = a.Y * b.X - a.X * b.Y;
        area += k;
        x += ( a.X + b.X ) * k;
        y += ( a.Y + b.Y ) * k;

        b = a;
    }
    area *= 3;

    return ( area == 0 ) ? Point.Empty : new Point( x /= area, y /= area );
}
JG in SD
  • 5,427
  • 3
  • 34
  • 46
Jer2654
  • 81
  • 1
  • 3
  • 1
    Looks like this is based on the same algorithm as used here: http://coding-experiments.blogspot.com/2009/09/xna-quest-for-centroid-of-polygon.html – RenniePet Nov 03 '13 at 05:53
5

You may check if .NET 4.5 DbSpatialServices functions, like DbSpatialServices.GetCentroid

ITmeze
  • 1,902
  • 2
  • 21
  • 32
0

For 3d Points, I created a method in C # I hope I can help you:

public static double[] GetCentroid(List<double[]> listOfPoints)
    {
        // centroid[0] = X
        // centroid[1] = Y
        // centroid[2] = Z
        double[] centroid = new double[3];

        // List iteration
        // Link reference:
        // https://en.wikipedia.org/wiki/Centroid
        foreach (double[] point in listOfPoints)
        {
            centroid[0] += point[0];
            centroid[1] += point[1];
            centroid[2] += point[2];
        }

        centroid[0] /= listOfPoints.Count;
        centroid[1] /= listOfPoints.Count;
        centroid[2] /= listOfPoints.Count;

        return centroid;
    }
-1

My implementation:

GpsCoordinates GetCentroid(ICollection<GpsCoordinates> polygonCorners)
{
    return new GpsCoordinates(polygonCorners.Average(x => x.Latitude), polygonCorners.Average(x => x.Longitude));
}

public readonly struct GpsCoordinates
{
    public GpsCoordinates(
        double latitude,
        double longitude
        )
    {
       
        Latitude = latitude;
        Longitude = longitude;
    }

    public double Latitude { get; }
    public double Longitude { get; }
}

The idea taken from here

xhafan
  • 2,140
  • 1
  • 26
  • 26