Auto exporting of multiple single photos to ortho

[mwillis]

Hi everyone, I'm trying to figure out how to do the following in PhotoScan with Python.

As an example, I have 20 photographs making up one 3D model.  I would like to automatically export one orthophoto of each individual image in the "Top XY" projection plane.  So that in the end, I have 20 separate orthophotos for each of the original photographs.

The number of photographs and the names of the photographs will change in each project.  The data is also planar without geographic coordinates.

Any help is greatly appreciated.

Thanks!

Mark

[Alexey Pasumansky]

Hello Mark,

Just to make sure, I've understood the task correctly: you have aligned set of cameras and mesh model has been already reconstructed. All you need is just to export set of orthophotos based on individual cameras in chunk using the same resolution and projection (TopXY), right?

What about export resolution? Are your projects scaled?

[mwillis]

Alexey, yes.  I have a set of fully aligned photos and mesh model built with texture. And yes, I need to export a set of orthophotos based on individual cameras in the chunk.  Resolution for each orthophoto and projection (TopXY) should be the same.

I would like to export at the same resolution as if I were exporting all the files as a single orthophoto.  The project are not scaled.

I am taking photographs of a large panel of prehistoric rock art with neutral light.  Then I'm taking detailed photographs of the rock art with light from various directions to bring out faint details.  My point is bring the neutral data into Photoshop as an orthophoto background and then have individual layers for each of the detailed photographs.  Does that make sense?

How I do it now:
I export neutral lighting 3D model (big) orthophoto (top XY) of entire panel
Next I  retexture the model with a single photograph and export each detailed photograph as single orthophoto, same resolution as above and with the same extents as the big orthophoto (same extent as step above), TopXY
Then I import these as layers into Photoshop.

The point is to be able to see all the detailed photographs with good lighting for tracing.

I can do this manually but very time consuming.

Thanks for all your help.

Mark

[mwillis]

Hi Alexey, do you think it is possible?  You don't need to code it for me but I'm not sure how to deal with variable names of the single photos.  Maybe it can't be done?

Thanks for time,

Mark

[Alexey Pasumansky]

Hello Mark,

To export orthophoto based on the single camera only, you need to disable all but one cameras in the project.
You can do it in the loop by enabling new camera and disabling the previous one, for example (providing that all the cameras have been initially disabled).

For TopXY projection you need to use corresponding matrix as projection argument (for unreferenced chunks it should be an identity matrix).

[mwillis]

Thanks for all your help Alexey.

[mwillis]

Hello again, does anyone have a script that will export each individual image in a chunk as its own orthophoto?

Thanks

[Alexey Pasumansky]

Hello Mark,

We have a sample code, but please note that export resolution is defined in the script body in degrees (for WGS84 export):

Code: [Select]

import os, 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 message



[mwillis]

Alexey, thank you for always being so helpful and for your excellent coding skills.

Mark

[Alexey Pasumansky]

Have posted more convenient script for batch orthophoto generation based on individual cameras:
http://wiki.agisoft.com/wiki/Export_Individual_Orthophotos.py

Any comments are appreciated.

[bisenberger]

Yep I did a save as.

[Alexey Pasumansky]

Hello bisenberger,

You need to copy only the part related to the script body, not the whole page. So actually "save as" for the link wouldn't work.

[bisenberger]

Hello Alexey,
I copied the script body, but it doesn't do anything when I run it. I've attached a rar copy of the script file I'm trying to use.

[bisenberger]

Oops!
I see it made a new command Export individual orthophotos in the pull down Custom.

[bisenberger]

It seems to be working, but I can't tell where it is saving them. Can there be an option to indicate what directory they are saved in?

Never mind I see it is saving them in the same folder as the project file (.psz). 

Thanks Alexey