Triangulating an AutoCAD sub-division mesh from a set of points using .NET

I'm just settling in after a week of very enjoyable paternity leave (thanks to all of you for your congratulatory messages :-). I have a few topics planned for the next few weeks, but in the meantime I'm going to post some code provided by our old friend Zeljko Gjuranic as a follow-up to this previous post.

Zeljko responded to the suggestion that – rather than creating the limited PolyFaceMesh object – his code should instead create the more modern SubDMesh (an object introduced in AutoCAD 2010). I've taken Zeljko's code and integrated it into the previous example, allowing creation of both PolyFaceMesh and SubDMesh objects from two separate commands. This allows us to compare the results. For even greater analysis capabilities, we'll see this code extended in a future post to add a depth to the SubDMesh and create a Solid3d from it, allowing more advanced volumetric analysis.

Here's Zeljko's C# code (edited for formatting and to run side-by-side with the previous implementation):

using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using System;
 
public class Triangulate
{
  public bool circum(
    double x1, double y1, double x2,
    double y2, double x3, double y3,
    ref double xc, ref double yc, ref double r)
  {
    // Calculation of circumscribed circle coordinates and
    // squared radius
 
    const double eps = 1e-6;
    const double big = 1e12;
    bool result = true;
    double m1, m2, mx1, mx2, my1, my2, dx, dy;
 
    if ((Math.Abs(y1 - y2) < eps) && (Math.Abs(y2 - y3) < eps))
    {
      result = false;
      xc = x1; yc = y1; r = big;
    }
    else
    {
      if (Math.Abs(y2 - y1) < eps)
      {
        m2 = -(x3 - x2) / (y3 - y2);
        mx2 = (x2 + x3) / 2;
        my2 = (y2 + y3) / 2;
        xc = (x2 + x1) / 2;
        yc = m2 * (xc - mx2) + my2;
      }
      else if (Math.Abs(y3 - y2) < eps)
      {
        m1 = -(x2 - x1) / (y2 - y1);
        mx1 = (x1 + x2) / 2;
        my1 = (y1 + y2) / 2;
        xc = (x3 + x2) / 2;
        yc = m1 * (xc - mx1) + my1;
      }
      else
      {
        m1 = -(x2 - x1) / (y2 - y1);
        m2 = -(x3 - x2) / (y3 - y2);
        if (Math.Abs(m1 - m2) < eps)
        {
          result = false;
          xc = x1;
          yc = y1;
          r = big;
        }
        else
        {
          mx1 = (x1 + x2) / 2;
          mx2 = (x2 + x3) / 2;
          my1 = (y1 + y2) / 2;
          my2 = (y2 + y3) / 2;
          xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
          yc = m1 * (xc - mx1) + my1;
        }
      }
    }
    dx = x2 - xc;
    dy = y2 - yc;
    r = dx * dx + dy * dy;
    return result;
  }
 
  private enum OutputObjectType
  {
    PolyFaceMesh = 1,
    SubDMesh = 2,
    All = 3
  }
 
  private void TriangulatePoints(
    OutputObjectType objType, int maxpoints
  )
  {
    Document doc =
      Application.DocumentManager.MdiActiveDocument;
    Database db = doc.Database;
    Editor ed = doc.Editor;
 
    bool createSubDMesh =
      (objType & OutputObjectType.SubDMesh) > 0;
    bool createPolyFaceMesh =
      (objType & OutputObjectType.PolyFaceMesh) > 0;
 
    TypedValue[] tvs = { new TypedValue(0, "POINT") };
    SelectionFilter sf =
      new SelectionFilter(tvs);
    PromptSelectionOptions pso = new PromptSelectionOptions();
    pso.MessageForAdding = "\nSelect points:";
    pso.AllowDuplicates = false;
    PromptSelectionResult psr = ed.GetSelection(pso, sf);
 
    if (psr.Status == PromptStatus.Error) return;
    if (psr.Status == PromptStatus.Cancel) return;
 
    SelectionSet ss = psr.Value;
    int npts = ss.Count;
    if (npts < 3)
    {
      ed.WriteMessage("Minimum of 3 points must be selected!");
      return;
    }
    if (npts > maxpoints)
    {
      ed.WriteMessage("Maximum number of points exceeded!");
      return;
    }
 
    int i, j, k, ntri, ned, nouted, status1 = 0, status2 = 0;
    bool status;
 
    // Point coordinates
 
    double[] ptx = new double[maxpoints + 3];
    double[] pty = new double[maxpoints + 3];
    double[] ptz = new double[maxpoints + 3];
 
    // Triangle definitions
 
    int[] pt1 = new int[maxpoints * 2 + 1];
    int[] pt2 = new int[maxpoints * 2 + 1];
    int[] pt3 = new int[maxpoints * 2 + 1];
 
    // Circumscribed circle
 
    double[] cex = new double[maxpoints * 2 + 1];
    double[] cey = new double[maxpoints * 2 + 1];
    double[] rad = new double[maxpoints * 2 + 1];
    double xmin, ymin, xmax, ymax, dx, dy, xmid, ymid;
    int[] ed1 = new int[maxpoints * 2 + 1];
    int[] ed2 = new int[maxpoints * 2 + 1];
 
    ObjectId[] idarray = ss.GetObjectIds();
    Transaction tr =
      db.TransactionManager.StartTransaction();
    using (tr)
    {
      DBPoint ent;
      k = 0;
      for (i = 0; i < npts; i++)
      {
        ent =
          (DBPoint)tr.GetObject(idarray[k], OpenMode.ForRead, false);
        ptx[i] = ent.Position[0];
        pty[i] = ent.Position[1];
        ptz[i] = ent.Position[2];
        for (j = 0; j < i; j++)
          if ((ptx[i] == ptx[j]) && (pty[i] == pty[j]))
          {
            i--; npts--; status2++;
          }
        k++;
      }
      tr.Commit();
    }
 
    if (status2 > 0)
      ed.WriteMessage(
        "\nIgnored {0} point(s) with same coordinates.",
        status2
      );
 
    // Supertriangle 
 
    xmin = ptx[0]; xmax = xmin;
    ymin = pty[0]; ymax = ymin;
    for (i = 0; i < npts; i++)
    {
      if (ptx[i] < xmin) xmin = ptx[i];
      if (ptx[i] > xmax) xmax = ptx[i];
      if (pty[i] < xmin) ymin = pty[i];
      if (pty[i] > xmin) ymax = pty[i];
    }
    dx = xmax - xmin; dy = ymax - ymin;
    xmid = (xmin + xmax) / 2; ymid = (ymin + ymax) / 2;
    i = npts;
    ptx[i] = xmid - (90 * (dx + dy)) - 100;
    pty[i] = ymid - (50 * (dx + dy)) - 100;
    ptz[i] = 0;
    pt1[0] = i;
    i++;
    ptx[i] = xmid + (90 * (dx + dy)) + 100;
    pty[i] = ymid - (50 * (dx + dy)) - 100;
    ptz[i] = 0;
    pt2[0] = i;
    i++;
    ptx[i] = xmid;
    pty[i] = ymid + 100 * (dx + dy + 1);
    ptz[i] = 0;
    pt3[0] = i;
    ntri = 1;
    circum(
      ptx[pt1[0]], pty[pt1[0]], ptx[pt2[0]],
      pty[pt2[0]], ptx[pt3[0]], pty[pt3[0]],
      ref cex[0], ref cey[0], ref rad[0]
    );
 
    // Main loop
 
    for (i = 0; i < npts; i++

)
    {
      ned = 0;
      xmin = ptx[i]; ymin = pty[i];
      j = 0;
      while (j < ntri)
      {
        dx = cex[j] - xmin; dy = cey[j] - ymin;
        if (((dx * dx) + (dy * dy)) < rad[j])
        {
          ed1[ned] = pt1[j]; ed2[ned] = pt2[j];
          ned++;
          ed1[ned] = pt2[j]; ed2[ned] = pt3[j];
          ned++;
          ed1[ned] = pt3[j]; ed2[ned] = pt1[j];
          ned++;
          ntri--;
          pt1[j] = pt1[ntri];
          pt2[j] = pt2[ntri];
          pt3[j] = pt3[ntri];
          cex[j] = cex[ntri];
          cey[j] = cey[ntri];
          rad[j] = rad[ntri];
          j--;
        }
        j++;
      }
 
      for (j = 0; j < ned - 1; j++)
        for (k = j + 1; k < ned; k++)
          if ((ed1[j] == ed2[k]) && (ed2[j] == ed1[k]))
          {
            ed1[j] = -1; ed2[j] = -1; ed1[k] = -1; ed2[k] = -1;
          }
 
      for (j = 0; j < ned; j++)
        if ((ed1[j] >= 0) && (ed2[j] >= 0))
        {
          pt1[ntri] = ed1[j]; pt2[ntri] = ed2[j]; pt3[ntri] = i;
          status =
            circum(
              ptx[pt1[ntri]], pty[pt1[ntri]], ptx[pt2[ntri]],
              pty[pt2[ntri]], ptx[pt3[ntri]], pty[pt3[ntri]],
              ref cex[ntri], ref cey[ntri], ref rad[ntri]
            );
          if (!status)
          {
            status1++;
          }
          ntri++;
        }
    }
 
    // Removal of outer triangles
 
    i = 0; nouted = 0;
    while (i < ntri)
    {
      if ((pt1[i] >= npts) || (pt2[i] >= npts) || (pt3[i] >= npts))
      {
        ntri--;
        pt1[i] = pt1[ntri];
        pt2[i] = pt2[ntri];
        pt3[i] = pt3[ntri];
        cex[i] = cex[ntri];
        cey[i] = cey[ntri];
        rad[i] = rad[ntri];
        i--;
      }
      i++;
    }
 
    tr = db.TransactionManager.StartTransaction();
    using (tr)
    {
      BlockTable bt =
        (BlockTable)tr.GetObject(
          db.BlockTableId,
          OpenMode.ForRead,
          false
        );
      BlockTableRecord btr =
        (BlockTableRecord)tr.GetObject(
          bt[BlockTableRecord.ModelSpace],
          OpenMode.ForWrite,
          false
        );
 
      if (createPolyFaceMesh)
      {
        PolyFaceMesh pfm = new PolyFaceMesh();
        btr.AppendEntity(pfm);
        tr.AddNewlyCreatedDBObject(pfm, true);
        for (i = 0; i < npts; i++)
        {
          PolyFaceMeshVertex vert =
       &

#160;    new PolyFaceMeshVertex(
              new Point3d(ptx[i], pty[i], ptz[i])
            );
          pfm.AppendVertex(vert);
          tr.AddNewlyCreatedDBObject(vert, true);
        }
        for (i = 0; i < ntri; i++)
        {
          FaceRecord face =
            new FaceRecord(
              (short)(pt1[i] + 1),
              (short)(pt2[i] + 1),
              (short)(pt3[i] + 1),
              0
            );
          pfm.AppendFaceRecord(face);
          tr.AddNewlyCreatedDBObject(face, true);
        }
      }
 
      if (createSubDMesh)
      {
        Point3dCollection vertarray = new Point3dCollection();
        Int32Collection facearray = new Int32Collection();
 
        for (i = 0; i < npts; i++)
          vertarray.Add(new Point3d(ptx[i], pty[i], ptz[i]));
 
        j = 0;
        for (i = 0; i < ntri; i++)
        {
          facearray.Add(3);
          facearray.Add(pt1[i]);
          facearray.Add(pt2[i]);
          facearray.Add(pt3[i]);
        }
 
        SubDMesh sdm = new SubDMesh();
        sdm.SetDatabaseDefaults();
        sdm.SetSubDMesh(vertarray, facearray, 0);
        btr.AppendEntity(sdm);
        tr.AddNewlyCreatedDBObject(sdm, true);
      }
 
      tr.Commit();
    }
    if (status1 > 0)
      ed.WriteMessage(
        "\nWarning! {0} thin triangle(s) found!" +
        " Wrong result possible!",
        status1
      );
    Application.UpdateScreen();
  }
 
  [CommandMethod("PFT")]
  public void PolyFaceTriangulate()
  {
    TriangulatePoints(OutputObjectType.PolyFaceMesh, 32767);
  }
 
  [CommandMethod("SDT")]
  public void SubDTriangulate()
  {
    TriangulatePoints(OutputObjectType.SubDMesh, 200000);
  }
}

The code defines two commands: PFT, which triangulates a PolyFaceMesh from a set of points (limiting that set to 32,767, the upper-limit of vertices for this object) and SDT, which triangulates a SubDMesh (and here we've set an arbitrary limit to 200,000, but we could very easily increase that – we're not limited by the SubDMesh object itself).

Let's see what happens when we run the PFT and SDT commands, using each to select an identical (but adjacent) set of points:

Here's are the two meshes in a 3D view with the conceptual visual style applied:

Just to highlight one major difference (other than capacity) between the PolyFaceMesh and SubDMesh objects, here's what we see when we increase the SubDMesh object's smoothing level:

In the next post (in this series, at least), we'll add code to create a Solid3d from the SubDMesh.

13 responses to “Triangulating an AutoCAD sub-division mesh from a set of points using .NET”

Jorge

May 5, 2014 at 4:14 am

Hi Kean,

I know this article is a bit dated but I'm doing some testing with .net and

I'm running your code (translated to VB.net) in AutoCAD 2014 but I keep getting an error on this line:

ent = DirectCast(tr.GetObject(idarray(k), OpenMode.ForRead, False), DBPoint)

This happens when the k value is equal to the ss.count (points selected) value.

Adding this line

If k >= ss.Count Then Exit For

seems to fix the problem and the mesh is drawn, but I'm not sure if something is missing.
Any ideas ?

1. Kean Walmsley
  
  May 5, 2014 at 6:10 am
  
  Hi Jorge,
  
  There may be a translation error somewhere: k shouldn't go beyond the size of ss, unless perhaps the "for" loop contains more than that first "if" statement...
  
  Regards,
  
  Kean
  
  1. Jorge
    
    May 13, 2014 at 4:23 am
    
    Thanks Kean,
    
    This is the translated part of the code that is producing the error:
    
    Dim idarray As ObjectId() = ss.GetObjectIds()
    
    Dim tr As Transaction = db.TransactionManager.StartTransaction()
    
    Using tr
    
    Dim ent As DBPoint
    
    k = 0
    
    For i = 0 To npts - 1
    
    'this is the line I added to fix the error
    If k >= ss.Count Then Exit For
    
    ent = DirectCast(tr.GetObject(idarray(k), OpenMode.ForRead, False), DBPoint)
    
    ptx(i) = ent.Position(0)
    
    pty(i) = ent.Position(1)
    
    ptz(i) = ent.Position(2)
    
    For j = 0 To i - 1
    
    If (ptx(i) = ptx(j)) AndAlso (pty(i) = pty(j)) Then
    
    i -= 1
    
    npts -= 1
    
    status2 += 1
    
    End If
    
    Next
    
    k += 1
    
    Next
    
    tr.Commit()
    
    End Using
    
    ''''''''''''''''''''''''''''''''''''
    
    npts is equal to ss.Count so it should the same as ss.GetObjectIds()
    
    Can you spot the translation problem ?
    
    1. Kean Walmsley
      
      May 13, 2014 at 9:15 am
      
      There doesn't appear to be a translation problem, from what I can tell. k should never be greater than or equal to ss.Count - just looking at the logic of the code, even if npts can reduce over time when there are duplicates - so I still don't quite get how this issue has come up. I expect the C# code would also fail for your particular test case, if there's a problem in your VB.NET translation.
      
      Kean
      
Peter Rimkus

June 8, 2015 at 4:01 pm

I'm also coming in here late, but there's a typo in the "supertriangle" creation:

Shouldn't these two lines :

if (pty[i] < xmin) ymin = pty[i];
if (pty[i] > xmin) ymax = pty[i];

...look like this? :

if (pty[i] < ymin) ymin = pty[i];
if (pty[i] > ymax) ymax = pty[i];

(comparing X's to X's and Y's to Y's??)

Pete

1. Kean Walmsley
  
  June 8, 2015 at 4:19 pm
  
  I would think so. Thanks, Peter!
  
  Kean
  
  1. jaison joseph
    
    October 2, 2016 at 6:32 am
    
    Hai
    I am new in civil 3d. I have points and 3d polylines through points. But 3dpolyline z value is 0 is there any possibilities to call points z value to this 3d polyline vertex.
    Regards
    jaison
    
    1. Kean Walmsley
      
      October 2, 2016 at 4:39 pm
      
      Please post your question to the AutoCAD .NET or Civil 3D developer forum.
      
      Kean
      
Ben

November 20, 2017 at 12:46 am

Can the Brep API be used to create the mesh given the set of input points?

1. Kean Walmsley
  
  November 20, 2017 at 5:08 am
  
  I haven't looked at it in quite a while. Please post to the online forums - someone there will be able to help.
  
  Kean
  
Ali Mobasshir

October 26, 2018 at 8:16 pm

Can you please tell how to use the code without database? I mean taking points one by one from loop?

1. Kean Walmsley
  
  October 26, 2018 at 8:51 pm
  
  Please post your support request to the online forums.
  
  Kean
  
2. cretavpatua197333
  
  May 2, 2019 at 12:49 am
  
  Wanna have fun with a hot girl?😏
  My Whatsapp +4-850-861-85-37 😘
  ⠀
  ⠀
  ⠀
  ⠀
  ⠀
  Message id:5874119