h = np.dot(vp1, [1 / f**2, 1 / f**2, 1])
vp2 = np.cross(h, l3)
The original code will cause 'AxisError: axisa: axis -1 is out of bounds for array of dimension 0' error, because h is a number and l3 is a vector.
I guess what you want to implement is:
h = [vp1[0]*(1 / f**2), vp1[1]*(1 / f**2), vp1[2]]
vp2 = np.cross(h, l3)
It works well on the test image.