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:
Now let’s compare them in a 3D view with the conceptual visual style applied:
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! :-)
Subscribe via RSS