Friday, June 20, 2014

Finding Polygon Intersection By Using MS SQL Server

When doing AutoCAD programming, especially when the target AutoCAD vertical is AutoCAD Map 3D/Civil 3D, we often need to find out if one polygon intersects with another polygon, and if yes, then we often need to get the derived polygon because of the intersection.

There seems there is no easy, ready to call API support in plain either AutoCAD, or AutoCAD Map/Civil.

There are quite a few well-known algorithms that can be found in the Internet, some even has code implementation, or available as compiled library.

However, if the AutoCAD working environment has a accessible MS SQL Server (2008 or later version), it would be very simple to take advantage of SQL Server's spatial capability to fulfill this task fairly simple and easy.

This article demonstrates how to use MS SQL Server's spatial functionality to calculate polygon intersection.

Since MS SQL Server's spatial capabilities are actually "add-in" to SQL Server. so the spatial capabilities can be used with or without SQL Server being involved. That is, our code can be set to direct reference to related DLLs and do the calculation; or our code can access SQL Server with necessary input and let the SQL Server do the calculation, then obtain the calculation result from the SQL Server.

Firstly, because I need to use 2 different ways to take advantage of SQL Server's spatial functionality, I decide to define an interface for AutoCAD to call. With this interface, if I later decide to use other polygon intersection algorithms to do the calculation, I can simply implement this interface, so that the code on AutoCAD side would not need change no matter what method is used to find polygon intersection. Here is the simple interface code:

using Autodesk.AutoCAD.Geometry;
namespace PolygonIntersection
{
    public interface IPolyIntersecCalculator
    {
        Point2dCollection GetPolygonIntersection(
            Point2dCollection polygonA, Point2dCollection polygonB);
    } 
}

This interface defines an method GetPolygonIntersection(), which effectively says: "Give me 2 collections of 2D points, representing 2 polygons, I'll return either a collection of 2D points, representing the intersection area, if the 2 polygons intersect to each other; or nothing if the 2 polygons do not intersect."

1. Use SQL Server to do the polygon intersection calculation

Now, let's see how to make SQL Server to do the calculation of finding polygon intersection. Here the code uses ADO.NET to directly access SQL Server to call a stored procedure, which is executed inside SQL Server. The input to the stored procedure call is 2 text values, representing 2 polygons (WKT - Well Known Text). Then SQL Server returns calculation result, which is also represented as WKT, back to the code in AutoCAD.

Here is the stored procedure that stored in a SQL Server database:

-- ============================================================================
-- Author:  Norman Yuan
-- Create date: 2014-05-12
-- Description: This stored procedure uses Sql Server's
--  spatial computing capability to find out polygon's intersection. 
--  There is no data involved in an data table.
--  Inputs are 2 polygons described as WKT (Well-known text);
--  Output is WKT of geometry, being POLYGON, POINT, or 
--  GEOMETRYCOLLECTION EMPTY
-- ============================================================================
ALTER PROCEDURE [dbo].[SpatialComputing_GetPolygonIntersection]
 @PolygonWKT1 nvarchar(max), 
 @PolygonWKT2 nvarchar(max)
AS
BEGIN
 
 SET NOCOUNT ON;
 
 --Test data
 /*
 --Intersecting polygons, returns 'POLYGON((...))'
 SET @PolygonWKT1='POLYGON((0 0, 2 0, 3 3, 0 2, 0 0))';
 SET @PolygonWKT2='POLYGON((1 1, 3 1, 3 3, 1 3, 1 1))';
 */
 
 /*
 --polygons that do not intersecting to each other,
 --but has overlapped point, returns 'POINT(...)'
 SET @PolygonWKT1='POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))';
 SET @PolygonWKT2='POLYGON((1 1, 2 1, 2 2, 1 2, 1 1))';
 */
 
 /*
 --polygons that do not intersecting to each other at all,
 --returns 'GEOMETRYCOLLECTION EMPTY'
 SET @PolygonWKT1='POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))';
 SET @PolygonWKT2='POLYGON((2 2, 3 2, 3 3, 2 3, 2 2))';
 */
 
 DECLARE @Polygon1 geometry
 SET @Polygon1 = geometry::STPolyFromText(@polygonWKT1,0)
 
 DECLARE @Polygon2 geometry
 SET @Polygon2 = geometry::STPolyFromText(@polygonWKT2,0)
 
 DECLARE @OutputPolygon geometry
 SET @OutputPolygon = @Polygon1.STIntersection(@Polygon2)
 
 DECLARE @OutputWKT geometry
 SET @OutputWKT = @OutputPolygon.ToString()
 
 SELECT @OutputWKT.ToString() AS OutputWKT
    
END

With the stored procedure created in a SQL Server database (and proper execution permission assigned), I can now write some simple ADO.NET code to access SQL Server and call this stored procedure easily.

Here is class PolyIntersecSqlServerBackEndCalculator, which implements IPolyIntersecCalculator:

using System;
using System.Text;
using System.Data;
using System.Data.SqlClient;
using Autodesk.AutoCAD.Geometry;
 
namespace PolygonIntersection
{
    public class PolyIntersecSqlServerBackEndCalculator : 
        IPolyIntersecCalculator
    {
        private string _connectionString;
        private string _procedureName;
        private string _inputParam1;
        private string _inputParam2;
        public PolyIntersecSqlServerBackEndCalculator(
            string connectionString, 
            string procedureName,
            string inputParamName1,
            string inputParamName2)
        {
            _connectionString = connectionString;
            _procedureName = procedureName;
            _inputParam1 = inputParamName1.StartsWith("@") ? 
                inputParamName1 : "@" + inputParamName1;
            _inputParam2 = inputParamName2.StartsWith("@") ? 
                inputParamName2 : "@" + inputParamName2;
        }
 
        public Point2dCollection GetPolygonIntersection(
            Point2dCollection polygonA, Point2dCollection polygonB)
        {
            return FindIntersectionPolygon(polygonA, polygonB);
        }
 
        #region private methods
 
        private Point2dCollection FindIntersectionPolygon(
            Point2dCollection polygonA, Point2dCollection polygonB)
        {
            Point2dCollection points = null;
 
            string wkt1 = GetPolygonWKTFromPoints(polygonA);
            string wkt2 = GetPolygonWKTFromPoints(polygonB);
 
            using (SqlConnection cn = new SqlConnection(_connectionString))
            {
                using (SqlCommand cmd = cn.CreateCommand())
                {
                    cmd.CommandType = CommandType.StoredProcedure;
                    cmd.CommandText = _procedureName;
 
                    cmd.Parameters.AddWithValue(_inputParam1, wkt1);
                    cmd.Parameters.AddWithValue(_inputParam2, wkt2);
 
                    try
                    {
                        cn.Open();
                        string outputWKT = (string)cmd.ExecuteScalar();
                        points = ParseOutputWKT(outputWKT);
                    }
                    catch
                    {
                        throw;
                    }
                    finally
                    {
                        if (cn.State == ConnectionState.Open) cn.Close();
                    }
                }
            }
 
            return points;
        }
 
        private string GetPolygonWKTFromPoints(Point2dCollection points)
        {
            StringBuilder coords = new StringBuilder();
            foreach (Point2d p in points)
            {
                coords.Append(p.X.ToString() + " " + p.Y.ToString() + ", ");
            }
 
            //Add the first point as the WKT's last point, indicaating
            //the geometry is a closed polygon
            coords.Append(
                points[0].X.ToString() + " " + points[0].Y.ToString());
 
            return "POLYGON((" + coords.ToString() + "))";
        }
 
        private Point2dCollection ParseOutputWKT(string wkt)
        {
            if (wkt.ToUpper().StartsWith("POLYGON"))
            {
                Point2dCollection points = new Point2dCollection();
 
                int pos = wkt.IndexOf("((");
                string temp = wkt.Substring(pos + 2).Replace("))""");
                string[] coords = temp.Trim().Split(',');
                foreach (var coord in coords)
                {
                    string[] val = coord.Trim().Split(' ');
 
                    Point2d p = new Point2d(
                        Convert.ToDouble(val[0]), Convert.ToDouble(val[1]));
                    points.Add(p);
                }
 
                return points;
            }
            else
            {
                return null;
            }
        }
 
        #endregion
    }
}

As you can see, the code does not do any calculation on polygon intersecting. It simply interprets the input and output to/from SQL Server from WKT into Point2dCollection. Naturally, if the SQL Server access is exposed through a service layer, I can simply swap the ADO.NET data access code with service access code. The point is that the polygon intersection calculation is done on the SQL Server side.

2. Use Microsoft.SqlServer.Types.dll to do polygon intersection calculation

We can directly call a .NET API exposed by Microsoft.SqlServer.Types.dll to calculate polygon intersection without SQL Server involved at all. This is because that SQL Server's spatial functionality is implemented independently from SQL Server. SQL Server uses it as add-in.

To use Microsoft.SqlServer.Types.dll, you need to add reference to it in the Visual Studio project:



More explanation later on setting reference to Microsoft.SqlServer.Types.dll and making it work outside SQL Server.

With Microsoft.SqlServer.Types namespace available, I created another class that implements IPolyIntersecCalculator, called PolyIntersecSqlServerSpatialLibCalculator:

using System;
using Microsoft.SqlServer.Types;
using Autodesk.AutoCAD.Geometry;
 
namespace PolygonIntersection
{
    public class PolyIntersecSqlServerSpatialLibCalculator :
        IPolyIntersecCalculator
    {
        public Point2dCollection GetPolygonIntersection(
            Point2dCollection polygonA, Point2dCollection polygonB)
        {
            return FindIntersectionPolygon(polygonA, polygonB);
        }
 
        #region private methods
 
        private Point2dCollection FindIntersectionPolygon(
            Point2dCollection polygonA, Point2dCollection polygonB)
        {
            Point2dCollection points = new Point2dCollection();
 
            SqlGeometry polygon1 = 
                CreateSqlPolygonGeometryFromPoint2ds(polygonA);
            polygon1.MakeValid();
 
            SqlGeometry polygon2 = 
                CreateSqlPolygonGeometryFromPoint2ds(polygonB);
            polygon2.MakeValid();
 
            //Debugging code
            string wktP1 = polygon1.ToString();
            string wktP2 = polygon2.ToString();
 
            SqlGeometry area = polygon1.STIntersection(polygon2);
            if (!area.IsNull)
            {
                if (area.STGeometryType().ToString().ToUpper() == "POLYGON")
                {
                    for (int i = 1; i <= area.STNumPoints(); i++)
                    {
                        SqlGeometry geoPoint = 
                            area.STPointN(Convert.ToInt32(i));
 
                        double x = geoPoint.STX.Value;
                        double y = geoPoint.STY.Value;
                        Point2d pt = new Point2d(x, y);
                        points.Add(pt);
 
                    }
                }
            }
 
            return points;
        }
 
        private static SqlGeometry CreateSqlPolygonGeometryFromPoint2ds(
            Point2dCollection points)
        {
            SqlGeometry geometry = null;
 
            SqlGeometryBuilder gb = new SqlGeometryBuilder();
            gb.SetSrid(0);
            gb.BeginGeometry(OpenGisGeometryType.Polygon);
 
            gb.BeginFigure(points[0].X, points[0].Y);
 
            for (int i = 1; i <= points.Count; i++)
            {
                if (i < points.Count)
                {
                    gb.AddLine(points[i].X, points[i].Y);
                }
                else
                {
                    gb.AddLine(points[0].X, points[0].Y);
                }
            }
 
            gb.EndFigure();
 
            gb.EndGeometry();
            geometry = gb.ConstructedGeometry;
 
            return geometry;
        }
 
        #endregion
    }
}

Following CommandClass MyCadCommands uses the 2 IPolyIntersecCalculator classes in action:

using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.Geometry;
using CadApp = Autodesk.AutoCAD.ApplicationServices.Application;
using CadDb = Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.GraphicsInterface;
 
[assemblyCommandClass(typeof(PolygonIntersection.MyCadCommands))]
 
namespace PolygonIntersection
{
    public class MyCadCommands
    {
        [CommandMethod("PolyIntersec")]
        public static void GetPolygonIntersection()
        {
            Document dwg = CadApp.DocumentManager.MdiActiveDocument;
            Editor ed = dwg.Editor;
 
            try
            {
                ObjectId polyId1 = PickPolygon(
                    ed, "Select a closed polyline:");
                if (polyId1.IsNull)
                {
                    ed.WriteMessage("\n*Cancel*");
                    return;
                }
 
                ObjectId polyId2 = PickPolygon(
                    ed, "Select another closed polyline:");
                if (polyId2.IsNull)
                {
                    ed.WriteMessage("\n*Cancel*");
                    return;
                }
 
                Point2dCollection polygonA = GetPolylineVertices(polyId1);
                Point2dCollection polygonB = GetPolylineVertices(polyId2);
 
                IPolyIntersecCalculator calculator = null;
 
                //Use SqlServerSpatial110.dll and 
                //its .NET wrapper Microsoft.SqlServer.Types.dll
                calculator = new PolyIntersecSqlServerSpatialLibCalculator();
 
                Point2dCollection interPolygon1 = 
                    calculator.GetPolygonIntersection(polygonA, polygonB);
 
                if (interPolygon1 != null)
                {
                    ShowIntersectionArea(interPolygon1, 
                        "Intersectiion area found by " +
                        "SqlServerSpatialLib calculator.");
                    ed.UpdateScreen();
                }
 
                //Use SqlServer DB backend
                string connection = "User ID=XXXXXX;" +
                    "Password=xxxx;Initial Catalog=" +
                    "[DatabaseName];Data Source=[ServerName]";
                string procedureName = 
                    "SpatialComputing_GetPolygonIntersection";
                string param1 = "PolygonWKT1";
                string param2 = "PolygonWKT2";
                calculator = new PolyIntersecSqlServerBackEndCalculator(
                    connection, procedureName, param1, param2);
 
                Point2dCollection interPolygon2 =
                    calculator.GetPolygonIntersection(polygonA, polygonB);
 
                if (interPolygon2 != null)
                {
                    ShowIntersectionArea(interPolygon2, 
                        "Intersection area found by " +
                        "SqlServerBackEnd calculator.");
                    ed.UpdateScreen();
                }
            }
            catch (System.Exception ex)
            {
                ed.WriteMessage("\nError: {0}\n{1}", ex.Message, ex.StackTrace);
            }
            finally
            {
                Autodesk.AutoCAD.Internal.Utils.PostCommandPrompt();
            }
        }
 
        #region private methods
 
        private static ObjectId PickPolygon(Editor ed, string message)
        {
            PromptEntityOptions opt = new PromptEntityOptions(
                "\n" + message);
            opt.SetRejectMessage("\nLnvalid pick: must be a polyline.");
            opt.AddAllowedClass(typeof(CadDb.Polyline), true);
            PromptEntityResult res = ed.GetEntity(opt);
            if (res.Status==PromptStatus.OK)
            {
                return res.ObjectId;
            }
            else
            {
                return ObjectId.Null;
            }
        }
 
        private static Point2dCollection GetPolylineVertices(
            ObjectId polylineId)
        {
            Point2dCollection points = new Point2dCollection();
 
            using (Transaction tran=
                polylineId.Database.TransactionManager.
                StartOpenCloseTransaction())
            {
                CadDb.Polyline pl = (CadDb.Polyline)tran.GetObject(
                    polylineId, OpenMode.ForRead);
 
                for (int i = 0; i < pl.NumberOfVertices; i++)
                {
                    points.Add(pl.GetPoint2dAt(i));
                }
                    
                tran.Commit();
            }
 
            return points;
        }
 
        private static void ShowIntersectionArea(
            Point2dCollection points, string msg)
        {
            Editor ed = CadApp.DocumentManager.MdiActiveDocument.Editor;
 
            using (CadDb.Polyline pl = CreatePolyline(points))
            {
                TransientManager ts = TransientManager.CurrentTransientManager;
 
                ts.AddTransient(
                    pl,
                    TransientDrawingMode.DirectTopmost,
                    128,
                    new IntegerCollection());
 
                ed.UpdateScreen();
 
                ed.GetString("\n" + msg + " Press any key to continue...");
 
                ts.EraseTransient(pl, new IntegerCollection());
            }
        }
 
        private static CadDb.Polyline CreatePolyline(Point2dCollection points)
        {
            CadDb.Polyline pl = new CadDb.Polyline(points.Count);
            for (int i=0; i<points.Count; i++)
            {
                pl.AddVertexAt(i, points[i], 0.0, 0.0, 0.0);
            }
 
            pl.Closed = true;
            pl.ColorIndex = 1;  //Show it in red
 
            return pl;
        }
 
        #endregion
    }
}

This video clip shows how the code works.

Now more about how to get Microsoft.SqlServer.Types.dll.

As a developer who does .NET programming, it is very likely the computer has a version of Visual Studio installed, and often also has a version of SQL Server (possible SQL Server Express) installed. In this case, Visual Studio's adding references dialog box should show Microsoft.SqlServer.Types in the list, as the picture above shows. However Microsoft.SqlServer.Types.dll is only a .NET assembly that wraps some native code that does the real heavy lifting work. That means, the second way of calculating polygon intersection works only if the computer has SQL Server (2008 or later) installed. Since Microsoft.SqlServer.Types.dll get installed into GAC due to SQL Server installation, in the Visual Studio project, the reference's "Copy Local" can be set to "False" (it is default to "False" anyway).

In reality, not all CAD computers would have a version of SQL Server installed for sure. In this case, we need to have a copy of Microsoft.SqlServer.Types.dll and the native DLL it wraps. What it is, then? It turned out that SqlServerSpatial.dll (SQL Server2008) or SqlServerSpatial110.dll (SQL Server 2012. I uses SqlServerSpatial110.dll here). To get both Microsoft.SqlServer.Types.dll and SqlServerSpatial110.dll installed to the computer, one can download Microsoft SQL Server 2012 Feature Pack, or find a copy of SqlServerSpatial11.dll in a computer that has SQL Server 2012 installed (in C:\Windows.System32 folder).

As for possible license issue of using these 2 DLLs, it is beyond the technical discussion and I'll leave it to user's discretion.

With these 2 DLLs available, make sure them to be build into the output of the project. That is, for the reference of Microsoft.SqlServer.Types.dll, its "Copy Local" should be set to "True"; and for the SqlServerSpatial110.dll, add it into the project and set "Copy to Output Directory" to "Copy always". See picture below:



So, after building the project, besides the AutoCAD add-in DLL in the bin folder, there are also both Microsoft.SqlServer.Types.dll and SqlServerSpatial110.dll. The 3 DLLs should go together to when used with other AutoCAD computers, no need to install SQL Server.

Download source code here.







No comments:

Followers

About Me

My photo
After graduating from university, I worked as civil engineer for more than 10 years. It was AutoCAD use that led me to the path of computer programming. Although I now do more generic business software development, such as enterprise system, timesheet, billing, web services..., AutoCAD related programming is always interesting me and I still get AutoCAD programming tasks assigned to me from time to time. So, AutoCAD goes, I go.