Agisoft Metashape
Agisoft Metashape => Python and Java API => Topic started by: Paulo on July 31, 2019, 01:08:41 AM
-
Hi all,
I am trying to write a script to export to a file the coordinates of all Tie points (after optimizing cameras with Estimate Tie point covariance option) along with the RMSe for X, Y and Z...
from each tie point covariance matrix, I extract diagonal values to calculate X, Y , Z RMSe as in:
chunk = Metashape.app.document.chunk
M = chunk.transform.matrix
s = chunk.transform.scale
points = chunk.point_cloud.points
npoints = len(points)
for point_id in range(npoints):
if not points[point_id].valid:
continue
if chunk.crs:
V = M * (points[point_id].coord)
V.size = 3
#print (V)
X = chunk.crs.project(V)
else:
V = M * (points[point_id].coord)
X = V
rmsXYZ = s * Metashape.Vector([math.sqrt(points[point_id].cov[0,0]), math.sqrt(points[point_id].cov[1,1]), math.sqrt(points[point_id].cov[2,2])])
However, I think this represents RMSe along internal X, Y , Z coordinate system....
How do I get the RMSe values corresponding to X, Y, Z of CRS projected coordinate system?
Any help will be greatly appreciated...
-
Any help on this topic?
-
Hello Paulo,
Do you intend to save the variance values as XYZ components for each tie point with this code?
-
Yes Alexey,
i would like to save in a TPs file, the X, Y, Z, and Error X, Error Y, Error Z for each point of course related to crs projected system as in
Coordenadas de puntos de enlace con sus errores RMS en X,Y,Z, error de reproyeccion promedio y numero de proyecciones
Sistema de coordenadas WGS 84 / UTM zone 16N + GGM2010 geoid height
Index X Y Z Error_X(cm) Error_Y(cm) Error_Z(cm) Error_total(cm) num_obs errrep(pix)
0 236376.294 2325813.957 8.117 2.649 3.897 5.733 7.421 2 0.508
Right now I am using this piece of code to transform error components from internal CS to projected CRS:
M = chunk.transform.matrix
s = chunk.transform.scale
points = chunk.point_cloud.points
npoints = len(points)
....
for point_id in range(npoints):
if not points[point_id].valid:
continue
if chunk.crs:
V = M * (points[point_id].coord)
V.size = 3
#print (V)
X = chunk.crs.project(V)
m = chunk.crs.localframe(V)
R = m.rotation() * M.rotation()
rmsXYZ = s * Metashape.Vector([math.sqrt(points[point_id].cov[0,0]), math.sqrt(points[point_id].cov[1,1]), math.sqrt(points[point_id].cov[2,2])])
rmsXYZ = R * rmsXYZ
rmsXYZ = Metashape.Vector([abs(rmsXYZ.x),abs(rmsXYZ.y),abs(rmsXYZ.z)])
rmsXYZ = 100 * rmsXYZ
I do appreciate your input,
-
Hello Paul,
The following script should save XYZ coordinates for the tie points, variance vector length and it's XYZ components:
import Metashape, math
path = Metashape.app.getSaveFileName("Specify the export file path:", filter = "Text file (*.txt);;All formats (*.*)")
file = open(path, "wt")
file.write("Id\tX\tY\tZ\tvar\tcov_x\tcov_y\tcov_z\n")
chunk = Metashape.app.document.chunk
T = chunk.transform.matrix
if chunk.transform.translation and chunk.transform.rotation and chunk.transform.scale:
T = chunk.crs.localframe(T.mulp(chunk.region.center)) * T
R = T.rotation() * T.scale()
for point in chunk.point_cloud.points:
if not point.valid:
continue
cov = point.cov
coord = point.coord
coord = T * coord
cov = R * cov * R.t()
u, s, v = cov.svd()
var = math.sqrt(sum(s)) #variance vector length
vect = (u.col(0) * var)
file.write(str(point.track_id))
file.write("\t{:.6f}\t{:.6f}\t{:.6f}\t{:.6f}".format(coord[0], coord[1], coord[2], var))
file.write("\t{:.6f}\t{:.6f}\t{:.6f}".format(vect.x, vect.y, vect.z))
file.write("\n")
file.close()
-
Thanks a thousand, Alexey.
I have adapted your code to my script....
the only detail is that it prints out the X,Y,Z coordinates of TiePts in LRS (East, North, Up) with origin at region center as seen in:
Id X Y Z var cov_x cov_y cov_z
21 -174.021457 -396.605524 640.282965 0.033728 -0.006520 -0.002776 0.032975
so I changed coord = T * point.coord
with V = chunk.transform.matrix * point.coord
V.size = 3
X, Y, Z = chunk.crs.project(V)
also I use directly point.id instead of point.track_id and absolute value of each vect component...
and resulting first lines of export file is (RMSe values in cm)
Coordenadas de puntos de enlace con sus errores RMS en X,Y,Z, numero de proyecciones y error de reproyeccion promedio
Sistema de coordenadas WGS 84 / UTM zone 14N + GGM2010 geoid height
Index X Y Z Error X(cm) Error Y(cm) Error Z(cm) Total(cm) num_obs errrep(pix)
0 346248.403 2852992.029 659.232 0.652 0.278 3.298 3.373 8 0.380
-
Hello Paul
I need to perform the same task
I wanted to ask if you can clarify your modification you pointed out in your last post
so I changed
Code: [Select]
coord = T * point.coord
with
Code: [Select]
V = chunk.transform.matrix * point.coord
V.size = 3
X, Y, Z = chunk.crs.project(V)
As you introduce V = chunk.transform.matrix * point.coord
V.size = 3
X, Y, Z = chunk.crs.project(V)
I was wondering if you still use
T = chunk.crs.localframe(T.mulp(chunk.region.center)) * T
THANKS
antoine
-
Hello Antoine,
the code change was to print out the coordinates of each tiepoint in chunk crs....
the line T = chunk.crs.localframe(T.mulp(chunk.region.center)) * T is still used to calculate the matrix R that lets you transform the covariance matrix from internal CS to LSE CS....
-
thank you Paul
-
@Paulo,
Is num_obs how many cameras (or how many observations-camera) are associated with the tiepoint?
What is errrep(pix) ?
How did you cacluate/extract num_obs errrep(pix)?
I wonder if num_obs can give a sense to the population that a variance is calculated for each tie point.
Thanks
-
Hello 3DWinter,
the num_obs for each Tie Point is the number of images where it is observed. reperr(pix) is the RMS reprojection error for the number of observations i.e. for a given camera with a tie point and its measurement proj on the camera then the following gives the square of reprojection error:
camera.error(points[point_id].coord, proj.coord).norm() **2
where points[point_id] is the tie point (point in sparse cloud with id = point_id)
For each proj whre the tie point is observed you sum these squared errors, divide by num_obs and extract square root to get RMS reperror for the point...
-
Paul,
I don't understand what is the camera object in your code. Is camera based on :
for camera in chunk.cameras:
camera.error()
If so, I don't understand how you are mixing camera from cameras iteration and point from chunk.point_cloud.points?
Also, you are using a norm() method, but I am also interested in the number of observations for a single tie-point (num_obs). Is the num_obs and attribute or method within point?
You are very knowledgeable about Metashape API , way outside what I have found in the Metashape API manual. Thanks for your patience to explain and provide code snippets.
-
3D Winter,
thank you for the kind words...
the point_ids list is built as:
point_cloud = chunk.point_cloud
projections = point_cloud.projections
points = point_cloud.points
npoints = len(points)
tracks = point_cloud.tracks
point_ids = [-1] * len(point_cloud.tracks)
for point_id in range(0, npoints):
point_ids[points[point_id].track_id] = point_id
then you iterate on chunk.cameras and then projections for each camera as:
points_errors = {}
for camera in chunk.cameras:
if camera.type == Metashape.Camera.Type.Keyframe:
continue # skipping Keyframes
if not camera.transform:
continue
all_proj += len(projections[camera])
for proj in projections[camera]:
track_id = proj.track_id
point_id = point_ids[track_id]
if point_id < 0:
continue
if not points[point_id].valid:
continue
dist = camera.error(points[point_id].coord, proj.coord).norm() **2
distn = camera.error(points[point_id].coord, proj.coord).norm()
if point_id in points_errors.keys():
point_id = int(point_id)
points_errors[point_id].x += dist # squared
points_errors[point_id].y += distn # norm
points_errors[point_id].z += 1
else:
points_errors[point_id] = PhotoScan.Vector([dist, distn, 1])
for a given point_id in point cloud:
num_obs = points_errors[point_id].z
RMS reperr = math.sqrt(points_errors[point_id].x / points_errors[point_id].z)
avg reperr = points_errors[point_id].y / points_errors[point_id].z
hope this is clearer,
-
Paul,
Thanks for the deeper explanation. It is impressive your mastery of the Metashape API despite the limited reading content in the manual and assuming you are work/worked for Metashape.
First, is a "Track" the line that connects two adjacent tie-points that RANSAC filters?
Second, could you give me some advice on how to better understand Metashape API functionality and be able to export the things I am interested in?
Third, I supposed a similar analysis can be done on the dense cloud? Do you have a blog/or book about your work with Metashape API?
Again thanks a million!
-
Hi 3dwinter,
my knowledge has come from a lot of searching, looking at forum posts and painfully trying to understand the concepts....a question of time and experience
From my undertanding tracks are all the matched feature points in the project (valid and invalid). The valid tracks are in fact the tie points (point_cloud.points). The invalid tracks are matched points that have been rejected by software according to some criteria after adjustment or tie points deleted after gradual selection. Valid tracks or tie points appear blue in image while invalid tracks appear greyish...
so len(chunk.point_cloud.points) < len((chunk.point_cloud.tracks) and points are a sub set of tracks...see attachment for better visual understanding...
Dense cloud is completely different to point cloud as they result from depth maps generation while tie points (point_cloud.points) result from image matching and adjustment....
2nd example shows a track id=2504 matched in 6 images not a tie point as it did not pass the adjustment criteria (colored grey in images) ... or a tie point deleted after gradual selection
-
Dear Paulo,
I'm involved in adjusting a script by adding the BA-derived tie points precisions.
In the previous discussion and posted code snippet you introduced the tranformation matrix T, which basically allows to transform the model (or "internal") coordinates into a local cartesian coordinates system. I was wondering why did you use this reference system for precisions computation insted of the geocentric (or the world reference) one? Am I missing some detail?
Many thanks.
Mauro
-
Hola Mauro,
the idea in calculating TPs variance in LRS or local East North Up system reference is that its axii are representative RMS TP east, North up, while using RMS in geocentric is not meaningful and using PCS (ie. UTM or other) is not possible as conversion from geocentric to PCS cannot be represented by a 4x4 transformation matrix which needed in the propagation of variance-covariance matrix from internal CS...