Kean Walmsley


  • About the Author
    Kean on Google+

August 2014

Sun Mon Tue Wed Thu Fri Sat
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            








« Batch-processing AutoCAD drawings from LISP without SDI (take 2) | Main | LINQ post over on The Building Coder »

July 06, 2009

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 =

            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:

Results of PFT and SDT commands in 2D wireframe plan view

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

Results of PFT and SDT commands in conceptual 3D view

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:

Comparing the standard PolyFaceMesh with a smoothed SubDMesh

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

TrackBack

TrackBack URL for this entry:
http://www.typepad.com/services/trackback/6a00d83452464869e2011571c82f65970b

Listed below are links to weblogs that reference Triangulating an AutoCAD sub-division mesh from a set of points using .NET:

blog comments powered by Disqus

10 Random Posts