Topic: centroid of closed polyline

[Chumplybum]

hi all, i'm trying to find the centroid of a closed polyline.

i've read a few lisp examples that can do this (including one on the swamp http://www.theswamp.org/index.php?topic=18725.0) but haven't been able to transform this code to .net. From the lisp examples it seems that the polyline is converted to a region and then the centriod is found from the region, however from what i can see there is no centroid property (although there is a centroid property under solid3dmassproperties).
So i couldn't convert lisp to .net and went off and found the formula for the centroid of a polyline (http://en.wikipedia.org/wiki/Centroid) as well as a vb.net sample (http://www.vb-helper.com/howto_net_polygon_centroid.html) after slightly modifying the code i've come up with the following...

Code: [Select]

        ' For polylines
        Dim pl As AADS.Polyline = TryCast(context.PickedObject, AADS.Polyline)

        If pl Is Nothing OrElse pl.Closed = False Then Return

        Dim RunningX As Long
        Dim RunningY As Long
        Dim second_factor As Long

        For Counter As Integer = 0 To pl.NumberOfVertices - 2
            second_factor = pl.GetPoint2dAt(Counter).X * pl.GetPoint2dAt(Counter + 1).Y - pl.GetPoint2dAt(Counter + 1).X * pl.GetPoint2dAt(Counter).Y

            RunningX += (pl.GetPoint2dAt(Counter).X + pl.GetPoint2dAt(Counter + 1).X) * second_factor
            RunningY += (pl.GetPoint2dAt(Counter).Y + pl.GetPoint2dAt(Counter + 1).Y) * second_factor
        Next


        Dim X, Y As Long

        ' Divide by 6 times the polygon's area.
        Dim polygon_area As Single = pl.Area
        X = RunningX / (6 * pl.Area)
        Y = RunningY / (6 * pl.Area)

        ' If the values are negative, the polygon is
        ' oriented counterclockwise. Reverse the signs.
        If X < 0 Then
            X = -X
            Y = -Y
        End If

        Dim CentroidPoint As AAG.Point3d = New AAG.Point3d(X, Y, pl.Elevation)


the code works for a rectangle that is drawn at 0,0... however, when the polyline moves away from 0,0 the code produces some crazy results and if the polyline is in a different quadrant about the usc it produces some more crazy results

does anyone know of a better / easier way to find the centroid, or a problem with my code???

any help would be great


cheers, Mark

[pkohut]

You have to do the calculations relative to the object you're trying to measure.  The code would
look something like this (sorry I don't have VB skills)

Code: [Select]

dim ptRef as Point2d
dim pt1 as Point2d
dim pt2 as Point2d
dim second_factor as double
dim drX as double
dim drY as double

For Counter as Integer = 0 To pl.NumberOfVertices - 2
  pt1 = pl.GetPoint2dAt(Counter)
  Pt2 = pl.GetPoint2dAt(counter + 1)
  if i == 0 then
      ptRef = pt1
  endif
  second_factor = (pt1.x - ptRef.x) * (pt2.y - ptRef.y) - (pt2.x - ptRef.x) * (pt1.y - ptRef.y);
  drX += ((pt1.x - ptRef.x) + (pt2.x - ptRef.x)) * b0;
  drY += ((pt1.y - ptRef.y) + (pt2.y - ptRef.y)) * b0;
Next

x = drX / (6 * pl.area) + ptRef.x
y = drY / (6 * pl.area) + ptRef.y


With this there is no need to reverse sign at the end.  Your example code is declaring some variables
as Long when they should be doubles (or whatever VB calls them).  Was your code returning
a float or double for those computations?  If they are then I would turn on the Option strict to get a
little tighter control over the VB.

Paul

[Cathy]

It looks to me like you are missing one side of the polygon -- the one from the last point to the first point.  I note that the sample code at http://www.vb-helper.com/howto_net_polygon_centroid.html makes a copy of the first point at the end, but I don't see where you do that.

[SomeCallMeDave]

You could also used TF.net  http://code.google.com/p/tf-net/ and let it do the work for you.


C# code, but I'm sure you can convert it

Code: [Select]


using System;
using System.Linq;
using System.Text;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Topology.Geometries;
using Topology.IO.Dwg;

namespace TF_Test
{

    public class TF_Test
    {
        [CommandMethod("TFCentroid")]
        public void FindCentroid()
        {
            {
                Document doc = Application.DocumentManager.MdiActiveDocument;
                Database db = doc.Database;
                Editor ed = doc.Editor;

                PromptEntityOptions opts = new PromptEntityOptions("\nSelect Entity to Query: ");
                opts.AllowNone = true;
                opts.AllowObjectOnLockedLayer = true;
                opts.SetRejectMessage("\nDon't pick that kind ");
                opts.AddAllowedClass(typeof(Polyline), false);


                //Select the object on-screen
                PromptEntityResult res = ed.GetEntity(opts);
                if (res.Status == PromptStatus.OK)
                {
                    ObjectId ObjID = res.ObjectId;
                    try
                    {
                        using (Transaction trans = db.TransactionManager.StartTransaction())
                        {
                            Entity ent = (Entity)trans.GetObject(ObjID, OpenMode.ForRead, false);
                            DwgReader reader = new DwgReader();
                            DwgWriter writer = new DwgWriter();
                            Geometry geometry = (Geometry)reader.ReadGeometry(ent);
                            IPoint centroid = geometry.Centroid;
                            ed.WriteMessage("\nCentroid is " + centroid.ToString());




                        }//using
                    }//end try
                    catch (System.Exception ex)
                    {
                        ed.WriteMessage("Error Reading Centroid: " + ex.Message);

                    }
                }// if

            }
        }



[Chumplybum]

thanks pkohut for your help / code...

i've amended my code to suit your changes (see below), however, the code only seems to produce an accurate result if the polyline is drawn counterclockwise (see attached images, the polyline has been mirrod to get counterclockwise vs clockwise)

Code: [Select]

       

imports AADS = AutoDesk.AutoCAD.DatabaseServices
imports AAG = Autodesk.AutoCAD.Geometry

Public Shared Function GetCentroid(ByVal polyline As AADS.Polyline) As AAG.Point3d

            Dim ReferencePoint As AAG.Point2d
            Dim RunningX As Double
            Dim RunningY As Double
            Dim second_factor As Double
            Dim Point1 As AAG.Point2d
            Dim Point2 As AAG.Point2d

            For Counter As Integer = 0 To polyline.NumberOfVertices - 2
                Point1 = polyline.GetPoint2dAt(Counter)
                Point2 = polyline.GetPoint2dAt(Counter + 1)

                If Counter = 0 Then ReferencePoint = Point1

                second_factor = (Point1.X - ReferencePoint.X) * (Point2.Y - ReferencePoint.Y) - (Point2.X - ReferencePoint.X) * (Point1.Y - ReferencePoint.Y)

                RunningX += ((Point1.X - ReferencePoint.X) + (Point2.X - ReferencePoint.X)) * second_factor
                RunningY += ((Point1.Y - ReferencePoint.Y) + (Point2.Y - ReferencePoint.Y)) * second_factor
            Next

            Dim X, Y As Double

            X = RunningX / (6 * polyline.Area) + ReferencePoint.X
            Y = RunningY / (6 * polyline.Area) + ReferencePoint.Y


            Return New AAG.Point3d(X, Y, polyline.Elevation)

        End Function


[pkohut]

The easiest way to fix this is to keep a running total of second_factor during the loop.  Then
once out of the loop if the summation is negative then you'll need to multiply your X and Y
calculations with -1.0.  FYI, the area of the polyline is the abs(summation) * 0.5   (straight
line segments only).

Cathy mentions possibly missing the final segment.  You'll need to check if the first and last
vertex are the same or not.  You already know the polyline is closed, so the final segment
may be interpolated across the final vertex to the first vertex, and you would miss that.

Paul

[pkohut]

Need to correct a sentence in my last post from:

Then
once out of the loop if the summation is negative then you'll need to multiply your X and Y
calculations with -1.0.

to:

Then once out of the loop if the summation is negative then you'll need to multiply
RunningX and RunningY with -1.0.  Here's an example:

dim dVal as double
dVal = 1.0
if summation < 0.0 then dVal = -1.0
    X = (RunningX * dVal) / (6 * polyline.Area) + ReferencePoint.X
    Y = (RunningY * dVal) / (6 * polyline.Area) + ReferencePoint.Y

Paul

[Chumplybum]

awesome... thanks pkohut, that fixed it up


cheers, Mark

[gile]

Hi,

It was quite difficult for me (I'm a complete beginer with .NET), but it was a fine exercice.
Here's a translation (C#) of the LISP routines I gave here.

All comments are welcome.

It works with arc polylines.

EDIT: added a testing command and a LISP function.
The "pline_centroid" command draws a point on the selected polyline centroid.
(pline_centroid (car (entsel))) returns the WCS coordinates of the centroid.

Code: [Select]

using System;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;

namespace PolylineCentroid
{
    // Main class
    public class Centroid
    {
        public Point3d GetCentroid(Polyline pl)
        {
            Point2d p0 = pl.GetPoint2dAt(0);
            Point2d cen = new Point2d(0.0, 0.0);
            double area = 0.0;
            double bulge = pl.GetBulgeAt(0);
            int last = pl.NumberOfVertices - 1;
            double tmpArea;
            Point2d tmpPoint;

            if (bulge != 0.0)
            {
                double[] datas = getArcGeom(pl, bulge, 0, 1);
                area = datas[0];
                cen = new Point2d(datas[1], datas[2]) * datas[0];
            }
            for (int i = 1; i < last; i++)
            {
                tmpArea = triangleAlgebricArea(p0, pl.GetPoint2dAt(i), pl.GetPoint2dAt(i + 1));
                tmpPoint = triangleCentroid(p0, pl.GetPoint2dAt(i), pl.GetPoint2dAt(i + 1));
                cen += (tmpPoint * tmpArea).GetAsVector();
                area += tmpArea;
                bulge = pl.GetBulgeAt(i);
                if (bulge != 0.0)
                {
                    double[] datas = getArcGeom(pl, bulge, i, i + 1);
                    area += datas[0];
                    cen += new Vector2d(datas[1], datas[2]) * datas[0];
                }
            }
            bulge = pl.GetBulgeAt(last);
            if (bulge != 0.0)
            {
                double[] datas = getArcGeom(pl, bulge, last, 0);
                area += datas[0];
                cen += new Vector2d(datas[1], datas[2]) * datas[0];
            }
            cen = cen.DivideBy(area);
            Point3d result = new Point3d(cen.X, cen.Y, pl.Elevation);
            return result.TransformBy(Matrix3d.PlaneToWorld(pl.Normal));
        }

        public double[] GetArcGeom(Polyline pl, double bulge, int index1, int index2)
        {
            CircularArc2d arc = (pl.GetArcSegment2dAt(index1));
            double arcRadius = arc.Radius;
            Point2d arcCenter = arc.Center;
            double arcAngle = 4.0 * Math.Atan(bulge);
            double tmpArea = arcAlgebricArea(arcRadius, arcAngle);
            Point2d tmpPoint = arcCentroid(pl.GetPoint2dAt(index1), pl.GetPoint2dAt(index2), arcCenter, tmpArea);
            return new double[3] { tmpArea, tmpPoint.X, tmpPoint.Y };
        }

        public Point2d TriangleCentroid(Point2d p0, Point2d p1, Point2d p2)
        {
            return (p0 + p1.GetAsVector() + p2.GetAsVector()) / 3.0;
        }

        public double TriangleAlgebricArea(Point2d p0, Point2d p1, Point2d p2)
        {
            return (((p1.X - p0.X) * (p2.Y - p0.Y)) - ((p2.X - p0.X) * (p1.Y - p0.Y))) / 2.0;
        }

        public Point2d ArcCentroid(Point2d start, Point2d end, Point2d cen, double tmpArea)
        {
            double chord = start.GetDistanceTo(end);
            double angle = angleFromTo(start, end);
            return polar2d(cen, angle - (Math.PI / 2.0), (chord * chord * chord) / (12.0 * tmpArea));
        }

        public double ArcAlgebricArea(double rad, double ang)
        {
            return rad * rad * (ang - Math.Sin(ang)) / 2.0;
        }

        public double AngleFromTo(Point2d p1, Point2d p2)
        {
            return (p2 - p1).Angle;
        }

        public Point2d Polar2d(Point2d org, double angle, double distance)
        {
            return new Point2d(org.X + distance, org.Y).RotateBy(angle, org);
        }
    }

    // Testing command
    public class MyClass
    {
        [CommandMethod("pline_centroid")]
        public void centroid()
        {
            Document doc = Application.DocumentManager.MdiActiveDocument;
            Database db = doc.Database;
            Editor ed = doc.Editor;

            PromptEntityOptions opts = new PromptEntityOptions("\nSelect a polyline: ");
            opts.AllowNone = true;
            opts.AllowObjectOnLockedLayer = true;
            opts.SetRejectMessage("\nEntité non valide.");
            opts.AddAllowedClass(typeof(Polyline), false);
            PromptEntityResult pline = ed.GetEntity(opts);
            if (pline.Status == PromptStatus.OK)
            {
                ObjectId ObjID = pline.ObjectId;
                try
                {
                    using (Transaction trans = db.TransactionManager.StartTransaction())
                    {
                        Polyline ent = (Polyline)trans.GetObject(ObjID, OpenMode.ForRead, false);
                        Centroid c = new Centroid();
                        Point3d cen = c.GetCentroid(ent);
                        DBPoint pt = new DBPoint(cen);
                        BlockTable bTable = (BlockTable)trans.GetObject(db.BlockTableId, OpenMode.ForRead);
                        BlockTableRecord mSpace = (BlockTableRecord)trans.GetObject(bTable[BlockTableRecord.ModelSpace], OpenMode.ForWrite);
                        mSpace.AppendEntity(pt);
                        trans.AddNewlyCreatedDBObject(pt, true);
                        trans.Commit();
                    }
                }
                catch (System.Exception ex)
                {
                    ed.WriteMessage("Error: " + ex.Message);
                }
            }
        }
    }

    // LISP function
    public class LispCentoid
    {
        [LispFunction("pline_centroid")]
        public static object centroid(ResultBuffer buff)
        {
            Document doc = Application.DocumentManager.MdiActiveDocument;
            Database db = doc.Database;
            Editor ed = doc.Editor;

            if (buff == null)
            {
                ed.WriteMessage("Error: none argument.\n");
                return false;
            }

            TypedValue[] args = buff.AsArray();
            if (args[0].TypeCode != (int)LispDataType.ObjectId)
            {
                ed.WriteMessage("Error: incorrect argument type: lentityp {0}.\n", args[0].Value);
                return false;
            }
            ObjectId objID = (ObjectId)args[0].Value;
            try
            {
                using (Transaction trans = db.TransactionManager.StartTransaction())
                {
                    Entity ent = (Entity)trans.GetObject(objID, OpenMode.ForRead, false);
                    string entType = ent.GetType().ToString().Substring(34);
                    if (entType == "Polyline")
                    {
                        Polyline pl = (Polyline)trans.GetObject(objID, OpenMode.ForRead, false);
                        Centroid c = new Centroid();
                        return c.GetCentroid(pl);
                    }
                    else
                    {
                        if (entType.Contains("DB")) entType = entType.Substring(2);
                        ed.WriteMessage("Error: incorrect entity type: " + entType + "\n");
                        return null;
                    }
                }
            }
            catch (System.Exception ex)
            {
                ed.WriteMessage("\nError: " + ex.Message + "\n");
                return false;
            }
        }
    }
}

[gile]

I edited the Centroid class code so that some methods, as angleFromTo, polar2d, triangleAlgebricArea, triangleCentroid and so on, which may be usefull are more easily readable.

[It's Alive!]

Excellent gile!

[gile]

Thanks Daniel..

[Chumplybum]

nice one, thanks gile

[pkohut]

Ya Gile looks good.  Since your a self admitted beginner at .Net
may I suggest making all those "int" doubles into doubles
instead of letting the compiler upcast them for you. 
return (p0 + p1.GetAsVector() + p2.GetAsVector()) / 3
to
return (p0 + p1.GetAsVector() + p2.GetAsVector()) / 3.0

if(bulge != 0)
to
if(bulge != 0.0)

The compiler will happily up cast for you, but that isn't something
you should let it do.  Make your intent known to the compiler
and it will help you when you make a mistake.  It also helps
you and others with a visual clue about your intent.

Later,
Paul

[gile]

Thanks Paul,

It makes sense, I edit the code according to your suggestion.