Agisoft Metashape
Agisoft Metashape => Python and Java API => Topic started by: mwillis on July 10, 2014, 11:31:32 AM
-
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
-
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?
-
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
-
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
-
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).
-
Thanks for all your help Alexey.
-
Hello again, does anyone have a script that will export each individual image in a chunk as its own orthophoto?
Thanks
-
Hello Mark,
We have a sample code, but please note that export resolution is defined in the script body in degrees (for WGS84 export):
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
-
Alexey, thank you for always being so helpful and for your excellent coding skills.
Mark
-
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.
-
Yep I did a save as.
-
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.
-
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.
-
Oops!
I see it made a new command Export individual orthophotos in the pull down Custom.
-
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
-
Hello bisenberger,
By default the orthophotos are saved to the folder with the project file.
But you can change the assignment of the path variable in def exp_ortho(self) function to:
path = PhotoScan.app.getExistingDirectory("Browse for the export folder:")
-
Hi Alexey,
I went ahead and downloaded the rest of the scripts posted on the wiki.
http://wiki.agisoft.com/wiki/Python
All are very useful.
Thanks,
Bill
-
Works nicely. I can export large areas with low RAM (8Gb) and then recombine in GlobalMapper
-
Dear Alexey,
I have an problem while running your script.
My Dataset aint got no camera position files, but an generic alignment and just markers for generating orthophotos.
If i run the script after proceeding every single step i`ve got the following error:
Traceback (most recent call last):
File "G:/2012-01_STC_Venezuela_Lake-Valencia/Thomas/prio5/ost/prj/try.py", line 83, in <lambda>
proc_exp = lambda : self.exp_ortho()
File "G:/2012-01_STC_Venezuela_Lake-Valencia/Thomas/prio5/ost/prj/try.py", line 272, in exp_ortho
v1 = PhotoScan.Vector((crd[0], crd[1], 0) )
TypeError: 'NoneType' object is not subscriptable
I would be delighted to read from you.
Best regards Tommek
-
Hello Tommek,
Could you please specify, what coordinate system is defined in the Reference pane settings dialog?
-
Hey Akexey,
thanks for the reply. It's WGS84/UTM zone 19N.
Regards
-
Hello Tommek,
Thanks, it seems that wrong criterion has been used when checking if the coordinate system is projected or geographic. I've modified the script (changed only one line):
...
#recalculating WGS84 resolution from degrees into meters if required
if chunk.crs:
if not ('PROJCS' in proj.wkt):
...
So please check if it works fine now.
-
It works. That ist amazing.
thanks a bunch!!!!! ;D
-
OK, nice to hear that.
Note that in the next version (that is currently available as pre-release), we've implemented similar feature to GUI interface for user convenience.
-
Hey Alexey,
I have one question in addition. Is it possible to include an all over color correction befor exporting the single orthos?
Blessings Tommek
-
Hello Tommek,
At the moment it is not possible, since color correction is only applicable to the images used for the orthophoto export, and if there's only one image per individual orthophoto, color correction couldn't be applied.
-
Hi Alexey,
In running Export Individual Orthophotos.py I am not getting a result from the script when I run with the option for a single selected image. In the console, although the process appears to be working it exits with code 0 and fails to produce an ortho. Just wondering if I might be executing something improperly. Thanks for any advice you might be able to offer.
Rob
-
Hello Rob,
Maybe you can send the screenshot of the script dialog and the full Console pane ouput related to the export operation?
But probably, you can update to the version 1.2 and use built-in functionality in Tools Menu -> Export -> Export Orthophotos to save individual images. This functionality can be used after the orthomosaic is built using the corresponding option in Workflow menu.
-
Thanks for the reply Alexey,
I'm attaching a screen shot of the console in addition to the startup and finish dialog from running the script. I'll update this evening and give it a shot.
Rob
--
-
Hello Rob,
I have two ideas of possible reasons causing failed export:
1) Export resolution is two high so the size of individual orthophotos is higher than 4 GB for GeoTIFFs.
2) The document has not been saved prior to the script start. Script exports the orthophotos to the folder containing the project file.
-
Thanks for the tips Alexey,
By the way I did upgrade to 1.2, however, I'm not sure how to export individual orthos. I did save the psz to pox but the Tools > Export > Export Orthos option appears greyed out. I thought perhaps I needed to build the orthomosaic first however, when I tried to do so I receive and empty extent error (without the setup boundary check box) and empty camera list (with it checked). (Console Log Attached) Don't know if you have any other pointers.
Thanks again for any advice.
Rob
--