Topic: Usable code for TIN volume calculation

[dgorsman]

I found a lot of useful information here: http://paulbourke.net/papers/triangulate/

My standard reference as well.  I had a similar requirement for our civil/structural department, since most of our work is in a very flat part of the world a few thousand points would be sufficient for most of our projects.  Got pretty far with it, including managing basic contour lines along with the triangles (yes, triangular 3D Face objects).  Other priorities took hold, and our parent company gave us access to Civil3D licenses, so development was stopped at that point.

For differential volumes, I was going to compare the points of the two triangulated areas.  The outer points for one area would be projected onto the triangles of the other, interpolating points where necessary (intersection of vertical line with 3-point plane just requires a little math) so there is a common bounding area.  Thats the point where I got stuck with the odd end-shape volume calculations.  Eventually I went to thinking about the typical areas and volumes involved, the equipment/methods involved in the field, and the additional percentage that gets added (15% is the number I usually see for estimates).  A reasonable estimate on the volume of that odd end piece will easily get lost in all that.

[WILL HATCH]

I suspected this might be a civil cut/fill calculation tool... That's great because it means we don't have to consider the cases of extremely complex cases such as when a surface folds back on itself.

I found a lot of useful information here: http://paulbourke.net/papers/triangulate/

My standard reference as well.  I had a similar requirement for our civil/structural department, since most of our work is in a very flat part of the world a few thousand points would be sufficient for most of our projects.  Got pretty far with it, including managing basic contour lines along with the triangles (yes, triangular 3D Face objects).  Other priorities took hold, and our parent company gave us access to Civil3D licenses, so development was stopped at that point.

For differential volumes, I was going to compare the points of the two triangulated areas.  The outer points for one area would be projected onto the triangles of the other, interpolating points where necessary (intersection of vertical line with 3-point plane just requires a little math) so there is a common bounding area.  Thats the point where I got stuck with the odd end-shape volume calculations.  Eventually I went to thinking about the typical areas and volumes involved, the equipment/methods involved in the field, and the additional percentage that gets added (15% is the number I usually see for estimates).  A reasonable estimate on the volume of that odd end piece will easily get lost in all that.

Are you able to share any of what you had previously achieved?  If you had gotten as far as truncating triangular pyramids then it would be very simple to extend that with something similar to the volume calculation method I posted.

[huiz]

To be honest I don't have any code written yet, just a bunch of urls I saved and about 10 evenings time spent reading and rereading. If I am not able to easily find or write code for accurate volume calculation, then there is no need to create a module. If it would take weeks of programming, then there is also no need to start creating such a module.

I'll try to put a a scheme for the steps to take next week, probably I just need a few hints for little things or optimalisation. To be continued

[huiz]

see this thread
http://www.theswamp.org/index.php?topic=15784.0

Wow, another evening of reading

[WILL HATCH]

I looked at that other thread a bit... Very long!

You mentioned you have it as far as a TIN, did you write a routine for that?  I don't have access to Civil3D here and don't want to redo anything you've already done.  If you have it as far as a tin I've got some ideas that we may be able to work with extending that triangular prism volume calculation.  Just need to work out the algorithm to properly define the prisms...

[WILL HATCH]

Huiz,

Sorry I've been pretty busy lately, but I took a quick look at the drawing you sent me and put this together to do a quick analysis of the problem.  It goes as far as collecting all the points and faces, doing some basic organization (in preparation for the final volume calculation later) and projects them onto a plane (in this case the x-y plane).  To continue from this, you wanted to find the volume contained in the overlapping parts of the 2 TINs.  I would look at all the faces on surface1, and connect them all to a point on the surface2, then look at all the faces on surface 2 and connect them to a point on surface 1.  Finally, if there are any points that have not been matched to a face they will need to be matched to a face on one of the triangular pyramids which we have created in the previous process.

If somebody can suggest a better way to approach this sort of problem I would love to hear it!

Code - C#: [Select]

1. public class Volumetry
2.     {
3.         [CommandMethod("Testing","TestVolume",0)]
4.         public void TestVolume()
5.         {
6.             Document doc = Application.DocumentManager.MdiActiveDocument;
7.             Editor ed = doc.Editor;
8.             Database db = doc.Database;
9.             using (Transaction tr = doc.TransactionManager.StartTransaction())
10.             {
11.  
12.                 BlockTableRecord current = (BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
13.                    
14.                 SelectionFilter sf = new SelectionFilter(new TypedValue[] { new TypedValue((int)DxfCode.Start, "3DFACE") });
15.                 PromptSelectionOptions pso = new PromptSelectionOptions();
16.                 pso.MessageForAdding = "\nSelect faces on surface 1";
17.                 PromptSelectionResult psr1 = ed.GetSelection(pso, sf);
18.                 if (psr1.Status != PromptStatus.OK) return;
19.                 pso.MessageForAdding = pso.MessageForAdding.Replace('1', '2');
20.                 PromptSelectionResult psr2 = ed.GetSelection(pso, sf);
21.                 if (psr2.Status != PromptStatus.OK) return;
22.                 HashSet<myFace> surface1;
23.                 HashSet<myFace> surface2;
24.                 Dictionary<Point3d, myPoint3d> surf1Points;
25.                 Dictionary<Point3d, myPoint3d> surf2Points;
26.                 Plane pl = new Plane(Point3d.Origin, Vector3d.ZAxis);//plane to do 2d projections onto
27.                 GetFacesPoints(psr1, pl, out surface1, out surf1Points);
28.                 GetFacesPoints(psr2, pl, out surface2, out surf2Points);
29.                 foreach (myFace face in surface1)
30.                 {
31.                     face.Face2d.ColorIndex = 3;
32.                     current.AppendEntity(face.Face2d);
33.                     tr.AddNewlyCreatedDBObject(face.Face2d, true);
34.                     DBText text = new DBText();
35.                     text.SetDatabaseDefaults();
36.                     text.Position = face.Center;
37.                     text.TextString = face.index.ToString();
38.                     current.AppendEntity(text);
39.                     tr.AddNewlyCreatedDBObject(text, true);
40.                 }
41.                 foreach (myFace face in surface2)
42.                 {
43.                     face.Face2d.ColorIndex = 5;
44.                     current.AppendEntity(face.Face2d);
45.                     tr.AddNewlyCreatedDBObject(face.Face2d, true);
46.                     DBText text = new DBText();
47.                     text.SetDatabaseDefaults();
48.                     text.Position = face.Center;
49.                     text.TextString = face.index.ToString();
50.                     current.AppendEntity(text);
51.                     tr.AddNewlyCreatedDBObject(text, true);
52.                 }
53.                 foreach (KeyValuePair<Point3d, myPoint3d> kvp in surf1Points)
54.                 {
55.                     Line newLine = new Line(Point3d.Origin, kvp.Value.point);
56.                     current.AppendEntity(newLine);
57.                     tr.AddNewlyCreatedDBObject(newLine, true);
58.                     DBText text = new DBText();
59.                     text.SetDatabaseDefaults();
60.                     string temp = null;
61.                     foreach (int num in kvp.Value.ownerFaces)
62.                     {
63.                         temp += string.Format("{0}, ", num);
64.                     }
65.                     text.TextString = kvp.Value.pointIndex.ToString() + "Faces: " + temp;
66.                     text.Position = kvp.Value.point;
67.                     current.AppendEntity(text);
68.                     tr.AddNewlyCreatedDBObject(text, true);
69.                 }
70.                 foreach (KeyValuePair<Point3d, myPoint3d> kvp in surf2Points)
71.                 {
72.                     Line newLine = new Line(Point3d.Origin, kvp.Value.point);
73.                     current.AppendEntity(newLine);
74.                     tr.AddNewlyCreatedDBObject(newLine, true);
75.                     DBText text = new DBText();
76.                     text.SetDatabaseDefaults();
77.                     string temp = null;
78.                     foreach (int num in kvp.Value.ownerFaces)
79.                     {
80.                         temp += string.Format("{0}, ", num);
81.                     }
82.                     text.TextString = kvp.Value.pointIndex.ToString() + "Faces: " + temp;
83.                     text.Position = kvp.Value.point;
84.                     current.AppendEntity(text);
85.                     tr.AddNewlyCreatedDBObject(text, true);
86.                 }
87.                 tr.Commit();
88.  
89.             }
90.         }
91.  
92.         private void GetFacesPoints(PromptSelectionResult psr, Plane workingPlane, out HashSet<myFace> surface, out Dictionary<Point3d, myPoint3d> surfPoints)
93.         {
94.             surface = new HashSet<myFace>();
95.             surfPoints = new Dictionary<Point3d, myPoint3d>();
96.             for (int faceIndex = 0; faceIndex < psr.Value.Count; faceIndex++)
97.             {
98.                 SelectedObject so = psr.Value[faceIndex];
99.                 surface.Add(new myFace(faceIndex, (Face)so.ObjectId.GetObject(OpenMode.ForRead), workingPlane));
100.             }
101.             int pointIndex = 0;
102.             foreach (myFace face in surface)
103.             {
104.                 Face thisFace = face.Face;
105.                 for (short i = 0; i < 3; i++)
106.                 {
107.                     Point3d thisVertex = thisFace.GetVertexAt(i);
108.                     if (surfPoints.ContainsKey(thisVertex))
109.                     {
110.                         surfPoints[thisVertex].ownerFaces.Add(face.index);
111.  
112.                     }
113.                     else
114.                     {
115.                         surfPoints.Add(thisVertex, new myPoint3d(pointIndex++, thisVertex, face.index, workingPlane));
116.  
117.                     }
118.                 }
119.                
120.  
121.  
122.             }
123.            
124.         }
125.         private class myPoint3d
126.         {
127.             public myPoint3d(int index, Point3d point, int faceIndex, Plane pl)
128.             {
129.                 this.pointIndex = index;
130.                 this.point = point;
131.                 this.ownerFaces = new HashSet<int>();
132.                 this.ownerFaces.Add(faceIndex);
133.                 this.matchedToFace = false;
134.                 this.point2d = pl.GetClosestPointTo(point).GetPoint();
135.             }
136.             public int pointIndex { get; set; }
137.             public Point3d point { get; set; }
138.             public Point3d point2d { get; set; }
139.             public HashSet<int> ownerFaces { get; set; }
140.             public bool matchedToFace { get; set; }
141.         }
142.         private class myFace
143.         {
144.             public myFace(int index, Face Face, Plane pl)
145.             {
146.                 this.index = index;
147.                 this.Face = Face;
148.                 this.vch = new bool[3] { false, false, false };
149.                 Point3d vertex0 = Face.GetVertexAt(0);
150.                 Point3d vertex1 = Face.GetVertexAt(1);
151.                 Point3d vertex2 = Face.GetVertexAt(2);
152.                 this.Center = new Point3d(
153.                     (vertex0.X + vertex0.X + vertex0.X) / 3,
154.                     (vertex1.Y + vertex1.Y + vertex1.Y) / 3,
155.                     (vertex2.Z + vertex2.Z + vertex2.Z) / 3
156.                     );
157.                 this.Face2d = new Face();
158.                 this.Face2d.SetVertexAt(0, pl.GetClosestPointTo(vertex0).GetPoint());
159.                 this.Face2d.SetVertexAt(1, pl.GetClosestPointTo(vertex1).GetPoint());
160.                 this.Face2d.SetVertexAt(2, pl.GetClosestPointTo(vertex2).GetPoint());
161.                 this.Face2d.SetVertexAt(3, pl.GetClosestPointTo(vertex2).GetPoint());
162.                 this.matchedToPoint = false;
163.             }
164.             public int index { get; set; }
165.             public Face Face { get; set; }
166.             public Face Face2d { get; set; }
167.             public Point3d Center { get; set; }
168.             public bool[] vch { get; set; }
169.             public bool matchedToPoint { get; set; }
170.         }
171.     }
172.  

[nobody]

Interesting.

[huiz]

Thanks Will, it is sure useful. I'll try to work out the exact problem and steps to find the solution. As far as I am now it is exactly how you describe, projecting the points of the surfaces to project onto the other.

I'm not sure about the lines which I get with a quick run Have to dig in deeper.

[WILL HATCH]

The lines were an intentional debugging portion, inserted on lines 55-57, 72-74.


TT:

You showed me once that using a HashSet is faster than a list, in this case I used a dictionary with the vertex in question as the key as a mechanism to see if it was being shared with another vertex:

Code - C#: [Select]

1. private void GetFacesPoints(PromptSelectionResult psr, Plane workingPlane, out HashSet<myFace> surface, out Dictionary<Point3d, myPoint3d> surfPoints)
2.         {
3.             surface = new HashSet<myFace>();
4.             surfPoints = new Dictionary<Point3d, myPoint3d>();
5.             for (int faceIndex = 0; faceIndex < psr.Value.Count; faceIndex++)
6.             {
7.                 SelectedObject so = psr.Value[faceIndex];
8.                 surface.Add(new myFace(faceIndex, (Face)so.ObjectId.GetObject(OpenMode.ForRead), workingPlane));
9.             }
10.             int pointIndex = 0;
11.             foreach (myFace face in surface)
12.             {
13.                 Face thisFace = face.Face;
14.                 for (short i = 0; i < 3; i++)
15.                 {
16.                     Point3d thisVertex = thisFace.GetVertexAt(i);
17.                     if (surfPoints.ContainsKey(thisVertex))
18.                     {
19.                         surfPoints[thisVertex].ownerFaces.Add(face.index);
20.  
21.                     }
22.                     else
23.                     {
24.                         surfPoints.Add(thisVertex, new myPoint3d(pointIndex++, thisVertex, face.index, workingPlane));
25.                     }
26.                 }
27.                
28.  
29.  
30.             }
31.            
32.         }

I've just started learning about what is really going on behind the scenes in the HashSet and particularly generating custom hash functions I can see that I could improve performance here by implementing a hash function based only on the vertex and not the extra stuff I store in the myPoint3d object I've created.  By making the following modification to myPoint3d:

****UPDATE: had to add override equals to make this work properly.  Still need to sort out accessing the existing item in the HashSet to add ownerFace 

Code - C#: [Select]

1.         private class myPoint3d
2.         {
3.             public myPoint3d(Point3d point)//int index, Point3d point, int faceIndex, Plane pl)
4.             {
5.                 //this.pointIndex = index;
6.                 this.point = point;
7.                 this.ownerFaces = new HashSet<int>();
8.                 //this.ownerFaces.Add(faceIndex);
9.                 this.matchedToFace = false;
10.                 //this.point2d = pl.GetClosestPointTo(point).GetPoint();
11.             }
12.             public int pointIndex { get; set; }
13.             public Point3d point { get; set; }
14.             public Point3d point2d { get; set; }
15.             public HashSet<int> ownerFaces { get; set; }
16.             public bool matchedToFace { get; set; }
17.             public override int GetHashCode()
18.             {
19.                 return this.point.GetHashCode();
20.             }
21.             public override bool Equals(object obj)
22.             {
23.                 return this.point.Equals(((myPoint3d)obj).point);
24.             }
25.         }

I think I can significantly improve the speed of this operation by changing my GetFacesPoint() routine to:

Code - C#: [Select]

1.         private void GetFacesPoints(PromptSelectionResult psr, Plane workingPlane, out HashSet<myFace> surface, out HashSet<myPoint3d> surfPoints)
2.         {
3.             surface = new HashSet<myFace>();
4.             surfPoints = new HashSet< myPoint3d>();
5.             for (int faceIndex = 0; faceIndex < psr.Value.Count; faceIndex++)
6.             {
7.                 SelectedObject so = psr.Value[faceIndex];
8.                 surface.Add(new myFace(faceIndex, (Face)so.ObjectId.GetObject(OpenMode.ForRead), workingPlane));
9.             }
10.             int pointIndex = 0;
11.             foreach (myFace face in surface)
12.             {
13.                 Face thisFace = face.Face;
14.                 for (short i = 0; i < 3; i++)
15.                 {
16.                     myPoint3d thisVertex = new myPoint3d(thisFace.GetVertexAt(i));
17.                     if (surfPoints.Add(thisVertex))
18.                     {
19.                         thisVertex.pointIndex = pointIndex++;
20.                         thisVertex.point2d = workingPlane.GetClosestPointTo(thisVertex.point).GetPoint();
21.                         thisVertex.ownerFaces.Add(face.index);
22.                     }
23.                     else
24.                         surfPoints.ElementAt(thisVertex.GetHashCode()).ownerFaces.Add(face.index); //****UPDATE: fails here
25.                 }
26.             }
27.         }

Any comments are greatly appreciated!

Cheers

[Jeff H]

A little info for Hashset vs list
http://stackoverflow.com/questions/150750/hashset-vs-list-performance

[WILL HATCH]

After I started down that path I've learned to the best of my knowledge that there is no way to access a member of the HashSet without enumerating over the whole set :realmad: I though I would be able to access the items with the hashcode