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) 999**B** 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.

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)

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).