Convex Hull

[(gile)]

Here're a C# and a F# implementation of the 'Graham's scan' algorythm to get the convex hull of a 2d points collection.

[EDIT: added a C# implementation using Linq]

C# (targeting NET Framework 2.0 which is the default for A2007 -> A2009).

Code: [Select]

using System.Collections.Generic;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
using AcAp = Autodesk.AutoCAD.ApplicationServices.Application;

namespace ConvexHull
{
    public class Commands
    {
        private Point2d _p0;

        private bool Clockwise(Point2d p1, Point2d p2, Point2d p3)
        {
            return ((p2.X - p1.X) * (p3.Y - p1.Y) - (p2.Y - p1.Y) * (p3.X - p1.X)) < 1e-9;
        }

        private int ComparePoints(Point2d p1, Point2d p2)
        {
            if (p1.IsEqualTo(p2)) return 0;
            double d1 = _p0.GetDistanceTo(p1);
            double d2 = _p0.GetDistanceTo(p2);
            if (d1 == 0.0) return -1;
            if (d2 == 0.0) return 1;
            double cos = (p2.X - _p0.X) / d2 - (p1.X - _p0.X) / d1;
            if (cos < -1e-9) return -1;
            if (cos > 1e-9) return 1;
            return d1.CompareTo(d2);
        }

        private List<Point2d> ConvexHull(List<Point2d> pts)
        {
            _p0 = pts[0];
            for (int i = 1; i < pts.Count; i++)
            {
                Point2d pt = pts[i];
                if (pt.Y < _p0.Y || (pt.Y == _p0.Y && pt.X < _p0.X))
                    _p0 = pt;
            }
            pts.Sort(ComparePoints);
            for (int i = 1; i < pts.Count - 1; i++)
            {
                while (i > 0 && Clockwise(pts[i - 1], pts[i], pts[i + 1]))
                {
                    pts.RemoveAt(i);
                    i--;
                }
            }
            return pts;
        }

        [CommandMethod("ch", CommandFlags.UsePickSet)]
        public void testCh()
        {
            Document doc = AcAp.DocumentManager.MdiActiveDocument;
            Database db = doc.Database;
            Editor ed = doc.Editor;
            TypedValue[] filter = new TypedValue[1] { new TypedValue(0, "POINT") };
            PromptSelectionResult psr = ed.GetSelection(new SelectionFilter(filter));
            if (psr.Status != PromptStatus.OK) return;
            using (Transaction tr = db.TransactionManager.StartTransaction())
            using (Polyline pline = new Polyline())
            {
                List<Point2d> pts = new List<Point2d>();
                foreach (SelectedObject so in psr.Value)
                {
                    DBPoint dbPt = (DBPoint)tr.GetObject(so.ObjectId, OpenMode.ForRead);
                    pts.Add(new Point2d(dbPt.Position.X, dbPt.Position.Y));
                }
                for (int i = 0; i < ConvexHull(pts).Count; i++)
                {
                    pline.AddVertexAt(i, pts[i], 0.0, 0.0, 0.0);
                }
                pline.Closed = true;
                pline.SetDatabaseDefaults();
                BlockTableRecord btr = (BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
                btr.AppendEntity(pline);
                tr.AddNewlyCreatedDBObject(pline, true);
                tr.Commit();
            }
        }
    }
}
C# (using Linq, NET Framework 3 or upper required)

Code: [Select]

using System;
using System.Collections.Generic;
using System.Linq;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
using AcAp = Autodesk.AutoCAD.ApplicationServices.Application;

namespace LinqConvexHull
{
    public class Commands
    {
        private Point2d _p0;

        private bool Clockwise(Point2d p1, Point2d p2, Point2d p3)
        {
            return ((p2.X - p1.X) * (p3.Y - p1.Y) - (p2.Y - p1.Y) * (p3.X - p1.X)) < 1e-9;
        }

        private double Cosine(Point2d pt)
        {
            double d = _p0.GetDistanceTo(pt);
            return d == 0.0 ? 1.0 : Math.Round((pt.X - _p0.X) / d, 9);
        }

        private List<Point2d> ConvexHull(List<Point2d> pts)
        {
            _p0 = pts.OrderBy(p => p.Y).ThenBy(p => p.X).First();
            pts = pts.OrderByDescending(p => Cosine(p)).ThenBy(p => _p0.GetDistanceTo(p)).ToList();
            for (int i = 1; i < pts.Count - 1; i++)
            {
                while (i > 0 && Clockwise(pts[i - 1], pts[i], pts[i + 1]))
                {
                    pts.RemoveAt(i);
                    i--;
                }
            }
            return pts;
        }

        [CommandMethod("Test")]
        public void Test()
        {
            Document doc = AcAp.DocumentManager.MdiActiveDocument;
            Database db = doc.Database;
            Editor ed = doc.Editor;
            TypedValue[] filter = new TypedValue[1] { new TypedValue(0, "POINT") };
            PromptSelectionResult psr = ed.GetSelection(new SelectionFilter(filter));
            if (psr.Status != PromptStatus.OK) return;
            using (Transaction tr = db.TransactionManager.StartTransaction())
            using (Polyline pline = new Polyline())
            {
                List<Point2d> pts = new List<Point2d>();
                foreach (SelectedObject so in psr.Value)
                {
                    DBPoint dbPt = (DBPoint)tr.GetObject(so.ObjectId, OpenMode.ForRead);
                    pts.Add(new Point2d(dbPt.Position.X, dbPt.Position.Y));
                }
                pts = ConvexHull(pts);
                for (int i = 0; i < pts.Count; i++)
                {
                    pline.AddVertexAt(i, pts[i], 0.0, 0.0, 0.0);
                }
                pline.Closed = true;
                pline.SetDatabaseDefaults();
                BlockTableRecord btr = (BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
                btr.AppendEntity(pline);
                tr.AddNewlyCreatedDBObject(pline, true);
                tr.Commit();
            }
        }
    }
}

F#

Code: [Select]

module ConvexHull

open Autodesk.AutoCAD.ApplicationServices
open Autodesk.AutoCAD.DatabaseServices
open Autodesk.AutoCAD.EditorInput
open Autodesk.AutoCAD.Geometry
open Autodesk.AutoCAD.Runtime

type AcAp = Autodesk.AutoCAD.ApplicationServices.Application

let clockwise (p1:Point2d) (p2:Point2d) (p3:Point2d) =
    (p2.X - p1.X) * (p3.Y - p1.Y) - (p2.Y - p1.Y) * (p3.X - p1.X) < 1e-9
   
let convexHull (pts : Point2d seq) =
    let rec fill acc pt =
        match acc with
        | a :: b :: _ when clockwise b a pt -> fill acc.Tail pt
        | _ -> pt :: acc
    let p0 = pts
             |> Seq.reduce (fun p1 p2 ->
                if p2.Y < p1.Y || (p1.Y = p2.Y && p2.X < p1.X) then p2 else p1)
    pts
    |> List.ofSeq
    |> List.sortBy (fun p ->
        let d = p0.GetDistanceTo(p)
        (Math.Round((p0.X - p.X) / d, 8), d))
    |> List.fold fill []
    |> List.rev

[<CommandMethod("ch")>]
let Test() =
    let doc = AcAp.DocumentManager.MdiActiveDocument
    let db = doc.Database
    let ed = doc.Editor
    let psr = ed.GetSelection(new SelectionFilter([| new TypedValue(0, "POINT") |]))
    if psr.Status = PromptStatus.OK then
        use tr = db.TransactionManager.StartTransaction()
        use pl = new Polyline()
        psr.Value
        |> Seq.cast<_>
        |> Seq.map (fun (so : SelectedObject) ->
            let pt = tr.GetObject(so.ObjectId, OpenMode.ForRead) :?> DBPoint
            new Point2d(pt.Position.X, pt.Position.Y))
        |> convexHull
        |> List.fold(fun i p ->  pl.AddVertexAt(i, p, 0.0, 0.0, 0.0); i + 1) 0
        |> ignore
        pl.Closed <- true
        pl.SetDatabaseDefaults()
        let btr = tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite) :?> BlockTableRecord
        btr.AppendEntity(pl) |> ignore
        tr.AddNewlyCreatedDBObject(pl, true)
        tr.Commit()
target audience:{advanced}

[rom1]

Hi,

Thanks, this function seems to be very interesting.

I guess it's useful for topology plan? There is no similar function natively in AutoCAD?

Thanks

[(gile)]

Hi,

There's no similar function in AutoCAD, maybe in MAP or Civil, I don't know.

I thaught it was intersting to show how to use a functional style with F# and an imperative* one with C# (or VB) to solve the same problem.

* even if the use of a delegate with the Sort method is more a functional feature.

[(gile)]

I added a C# code using Linq: a more 'functional style' with C# (can be done with VB too).

[fixo]

Another convexhull code, not sure about if this
working good

Code: [Select]

        // See for more
       // https://github.com/lein/ConvexHull
       // http://www.chrisharrison.net/projects/convexHull/jarvis.cpp
        [CommandMethod("cxx")]
        public void testConvexHull()
        {
            Document doc = Autodesk.AutoCAD.ApplicationServices.Application.DocumentManager.MdiActiveDocument;
            Database db = doc.Database;
            Editor ed = doc.Editor;
            testConvexHull(db, ed);
        }

        public static List<Point2d> ConvexHull(List<Point2d> points) 
        {     
            if (points.Count < 3) 
            {       
                throw new ArgumentException("At least 3 points reqired", "points"); 
 
            }       
            List<Point2d> convex = new List<Point2d>();   
            // get leftmost point   
            Point2d ptOnHull = points.Where(p => p.X == points.Min(min => min.X)).First();     
            Point2d ptNext;     
            do         
            {
                convex.Add(ptOnHull);     
                ptNext = points[0];     
                for (int i = 1; i < points.Count; i++)       
                {             
                    if ((ptOnHull == ptNext)
                        || (CCW(ptOnHull, ptNext, points[i]) == -1))    // check direction
                    {                     
                        ptNext = points[i];     
                    }           
                }           
                ptOnHull = ptNext; 
            }
            while (ptNext != convex[0]);
            return convex;     
        }   

        private static int CCW(Point2d p1, Point2d p2, Point2d p)     
        {           
            // Determinant   
            int ptDir = Convert.ToInt32((p2.X - p1.X) * (p.Y - p1.Y) - (p.X - p1.X) * (p2.Y - p1.Y));
            if (ptDir > 0)       
                return -1;  // clockwise   
            if (ptDir < 0)           
                return 1; // counterclockwise 
            return 0; //  collinear direction   
        }

        //--------------------------------------------------------------------------//

        public void testConvexHull(Database db, Editor ed)
        {
            List<Point2d> points = new List<Point2d>();
            List<Point2d> convex = new List<Point2d>();
            Point2d pt;
            Matrix3d ucs = ed.CurrentUserCoordinateSystem;
            using (Transaction tr = db.TransactionManager.StartTransaction())
            {

                BlockTableRecord btr = (BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
                TypedValue[] values = { new TypedValue(0, "point,line,lwpolyline,circle") };// add other types you need here
                SelectionFilter filter = new SelectionFilter(values);

                PromptSelectionResult res = ed.GetSelection(filter);

                if (res.Status != PromptStatus.OK) return;

                SelectionSet sset = res.Value;
                //   iterate through selected objects
                foreach (SelectedObject selobj in sset)
                {
                    Entity pEnt = (Entity)tr.GetObject(selobj.ObjectId, OpenMode.ForRead, false);

                    if (pEnt == null) return;
                    //----------------------------DBPOint----------------------------//
                    if (pEnt is DBPoint)
                    {
                        DBPoint dp = pEnt as DBPoint;
                        if (dp != null)
                        {

                            Point3d pp = dp.Position.TransformBy(ucs);
                            pt = new Point2d(pp.X, pp.Y);
                            points.Add(pt);
                        }
                    }
                    //----------------------------Line----------------------------//
                    if (pEnt is Line)
                    {
                        Line ln = pEnt as Line;
                        if (ln != null)
                        {

                            Point3d pp = ln.StartPoint.TransformBy(ucs);
                            pt = new Point2d(pp.X, pp.Y);
                            points.Add(pt);
                            pp = ln.EndPoint.TransformBy(ucs);
                            pt = new Point2d(pp.X, pp.Y);
                            points.Add(pt);
                        }
                    }

                    //---------------------------Polyline----------------------------//
                    if (pEnt is Polyline)
                    {
                        Polyline poly = pEnt as Polyline;
                        if (poly != null)
                        {
                            for (int m = 0; m < poly.NumberOfVertices; m++)
                            {
                                Point3d pp = poly.GetPoint3dAt(m).TransformBy(ucs);
                                pt = new Point2d(pp.X, pp.Y);
                                points.Add(pt);
                            }

                        }
                    }
                    //---------------------------Circle----------------------------//
                    if (pEnt is Circle)
                    {
                        Circle circ = pEnt as Circle;
                        if (circ != null)
                        {
                            int num = 64;
                            double dp = (circ.EndParam - circ.StartParam) / num;
                            for (int a = 0; a < num; a++)
                            {
                                Point3d pp = circ.GetPointAtParameter(dp * a).TransformBy(ucs);
                                pt = new Point2d(pp.X, pp.Y);
                                points.Add(pt);
                            }

                        }
                    }
                }
                convex=ConvexHull(points);
                ed.WriteMessage("\nAll Points Count:\t{0}\n", points.Count.ToString());
                ed.WriteMessage("\nConvexHull points Count:\t{0}\n", convex.Count.ToString());
                //-------------------------------------------------------//
                // Create a polyline by boundary points
                Polyline pline = new Polyline();
                pline.SetDatabaseDefaults();

                for (int cnt = 0; cnt < convex.Count; cnt++)
                {

                    pline.AddVertexAt(cnt, convex[cnt], 0, 0, 0);

                }

                // close boundary
                pline.Closed = true;
                pline.ColorIndex = 121;

                // Add the new object to the block table record and the transaction
                btr.AppendEntity(pline);
                tr.AddNewlyCreatedDBObject(pline, true);
                tr.Commit();

            }
        }