Find closest point in 2D mashed array

Question

To give y'all some context, I'm doing this inversion technique where I am trying to reproduce a profile using the integrated values. To do that I need to find the value within an array along a certain line(s). To exemplify my issue I have the following code:

fig, ax = plt.subplots(1, figsize = (10,10))


#Create the grid (different grid spacing):
X = np.arange(0,10.01,0.25)
Y = np.arange(0,10.01,1.00)

#Create the 2D array to be plotted
Z = []
for i in range(np.size(X)):
    Zaux = []
    for j in range(np.size(Y)):
        Zaux.append(i*j + j)
        
        ax.scatter(X[i],Y[j], color = 'red', s = 0.25)
    Z.append(Zaux)

#Mesh the 1D grids:
Ymesh, Xmesh = np.meshgrid(Y, X)

#Plot the color plot:
ax.pcolor(Y,X,  Z, cmap='viridis', vmin=np.nanmin(Z), vmax=np.nanmax(Z))

#Plot the points in the grid of the color plot:
for i in range(np.size(X)):
    for j in range(np.size(Y)):
        ax.scatter(Y[j],X[i], color = 'red', s = 3)

#Create a set of lines:

for i in np.linspace(0,2,5):
    X_line = np.linspace(0,10,256)
    Y_line = i*X_line*3.1415-4

    #Plot each line:
    ax.plot(X_line,Y_line, color = 'blue')

ax.set_xlim(0,10)
ax.set_ylim(0,10)

plt.show()

That outputs this graph:

I need to find the closest points in Z that are being crossed by each of the lines. The idea is to integrate the values in Z that are crossed by the blue lines and plot that as a function of slope of the lines. Anyone has a good solution for it? I've tried a set of for loops, but I think it's kind of clunky.

Anyway, thanks for your time...

Answer 1

I am not sure about the closest points thing. That seems "clunky" too. What if it passes exactly in the middle between two points? Also I already had written code that weighs the four neighbor pixels by their closeness for an other project so I am going with that. Also I take the liberty of not rescaling the picture.

i,j = np.meshgrid(np.arange(41),np.arange(11))
Z = i*j + j

class Image_knn():
    def fit(self, image):
        self.image = image.astype('float')

    def predict(self, x, y):
        image = self.image
        weights_x = [1-(x % 1), x % 1]
        weights_y = [1-(y % 1), y % 1]
        start_x = np.floor(x).astype('int')
        start_y = np.floor(y).astype('int')
        return sum([image[np.clip(np.floor(start_x + x), 0, image.shape[0]-1).astype('int'),
                          np.clip(np.floor(start_y + y), 0, image.shape[1]-1).astype('int')] * weights_x[x]*weights_y[y] 
                    for x,y in itertools.product(range(2),range(2))])

And a little sanity check it returns the picture if we give it it's coordinates.

image_model = Image_knn()
image_model.fit(Z)
assert np.allclose(image_model.predict(*np.where(np.ones(Z.shape, dtype='bool'))).reshape((11,41)), Z)

I generate m=100 lines and scale the points on them so that they are evenly spaced. Here is a plot of every 10th of them.

n = 1000
m = 100
slopes = np.linspace(1e-10,10,m)
t, slope = np.meshgrid(np.linspace(0,1,n), slopes)
x_max, y_max = Z.shape[0]-1, Z.shape[1]-1
lines_x = t
lines_y = t*slope
scales = np.broadcast_to(np.stack([x_max/lines_x[:,-1], y_max/lines_y[:,-1]]).min(axis=0), (n,m)).T
lines_x *= scales
lines_y *= scales

And finally I can get the "points" consisting of slope and "integral" and draw it. You probably should take a closer look at the "integral" it's just a ruff guess of mine.

%%timeit
points = np.array([(slope, np.mean(image_model.predict(lines_x[i],lines_y[i]))
                    *np.linalg.norm(np.array((lines_x[i,-1],lines_y[i,-1]))))
                   for i,slope in enumerate(slopes)])

plt.scatter(points[:,0],points[:,1])

Notice the %%timeit in the last block. This takes ~38.3 ms on my machine and therefore wasn't optimized. As Donald Knuth puts it "premature optimization is the root of all evil". If you were to optimize this you would remove the for loop, shove all the coordinates for line points in the model at once by reshaping and reshaping back and then organize them with the slopes. But I saw no reason to put myself threw that for a few ms.

And finally we get a nice cusp as a reward. Notice that it makes sense that the maximum is at 4 since the diagonal is at a slope of 4 for our 40 by 10 picture. The intuition for the cusp is a bit harder to explain but I guess you probably have that already. For the length it comes down to the function (x,y) -> sqrt(x^2+y^2) having different directional differentials when going up and when going left on the rectangle.