Author Topic: Can I get an assist from Math Geniuses  (Read 10065 times)

0 Members and 1 Guest are viewing this topic.

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Can I get an assist from Math Geniuses
« on: March 20, 2013, 09:41:07 PM »
Posted in autolisp forum but someone suggested this would be better location. I agree.

I need to get estimated normal values for points extracted from a point cloud.  Specifically estimated from covariance analysis. The point cloud would be reduced by preprocessing methods as well as isolated to individual buildings.

What I'd like to do is demonstrated here:

http://pointclouds.org/documentation/tutorials/normal_estimation.php#normal-estimation

This is beyond me and my skills so I'm begging for anyone's help.

Thanks!

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #1 on: March 20, 2013, 10:46:52 PM »
I can help you next week if you don't manage anything before then.

Basically you need to figure out what that summation equation really means and do it for every point in the cloud.

It looks complicated, but really it's not too bad.  Very similar to doing audio effects on a DSP.

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #2 on: March 21, 2013, 03:20:03 PM »
Thanks Will! I've been trying to find an alternative, but no luck. Will let you know.  I'd be happy to compensate you for your efforts. - Alien
« Last Edit: March 21, 2013, 10:14:41 PM by Alien »

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #3 on: March 24, 2013, 10:10:47 PM »
after a quick browse through AcPointCloud.h it appears there may be a function in place already for getting the normal vector.  Not quite sure yet if that's wuite what you're needing.

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #4 on: March 25, 2013, 02:58:31 AM »
It could be.  The software I'm using is a piece of a beta project which uses xyz and n (n being normal vector info) to construct a building best as possible using aerial LIDAR data. That page was the one the gentleman who sent me the software directed me to.  I shot him an email to see if they developed an automated way to retrieve the "n" data.  I believe I downloaded the C++ (cpp) files and actual installation files for the program on that page. I'm trying to dig through and figure it out now.

I'm new to C#.Net but if accessing point cloud normal data is not different than accessing data for text or polylines, I should have no problem getting the info.  The important part really would be to insure the vector is aiming toward the centroid (inside) of the points... which I think is required to construct the model properly.

Thanks for the advice. 
« Last Edit: March 28, 2013, 10:29:16 PM by Alien »

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #5 on: March 25, 2013, 01:35:05 PM »
That said I'll look into what needs to be done to implement the needed Point Cloud Library materials.  Can you give me a point cloud to play with?

Thanks.

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #6 on: March 26, 2013, 09:47:28 AM »
Definitely! I'm at work but the data I have is at home.  I'll upload this evening.

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #7 on: March 27, 2013, 12:52:40 AM »
point cloud attached.

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #8 on: March 28, 2013, 01:54:05 AM »
Downloading the point cloud library stuff... damn it's almost 400 MB!

Sooo the point cloud you've uploaded is coming into my AutoCAD as an xref which is in your MyDocuments folder...
I have vanilla AutoCAD 2013, I've never worked with point clouds before so I don't know if this is something I need an object enabler to view
« Last Edit: March 28, 2013, 01:58:29 AM by WILL HATCH »

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #9 on: March 28, 2013, 03:00:34 AM »
Crud! Sorry Will!  I believe point clouds are civil 3d / map objects... not sure if they created an object enabler for it.... without that, I don't think you can access the methods for it...  I'll try to see if I can't write something noobish this weekend that might access it ... if you don't mind I might be asking you for a little guidance if needed... is that okay? Hope I'm not waisting your time... I enjoy this stuff so it's basically a hobby for me

I did find they created the .exe to complete the normal vectors on that site I sent. I had to download a trillion things to get it... Unfortunately I have to convert the .las file (LIDAR file) to a .pcd file first (their own made up filetype), using another .exe they developed lol... the input needed though is technical, consisting of coordinates and related to the lidar flight paths, etc itself... a photogrammatrist would probably have no problems... unfortunately I am just a CAD geek... so I'm stuck having to dig in and learn about that lol... I'm scared that when I figure out how to convert to .pcd I'll have to do something else.... it's like a black hole......eventually it's gonna smoosh me lol

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #10 on: March 28, 2013, 03:02:58 AM »
After digging a little more in your drawing, it's almost entirely made up of unused block references.  The point cloud is just xref'd in.  Can you bind it, purge, and then upload again?

Playing with getting point cloud library set up... This is some very cool software!

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #11 on: March 28, 2013, 10:18:58 PM »
Agreed! You can download almost any area/city through USGS website (.las file)...

I'm not sure how to bind my pointcloud...should be in my drawing (though it's a civil 3d drawing...with a bunch of stupid styles with blocks). I believe you are using the pointcloud part of autocad (pointcloudattach)... which I think is different than civil 3d pointclouds (see my attached)... but won't matter in the end... as long as there's a way to find the vectors...

If you're using the autocad version of attaching point clouds just need to type pointcloudattach and add this database (might have to dropdown your extensions to pick the right type).

Let me know. Thanks will!
« Last Edit: March 28, 2013, 10:27:25 PM by Alien »

nobody

  • Swamp Rat
  • Posts: 861
  • .net stuff
Re: Can I get an assist from Math Geniuses
« Reply #12 on: March 28, 2013, 10:23:47 PM »
I have developed one thought concerning this matter. I'm wondering if the centroid is totally necessary...I don't know.

If I simply pick a point in the middle somewhere, find the lowest point in the pointcloud (or just regular autocad points for that matter), go down an additional 10 meters or whatever, and find the vector from every point to it I think it would accomplish what I need...though I'm not sure because the vector needed is for "surface normals".... I don't understand that equation I posted but I the way the vectors are generated are actually by connecting 3 points, finding the vector perpendicular to the surface generated from those points toward the centroid, and then moving on... this would mean tons and tons and tons of calculations for even a small file. What do you think?

If I can't get  the pointcloud centroid to work I'll have to try the first option I mentioned... don't think i could pull off the second.
« Last Edit: March 28, 2013, 10:35:32 PM by Alien »

WILL HATCH

  • Bull Frog
  • Posts: 450
Re: Can I get an assist from Math Geniuses
« Reply #13 on: April 12, 2013, 11:03:29 AM »
I did some playing with the kd tree software gile recently posted (again MANY kudos gile)and managed to get something of a result.  It is quirky, and is by no means complete but could be something to get started.  Some of the points do not get a normal added.  I think this is due to a weak understanding of how to use the tree on my behalf.

Code - C#: [Select]
  1.         [CommandMethod("MyPointCloud")]
  2.         public void GeneratePoints()
  3.         {
  4.             Document doc = Application.DocumentManager.MdiActiveDocument;
  5.             Editor ed = doc.Editor;
  6.             Database db = doc.Database;
  7.             using (Transaction tr = doc.TransactionManager.StartTransaction())
  8.             {
  9.                 BlockTable bt = (BlockTable)tr.GetObject(db.BlockTableId, OpenMode.ForRead);
  10.                 BlockTableRecord btr = (BlockTableRecord)tr.GetObject(bt[BlockTableRecord.ModelSpace], OpenMode.ForWrite);
  11.                 Random rnd = new Random();
  12.                 for (int i = 0; i < 1000; i++)
  13.                 {
  14.                     DBPoint p = new DBPoint(new Point3d(5 * rnd.NextDouble(), 100 * rnd.NextDouble(), 100 * rnd.NextDouble()));
  15.                     btr.AppendEntity(p);
  16.                     tr.AddNewlyCreatedDBObject(p, true);
  17.                    
  18.                 }
  19.                 tr.Commit();
  20.             }
  21.         }
  22.         [CommandMethod("GetNormals")]
  23.         public void GetNormals()
  24.         {
  25.             Document doc = Application.DocumentManager.MdiActiveDocument;
  26.             Editor ed = doc.Editor;
  27.             Database db = doc.Database;
  28.             RXClass rxc = RXClass.GetClass(typeof(DBPoint));
  29.            
  30.  
  31.  
  32.                 using (BlockTableRecord btr = (BlockTableRecord)db.CurrentSpaceId.Open(OpenMode.ForWrite))
  33.                 {
  34.                     ObjectId[] ids = btr.Cast<ObjectId>().Where(id => id.ObjectClass == rxc).ToArray();
  35.                     int len = ids.Length;
  36.                     Point3d[] pts = new Point3d[len];
  37.                     for (int i = 0; i < len; i++)
  38.                     {
  39.                         using (DBPoint pt = (DBPoint)ids[i].Open(OpenMode.ForRead))
  40.                         {
  41.                             pts[i] = pt.Position;
  42.                         }
  43.                     }
  44.                     Point3dTree tree = new Point3dTree(pts, false);
  45.                     Dictionary<Point3d, List<Point3d>> neighbors = new Dictionary<Point3d, List<Point3d>>();
  46.                     getNeighbors(tree, tree.Root, 15, neighbors);
  47.                     List<Vector3d> vectors = new List<Vector3d>();
  48.                     foreach (var set in neighbors)
  49.                     {
  50.                         Vector3d ave = new Vector3d(0, 0, 0);
  51.                         if (set.Value.Count < 2) continue;
  52.                         Point3d pt = Nearest(set.Key,set.Value);
  53.                         int count = set.Value.Count;
  54.                         set.Value.Remove(pt);
  55.                         for (int i = set.Value.Count-1; i > -1; i--)
  56.                         {
  57.                             vectors.Add(set.Key.GetVectorTo(pt));
  58.                             pt = Nearest(pt, set.Value);
  59.                             set.Value.Remove(pt);
  60.                         }
  61.                         Vector3d tmp;
  62.                         for (int i = 0; i < vectors.Count-1; i++)
  63.                         {
  64.                             tmp = vectors[i].CrossProduct(vectors[i + 1]);
  65.                             if (ave.DotProduct(tmp) < 0)
  66.                                 ave -= tmp;
  67.                             else
  68.                                 ave += tmp;
  69.                         }
  70.                         tmp = vectors[0].CrossProduct(vectors[vectors.Count - 1]);
  71.                         if (ave.DotProduct(tmp) < 0)
  72.                             ave -= tmp;
  73.                         else
  74.                             ave += tmp;
  75.                         vectors.Clear();
  76.                         ave = ave.DivideBy(ave.Length / 10);
  77.                         using (Line line = new Line(set.Key,set.Key+ave))
  78.                         {
  79.                             btr.AppendEntity(line);
  80.                         }
  81.                     }
  82.                
  83.                
  84.             }
  85.         }
  86.         private Point3d Nearest(Point3d pt, List<Point3d> pts)
  87.         {
  88.             double best = double.MaxValue;
  89.             Point3d ret = new Point3d();
  90.             foreach (Point3d pnt in pts)
  91.             {
  92.                 double cur =pt.DistanceTo(pnt);
  93.                 if (cur<best)
  94.                 {
  95.                     best = cur;
  96.                     ret = pnt;
  97.                 }
  98.             }
  99.             return ret;
  100.         }
  101.         private void getNeighbors(Point3dTree tree, KdTreeNode<Point3d> node, double dist, Dictionary<Point3d, List<Point3d>> neighbors)
  102.         {
  103.             if (node == null) return;
  104.             node.Processed = true;
  105.             if (neighbors.ContainsKey(node.Value))//this is just a sanity check, MUST be removed for processing large items once method is proven
  106.             {
  107.                 Application.DocumentManager.MdiActiveDocument.Editor.WriteMessage("\nPoint processed again");
  108.             }
  109.             else
  110.             {
  111.                 List<Point3d> pts = new List<Point3d>();
  112.                 foreach (Point3d pt in tree.NearestNeighbours(node, dist))
  113.                 {
  114.                     pts.Add(pt);
  115.                 }
  116.                 neighbors.Add(node.Value, pts);
  117.             }
  118.             getNeighbors(tree, node.LeftChild, dist, neighbors);
  119.             getNeighbors(tree, node.RightChild, dist, neighbors);
  120.            
  121.         }
« Last Edit: April 12, 2013, 02:20:39 PM by WILL HATCH »

gile

  • Gator
  • Posts: 2507
  • Marseille, France
Re: Can I get an assist from Math Geniuses
« Reply #14 on: April 12, 2013, 12:38:17 PM »
Hi WILL HATCH,

It seems you haven't uderstand the advantages of a Kd Tree structure when searching neighbours.
The interest of Kd trees is to avoid searching points in portions of the space (whole sub-trees). At each recursive loop the coordinate of the current node corresponding to the current depth is compared to this of the point which neighbours are searched to see if the search will continue only in the left child node, only in the right one or in both.
If you search systematically in both, you don't need to use a kd tree, a linear search in a classic collection should be faster.
See my last attempt with some more methods.
Speaking English as a French Frog