export individual orthorectified images (corresponding to raw data)?

[MartinTolnai]

Hi!

Is there a way to export every single images i use in a project as individual orthorectified images (not as a mosaic) ?

[Alexey Pasumansky]

Hello Martin,

You can disable all but one photo and use Export Orthophoto option. Such operation can be automated using the Python script, so if you are interested, please, let me know.

[MartinTolnai]

Thank you for the helpfull answer Alexey!
This is the solution I'm looking for. Sure I'm intrested. Is such script existing allready?

[Alexey Pasumansky]

Hello Martin,

We have the script prototype for such task. Please note that it works for WGS84 projection and export resolution is defined in the script body. Also the exported files are generated in the project folder.

Code: [Select]

#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 not photo.transform:
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 message
If anything is not working or any modifications are required, just let me know.

[MartinTolnai]

Thank you very much!
Hopefully I can make the needed modifications.

[aggieair]

This thread is old, but I wanted to reference the script Alexey mentioned above.

We used the above script to at least get the four corner (x,y) of the refortified frames.  Using those we imported it into ArcGIS and developed the polygons that way.  However,


   #We checked the coordinates of this polygon and it is off from Agisoft's output, the corners are 6-10 meters inside the polygon

   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])

   #Why does he do this?  Why divide by 20?
   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.

[Alexey Pasumansky]

Hello aggieair,

I don't remember why it was done, but I assume that the size of the region has been extended by 10% by some reason.

[htinez_ojenra]

Just want to say "Thanks for this topic answered prayer "

[Alexey Pasumansky]

Note that in the version 1.2 there's Export Orthophotos option in the tools menu that allows to do the same from the GUI in a more convenient way.