Topics

1

Python and Java API / Export Alignement result to another software

« on: December 11, 2014, 04:13:42 PM »

Hi Everyone,

I use two photogrammetric(/structure from motion) softwares for my scientific research, micmac (command line tools) and photoscan. I like both softwares and I feel very insteresting to compare them in order to  improve my understanding of the fascinating computer vision algorithms. I recently exported the result of the alignement of photoscan in micmac by means of the following script. I would like to share it just in case someone else is willing to use it. In addition, I would thank Alexey for his help.

Jo

Code: [Select]

# compatibility PhotoScan Pro 1.1.0
# saves multiple files with tie-points per photo pair -- Alexey Pasumansky (AgiSoft LLC) (great thanks Alexey!), Jonathan Lisein and Samuel Quevauviller (Ulg), 12/2014
# the aim of the script is to use the alignment result (tie point, camera orientation and calibration) of photoscan in the photogrammetric suite Micmac

import time
import PhotoScan
import os

doc = PhotoScan.app.document
chunk = doc.chunk

print("Export of tie points and camera orientation from Photoscan 1.1.0 to Mimac software : started...")
path = PhotoScan.app.getExistingDirectory("Specify the output folder for Homol and Orientation database :")  #Opens save-to dialog box
Ori = PhotoScan.app.getString("How do you want to name the orientation database? ex : PSmax1000TP", "PS")


t0 = time.time()

pathHomol = os.path.join(path+"/Homol/")
if not os.path.exists(pathHomol):
os.mkdir(pathHomol)

OriBD = "Ori-"+Ori
pathOri = path+ "/" + OriBD
if not os.path.exists(pathOri):
os.mkdir(pathOri)

point_cloud = chunk.point_cloud
point_proj = point_cloud.projections

flag = False

# loop on every camera of the image block
for i in range(0, len(chunk.cameras)):

photo1 = chunk.cameras[i]
#camera.transform is the rotation matrix + translation matrix(= camera center). The following code ligne is a way to test if the camera is aligned.
if not photo1.transform:
continue

#---------------------------export of camera external orientation
# matrice de rotation
rot = photo1.transform
# creation of orientation file for the camera (micmac format)
oriFile = pathOri + "/Orientation-" + photo1.label + ".xml"
file = open(oriFile, "w")
file.write("<?xml version=\"1.0\" ?>\n")
file.write("<ExportAPERO>\n")
file.write(" <OrientationConique>\n")
file.write(" <OrIntImaM2C>\n")
# all these tags are useless in this situation but are still mandatory for compliance with micmac convention
file.write(" <I00>0 0</I00>\n")
file.write(" <V10>1 0</V10>\n")
file.write(" <V01>0 1</V01>\n")
file.write(" </OrIntImaM2C>\n")
file.write(" <TypeProj>eProjStenope</TypeProj>\n")
# get the focal on the camera.  In micmac, the calibration file in named AutoCal + the focal in 1/10 mm
foc = 10*int(photo1.photo.meta["Exif/FocalLength"])
file.write(" <FileInterne>" + OriBD +"/AutoCal" +str(foc) +".xml</FileInterne>\n")
file.write(" <RelativeNameFI>true</RelativeNameFI>\n")
file.write(" <Externe>\n")
# AltiSol (flight altitude) and Profondeur (english: depth) are redundant info which do not need to be accurate here
tiep_proj = point_cloud.projections[photo1][0]
XYZ = point_cloud.points[tiep_proj.track_id].coord
file.write(" <AltiSol>"+ str(photo1.center[2] - XYZ[2]) +"</AltiSol>\n")
file.write(" <Profondeur>"+ str(photo1.center[2] - XYZ[2]) +" </Profondeur>\n")
file.write(" <Time>-1.00000000000000002e+30</Time>\n")
file.write(" <KnownConv>eConvApero_DistM2C</KnownConv>\n")
# camera position center
center = photo1.center
file.write(" <Centre>"+ str(center[0]) + " " + str(center[1]) + " " + str(center[2]) + "</Centre>\n")
file.write(" <ParamRotation>\n")
file.write(" <CodageMatr>\n")
# rotation matrix.
file.write(" <L1>" + str(rot[0, 0]) + " " + str(rot[0, 1]) + " " + str(rot[0, 2]) + "</L1>\n")
file.write(" <L2>" + str(rot[1, 0]) + " " + str(rot[1, 1]) + " " + str(rot[1, 2]) + "</L2>\n")
file.write(" <L3>" + str(rot[2, 0]) + " " + str(rot[2, 1]) + " " + str(rot[2, 2]) + "</L3>\n")
file.write(" </CodageMatr>\n")
file.write(" </ParamRotation>\n")
file.write(" </Externe>\n")
file.write(" <Verif>\n")
file.write(" <Tol>0.00100000000000000002</Tol>\n")
file.write(" <ShowMes>true</ShowMes>\n")

# add 3 tie points for verification purpose, required in micmac
for i in (0, 1, 2):
tiep_proj = point_cloud.projections[photo1][i]
XYZ = point_cloud.points[tiep_proj.track_id].coord
file.write(" <Appuis>\n")
file.write(" <Num>"+str(i)+"</Num>\n")
# image measurement of the tie point : U, V
file.write("        <Im>" + str(tiep_proj.coord[0]) + " " + str(tiep_proj.coord[1]) + "</Im>\n")
# model measurement of the tie point : X Y Z
file.write(" <Ter>" + str(XYZ[0]) + " " + str(XYZ[1]) + " " + str(XYZ[2]) + "</Ter>\n")
file.write(" </Appuis>\n")
file.write(" </Verif>\n")
file.write(" <ConvOri>\n")
file.write(" <KnownConv>eConvApero_DistM2C</KnownConv>\n")
file.write(" </ConvOri>\n")
file.write("</OrientationConique>\n")
file.write("</ExportAPERO>\n")
file.close()
print("End of orientation export for camera "+str(photo1.label))

# create the folder for storing the set of tie points for this image
pathPastis1=os.path.join(pathHomol+"Pastis"+photo1.label)                   
if not os.path.exists(pathPastis1):
os.mkdir(pathPastis1)

# loop on all the remainder cameras of the image block (=chunk)
for j in range(i + 1, len(chunk.cameras)):

photo2 = chunk.cameras[j]
# is the camera aligned?
if not photo2.transform:
continue

pathPastis2=os.path.join(pathHomol+"Pastis"+photo2.label)                   
if not os.path.exists(pathPastis2):
os.mkdir(pathPastis2)

# associative arrays = dictionary object class in python
matches1 = dict()
matches2 = dict()

for proj in point_proj[photo1]:
# the dictionnary index for this projection (U, V) is the track index of the tie point (X,Y,Z)
matches1[proj.track_id] = proj.coord
for proj in point_proj[photo2]:
matches2[proj.track_id] = proj.coord

# check if there are tie points which are shared by the camera pair
if len(set(matches1.keys()).intersection(set(matches2.keys()))):

# Tie point of the camera pair are stored in two (redundant) txt files:
fileHomol12 = open(pathPastis1 + "\\" + photo2.label + ".txt", "wt")
fileHomol21 = open(pathPastis2 + "\\" + photo1.label + ".txt", "wt")

# In photoscan, there are a few tie point duplicate in the sparce point cloud. Duplicates are not allowed in micmac. These duplicate have same value for model measurements X,Y,Z but different track_id.
# Dictionaries do not contain duplicate keys. The following code remove the tie point duplicates
list=[]
for key in matches1.keys():
u,v=matches1[key]
list.append(str(u)+";"+str(v)+"#"+str(key))
list.sort()
matches11=dict()
for i in range(len(list)-1):
ss1=list[i].split("#")
ss2=list[i+1].split("#")
if ss1[0]!=ss2[0]:
ss3=ss1[0].split(";")
matches11[int(ss1[1])]=[ss3[0],ss3[1]]
# same removal of duplicates but on the second matches list
list=[]
for key in matches2.keys():
u,v=matches2[key]
list.append(str(u)+";"+str(v)+"#"+str(key))
list.sort()
matches21=dict()
for i in range(len(list)-1):
ss1=list[i].split("#")
ss2=list[i+1].split("#")
if ss1[0]!=ss2[0]:
ss3=ss1[0].split(";")
matches21[int(ss1[1])]=[ss3[0],ss3[1]]

# loop on every tie points shared by the two camera
for tiepoint in set(matches11.keys()).intersection(set(matches21.keys())):

u1, v1 = matches1[tiepoint]
u2, v2 = matches2[tiepoint]

# u1, v1, u2, v2
fileHomol12.write("{:.2f}".format(u1) + " {:.2f}".format(v1) + " {:.2f}".format(u2) + " {:.2f}\n".format(v2))
#  u2, v2, u1, v1,
fileHomol21.write("{:.2f}".format(u2) + " {:.2f}".format(v2) + " {:.2f}".format(u1) + " {:.2f}\n".format(v1))

fileHomol12.close()
fileHomol21.close()

for calib in chunk.sensors:

# generation of the calibration file / distorsion model "RadialBasic"
oriFile = pathOri + "/AutoCal" + str(int(10*calib.focal_length))+ ".xml"
file = open(oriFile, "w")
file.write("<?xml version=\"1.0\" ?>\n")
file.write("<ExportAPERO>\n")
file.write(" <CalibrationInternConique>\n")
file.write(" <KnownConv>eConvApero_DistM2C</KnownConv>\n")
file.write(" <PP>"+str(calib.calibration.cx) + "  " +  str(calib.calibration.cy)   +"</PP>\n")
file.write(" <F>"+ str(calib.calibration.fx) + "</F>\n")
file.write(" <SzIm>" + str(calib.calibration.width) + "  " +  str(calib.calibration.height) +"</SzIm>\n")
file.write(" <CalibDistortion>\n")
file.write(" <ModRad>\n")
file.write(" <CDist>" + str(calib.calibration.cx) + "  " +  str(calib.calibration.cy) +"</CDist>\n")
# set additionnal parameters of the calibration to 0 because micmac convention are different and I was lazy to convert the calibration model.
file.write(" <CoeffDist>0</CoeffDist>\n")
file.write(" <CoeffDist>0</CoeffDist>\n")
file.write(" <PPaEqPPs>true</PPaEqPPs>\n")
file.write(" </ModRad>\n")
file.write(" </CalibDistortion>\n")
file.write(" </CalibrationInternConique>\n")
file.write("</ExportAPERO>\n")
file.close()


t1 = time.time()
t1 -= t0
t1 = float(t1)

print("Script finished in " + "{:.2f}".format(t1) + " seconds.")