Forum

Author Topic: Exporting TPs with RMS error in xyz  (Read 8352 times)

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Exporting TPs with RMS error in xyz
« 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:
Code: [Select]
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...
« Last Edit: August 01, 2019, 08:08:56 PM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #1 on: August 03, 2019, 05:42:36 PM »
Any help on this topic?

Best Regards,
Paul Pelletier,
Surveyor

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #2 on: August 05, 2019, 06:57:47 PM »
Hello Paulo,

Do you intend to save the variance values as XYZ components for each tie point with this code?
Best regards,
Alexey Pasumansky,
Agisoft LLC

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #3 on: August 05, 2019, 07:11:49 PM »
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
Code: [Select]
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:
Code: [Select]
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,
« Last Edit: August 06, 2019, 03:36:47 AM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #4 on: August 18, 2019, 07:55:45 PM »
Hello Paul,

The following script should save XYZ coordinates for the tie points, variance vector length and it's XYZ components:
Code: [Select]
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()
Best regards,
Alexey Pasumansky,
Agisoft LLC

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #5 on: August 19, 2019, 03:37:12 AM »
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:

Code: [Select]
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
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)

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)
Code: [Select]
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

« Last Edit: August 20, 2019, 08:04:01 AM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor

antoine billault

  • Newbie
  • *
  • Posts: 9
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #6 on: July 26, 2021, 10:21:10 AM »
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
Quote
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
Code: [Select]
V = chunk.transform.matrix * point.coord
V.size = 3
X, Y, Z = chunk.crs.project(V)

I was wondering if  you still use

Code: [Select]
T = chunk.crs.localframe(T.mulp(chunk.region.center)) * T
THANKS

antoine

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #7 on: July 26, 2021, 11:06:37 AM »
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....

Best Regards,
Paul Pelletier,
Surveyor

antoine billault

  • Newbie
  • *
  • Posts: 9
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #8 on: July 27, 2021, 01:39:35 AM »
thank you Paul

3DWinter

  • Full Member
  • ***
  • Posts: 103
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #9 on: June 06, 2022, 11:07:28 PM »
@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

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #10 on: June 07, 2022, 12:10:47 AM »
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:
Code: [Select]
camera.error(points[point_id].coord, proj.coord).norm() **2where 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...
Best Regards,
Paul Pelletier,
Surveyor

3DWinter

  • Full Member
  • ***
  • Posts: 103
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #11 on: June 08, 2022, 02:19:19 AM »
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.

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #12 on: June 08, 2022, 04:59:06 AM »
3D Winter,

thank you for the kind words...

the point_ids list is built as:
Code: [Select]
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:
Code: [Select]
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,

« Last Edit: June 08, 2022, 05:51:41 AM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor

3DWinter

  • Full Member
  • ***
  • Posts: 103
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #13 on: June 09, 2022, 09:10:14 AM »
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!

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: Exporting TPs with RMS error in xyz
« Reply #14 on: June 09, 2022, 03:45:04 PM »
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
« Last Edit: June 10, 2022, 04:17:44 PM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor