Q:

efficiently apply a function on every pixel based on its location

i want to efficiently iterate over every pixel of a given image, and based on the location of the pixel, apply a function and put the output in this location.
Here's what i've tried:
def point_gaussian_value(p1, p2, sigma=0.8):
    x1, y1 = p1
    x2, y2 = p2

    return np.exp(-1 * (np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) / sigma ** 2))

point_annotation = np.zeros_like(sample).astype(np.float32)
for j in range(sample.shape[0]):
    for k in range(sample.shape[1]):

        value = point_gaussian_value(p1=(j, k), p2=(row[header[i]], row[header[i + 1]]))
        point_annotation[j, k] = point_gaussian_value(p1=(j, k), p2=(20, 20))

Though it's naive, and extremely unefficient, specially for large images. I know that i can vectorize a function, such that all calculations will be done concurrently and efficiently, but i'm not sure if its possible/how to vectorize a function that gets a location of a pixel for every pixel.

A:

Given that what the function will be receiving are the indices of both the rows and columns, you can directly compute all the values in a vectorized manner by using np.aranges as x1 and y1 coordinates:
rows = np.arange(sample.shape[0])
cols = np.arange(sample.shape[1])
x2, y2 = (row[header[i]], row[header[i + 1]])
out = np.exp(-1 * (np.sqrt((x2 - rows[:,None])**2 + (y2 - cols)**2) / 0.8**2))

We can also speed up the computation by using the numexpr module:
import numexpr as ne
rows = np.arange(sample.shape[0])
cols = np.arange(sample.shape[1])
x2, y2 = (row[header[i]], row[header[i + 1]])
out = ne.evaluate('exp(-1 * sqrt((x2- rows2D)**2 + (y2- cols)**2) / 0.8**2)', 
                 {'rows2D': rows[:,None]})

Quick check and timings:
n_rows = 1000
n_cols = 1000
rows = np.arange(n_rows)
cols = np.arange(n_cols)
x2, y2 = 10, 12
out = ne.evaluate('exp(-1*sqrt((x2- rows2D)**2 + (y2- cols)**2) / 0.8**2)', 
                 {'rows2D': rows[:,None]})

np.allclose(out, point_annotation, rtol=1e-5)
# True

def current_approach(n_rows, n_cols, x2, y2):
    point_annotation = np.zeros((n_rows,n_cols)).astype(np.float32)
    for j in range(n_rows):
        for k in range(n_cols):
            value = point_gaussian_value(p1=(j, k), p2=(x2,y2))
            point_annotation[j, k] = point_gaussian_value(p1=(j, k), p2=(x2,y2))

def approach_1(n_rows, n_cols, x2, y2):
    rows = np.arange(n_rows)
    cols = np.arange(n_cols)
    out = np.exp(-1 * (np.sqrt((x2- rows[:,None])**2 + (y2- cols)**2) / 0.8**2))

def approach_2(n_rows, n_cols, x2, y2):
    rows = np.arange(n_rows)
    cols = np.arange(n_cols)
    out = ne.evaluate('exp(-1*sqrt((x2- rows2D)**2 + (y2- cols)**2) / 0.8**2)', 
                                   {'rows2D': rows[:,None]})

%timeit current_approach(n_rows, n_cols, x2, y2)
# 10.7 s ± 471 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit approach_1(n_rows, n_cols, x2, y2)
# 59.3 ms ± 426 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit approach_2(n_rows, n_cols, x2, y2)  
# 965 µs ± 43.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Up to a 11000x speedup with the second approach!

