generating a multidimensional numpy array such that a[x,y] = [x,y] (calculating essential matrix from dense optical flow)

I need to take the output flow (shaped (cols,rows,2)) of a dense optical flow calculated through this:

flow = cv2.calcOpticalFlowFarneback(im0, im1, None, 0.5, 3, 15, 3, 5, 1.2, 0)

and generate 2 one-dimensional lists of points old_pts and new_pts such that:

• old_pts is a list of pixel coordinates [[y,x],...]
• new_pts is a list of corrected pixel coordinates: [[y+flow[y,x][1],x+flow[y,x][0]],...]

it is crucial that the 2 lists are ordered in the same way so that old_pts[i] is referring to same pixel as new_pts[i]
i.e. old_pts1 = [0,1] => new_pts1 = [0+flow0,1,1+flow[0,1][0]]

opencv docs reference

to feed them into this:
E,_ = cv2.findEssentialMat(old_pts,new_pts,ch.K)

i've got the following working solution (but it's terribly slow [~15" on a 4k image]):

    old_pts = []
    new_pts = []
    for y in range(cols):
        for x in range(rows):
            old_pts.append([y,x])
            new_pts.append([y+flow[y,x][1],x+flow[y,x][0]])
    
    old_pts = np.array(old_pts)
    new_pts = np.array(new_pts)

is there a better (numpythonic) way to generate the two point lists?

what i think it should be done is:

• generating and filling an image-shaped array "a" such that a[x,y] = [x,y]
• flipping flow axis 2
• adding a to the flipped flow
• reshape and transpose a and the result to generate the two points lists

Example of what i need:

this is `flow` (in real life is 3000x4000):  
| y,x |         0        |       1        |        2       |
|-----|------------------|----------------|----------------|
|   0 | [[[-0.81,-0.55], | [0.73,0.83],   | [-0.3,-0.86]], |
|   1 | [[-0.33,-0.71],  | [0.86,-0.27],  | [0.11,-0.03]], |
|   2 | [[0.46,-0.51],   | [-0.35,-0.88], | [0.4,-0.7]]]   |

this is what i need:  
| old_pts |        new_pts      |
|---------|---------------------|
| [[0,0], | [[0+-0.55,0+-0.81], |
| [1,0],  | [1+-0.71,0+-0.33],  |
| [2,0],  | [2+-0.51,0+0.46],   |
| [0,1],  | [0+0.83,1+0.73],    |
| [1,1],  | [1+-0.27,1+0.86],   |
| [2,1],  | [2+-0.88,1+-0.35],  |
| [0,2],  | [0+-0.86,2+-0.3],   |
| [1,2],  | [1+-0.03,2+0.11],   |
| [2,2]]  | [2+-0.7,2+0.4]]     |

Found a much faster solution:

    xes = np.tile(np.arange(im0.shape[1]),(im0.shape[0],1))
    yes = np.tile(np.arange(im0.shape[0])[:,None],(1,im0.shape[1]))

    nxes = xes + flow[:,:,0]
    nyes = yes + flow[:,:,1]

    xy  = np.stack((yes,xes),axis = 2)
    nxy = np.stack((nyes,nxes),axis = 2)

    old_pts = np.reshape(xy,(-1, xy.shape[-1]))
    new_pts = np.reshape(nxy,(-1, nxy.shape[-1]))