12746
« on: September 23, 2014, 03:06:07 PM »
Добрый день,
Ниже пример такого скрипта. Обратите внимание, что разрешение экспорта указано константами в теле скрипта в градусах, т.е. применимо к системе WGS84. При использовании спроецированных СК потребуется указание разрешения экспорта в метрах.
#Batch export of orthophotos based on each photo
#compatibility Agisoft PhotoScan Pro 1.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 message
Show Posts
