Count the number of neighbors

I have a numpy array called contour of $$\N\$$ coordinates (x,y) of dimension (N, 2).

For each point in this table, I would like to create a square centered at that point and perform a test on the square formed.

My current code is slow, probably due to the for loop. Can the execution speed be improved?

My code:

def neighbour(point , mask  , n ): # Create the square around this point and count the number of neighbors.
mask = mask[point[0] - int(n/2) : point[0] + int(n/2) + 1,point[1] - int(n/2):point[1] + int(n/2) + 1]

def max_neighbour(contour , mask=maske , n=ne): # Find the point with as many neighbors as possible
t = np.zeros(len(contour)) # contour is the numpy array of dimension (2,N)
for i in range(len(contour)):
return contour[np.argmax(t)] # t contains the list of neighbors for each contour point.

• Take a look at scipy.spatial.cKDTree. query_ball_tree() with r=n and p=1 should return a list of points with a manhatten distance <= n. – RootTwo Aug 7 at 1:46

Okay, so there are a lot of problems here. I think the best way to address them is in stages. I want to show you how to go through your own code critically and address problems in it as you write.

Before I begin, for the benefit of someone who will review your code, it is helpful to include all the imports you used so that the person doing the review can just copy your code into their environment and run it if they need to.

So a good first inclusion would be:

import numpy as np

My first stylistic point is that you really should conform to PEP8, it creates a standard that all Python programmers should follow and while some of its requirements can be a bit frustrating at times it is, for the most part, a very good style guideline.

def count_neighbours(point, mask, n):
# Create the square around this point and count the number of neighbors.
# Change 1: changed name from 'neighbour' to 'count neighbours'.
#     It is often helpful to have expressive names for functions and
#     variables. Since functions DO something I find verb-object phrases
#     very intuitive when I'm not working on a project where function
#     "what does this function do?", answer: oh it "counts neighbours".
# Change 2: fixed spacing in arguments. Look into PEP8 and a linter that
#     will tell you when you place spaces in poor or unconventional
#     locations.
mask = mask[point[0] - int(n/2) : point[0] + int(n/2) + 1,point[1] - int(n/2):point[1] + int(n/2) + 1]
# Change 3: made indent of count_neighbours function a multiple of four.


mask = mask[point[0] - int(n/2) : point[0] + int(n/2) + 1,point[1] - int(n/2):point[1]

There's no way around the fact that this is quite ungainly. What can we do about it? The biggest problem, aesthetics aside, is that at first glance I have no idea what it's actually supposed to do. Here I approve very much that you added a comment to give some explantion as to what this function does, that comment is the only thing that made figuring out the intention easier.

Python actually attaches a string encountered directly after a function defintion as a docstring for that function. You can find style guidelines for writing doctrings here. Let's make that comment more explicit and put it in a docstring.

def count_neighbours(point, mask, n):
"""
Count the neighbours of a point, where neighbours are other points in
mask that are within a square of side length n, centered on point.

If a copy of point appears in mask it is counted as a neighbour.

Arguments:
point -- the point on which the square will be centered.
mask -- the set of points to check.
n -- the side length of the square.
"""
mask = mask[point[0] - int(n/2) : point[0] + int(n/2) + 1,point[1] - int(n/2):point[1] + int(n/2) + 1]


Here I have guessed at what is supposed to be contained in mask and written it in the docstring. I have made this assumption based on what the function is supposed to do, but based on the call in your other function: neighbour(contour[i], mask, n) where you have said that contour is in fact your list of points I suspect there may be some confusion at the variable naming stage.

I am going to give mask the revised name all_points in order to make it clearer what the argument is, and update the docstring accordingly. Now, we've finally arrived back at this
all_points = all_points[point[0] - int(n/2) : point[0] + int(n/2) + 1,point[1] - int(n/2):point[1].
How can I make this easier to follow? One obvious idea is that, since we seem to be concerned with a plane we can name some variables of interest:

px = point[0]  # x-value (first coordinate) of the point p
py = point[1]  # y-value (second coordinate) of the point p
offset = n/2  # half of square side length


And now finally, after some cleaning up, we actually realise that there is a deeper problem here: the confusion expression was not only poorly written, but logically flawed! Why are we slicing using the nearest integer to the boundaries of the square as indicies? This makes no sense. We are interested in comparing the bounds px-offset and px+offset to the contents of all_points, not to the indicies of all_points.

So, what we really want to know which points in all_points satisfy these inequalities:

neighbours_x = all_points[:, 0] >= px-offset
# is a point >= the lower bound in x?
neighbours_x &= all_points[:, 0] <= px+offset
# is a point <= the upper bound in x?
neighbours_y = all_points[:, 1] >= py-offset  # lower bound in y
neighbours_y &= all_points[:, 1] <= py+offset  # upper bound in y
neighbours = neighbours_x & neighbours_y  # neighbours in both axes


Go through each step of this process and print the output at each step using the example given at the end and see what each step does.

Finally we want to return the number of neighbours. Both of the below expressions work. The former is faster, but the latter might give you more intuition about what we're doing:

return neighbours.sum()
return all_points[neighbours].shape[0]


I also recommend you print all_points[neighbours] and see what the output is, again using the example at the end.

import numpy as np
from itertools import product

def count_neighbours(point, all_points, n):
"""
Count the neighbours of a point, where neighbours are other points in
all_points that are within a square of side length n, centered on point.

Arguments:
point -- the point on which the square will be centered.
all_points -- the set of points to check.
n -- the side length of the square.
"""
px = point[0]  # x-value (first coordinate) of the point p
py = point[1]  # y-value (second coordinate) of the point p
offset = n/2  # half of square side length
neighbours_x = all_points[:, 0] >= px-offset
# is a point >= the lower bound in x?
neighbours_x &= (all_points[:, 0] <= px+offset)
# is a point <= the upper bound in x?
neighbours_y = all_points[:, 1] >= py-offset  # lower bound in y
neighbours_y &= all_points[:, 1] <= py+offset  # upper bound in y
neighbours = neighbours_x & neighbours_y  # neighbours in both axes

return neighbours.sum()

# example 1
X = np.array(list(product([-2,-1,0,1,2], [-2,-1,0,1,2])))
count_neighbours(X[0], X, 2)

# example 2
X = np.array(list(product([-2.1,-1,0.3,1.7,2], [-2.2,-1.8,0.1,1,2])))
count_neighbours(X[0], X, 4.7)


I hope this helps you identify logical errors in future. I think you should be able to review and refactor the other function on your own. Feel free to edit your question with the updated version later if you still want a second opinion one it is clearer what the function is doing and then we can address performance issues.

• Thank you for your answer, I agree with you. I deliberately neglected the presentation because it seemed pretty basic to me. I tried to follow your code to execute it on all my points and not just one but I can't do it. So for the moment the subject is not closed because I have no solution to propose. Of course I'll try to look for it myself and avoid bothering you again. Small detail for the return of your function, you wanted to write 'neighbours.sum()' – Sword Fish Aug 7 at 9:18
• Corrected the error. In order to run the function on all points you must iterate through all_points and call the count_neighbours() on each point. – JPI Aug 7 at 12:05
• Ok i tried with your version and it's 6 time longer than mine.. – Sword Fish Aug 7 at 13:10
• Okay, we can look at some other options. Most of the inefficiency in my suggestion come from the fact that we must perform 4 checks on all_points and then a boolean operation on an array of the same size. Before we try and figure out a one-pass way to do it ourselves, have you tried user @RootTwo's suggestion of using scipy.spatial.cKDTree.query_ball_tree()? – JPI Aug 9 at 20:15
• Not yet, i try to understand how to use it in my problem ! But thank you again for you time ! – Sword Fish Aug 12 at 7:48