Transformation matrix for image reprojection and image real position

[Gall]

Hi!

I'm currently working on a full toolchain for UAV imagery and I have a couple of questions. All my data is georefenced and we plan on using the automatic GCP detection.

First point, I would like to export the estimated position of all images in world space (in the same coordinate system as the input). I understand that I can get the estimated camera position and lot of important information with the Omega Phi Kappa file, which is great, but I didn't find anything for the projected image.
The goal here would be to calculate the real position, in world space, of the four corners for each raw image. Is there a way to do that easily with a python script?

Second point, where I'm quite lost at, I need to extract all necessary informations to be able to reproject all pixels of the raw image to a reference plane later on. To do this, I think I need to export the corresponding part of the point cloud for each image (can I do this directly in python?) and calculate the final transformation matrix which I need to apply to each point (I don't know how the full equation would be...). Am I wrong?

Thanks for the help!

Edit: Actually I don't need to correct perspective for the projection of the raw images.

[Alexey Pasumansky]

Hello Gall,

If you wish to export multiple orthophotos based on the individual images, so can use custom Python scripts for this purpose. Also using Python script you can extract projected image corners.

Do you think these two option will be applicable for your tasks? In this case I can post a sample script here (although I think they could be already posted somewhere in Python Scripting section).

[Gall]

Hello Alexey, thanks for your answer.

I didn't think about exporting multiple orthophotos but this could be a good idea! I plan on doing all processes by python script so these options would be perfect.

Thanks.

[Alexey Pasumansky]

Hello Gall,

Here are script samples for PhotoScan Pro version 1.0.4.

This one should export the corners of the projected image (requires reconstructed model):

Code: [Select]

#compatibility - PhotoScan Professional v 1.0.4

import time
import PhotoScan

def cross(a, b):
result = PhotoScan.Vector([a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y *b.x])
return result

print("Script started")

cam_index = PhotoScan.app.getInt("Input camera index (starting from zero): ")
save_path = PhotoScan.app.getSaveFileName("Specify output file:")

t0 = time.time()
file = open(save_path, "wt")

doc = PhotoScan.app.document
chunk = doc.activeChunk
model = chunk.model
faces = model.faces
vertices = model.vertices

camera = chunk.photos[cam_index]
print(camera) #camera label

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

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)

file.write("{:>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")
#face.selected = True
break #finish when the first intersection is found


file.close()
t1 = time.time()
t1 -= t0
t1 = float(t1)
print("Script finished in " + "{:.2f}".format(t1) + " seconds.")

And this one can be used for individual orthophotos generation (by disabling all but one camera):

Code: [Select]

#Batch export of orthophotos based on each photo

#compatibility Agisoft PhotoScan Pro 1.0.0
#no arguments required

import os
import time
import PhotoScan

doc = PhotoScan.app.document

PhotoScan.app.messageBox("Prossing started.\nPress OK.")  #information message

def intersect(p0, pn, l0, l):
d = ((p0 - l0) * pn) / (l * pn)
return d * l + l0

def surf_height(chunk, photo):
points_h = list()
point_cloud = chunk.point_cloud
num = len(point_cloud.projections[photo])
num_valid = 0

for i in range (0, num):
x = point_cloud.projections[photo][i].index

if (point_cloud.points[x].valid == False):
continue

v = PhotoScan.Vector( (point_cloud.points[x].coord[0], point_cloud.points[x].coord[1], point_cloud.points[x].coord[2], 1) )
vt = chunk.transform * v
vt.size = 3
vt = chunk.projection.project(vt)
points_h.append(vt[2])
num_valid += 1

points_h.sort()
height = points_h[int(num_valid/2)]
return height


chunk = doc.activeChunk

proj = PhotoScan.GeoProjection() 
proj.init("EPSG::4326")

path = doc.path.rsplit("\\", 1)[0]

processed = 0

t0 = time.time()

for i in range (0, len(chunk.photos)):
photo = chunk.photos[i]
photo.enabled = False

PhotoScan.app.update()

for i in range (0, len(chunk.photos)):


photo = chunk.photos[i]

if (photo.transform == None):
continue

x0 = PhotoScan.Vector((0.0,0.0,0.0))
x1 = PhotoScan.Vector((0.0,0.0,0.0))
x2 = PhotoScan.Vector((0.0,0.0,0.0))
x3 = PhotoScan.Vector((0.0,0.0,0.0))

width = photo.width
height = photo.height

# vectors corresponding to photo corners

v0 = PhotoScan.Vector(( -photo.calibration.cx / photo.calibration.fx, -photo.calibration.cy / photo.calibration.fy, 1))
v1 = PhotoScan.Vector(( (width - photo.calibration.cx) / photo.calibration.fx, -photo.calibration.cy / photo.calibration.fy, 1))
v2 = PhotoScan.Vector(( -photo.calibration.cx / photo.calibration.fx, (height - photo.calibration.cy) / photo.calibration.fy, 1))
v3 = PhotoScan.Vector(( (width - photo.calibration.cx) / photo.calibration.fx, (height - photo.calibration.cy) / photo.calibration.fy, 1))
vc = photo.center

v0.size = 4
v1.size = 4
v2.size = 4
v3.size = 4
vc.size = 4
v0[3] = v1[3] = v2[3] = v3[3] = 0
vc[3] = 1

v0_gc = chunk.transform * photo.transform * v0
v1_gc = chunk.transform * photo.transform * v1
v2_gc = chunk.transform * photo.transform * v2
v3_gc = chunk.transform * photo.transform * v3
vc_gc = chunk.transform * vc

v0_gc.size = 3
v1_gc.size = 3
v2_gc.size = 3
v3_gc.size = 3
vc_gc.size = 3

# surface normal

cen_p = photo.center
cen_p.size = 4
cen_p[3] = 1
cen_t = chunk.transform * cen_p
cen_t.size = 3
cen_t = chunk.projection.project(cen_t)

h = surf_height(chunk, photo)

vloc = PhotoScan.Vector((cen_t[0], cen_t[1], h))
vloc_h = PhotoScan.Vector((cen_t[0], cen_t[1], h))
vloc_h[2] += 1

vloc_gc = chunk.projection.unproject(vloc)
vloc_h_gc = chunk.projection.unproject(vloc_h)
surf_n =  vloc_h_gc - vloc_gc

surf_n.normalize()
v0_gc.normalize()
v1_gc.normalize()
v2_gc.normalize()
v3_gc.normalize()

#intersection with the surface

x0 = intersect(vloc_gc, surf_n, vc_gc, v0_gc)
x1 = intersect(vloc_gc, surf_n, vc_gc, v1_gc)
x2 = intersect(vloc_gc, surf_n, vc_gc, v2_gc)
x3 = intersect(vloc_gc, surf_n, vc_gc, v3_gc)

x0 = chunk.projection.project(x0)
x1 = chunk.projection.project(x1)
x2 = chunk.projection.project(x2)
x3 = chunk.projection.project(x3)

x_0 = min(x0[0],  x1[0], x2[0], x3[0])
x_1 = max(x0[0],  x1[0], x2[0], x3[0])
y_0 = min(x0[1],  x1[1], x2[1], x3[1])
y_1 = max(x0[1],  x1[1], x2[1], x3[1])

x_0 -= (x_1 - x_0) / 20.
x_1 += (x_1 - x_0) / 20.
y_0 -= (y_1 - y_0) / 20.
y_1 += (y_1 - y_0) / 20.

reg = (x_0, y_0, x_1, y_1)

photo.enabled = True
PhotoScan.app.update()
p_name = photo.path.rsplit("/", 1)[1].rsplit(".",1)[0]
p_name = "ortho_" + p_name

proj = chunk.projection   ##if chunk projection units are in meters - dx and dy should be specified in meters:

if chunk.exportOrthophoto(path + "\\" + p_name + ".tif", format = "tif", blending = "average", color_correction = False, projection = proj, region = reg, dx = 2.93022e-06, dy = 2.12914e-06, write_world = True):
processed +=1
photo.enabled = False

for i in range (0, len(chunk.photos)):
photo = chunk.photos[i]
photo.enabled = True

PhotoScan.app.update()

t1 = time.time()

t1 -= t0
t1 = int(t1)

PhotoScan.app.messageBox("Processing finished.\nProcessed "+ str(processed) +" images to orthophotos.\nProcessing time: "+ str(t1)  +" seconds.\nPress OK.")  #information messageNote that in the latter script export resolution is defined in the script body, but you can modify the script to use default resolution.

Hope both of them will work fine for you, but in case of any problems or questions please let me know.

[Gall]

Thank you for these samples, I looked a bit at the extraction of the corners and it seems to work, even though I do not have data to verify the result... Therefore, half of the tasks do not necessarily require the dense point cloud to be generated and I noticed that the sparse point cloud does not contain all the corners point. So, is the dense point cloud garanteed to have all the corners for all images?
In case I don't need to know the altitude, is it possible to retrieve the coordinates in a faster way with something like this?

(Obvously, this code gives wrong result!)

Code: [Select]

currentChunk  = PhotoScan.app.document.activeChunk
camera = currentChunk.cameras[0]
print('\n{}'.format(currentChunk.projection.authority))
print('Image: {}'.format(camera.label))

# 3D coordinates to pixel coordinates
def point3D_to_pixel(x, y, z):
point3D = PhotoScan.Vector([x, y, z])
print('\nDefault 3D: {}'.format(point3D))

point_geocentric = currentChunk.projection.unproject(point3D)
point_internal = currentChunk.transform.inv().mulp(point_geocentric)

return camera.project(point_internal)

# pixel coordinates to 3D coordinates
def pixel_to_point3D(x, y):
point = PhotoScan.Vector([x, y])
point = camera.calibration.unproject(point)
# point = camera.transform.mulv(point)
point = currentChunk.transform.mulp(point)

return currentChunk.projection.project(point)

# EPSG::4326
# Just a test with double conversion-----
imgx, imgy = point3D_to_pixel(3.2881240439016697, 50.4191454593172070, 36.6314926855607150)
print('Pixel: {}, {}'.format(imgx, imgy))

point3D = pixel_to_point3D(imgx, imgy)
print('3D: {}, {}, {}'.format(*point3D))
#----------------------------------------------------

corners = [(0, 0), (camera.width - 1, 0), (camera.width - 1, camera.height - 1), (0, camera.height - 1)]
print()
for c in corners:
point3D = pixel_to_point3D(*c)
print('{} 3D: {}, {}, {}'.format(c, *point3D))

[Alexey Pasumansky]

Hello Gall,

Dense cloud is not used in the scripts that were provided earlier. However, polygonal model should be reconstructed for orthophoto generation and intersecting camera rays with the surface model.

In principle it is possible, that some of the image corners cannot be projected on the reconstructed model, for example, if the reconstruction volume was shrink or shots are very oblique.

[Gall]

Hello Alexey,

this is probably what my problem is, the shots are obliques... One of the goal of this task is to find out if a point (in WGS84) is contained in the image afterward in another processing chain. I'am not an expert in image processing but I was thinking of maybe calcutating the projected real coordinates of the middle point of the image (or principal point) and extrapolate the bounding box coordinates of the image. But, as the projected image on the ground will not a simple rectangle, I am not sure if I can do that... What do you think?

I also ran some tests with the export othophoto sample and I managed to get it to work, even though, after talking with my coworkers, it is not exactly what they need for the projection of the raw image... But I noticed that orthophotos are exported as top XY when I select geographic system (even the full orthophoto) and, after some tests, I don't get what I am supposed to provide in the projection parameter of the exportOrthophoto function to have bottom XY. Some help would be appreciated here!

[Alexey Pasumansky]

Hello Gall,

Projection matrices for the planar views are mentioned in the following post: http://www.agisoft.com/forum/index.php?topic=1322.msg15758#msg15758
however, you'll need to apply additional transformation if gereferenced chunks are used.

The problem with bounding box estimation is in fact that 2D point from single photo (for example, image center) is a ray in 3D space and not the 3D point.
In principle, it is possible to extend the bounding box, but if the images are very oblique (near horizontal), the corners of the images may be on the sky and could not be estimated as 3D point anyhow.

[Gall]

I saw this post and I tried to provide the matrice product of BottomXY with chunk.transform but no orthophoto is exported then. I just don't know which matrices are involved for the resulting projection matrix.

Well, the shots I deal with should be max 45° ~ 50° from nadir so I might give it a try, I have no other leads right now...

[jkova96]

Hi Alexey,

Could you please provide code "transformation" so script will be runable to make job in version 2.1.

Could you also please provide me information what script exactly does, does it specify image footprints corners into txt file or?


Thank you!

J.K.