# Iterate over million points to get unique closest point

Let's say we've two parallel galaxies with same amount of stars. What I want to do is to find the nearest neighbour of GalaxyA in GalaxyB. But if that particular neighbour is shared by other star(closer to another GalaxyA star) then look another nearest neighbour. This way we will get unique nearest neighbour for each star. I've implemented this logic in using 4-5 different algorithms. So far below one is fastest one:

node = hou.pwd()
geo = node.geometry()
pts = geo.points()
targetpts = node.inputs()[1].geometry().points()

if len(targetpts) >= len(pts):
from operator import itemgetter

# setup targetpts list
targetptslist = []
targetptslist = [(n, targetpt.position()) for n, targetpt in enumerate(targetpts)]

# get the distance to every point in target geo
for pt in pts:
neardistlist = []
p1 = pt.position()
neardistlist = [(targetptslist[i][0], (p1 - targetptslist[i][1]).length()) for i in range(len(targetptslist))]

# sort the list by min distance
neardistlist.sort(key = itemgetter(1))

# check the neardistlist to see if this point has already been taken then remove this from the targetptslist
nearestpt = (neardistlist[0][0])
for j in range(len(targetptslist)):
ptn = targetptslist[j][0]
if ptn == nearestpt:
del targetptslist[j]
break
if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

# update 'uniqueNeighbour' attribute
pt.setAttribValue('uniqueNeighbour', nearestpt)

if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

else:
raise hou.NodeError('Target points must be equal or more than source points!')


Where pts are GalaxyA stars and targetpts are GalaxyB stars. hou is software dependent module.

More readable code(same as above) without list comprehension:

node = hou.pwd()
geo = node.geometry()
pts = geo.points()
targetpts = node.inputs()[1].geometry().points()

if len(targetpts) >= len(pts):
from operator import itemgetter

# setup targetpts list
targetptslist = []
#targetptslist = [(n, targetpt.position()) for n, targetpt in enumerate(targetpts)] # short n fast version of below loop
for n, targetpt in enumerate(targetpts):
targetptinfo = (n, targetpt.position())
targetptslist.append(targetptinfo)
if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

# get the distance to every point in target geo
for pt in pts:
neardistlist = []
p1 = pt.position()
#neardistlist = [(targetptslist[i][0], (p1 - targetptslist[i][1]).length()) for i in range(len(targetptslist))] # short n fast version of below loop
for i in range(len(targetptslist)):
tptinfo = targetptslist[i]
p2 = tptinfo[1]
distance = (p1 - p2).length()
targetinfo = (tptinfo[0], distance)
neardistlist.append(targetinfo)
if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

# sort the list by min distance
#neardistlist.sort(key = lambda ptdist: ptdist[1])
neardistlist.sort(key = itemgetter(1)) # faster than lambda sorting

# check the neardistlist to see if this point has already been taken then remove this from the targetptslist
nearestpt = (neardistlist[0][0])
for j in range(len(targetptslist)):
ptn = targetptslist[j][0]
if ptn == nearestpt:
del targetptslist[j]
break
if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

# update 'uniqueNeighbour' attribute
pt.setAttribValue('uniqueNeighbour', nearestpt)

if hou.updateProgressAndCheckForInterrupt(): break # respect keyboard interruption

else:
raise hou.NodeError('Target points must be equal or more than source points!')


My question is: Can we optimize it more so it will be faster? Current version takes around 90 seconds for only 6000 stars, and I'm talking about millions. Using current algorithm forget about millions, it'll take just forever.

EDIT: This question is to refine algorithm. I need better way to sort and iterate over.

current algorithm: step 1 (if len(targetpts) >= len(pts)): Only run if GalaxyB have same or more number of stars. Initialize uniqueNeighbour to -1.

step 2 (targetptslist = [(n, targetpt.position()) for n, targetpt in enumerate(targetpts)]): Store GalaxyB star position in a list.

step3 (for l, pt in enumerate(pts)): Loop over GalaxyA stars to store each GalaxyB star distance in a list(neardistlist). Then sort it by minimum distance.

step 4 (for j in range(len(targetptslist))): Loop over targetptslist(outcome of step 2) to compare each element to currently added element in neardistlist(which is nearestpt) if matches then delete it from targetptslist and break out of loop

step 5: Update uniqueNeighbour variable.

EDIT 02: Now let's assume there're 5 stars in both galaxies. GalaxyA neighbour list(let's call it closestNeighbourList) in GalaxyB starting from closest to fartherest...

star 0: [4, 1, 0, 3, 2]

star 1: [2, 0, 4, 3, 1]

star 2: [2, 1, 3, 0, 4]

star 3: [0, 3, 2, 4, 1]

star 4: [0, 3, 1, 4, 2]

from these lists default first neighbour list:

star 0: [4]

star 1: [2]

star 2: [2]

star 3: [0]

star 4: [0]

Since I want unique neighbour in GalaxyB, but here you see a problem, star 2 in GalaxyB is closest to two GalaxyA stars 1 & 2 and same for star 0. So these are not unique neighbour as shared by two or more stars. So to overcome this problem what I'm doing first storing all neighbours in a list(closestNeighbourList) then looping over those to check if it's closer to(shared by) some other star then go to next neighbour and when it's find unique neighbour break out of the loop.

after applying unique neighbour algorithm first neighbour list will be:

star 0: [4]

star 1: [2]

star 2: [1]

star 3: [3]

star 4: [0]

• What imports have you used? You're making it needlessly hard on the reviewers by not including those.
– Mast
Sep 2, 2016 at 10:43
• Is the indentation correct? The StackExchange software converts tabs to four spaces, and I rather get the impression that it's messed up most of the contents of the main if block. Sep 2, 2016 at 11:08
• I've fixed the indentation. I'll clean the code bit more so it doesn't need any dependency. For time being you can consider pts and targetpts list of tuples with two elements where first one is point number and second one is point position. Sep 2, 2016 at 11:45
• A few pointers that might help: en.wikipedia.org/wiki/Computational_geometry en.wikipedia.org/wiki/Closest_pair_of_points_problem en.wikipedia.org/wiki/Rotating_calipers#Applications - those non obvious techniques might help you to get new ideas even if your problems seems inherently O(n^2). Sep 2, 2016 at 13:01
• Just wondering if the question has anything to do with K -nearest neighbour algorithm. If it does please indicate in the title and the question's body Sep 2, 2016 at 13:36

The question doesn't make it clear exactly what your problem is, and so it is hard for us to help you. You are trying to find a matching in a weighted bipartite graph, but what are the conditions on that matching?

• Are you looking for a stable matching? That is, one where there's no two items that could both be brought closer in the matching by pairing them.

• Are you looking for a locally minimum matching? That is, one where the sum of distances between the pairs in the matching cannot be reduced by swapping two of the assignments.

• Are you looking for a globally minimum matching? That is, one where the sum of distances between the pairs in the matching is minimized among all matchings.

The algorithm in your post does not find any of these types of matching. So if we are to believe your code, then any matching will do, in which case why not use the built-in function zip?

• Please check "EDIT 02". Sep 4, 2016 at 9:59
• @PradeepBarua: I understand what your implementation does, what I don't understand is what conditions you need the solution to satisfy. Sep 4, 2016 at 10:15
• I'm trying to make it faster. So far I tried 'map' function and list comprehension which you can see in my code and that makes it little bit faster but not fast enough. In EDIT 03 you can see my old implementation which is even slower. Sep 4, 2016 at 10:26
• @PradeepBarua: I still don't understand what conditions you need the solution to satisfy. Are you trying to find a stable matching, a minimum matching, or what? Until I understand what the program is supposed to do, I can't help you make it faster. Sep 4, 2016 at 11:51
• @PradeepBarua: Can you update the question, please, so that it says that you are trying to find a stable matching? Sep 4, 2016 at 14:49