Export measurement precision values as scalar field into point cloud

[Daniele]

Hello Alexy or everyone could help me :-),
I am totally new to this forum so I apologize if this thread has been already addressed, but I can't find a solution to my issue.
I am using metashape to build  dense point clouds of underwater habitats and perform in CloudCompare a change detection analysis. For this task, I would like to use the  M3C3 plugin (https://www.cloudcompare.org/doc/wiki/index.php?title=M3C2_(plugin)). However, to fill the Precision maps tab to enables the calculation of detectable change using measurement precision values stored in scalar fields of point clouds, I would need such information stored in the exported cloud from metashape. If I well understood such precision information could be assumed to be the reprojection errors, as stated by Maria in this other topic: https://www.agisoft.com/forum/index.php?topic=3771.msg20283#msg20283
However, because I am using Metashape 1.7.2  I suppose that the same script needs to be updated, am I right? If a new version is needed, please could you help me because I am totally foreign to python scripting
Finally, do you suggest using sparse point cloud or dense point cloud for such computation?  I read that the dense matching process does not optimise any aspects of the image network and, therefore, does not affect the underlying precision estimates, indeed tie point precision could be used to represent the main measurement contribution to surface model precision. Could you confirm such information?

Many thanks for your time and support in advance.

All the best,
Daniele

[Daniele]

Hi Alexey,
In the meanwhile, I run the script and a text file is generated so the script still works also for methashape. However, I checked that the reproject error is not my target because I would need precision estimates (sigma) of each point for X,Y and Z components (sigmaX, sigmaY and sigmaZ). If you have any suggestions to achieve this result I will be grateful for your support.

Thanks again
Daniele

[Alexey Pasumansky]

Hello Daniele,

Please check the following script, which should export Sigma values:

Code: [Select]

import Metashape, math

def rotation_dx(alpha):
sina = math.sin(alpha)
cosa = math.cos(alpha)
return Metashape.Matrix([[0, 0, 0], [0, -sina, -cosa], [0, cosa, -sina]])

def rotation_dy(alpha):
sina = math.sin(alpha)
cosa = math.cos(alpha)
return Metashape.Matrix([[-sina, 0, cosa], [0, 0, 0], [-cosa, 0, -sina]])

def rotation_dz(alpha):
sina = math.sin(alpha)
cosa = math.cos(alpha)
return Metashape.Matrix([[-sina, -cosa, 0], [cosa, -sina, 0], [0, 0, 0]])

chunk = Metashape.app.document.chunk

path = Metashape.app.getSaveFileName("Specify the export file path:", filter = "Text file (*.txt);;All formats (*.*)")

file = open(path, "wt")

file.write("#Cameras\n")
file.write("#Label\tX\tY\tZ\tYaw\tPitch\tRoll\tSigmaX\tSigmaY\tSigmaZ\tSigmaYaw\tSigmaPitch\tSigmaRoll\n")
for camera in chunk.cameras:
if not camera.transform:
continue
T = chunk.transform.matrix

if chunk.transform.translation and chunk.transform.rotation and chunk.transform.scale:
T = chunk.crs.localframe(T.mulp(camera.center)) * T

R = T.rotation() * T.scale()

transform = T * camera.transform
coord = chunk.crs.project(chunk.transform.matrix.mulp(camera.center))
ypr = Metashape.Utils.mat2ypr(transform.rotation() * Metashape.Matrix.Diag((1, -1, -1)))
opk = Metashape.Utils.mat2opk(transform.rotation() * Metashape.Matrix.Diag((1, -1, -1)))
covO = camera.location_covariance
covR = camera.rotation_covariance
if not covO or not covR:
continue
covO = R * covO * R.t()

R0 = transform.rotation()

Rz = Metashape.Utils.ypr2mat(Metashape.Vector([ypr[0], 0, 0]))
Rx = Metashape.Utils.ypr2mat(Metashape.Vector([0, ypr[1], 0]))
Ry = Metashape.Utils.ypr2mat(Metashape.Vector([0, 0, ypr[2]]))

J1 = Metashape.Matrix([
list((-1) * rotation_dz(- ypr[0] * math.pi / 180) * Rx * Ry),
list(Rz * rotation_dx(ypr[1] * math.pi / 180) * Ry),
list(Rz * Rx * rotation_dy(ypr[2] * math.pi / 180))])

J2 = Metashape.Matrix([
list(R0 * rotation_dx(0) * Metashape.Matrix.Diag((1, -1, -1))),
list(R0 * rotation_dy(0) * Metashape.Matrix.Diag((1, -1, -1))),
list(R0 * rotation_dz(0) * Metashape.Matrix.Diag((1, -1, -1)))
])

J = (J2 * J2.t()).inv() * (J2 * J1.t())

J = J * (180 / math.pi)

covR = J * covR * J.t()

file.write(camera.label)
file.write("\t{:.6f}\t{:.6f}\t{:.6f}".format(coord[0], coord[1], coord[2]))
file.write("\t{:.6f}\t{:.6f}\t{:.6f}".format(ypr[0], ypr[1], ypr[2]))
file.write("\t{:.6e}\t{:.6e}\t{:.6e}".format(math.sqrt(covO[0, 0]), math.sqrt(covO[1, 1]), math.sqrt(covO[2, 2])))
file.write("\t{:.6e}\t{:.6e}\t{:.6e}".format(math.sqrt(covR[1, 1]), math.sqrt(covR[2, 2]), math.sqrt(covR[0, 0])))
file.write("\n")

file.write("#Markers\n")
file.write("#Label\tX\tY\tZ\tSigmaX\tSigmaY\tSigmaZ\n")
for marker in chunk.markers:
if not marker.position:
continue
T = chunk.transform.matrix

if chunk.transform.translation and chunk.transform.rotation and chunk.transform.scale:
T = chunk.crs.localframe(T.mulp(marker.position)) * T

R = T.rotation() * T.scale()

coord = chunk.crs.project(chunk.transform.matrix.mulp(marker.position))
covO = marker.position_covariance
if not covO:
continue
covO = R * covO * R.t()

file.write(marker.label)
file.write("\t{:.6f}\t{:.6f}\t{:.6f}".format(coord[0], coord[1], coord[2]))
file.write("\t{:.6e}\t{:.6e}\t{:.6e}".format(math.sqrt(covO[0, 0]), math.sqrt(covO[1, 1]), math.sqrt(covO[2, 2])))
file.write("\n")

print("Script finished.")
file.close()


[Daniele]

Many Thanks, Alexey for your kind reply and great support!
I run the script provided and in the output, I found the sigma values associated with camera positions and markers. However, I was looking for different results since in the M3C3 plugin of CloudComapre the required fields are the sigma values concerning sparse point cloud (to build precision map). I apologize because probably my first request was not very clear.
I found a very useful discussion and python script concerning this aspect in the valuable paper of  James et al., (2017) available here https://onlinelibrary.wiley.com/doi/10.1002/esp.4125. From this paper, I read that some other software such as
‘Vision Measurement System’ (VMS; http://www. geomsoft.com)  provides point precision as standard output, so I suppose that this feature could be very useful if integrated into the next updates of Metashape. What do you think about this?

Thanks again for your time.

All the best,
Daniele

[Alexey Pasumansky]

Hello Daniele,

The covariance matrix for the points of the sparse point cloud can be accessed by:

Code: [Select]

for point in chunk.point_cloud.points:
   cov_matrix = point.covWith that you should be able to create your custom tie point cloud exporter.

You may adapt the following script for your needs, if the project is georeferenced:

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()
print("Script finished.")

[Daniele]

Hello Daniele,

The covariance matrix for the points of the sparse point cloud can be accessed by:

Code: [Select]

for point in chunk.point_cloud.points:
   cov_matrix = point.covWith that you should be able to create your custom tie point cloud exporter.

You may adapt the following script for your needs, if the project is georeferenced:

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()
print("Script finished.")

Many thanks again Alexey!!
The covariance matrix you cited in the first part of the script is the same that can be also accessed through the camera optimization panel, isn't it?
Regarding the georeferencing of the project, is it sufficient to have GCP with coordinates or you meant that also the cameras should have coordinates? This aspect for me is a key point because underwater I can't use GPS for direct georeferencing the images so I have only some GCP distributed inside the mapping area.

Best regards,
Daniele