I have found a way but it involves looping over all the faces (worked okay in a model of 10 images).

In our own setup we used to do this by simply solving the linear system when camera projection is know and the height. I can see why this is smarter as it dont need the height as its found by the model here. Its abit shower but more robust, so i think i am happy with this for now.

Code was adapted from one of your other posts.

steps = [(0,0), (camera.width,0),(camera.width,camera.height),(0,camera.height)]

for x,y in steps:

point = PhotoScan.Vector([x, y])

point = camera.calibration.unproject(point)

point = camera.transform.mulv(point)

vect = point

p = PhotoScan.Vector(camera.center)

for face in faces:

v = face.vertices

E1 = PhotoScan.Vector(vertices[v[1]].coord - vertices[v[0]].coord)

E2 = PhotoScan.Vector(vertices[v[2]].coord - vertices[v[0]].coord)

D = PhotoScan.Vector(vect)

T = PhotoScan.Vector(p - vertices[v[0]].coord)

P = cross(D, E2)

Q = cross(T, E1)

result = PhotoScan.Vector([Q * E2, P * T, Q * D]) / (P * E1)

if (0 < result[1]) and (0 < result[2]) and (result[1] + result[2] <= 1):

t = (1 - result[1] - result[2]) * vertices[v[0]].coord

u = result[1] * vertices[v[1]].coord

v_ = result[2] * vertices[v[2]].coord

res = chunk.transform.mulp(u + v_ + t)

res = chunk.crs.project(res)

print("{:>04d}".format(x + 1) + "\t" + "{:04d}".format(y + 1) + "\t" + "{:.4f}".format(res[0]) + "\t" + "{:.4f}".format(res[1]) + "\t" + "{:.4f}".format(res[2]) + "\n")

break #finish when the first intersection is found