Topic: Can I get an assist from Math Geniuses

[WILL HATCH]

Thank you gile! That is the response I was hoping for.  I was struggling to find a good way to handle that situation without a strong understanding on using the tree. Being very tight on time to work on this, that was my 1 hour attempt.  With your input it will be easy to finish this.

[WILL HATCH]

Does this pretty well point you in the right direction now?  This generates normals using the latest version of the kd tree and orients them towards the viewport like the point cloud library does.  Not sure why, but this works on my laptop for 100k points, but at 1M it shat a brick and threw an exception

Code - C#: [Select]

1.         [CommandMethod("GetNormals")]
2.         public void GetNormals()
3.         {
4.             Document doc = Application.DocumentManager.MdiActiveDocument;
5.             Editor ed = doc.Editor;
6.             Database db = doc.Database;
7.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
8.  
9.  
10.  
11.             using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
12.             {
13.                 ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
14.                 int len = ids.Length;
15.                 HashSet<Point3d> pts = new HashSet<Point3d>();
16.                 for (int i = 0; i < len; i++)
17.                 {
18.                     using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
19.                     {
20.                         pts.Add(pt.Position);
21.                     }
22.                 }
23.                 Point3dTree tree = new Point3dTree(pts, false);
24.                 ConcurrentBag<Tuple<Point3d, Vector3d>> normals = new ConcurrentBag<Tuple<Point3d, Vector3d>>();
25.                 Vector3d camera = ed.GetCurrentView().ViewDirection;
26.                 Parallel.ForEach(pts, pt =>
27.                 {
28.                     List<Point3d> neighborinos;
29.                     neighborinos = tree.NearestNeighbours(pt, (int)5);
30.                     Point3d npt = pt.Nearest(neighborinos), lpt;
31.                     Vector3d ave = new Vector3d(0, 0, 0), tmp;
32.                     int count = neighborinos.Count;
33.                     while (neighborinos.Count > 1)
34.                     {
35.                         lpt = npt;
36.                         neighborinos.Remove(npt);
37.                         npt = npt.Nearest(neighborinos);
38.                         tmp = pt.to(lpt).cross(pt.to(npt));
39.                         if (ave.dot(tmp) < 0)
40.                             ave -= tmp;
41.                         else
42.                             ave += tmp;
43.                     }
44.                     if (ave.dot(camera) < 0)
45.                     {
46.                         ave *= -1;
47.                     }
48.                     ave = ave / ave.Length;
49.                     normals.Add(new Tuple<Point3d, Vector3d>(pt, ave));
50.  
51.                 });
52.                 foreach (var item in normals)
53.                 {
54.                     using (Line line = new Line(item.Item1, item.Item1 + 200 * item.Item2))
55.                     {
56.                         btr.AppendEntity(line);
57.                     }
58.                 }
59.             }
60.         }

The routine requires the following extension methods to avoid making any api calls during parallel execution

Code - C#: [Select]

1.     public static class NormalExtensions
2.     {
3.         public static Vector3d to(this Point3d a, Point3d b)
4.         {
5.             return new Vector3d(b.X - a.X, b.Y - a.Y, b.Z - a.Z);
6.         }
7.         public static Vector3d cross(this Vector3d a, Vector3d b)
8.         {
9.             return new Vector3d(a.Y * b.Z - a.Z * b.Y, a.Z * b.X - a.X * b.Z, a.X * b.Y - a.Y * b.X);
10.         }
11.         public static double dot(this Vector3d v1, Vector3d v2)
12.         {
13.             return v1.X * v2.X + v1.Y * v2.Y + v1.Z * v2.Z;
14.         }
15.         public static Point3d Nearest(this Point3d pt, List<Point3d> pts)
16.         {
17.             double best = double.MaxValue;
18.             Point3d ret = new Point3d();
19.             foreach (Point3d pnt in pts)
20.             {
21.                 double cur = pt.DistanceTo(pnt);
22.                 if (cur < best)
23.                 {
24.                     best = cur;
25.                     ret = pnt;
26.                 }
27.             }
28.             return ret;
29.         }
30.     }

[gile]

Hi,

You do not need these extension methods, you can use the Vector3d/Point3d ones even in parallel execution because Vector3d, Point3d, ... are immutable strucures and the members of these structures are independent from AutoCAD.

[WILL HATCH]

I was getting a fatal error when calling Vector3d.CrossProduct(Vector3d) and creating the extension fixed that problem.

[gile]

Hi,

From your code, this works fine for me using the Vector3d built-in methods

Code - C#: [Select]

1.         [CommandMethod("GetNormals")]
2.         public void GetNormals()
3.         {
4.             Document doc = Application.DocumentManager.MdiActiveDocument;
5.             Editor ed = doc.Editor;
6.             Database db = doc.Database;
7.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
8.  
9.             using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
10.             {
11.                 ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
12.                 int len = ids.Length;
13.                 HashSet<Point3d> pts = new HashSet<Point3d>();
14.                 for (int i = 0; i < len; i++)
15.                 {
16.                     using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
17.                     {
18.                         pts.Add(pt.Position);
19.                     }
20.                 }
21.                 Point3dTree tree = new Point3dTree(pts, false);
22.                 Vector3d camera = ed.GetCurrentView().ViewDirection;
23.  
24.                 ParallelQuery<Tuple<Point3d, Vector3d>> normals =
25.                     pts.AsParallel()
26.                     .Select(pt => new Tuple<Point3d, Vector3d>(pt, GetNormal(pt, tree, camera)));
27.  
28.                 foreach (var item in normals)
29.                 {
30.                     using (Line line = new Line(item.Item1, item.Item1 + 200 * item.Item2))
31.                     {
32.                         btr.AppendEntity(line);
33.                     }
34.                 }
35.             }
36.         }
37.  
38.         private Vector3d GetNormal(Point3d pt, Point3dTree tree, Vector3d camera)
39.         {
40.             Vector3d normal = new Vector3d();
41.             Vector3d[] vecs = tree.NearestNeighbours(pt, (int)5)
42.                 .Cast<Point3d>()
43.                 .Select(p => p - pt)
44.                 .ToArray();
45.             for (int i = 0; i < vecs.Length - 1; i++)
46.             {
47.                 Vector3d tmp = vecs[i].CrossProduct(vecs[i + 1]);
48.                 normal = normal.DotProduct(tmp) < 0 ? normal - tmp : normal + tmp;
49.             }
50.             normal = normal / normal.Length;
51.             return normal.DotProduct(camera) < 0 ? normal.Negate() : normal;
52.         }

[WILL HATCH]

Thanks gile!
Something I think is important to note here is that when iterating through the neighbors it should be done by grabbing the nearest in the group as I have done.  Reason being is this approximates creating triangular faces and averaging them together weighted by face area.  This is because vector cross product a x b = |a|*|b|*sin(angle between).

Also, consider the array of vectors {a:(1,0,0),b:(-1,0,0),c:(0,1,0),d:(0,-1,0)}
If they were evaluated linearly the only non-zero result would be b x c

[gile]

Hi,
I didn't care about the normal calculation algorithm. I only wanted to confirm and show an example of what I said about using Point3d and Vector3d structures in a parallel execution.

[WILL HATCH]

That was more for anybody using this code for normal calculations.  Again, thank you for the clarification.  Your code is much cleaner than mine and you've shown me some new methods to investigate.

[nobody]

Thanks will! Can you please tell me what files you referenced?  Really appreciate your work with.

[WILL HATCH]

Have many in the class there (rightly named random experiments) but I think the ones you need are:

using System;
using System.Collections.Generic;
using System.Collections.Concurrent;
using System.Linq;
using System.Threading.Tasks;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;

and of course gile's kdtree.

keep me posted on how this goes plugging into other application, we would have use for this here too.

[nobody]

Will do!  I am getting one error when piecing everything together... I probably did it wrong:

Error   3   Cannot implicitly convert type 'Autodesk.AutoCAD.Geometry.Point3dCollection' to 'System.Collections.Generic.List<Autodesk.AutoCAD.Geometry.Point3d>'   C:\Users\TheCadExpert\Desktop\TestNormalVectors\NormalVectors\NormalVectors\myCommands.cs   49   36   NormalVectors

[WILL HATCH]

guess I should have mentioned I modded gile's tree there to return a List...

Code - C#: [Select]

1.         public List<Point3d> NearestNeighbours(Point3d point, int number)
2.         {
3.             List<Tuple<double, Point3d>> pairs = new List<Tuple<double, Point3d>>(number);
4.             this.center = point;
5.             if (ignoreZ)
6.             {
7.                 this.center2d = point.Flatten();
8.                 this.distanceTo = Distance2dTo;
9.             }
10.             else
11.             {
12.                 this.distanceTo = Distance3dTo;
13.             }
14.             GetNNeighbours(number, double.MaxValue, this.Root, pairs);
15.             List<Point3d> points = new List<Point3d>();
16.             for (int i = 0; i < pairs.Count; i++)
17.             {
18.                 points.Add(pairs[i].Item2);
19.             }
20.             return points;
21.         }

[gile]

Hi,

You can try this way, without changing anything to the Point3dTree class.
Note I changed the GetNormal() method to follow WILL HATCH advice in reply #20.

Code - C#: [Select]

1.         [CommandMethod("GetNormals")]
2.         public void GetNormals()
3.         {
4.             Document doc = Application.DocumentManager.MdiActiveDocument;
5.             Editor ed = doc.Editor;
6.             Database db = doc.Database;
7.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
8.  
9.             using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
10.             {
11.                 ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
12.                 int len = ids.Length;
13.                 HashSet<Point3d> pts = new HashSet<Point3d>();
14.                 for (int i = 0; i < len; i++)
15.                 {
16.                     using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
17.                     {
18.                         pts.Add(pt.Position);
19.                     }
20.                 }
21.                 Point3dTree tree = new Point3dTree(pts, false);
22.                 Vector3d camera = ed.GetCurrentView().ViewDirection;
23.  
24.                 ParallelQuery<Tuple<Point3d, Vector3d>> normals =
25.                     pts.AsParallel()
26.                     .Select(pt => new Tuple<Point3d, Vector3d>(pt, GetNormal(pt, tree, camera)));
27.  
28.                 foreach (var item in normals)
29.                 {
30.                     using (Line line = new Line(item.Item1, item.Item1 + 200 * item.Item2))
31.                     {
32.                         btr.AppendEntity(line);
33.                     }
34.                 }
35.             }
36.         }
37.  
38.         private Vector3d GetNormal(Point3d pt, Point3dTree tree, Vector3d camera)
39.         {
40.             Point3dCollection neighbors = tree.NearestNeighbours(pt, 5);
41.             Point3d npt = Nearest(pt, neighbors), lpt;
42.             Vector3d normal = new Vector3d(), tmp;
43.             while (neighbors.Count > 1)
44.             {
45.                 lpt = npt;
46.                 neighbors.Remove(npt);
47.                 npt = Nearest(npt, neighbors);
48.                 tmp = (lpt - pt).CrossProduct(npt - pt);
49.                 if (normal.DotProduct(tmp) < 0)
50.                     normal -= tmp;
51.                 else
52.                     normal += tmp;
53.             }
54.             if (normal.DotProduct(camera) < 0)
55.             {
56.                 normal = normal.Negate();
57.             }
58.             return normal / normal.Length;
59.         }
60.  
61.         private Point3d Nearest(Point3d pt, Point3dCollection pts)
62.         {
63.             double best = double.MaxValue;
64.             Point3d ret = new Point3d();
65.             foreach (Point3d pnt in pts)
66.             {
67.                 double cur = pt.DistanceTo(pnt);
68.                 if (cur < best)
69.                 {
70.                     best = cur;
71.                     ret = pnt;
72.                 }
73.             }
74.             return ret;
75.         }

[WILL HATCH]

gile,

I haven't gotten into the world of linq yet, but this has been a good learning experience for me! Not understanding how ParallelQuery works I thought to do the following with the hope of catching up to your method in execution time by not waiting for the vector calculations to finish before starting to draw lines.  But eventually discovered that ParallelQuery is doing the same thing, only slightly faster.  Thanks for the great examples!

Code - C#: [Select]

1.         [CommandMethod("cGetNormals")]
2.         public void GetNormals2()
3.         {
4.             Document doc = Application.DocumentManager.MdiActiveDocument;
5.             Editor ed = doc.Editor;
6.             Database db = doc.Database;
7.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
8.  
9.  
10.             System.Diagnostics.Stopwatch tmr = new System.Diagnostics.Stopwatch();
11.             tmr.Start();
12.             using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
13.             {
14.                 ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
15.                 int len = ids.Length;
16.                 HashSet<Point3d> pts = new HashSet<Point3d>();
17.                 for (int i = 0; i < len; i++)
18.                 {
19.                     using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
20.                     {
21.                         pts.Add(pt.Position);
22.                     }
23.                 }
24.                 ConcurrentQueue<Tuple<Point3d, Vector3d>> normals = new ConcurrentQueue<Tuple<Point3d, Vector3d>>();
25.                 Vector3d camera = ed.GetCurrentView().ViewDirection;
26.                 workerThread wt = new workerThread(pts, normals, camera);
27.                 Thread worker = new Thread(new ThreadStart(wt.doWork));
28.                 worker.Start();
29.                 Tuple<Point3d, Vector3d> point;
30.                 bool dequeued = false;
31.                 while (!worker.IsAlive) ;
32.                 while (worker.IsAlive || dequeued)
33.                 {
34.                     if (dequeued = normals.TryDequeue(out point))
35.                     {
36.                         using (Line line = new Line(point.Item1, point.Item1 + 200 * point.Item2))
37.                         {
38.                             btr.AppendEntity(line);
39.                         }
40.                     }
41.                 }
42.             }
43.             tmr.Stop();
44.             ed.WriteMessage("\nElapsed: {0}", tmr.ElapsedTicks);
45.         }

Code - C#: [Select]

1.     public class workerThread
2.     {
3.         HashSet<Point3d> pts;
4.         ConcurrentQueue<Tuple<Point3d, Vector3d>> normals;
5.         Point3dTree tree;
6.         Vector3d camera;
7.         public workerThread(HashSet<Point3d> p, ConcurrentQueue<Tuple<Point3d, Vector3d>> n, Vector3d c)
8.         {
9.             pts = p;
10.             normals = n;
11.             tree = new Point3dTree(pts, false);
12.             camera = c;
13.         }
14.         public void doWork()
15.         {
16.             Parallel.ForEach(pts, pt =>
17.             {
18.                 List<Point3d> neighborinos = tree.NearestNeighbours(pt, (int)5);
19.                 Point3d npt = pt.Nearest(neighborinos), lpt;
20.                 Vector3d ave = new Vector3d(0, 0, 0), tmp, lvt = pt.to(npt);
21.  
22.                 while (neighborinos.Count > 1)
23.                 {
24.                     lpt = npt;
25.                     neighborinos.Remove(npt);
26.                     npt = npt.Nearest(neighborinos);
27.                     tmp = lvt.cross(lvt = pt.GetVectorTo(npt));
28.                     if (ave.dot(tmp) < 0)
29.                         ave -= tmp;
30.                     else
31.                         ave += tmp;
32.                 }
33.  
34.                 if (ave.dot(camera) < 0)
35.                 {
36.                     ave *= -1;
37.                 }
38.                 ave = ave / ave.Length;
39.                 normals.Enqueue(new Tuple<Point3d, Vector3d>(pt, ave));
40.  
41.             });
42.         }
43.     }

[WILL HATCH]

Also, out of curiosity I compared the extensions I created to the API methods.  I suspected mine would be slightly faster since I was avoiding wrapper creation through calling the native code.  Any thoughts on why this is not the case for the dot product?

(this is working on a 100k point cloud)

My VectorTo Elapsed: 7118
Ac VectorTo Elapsed: 8509
My Dot Elapsed: 4590
Ac Dot Elapsed: 2371
My Cross Elapsed: 9527
Ac Cross Elapsed: 10009

Code - C#: [Select]

1.         [CommandMethod("VectorSpeed")]
2.         public void MeasureSpeed()
3.         {
4.             Document doc = Application.DocumentManager.MdiActiveDocument;
5.             Editor ed = doc.Editor;
6.             Database db = doc.Database;
7.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
8.  
9.  
10.             System.Diagnostics.Stopwatch tmr = new System.Diagnostics.Stopwatch();
11.  
12.             using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
13.             {
14.                 ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
15.                 int len = ids.Length;
16.                 HashSet<Point3d> pts = new HashSet<Point3d>();
17.                 for (int i = 0; i < len; i++)
18.                 {
19.                     using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
20.                     {
21.                         pts.Add(pt.Position);
22.                     }
23.                 }
24.                 Point3d[] points = pts.ToArray();
25.                 HashSet<Vector3d> vecs = new HashSet<Vector3d>();
26.                 for (int i = 0; i < points.Length - 1; i++)
27.                 {
28.                     vecs.Add(points[i].to(points[i + 1]));
29.                 }
30.                 Vector3d[] vtrs = vecs.ToArray();
31.                 tmr.Start();
32.                 for (int i = 0; i < points.Length - 1; i++)
33.                 {
34.                     points[i].to(points[i + 1]);
35.                 }
36.                 tmr.Stop();
37.                 ed.WriteMessage("\nMy VectorTo Elapsed: {0}", tmr.ElapsedTicks);
38.                 tmr.Restart();
39.                 for (int i = 0; i < vtrs.Length - 1; i++)
40.                 {
41.                     vtrs[i].cross(vtrs[i + 1]);
42.                 }
43.                 tmr.Stop();
44.                 ed.WriteMessage("\nMy Cross Elapsed: {0}", tmr.ElapsedTicks);
45.                 tmr.Restart();
46.                 for (int i = 0; i < vtrs.Length - 1; i++)
47.                 {
48.                     vtrs[i].dot(vtrs[i + 1]);
49.                 }
50.                 tmr.Stop();
51.                 ed.WriteMessage("\nMy Dot Elapsed: {0}", tmr.ElapsedTicks);
52.                
53.                
54.                 tmr.Start();
55.                 for (int i = 0; i < points.Length - 1; i++)
56.                 {
57.                     points[i].GetVectorTo(points[i + 1]);
58.                 }
59.                 tmr.Stop();
60.                 ed.WriteMessage("\nAc VectorTo Elapsed: {0}", tmr.ElapsedTicks);
61.                 vtrs = vecs.ToArray();
62.  
63.                 tmr.Restart();
64.                 for (int i = 0; i < vtrs.Length - 1; i++)
65.                 {
66.                     vtrs[i].CrossProduct(vtrs[i + 1]);
67.                 }
68.                 tmr.Stop();
69.                 ed.WriteMessage("\nAc Cross Elapsed: {0}", tmr.ElapsedTicks);
70.                 tmr.Restart();
71.                 for (int i = 0; i < vtrs.Length - 1; i++)
72.                 {
73.                     vtrs[i].DotProduct(vtrs[i + 1]);
74.                 }
75.                 tmr.Stop();
76.                 ed.WriteMessage("\nAc Dot Elapsed: {0}", tmr.ElapsedTicks);
77.  
78.             }
79.  
80.         }