December 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      










« Reminder: free API webcast series on using WPF in your AutoCAD applications | Main | Adding a dynamic distance property to an AutoCAD object using .NET »

July 13, 2009

Triangulating an AutoCAD 3D solid from a set of points using .NET

As promised in this previous post, we’re now going to an update to Zeljko’s triangulation code which creates Solid3d objects with a user-specified depth from our SubDMesh.

Here’s the C# code with the new/modified lines in red (here is the .cs source file). I’ve modified the code to reformat it and to provide the various options side-by-side in separate commands

    1 using Autodesk.AutoCAD.ApplicationServices;

    2 using Autodesk.AutoCAD.DatabaseServices;

    3 using Autodesk.AutoCAD.Runtime;

    4 using Autodesk.AutoCAD.EditorInput;

    5 using Autodesk.AutoCAD.Geometry;

    6 using System;

    7 

    8 public class Triangulate

    9 {

   10   public bool circum(

   11     double x1, double y1, double x2,

   12     double y2, double x3, double y3,

   13     ref double xc, ref double yc, ref double r)

   14   {

   15     // Calculation of circumscribed circle coordinates and

   16     // squared radius

   17 

   18     const double eps = 1e-6;

   19     const double big = 1e12;

   20     bool result = true;

   21     double m1, m2, mx1, mx2, my1, my2, dx, dy;

   22 

   23     if ((Math.Abs(y1 - y2) < eps) && (Math.Abs(y2 - y3) < eps))

   24     {

   25       result = false;

   26       xc = x1; yc = y1; r = big;

   27     }

   28     else

   29     {

   30       if (Math.Abs(y2 - y1) < eps)

   31       {

   32         m2 = -(x3 - x2) / (y3 - y2);

   33         mx2 = (x2 + x3) / 2;

   34         my2 = (y2 + y3) / 2;

   35         xc = (x2 + x1) / 2;

   36         yc = m2 * (xc - mx2) + my2;

   37       }

   38       else if (Math.Abs(y3 - y2) < eps)

   39       {

   40         m1 = -(x2 - x1) / (y2 - y1);

   41         mx1 = (x1 + x2) / 2;

   42         my1 = (y1 + y2) / 2;

   43         xc = (x3 + x2) / 2;

   44         yc = m1 * (xc - mx1) + my1;

   45       }

   46       else

   47       {

   48         m1 = -(x2 - x1) / (y2 - y1);

   49         m2 = -(x3 - x2) / (y3 - y2);

   50         if (Math.Abs(m1 - m2) < eps)

   51         {

   52           result = false;

   53           xc = x1;

   54           yc = y1;

   55           r = big;

   56         }

   57         else

   58         {

   59           mx1 = (x1 + x2) / 2;

   60           mx2 = (x2 + x3) / 2;

   61           my1 = (y1 + y2) / 2;

   62           my2 = (y2 + y3) / 2;

   63           xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);

   64           yc = m1 * (xc - mx1) + my1;

   65         }

   66       }

   67     }

   68     dx = x2 - xc;

   69     dy = y2 - yc;

   70     r = dx * dx + dy * dy;

   71     return result;

   72   }

   73 

   74   private enum OutputObjectType

   75   {

   76     PolyFaceMesh = 1,

   77     Solid3d = 2,

   78     SubDMesh = 4,

   79     All = 7

   80   }

   81 

   82   private void TriangulatePoints(

   83     OutputObjectType objType, int maxpoints

   84   )

   85   {

   86     Document doc =

   87       Application.DocumentManager.MdiActiveDocument;

   88     Database db = doc.Database;

   89     Editor ed = doc.Editor;

   90 

   91     bool createSubDMesh =

   92       (objType & OutputObjectType.SubDMesh) > 0;

   93     bool createPolyFaceMesh =

   94       (objType & OutputObjectType.PolyFaceMesh) > 0;

   95     bool createSolid3d =

   96       (objType & OutputObjectType.Solid3d) > 0;

   97 

   98     TypedValue[] tvs = { new TypedValue(0, "POINT") };

   99     SelectionFilter sf =

  100       new SelectionFilter(tvs);

  101     PromptSelectionOptions pso = new PromptSelectionOptions();

  102     pso.MessageForAdding = "\nSelect points:";

  103     pso.AllowDuplicates = false;

  104     PromptSelectionResult psr = ed.GetSelection(pso, sf);

  105 

  106     if (psr.Status == PromptStatus.Error) return;

  107     if (psr.Status == PromptStatus.Cancel) return;

  108 

  109     SelectionSet ss = psr.Value;

  110     int npts = ss.Count;

  111     if (npts < 3)

  112     {

  113       ed.WriteMessage("Minimum of 3 points must be selected!");

  114       return;

  115     }

  116     if (npts > maxpoints)

  117     {

  118       ed.WriteMessage("Maximum number of points exceeded!");

  119       return;

  120     }

  121 

  122     double zref = 0.0;

  123     if (createSolid3d)

  124     {

  125       PromptDoubleResult ps =

  126         ed.GetDouble(

  127           "\nEnter Z coordinate of reference plane:"

  128         );

  129       if (ps.Status != PromptStatus.OK) return;

  130       zref = ps.Value;

  131     }

  132 

  133     int i, j, k, ntri, ned, nouted, status1 = 0, status2 = 0;

  134     bool status;

  135 

  136     // Point coordinates

  137 

  138     double[] ptx = new double[maxpoints + 3];

  139     double[] pty = new double[maxpoints + 3];

  140     double[] ptz = new double[maxpoints + 3];

  141 

  142     // Triangle definitions

  143 

  144     int[] pt1 = new int[maxpoints * 2 + 1];

  145     int[] pt2 = new int[maxpoints * 2 + 1];

  146     int[] pt3 = new int[maxpoints * 2 + 1];

  147 

  148     // Circumscribed circle

  149 

  150     double[] cex = new double[maxpoints * 2 + 1];

  151     double[] cey = new double[maxpoints * 2 + 1];

  152     double[] rad = new double[maxpoints * 2 + 1];

  153     double xmin, ymin, xmax, ymax, dx, dy, xmid, ymid;

  154     int[] ed1 = new int[maxpoints * 2 + 1];

  155     int[] ed2 = new int[maxpoints * 2 + 1];

  156     int[] outed1 = null;

  157     if (createSolid3d)

  158       outed1 = new int[maxpoints + 1];

  159 

  160     ObjectId[] idarray = ss.GetObjectIds();

  161     Transaction tr =

  162       db.TransactionManager.StartTransaction();

  163     using (tr)

  164     {

  165       DBPoint ent;

  166       k = 0;

  167       for (i = 0; i < npts; i++)

  168       {

  169         ent =

  170           (DBPoint)tr.GetObject(

  171             idarray[k], OpenMode.ForRead, false

  172           );

  173         ptx[i] = ent.Position[0];

  174         pty[i] = ent.Position[1];

  175         ptz[i] = ent.Position[2];

  176         for (j = 0; j < i; j++)

  177           if ((ptx[i] == ptx[j]) && (pty[i] == pty[j]))

  178           {

  179             i--; npts--; status2++;

  180           }

  181         k++;

  182       }

  183       tr.Commit();

  184     }

  185 

  186     if (status2 > 0)

  187       ed.WriteMessage(

  188         "\nIgnored {0} point(s) with same coordinates.",

  189         status2

  190       );

  191 

  192     // Supertriangle

  193 

  194     xmin = ptx[0]; xmax = xmin;

  195     ymin = pty[0]; ymax = ymin;

  196     for (i = 0; i < npts; i++)

  197     {

  198       if (ptx[i] < xmin) xmin = ptx[i];

  199       if (ptx[i] > xmax) xmax = ptx[i];

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

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

  202     }

  203     dx = xmax - xmin; dy = ymax - ymin;

  204     xmid = (xmin + xmax) / 2; ymid = (ymin + ymax) / 2;

  205     i = npts;

  206     ptx[i] = xmid - (90 * (dx + dy)) - 100;

  207     pty[i] = ymid - (50 * (dx + dy)) - 100;

  208     ptz[i] = 0;

  209     pt1[0] = i;

  210     i++;

  211     ptx[i] = xmid + (90 * (dx + dy)) + 100;

  212     pty[i] = ymid - (50 * (dx + dy)) - 100;

  213     ptz[i] = 0;

  214     pt2[0] = i;

  215     i++;

  216     ptx[i] = xmid;

  217     pty[i] = ymid + 100 * (dx + dy + 1);

  218     ptz[i] = 0;

  219     pt3[0] = i;

  220     ntri = 1;

  221     circum(

  222       ptx[pt1[0]], pty[pt1[0]], ptx[pt2[0]],

  223       pty[pt2[0]], ptx[pt3[0]], pty[pt3[0]],

  224       ref cex[0], ref cey[0], ref rad[0]

  225     );

  226 

  227     // Main loop

  228 

  229     for (i = 0; i < npts; i++)

  230     {

  231       ned = 0;

  232       xmin = ptx[i]; ymin = pty[i];

  233       j = 0;

  234       while (j < ntri)

  235       {

  236         dx = cex[j] - xmin; dy = cey[j] - ymin;

  237         if (((dx * dx) + (dy * dy)) < rad[j])

  238         {

  239           ed1[ned] = pt1[j]; ed2[ned] = pt2[j];

  240           ned++;

  241           ed1[ned] = pt2[j]; ed2[ned] = pt3[j];

  242           ned++;

  243           ed1[ned] = pt3[j]; ed2[ned] = pt1[j];

  244           ned++;

  245           ntri--;

  246           pt1[j] = pt1[ntri];

  247           pt2[j] = pt2[ntri];

  248           pt3[j] = pt3[ntri];

  249           cex[j] = cex[ntri];

  250           cey[j] = cey[ntri];

  251           rad[j] = rad[ntri];

  252           j--;

  253         }

  254         j++;

  255       }

  256 

  257       for (j = 0; j < ned - 1; j++)

  258         for (k = j + 1; k < ned; k++)

  259           if ((ed1[j] == ed2[k]) && (ed2[j] == ed1[k]))

  260           {

  261             ed1[j] = -1; ed2[j] = -1; ed1[k] = -1; ed2[k] = -1;

  262           }

  263 

  264       for (j = 0; j < ned; j++)

  265         if ((ed1[j] >= 0) && (ed2[j] >= 0))

  266         {

  267           pt1[ntri] = ed1[j];

  268           pt2[ntri] = ed2[j];

  269           pt3[ntri] = i;

  270           status =

  271             circum(

  272               ptx[pt1[ntri]], pty[pt1[ntri]], ptx[pt2[ntri]],

  273               pty[pt2[ntri]], ptx[pt3[ntri]], pty[pt3[ntri]],

  274               ref cex[ntri], ref cey[ntri], ref rad[ntri]

  275             );

  276           if (!status)

  277           {

  278             status1++;

  279           }

  280           ntri++;

  281         }

  282     }

  283 

  284     // Removal of outer triangles

  285 

  286     i = 0; nouted = 0;

  287     while (i < ntri)

  288     {

  289       if ((pt1[i] >= npts) ||

  290           (pt2[i] >= npts) ||

  291           (pt3[i] >= npts))

  292       {

  293         if (createSolid3d)

  294         {

  295           if ((pt1[i] >= npts) &&

  296               (pt2[i] < npts) &&

  297               (pt3[i] < npts))

  298           {

  299             ed1[nouted] = pt2[i];

  300             ed2[nouted] = pt3[i];

  301             nouted++;

  302           }

  303           if ((pt2[i] >= npts) &&

  304               (pt1[i] < npts) &&

  305               (pt3[i] < npts))

  306           {

  307             ed1[nouted] = pt3[i];

  308             ed2[nouted] = pt1[i];

  309             nouted++;

  310           }

  311           if ((pt3[i] >= npts) &&

  312               (pt1[i] < npts) &&

  313               (pt2[i] < npts))

  314           {

  315             ed1[nouted] = pt1[i];

  316             ed2[nouted] = pt2[i];

  317             nouted++;

  318           }

  319         }

  320         ntri--;

  321         pt1[i] = pt1[ntri];

  322         pt2[i] = pt2[ntri];

  323         pt3[i] = pt3[ntri];

  324         cex[i] = cex[ntri];

  325         cey[i] = cey[ntri];

  326         rad[i] = rad[ntri];

  327         i--;

  328       }

  329       i++;

  330     }

  331 

  332     if (createSolid3d)

  333     {

  334       outed1[0] = 0;

  335       for (i = 1; i < nouted; i++)

  336         for (j = 1; j < nouted; j++)

  337           if (ed2[outed1[i - 1]] == ed1[j])

  338           {

  339             outed1[i] = j;

  340             j = nouted;

  341           }

  342       outed1[nouted] = 0;

  343     }

  344 

  345     tr = db.TransactionManager.StartTransaction();

  346     using (tr)

  347     {

  348       BlockTable bt =

  349         (BlockTable)tr.GetObject(

  350           db.BlockTableId,

  351           OpenMode.ForRead,

  352           false

  353         );

  354       BlockTableRecord btr =

  355         (BlockTableRecord)tr.GetObject(

  356           bt[BlockTableRecord.ModelSpace],

  357           OpenMode.ForWrite,

  358           false

  359         );

  360 

  361       if (createPolyFaceMesh)

  362       {

  363         PolyFaceMesh pfm = new PolyFaceMesh();

  364         btr.AppendEntity(pfm);

  365         tr.AddNewlyCreatedDBObject(pfm, true);

  366         for (i = 0; i < npts; i++)

  367         {

  368           PolyFaceMeshVertex vert =

  369             new PolyFaceMeshVertex(

  370               new Point3d(ptx[i], pty[i], ptz[i])

  371             );

  372           pfm.AppendVertex(vert);

  373           tr.AddNewlyCreatedDBObject(vert, true);

  374         }

  375         for (i = 0; i < ntri; i++)

  376         {

  377           FaceRecord face =

  378             new FaceRecord(

  379               (short)(pt1[i] + 1),

  380               (short)(pt2[i] + 1),

  381               (short)(pt3[i] + 1),

  382               0

  383             );

  384           pfm.AppendFaceRecord(face);

  385           tr.AddNewlyCreatedDBObject(face, true);

  386         }

  387       }

  388 

  389       if (createSubDMesh || createSolid3d)

  390       {

  391         Point3dCollection vertarray = new Point3dCollection();

  392         Int32Collection facearray = new Int32Collection();

  393 

  394         for (i = 0; i < npts; i++)

  395           vertarray.Add(new Point3d(ptx[i], pty[i], ptz[i]));

  396 

  397         if (createSolid3d)

  398           for (i = 0; i < nouted; i++)

  399             vertarray.Add(

  400               new Point3d(

  401                 ptx[ed1[outed1[i]]], pty[ed1[outed1[i]]], zref

  402               )

  403             );

  404 

  405         j = 0;

  406         for (i = 0; i < ntri; i++)

  407         {

  408           facearray.Add(3);

  409           facearray.Add(pt1[i]);

  410           facearray.Add(pt2[i]);

  411           facearray.Add(pt3[i]);

  412         }

  413 

  414         if (createSolid3d)

  415         {

  416           for (i = 0; i < nouted; i++)

  417           {

  418             facearray.Add(4);

  419             k = outed1[i];

  420             facearray.Add(ed1[k]);

  421             facearray.Add(ed2[k]);

  422             if (i == nouted - 1)

  423             {

  424               facearray.Add(npts);

  425             }

  426             else

  427             {

  428               facearray.Add(npts + i + 1);

  429             }

  430             facearray.Add(npts + i);

  431           }

  432           facearray.Add(nouted);

  433           for (i = 0; i < nouted; i++)

  434             facearray.Add(npts + i);

  435         }

  436 

  437         SubDMesh sdm = new SubDMesh();

  438         sdm.SetDatabaseDefaults();

  439         sdm.SetSubDMesh(vertarray, facearray, 0);

  440         btr.AppendEntity(sdm);

  441         tr.AddNewlyCreatedDBObject(sdm, true);

  442 

  443         if (createSolid3d)

  444         {

  445           Solid3d sol = null;

  446           try

  447           {

  448             sol = sdm.ConvertToSolid(false, false);

  449             btr.AppendEntity(sol);

  450             tr.AddNewlyCreatedDBObject(sol, true);

  451           }

  452           catch

  453           {

  454             ed.WriteMessage(

  455               "\nMesh was too complex to turn into a solid."

  456             );

  457           }

  458           if (!createSubDMesh)

  459             sdm.Erase();

  460         }

  461       }

  462 

  463       tr.Commit();

  464     }

  465     if (status1 > 0)

  466       ed.WriteMessage(

  467         "\nWarning! {0} thin triangle(s) found!" +

  468         " Wrong result possible!",

  469         status1

  470       );

  471     Application.UpdateScreen();

  472   }

  473 

  474   [CommandMethod("PFT")]

  475   public void PolyFaceTriangulate()

  476   {

  477     TriangulatePoints(OutputObjectType.PolyFaceMesh, 32767);

  478   }

  479 

  480   [CommandMethod("SDT")]

  481   public void SubDTriangulate()

  482   {

  483     TriangulatePoints(OutputObjectType.SubDMesh, 200000);

  484   }

  485 

  486   [CommandMethod("S3T")]

  487   public void Solid3dTriangulate()

  488   {

  489     TriangulatePoints(OutputObjectType.Solid3d, 200000);

  490   }

  491 }

Let’s now run the SDT command to create a SubDMesh on the left, and the S3T command to create a Solid3d with a particular depth on the right:

Results of SDT and S3T commands in 2D wireframe plan

Now let’s compare them in a 3D view with the conceptual visual style applied:

Results of SDT and S3T commands in conceptual 3D view

We can then use MASSPROP to determine the 3D solid’s volume (which I’m sure would be useful when analysing land surveys, for instance).

A quick word on performance: the Solid3d object is really not designed to represent large volumes derived from point clouds. This approach can work for modest datasets, but you should not expect it to work acceptably with many thousands of points.

Zeljko very cleverly implemented an alternative approach using (unsupported) P/Invoke calls to acdbEntMake() rather than creating the initial SubDMesh and converting it. This resulted in an improvement of nearly 50% (he ran it on a dataset of 5,000 points in 2 minutes, vs. 3:45 minutes with the above code), but ultimately the fundamental bottleneck is elsewhere: Solid3d objects are intended to be precise solid modelling primitives rather than containers for point clouds.

Thanks again for a fascinating series, Zeljko! :-)

TrackBack

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

Listed below are links to weblogs that reference Triangulating an AutoCAD 3D solid from a set of points using .NET:

blog comments powered by Disqus

Feed/Share

10 Random Posts