find all points in the grid-CodePudding

I have 1 million random points (float x, float y)，they are in the region of a rectangle. Now I want to divide the rectangle into many small rectangles, and the coordinates of the center point of each small rectangle can be obtained, like [(0.5,0.5),(0.5,1.5)]. The problem is how to know which small rectangle each random point belongs to in a very fast way!! Here is my code:

## generate the center point of each small rectangle. 
def gen_center_coor(row = None, col = None,bias = None):
    spatial = np.zeros((row*col,2))
    a,b=np.meshgrid(np.arange(row),np.arange(col),indexing = 'ij')
    spatial[:,0], spatial[:,1] = a.reshape(-1),b.reshape(-1) 
    spatial  = bias
    return spatial
centerdot = gen_center_coor(row =320, col = 320,bias = 0.5)


## generate the number of points per small rectangle (randomnums) 
## and random points (randdot) 

randomnums = np.random.randint(low=5,high=15,size=320 * 320)
randrow = np.random.uniform(low=-0.5, high=0.5, size=sum(randomnums))
randcol = np.random.uniform(low=-0.5, high=0.5, size=sum(randomnums))
randdot_row = np.repeat(centerdot[:, 0], randomnums)   randrow
randdot_col = np.repeat(centerdot[:, 1], randomnums)   randcol
randdot = np.hstack((randdot_row.reshape(-1,1), randdot_col.reshape(-1, 1)))

## in fact, randomnums: just to verify your code
## you only have the randdot
##  here is my main function, but it is too slow for me 
for i in range(100):
    gdnet = np.repeat(centerdot[i],randdot.shape[0]).reshape(randdot.shape[0],2)
    np.where(np.abs(randdot - gdnet).max(-1)<0.5)[0]

CodePudding user response：

With three assumptions:

The rectangles are aligned with the X/Y axes. If they aren't, you can first apply a rotation to the point coordinates first.
The large rectangle spans from (0, 0) to (Lx, Ly) - if not, again, transform your data so that it does.
All the small rectangles are identical in size - in other words, you're splitting your large rectangle into a uniform grid.

the fastest way I know of would be:

import numpy as np
# prep some mock data first
N = 10_000_000
Lx = 1
Ly = 1
points_x = np.random.uniform(0, Lx, N)
points_y = np.random.uniform(0, Ly, N)

N_grid = 1000    # splitting into 1000 uniform rectangles

# and now for the main dish:

dx = Lx / N_grid  # grid size
dy = Ly / N_grid

points_x_indices = (points_x / dx).astype(int)
points_y_indices = (points_y / dy).astype(int)