I have two sets of 3D points in different co ordinate sytems and orientations they contain a number of common points how can I find the best fit / least squares 3D transformation that will translate one set of data onto the other set of data using the common points.
You are after a variant (more or less) of this?
http://nghiaho.com/?page_id=671
What language you speak?
Hi Peter
That is exactly what I am after and wondered if it can be done in Grasshopper or Rhino VB. English is my language so any information you may have will be gratfully received.
Roger
Er … this is not the thing that I had in mind: C#,(my game) P or that dead (but alive) VB?
Anyway:
Imagine a point cloud that has a subset A that we assume identical(x) to a subset B of another point cloud. (x) If not no harm done.
If we can find (fast) 3 pts in A that define a triad same as a triad from 3 pts in B … case closed [Same means: the obvious (triangle to triangle equality) AND not the obvious (same topology of points up/below the triangle planes)].
Why the triads can cut the mustard? it’s a classic Matrix game that’s why: see attached.
Transformations_3PtsTo3Pts_V1.3dm (2.0 MB)
Transformations_3PtsTo3Pts_V1A.gh (140.9 KB)
So the big deal is to find these 2 triads.Sounds peanuts … but it’s not if you have 1Z * 2 points to search. I have stuff that does that almost in real-time but it’s strictly internal to the practice (since is a weapon for reverse engineering matters) .
I’ll try (I can’t promise any ETA) to write something else that can solve the puzzle … hopefully in less than a year and a half.
Did some very stupid thing: 500 milliseconds for 1000 random points in setA, setB [note: setB is shuffled randomly - FisherYates and the likes - so the indices differ 100%]. This means that it takes a year and a half to finish if we have 1M++ points and a decade if we have 1Z++ points.
Moral: pathetic
Hi Peter,
You are far more advanced with Grasshopper than myself, I am just a self taught Rhino VB programmer so this is all new to me. I would like to see the GH code for the finding of three points if you wouldn’t mind as partly shown in the last screenshot. It is all new information and logic to myself but I think I am starting to understand it.
Many thanks for your help to date it has really opened my eyes to a new approach of solving problems.
Hi Roger
Well … as I said I can’t post the real thing here since is strictly internal. So … the challenging part is to mastermind something that I can post … without using the restricted nuts and bits. See for instance results for 2.3 times the points at 1/10th of the time (achieved via the crap test C# captured above [that is deleted from my workstation anyway]) on an ancient I5 that I always use for testing performance. Compare the Method calls as well. Time drops (around 20-30 milliseconds) for 3K points on an I9.
So for the moment my help is 100% zero … but I’ll try to write something that has a performance that allows some real-life usage: because solving any problem is nothing … solving it fast (ideally in real-time) is everything.
In the mean time: how many are your common points? If they are, say, beyond the 10K mark… well …
BTW: The solution to similar problems IS NOT the obvious one: i.e. a triple nested loop on setB: for i = 0;i<setB.Count;i++, for j = i+1;j<setB.Count;j++, for k = j+1; k<setB.Count;k++. Why? well … because if you have 10K points you’ll need (worst case scenario) 999B comparisons … BEFORE counting the points that are above/below your candidate triads … meaning: adios amigos.
This is probably a naive question, but if the points are the same and only the coordinate systems are different, then as @PeterFotiadis said, we need to identify three points that are the unequivocally same in both systems. So what’s preventing you from finding say the first, second and third closest point to the centroid, assuming there’s no degeneracy?
Not naive at all IF YOU ARE 100% sure that the points are the same (i.e. an Intersection of LIDAR setA to setB). But in fact they are never the same: they are “likely” to be the same within some tolerance (if we are lucky, that is). On the other hand 999.998 points on setA that yield a centroid and 999.997 points on setB that yield another centroid => Armageddon. On the other hand a laser shooting aspects of a big object from different positions … well you understand the issue. On the other hand even if the Point3dList Methods are more or less fast … 1++M points are not 1K points. On the other hand if you do maths without rounding the values … well try it by using a Point3dList against an OEM random collection (skip David’s “as even as possible” distribution approach and do it the classic way) and then filter the coordinates by a user controlled variable (p = new Point3d (Math.Round(p.X, decimals) , …, …)).
Maybe @RogerD is lucky and he knows for sure they’re the same points, otherwise, point well taken.
I’ll post here a demo on Karma matters soon (Karma is the biggest thing in reverse engineering).
Indeed luck is the 1M question here. See attached: it varies some random points (if U != 0) in set B when transforming A to B … but depending on where the points are in the collection the triad in B may or may not be the same (3 first closest pts with regard the centroids etc etc). Note: we shooting (well … kinda) points in a 50x50x50 m cube and the distortion value is 30 cm for the snapshots below. Imagine doing RE in some big and complex refinery plant and the likes.
But this talking is 100% academic anyway since the fellas doing the LIDAR job they also take care of matching the partial sets as well … so why bother?
Points_Karma_V1.gh (117.6 KB)
I hear you: but what would be CAD without Karma? That indeed is the core of the matter.
1 Like
Hi Peter,
If I could accommodate up to 1000 common points it would meet my needs.
Well … I added an angle test as well (better safe than sorry).
Anyway you can try the attached (IF YOU ARE SURE that the common points are the same … if not the karma info panel tells you the truth out there). Just feed the 2d C# with your common sets and pray. If the sets differ in count the C# stops. If the sets contain only Point3d.Unset the C# stops as well.
Feed the triads in the other C# provided some replies above and get the transformation (4x4 classic matrix etc etc).
Points_Karma_V1A.gh (124.3 KB)
Update: Added some visual stuff:
Points_Karma_V1B.gh (120.7 KB)
Update (of the Update): Added the classic up/down test
Points_Karma_V1C.gh (129.6 KB)
1 Like
If all the points are subject to noise, the centroid should be the most stable value, since in a large data set, the sum of all the random displacement vectors would cancel out. Centroid to centroid would give the translation distance between the two sets, and from what I’ve been reading above, translating both set’s centroids to the origin is the first step in the SVD algorithm. As your example shows, if things get too noisy, you can’t just look for three points with any assurance they’re the same.
Well … use the first C# attached that does the trans matrix using the 2 pt triads and forget SVD.
That said a stable centroid doesn’t mean that the first 3 closest pts yield triads that match. Kinda buying a Harley Davidson: it’s a proper bike but only when viewed from a safe (1+M miles) distance.
BTW: In the real world you are searching for multiple triads (in setB) that match “as best as possible” the triad in setA. Then you evaluate each one and the candidate with the best score is the best fit triad (kinda like Kangaroo works: attempts to find some sort of optimum solution). If no solution is found then change triadA and retry (this means that the whole procedure MUST be as fast as possible AND NOT as accurate as possible … because the amount of repetitions … blah, blah). A vast variety of tricks is used in order to speed up a task that appears impossible on first sight.
BTW: Here’s the ultimate challenge for you: Imagine a random points cloud. We are after finding the best pairs. This means pairs where their distances deviate as little as possible VS the average distance AND no pair has a distance greater than a given limit (the max pt-pt distance in the cloud [find this via a double loop: for ( i=0; i<cloud.Count; i++) , for (j= i+1; j< cloud.Count; j++) … etc etc]). The classic Steiner graph puzzle (for wireless networks, war games and the likes that is). What could be the way to solve this? (NOT iterative proximity).