1
Python and Java API / Re: Reprojection Errors (script x chunk info)
« on: December 09, 2019, 02:27:45 PM »
Thank you Alexey!
Looking forward for your test results.
Looking forward for your test results.
This section allows you to view all posts made by this member. Note that you can only see posts made in areas you currently have access to.
def getPointErrors(path, i=[0]): # função: erro por ponto
i[0] += 1
f = open(path + "/_1_PONTOS_" + str(i) + ".txt", "wt")
f.write("#N_idx\tCoordX_m\tCoordY_m\tCoordZ_m\tErro_pix\t\n")
point_cloud = chunk.point_cloud
points = point_cloud.points
npoints = len(points)
projections = chunk.point_cloud.projections
M = chunk.transform.matrix # matriz de transformação
point_ids = [-1] * len(point_cloud.tracks)
point_errors = {}
for point_id in range(0, npoints):
point_ids[points[point_id].track_id] = point_id
for cam in chunk.cameras:
if not cam.transform:
continue
for proj in projections[cam]:
track_id = proj.track_id
point_id = point_ids[track_id]
if point_id < 0:
continue
point = points[point_id]
if not point.valid:
continue
error = cam.error(point.coord, proj.coord).norm() ** 2
if point_id not in point_errors.keys():
point_errors[point_id] = [error]
else:
point_errors[point_id].append(error)
for point_id in range(npoints):
if not points[point_id].valid:
continue
if chunk.crs:
w = M * points[point_id].coord
w.size = 3
X, Y, Z = chunk.crs.project(w)
errpto = sum(point_errors[point_id]) / len(point_errors[point_id])
errpto = errpto ** (5/10) # sqrt
f.write("{:6d}\t{:.6f}\t{:.6f}\t"
"{:.6f}\t{:.6f}\n".format(point_id, X, Y, Z, errpto))
f.close()
Hello Luiz,Yes the code you provided calculates precisely RMS and max reprojection for the chunk.
Can you please check, if the following script gives you the correct values for max and RMS reprojection errors:
I think that the different comes from the confusion - max. error is calculated per projection and not per point.I'm guessing you are quite right about the difference using points vs projections on the codes. Ok then, max reprojection error is the max error value within all projections and RMS is the reprojection error averaged over all valid tie points (and their projections) in chunk. Then with the code I'm using I get the reprojection error for each valid tie point (where every single one has at least 2 or more projections and respectively reprojections errors averaged into one value). Is this right?
def tiepointsRMS(path, i=[0]):
i[0] += 1
"""
Generates a TXT file with point index of Sparse Cloud, associated
coordinates and calculated reprojection error (pIndex, X, Y, Z, error)
"""
point_cloud = chunk.point_cloud
projections = point_cloud.projections
points = point_cloud.points
npoints = len(points)
M = chunk.transform.matrix
file = open(path + str(i) + ".txt", "wt")
print("Tiepoints reprojection error calculations started...")
t0 = time.time()
points_errors = {}
for photo in chunk.cameras:
if not photo.transform:
continue
T = photo.transform.inv()
calib = photo.sensor.calibration
pIndex = 0
for proj in projections[photo]:
track_id = proj.track_id
while pIndex < npoints and points[pIndex].track_id < track_id:
pIndex += 1
if pIndex < npoints and points[pIndex].track_id == track_id:
if not points[pIndex].valid:
continue
coord = T * points[pIndex].coord
coord.size = 3
dist = calib.error(coord, proj.coord).norm() ** 2
v = M * points[pIndex].coord
v.size = 3
if pIndex in points_errors.keys():
pIndex = int(pIndex)
points_errors[pIndex].x += dist
points_errors[pIndex].y += 1
else:
points_errors[pIndex] = PS.Vector([dist, 1])
for pIndex in range(npoints):
if not points[pIndex].valid:
continue
if chunk.crs:
w = M * points[pIndex].coord
w.size = 3
X, Y, Z = chunk.crs.project(w)
else:
X, Y, Z, w = M * points[pIndex].coord
error = math.sqrt(points_errors[pIndex].x /
points_errors[pIndex].y)
file.write("{:6d}\t{:.6f}\t{:.6f}\t"
"{:.6f}\t{:.6f}\n".format(pIndex, X, Y, Z, error))
t1 = time.time()
file.close()
print("Script finished in " + str(int(t1-t0)) + " seconds.")
try this workflow https://de.mathworks.com/help/vision/ug/single-camera-calibrator-app.html … i use it for every lens …concerning color … there is a bunch of LUTs out there for any camera you own.
regards nj
You should multiply it by the GSD.
import Metashape as PS
import math
doc = PS.app.document
chunk = doc.chunk
def getGSD(chunk):
"""
Based on:
https://www.agisoft.com/forum/index.php?topic=1881.msg10014#msg10014
"""
# Get average camera altitude from estimated positions (not source)
N_cameras = 0
CamAlt = 0
for idx, camera in enumerate(chunk.cameras):
if camera.transform:
N_cameras += 1
mmulp = chunk.transform.matrix.mulp
est = chunk.crs.project(mmulp(camera.center))
CamAlt += est.z
CamAlt /= N_cameras
# Get focal length in mm (imported from EXIF?)
FocalLength = chunk.sensors[0].focal_length
# Get focal length in pix values after calib (current value)
FocalPix = chunk.sensors[0].calibration.f
# Get the average ground altitude baed on every point on point cloud
GroundAlt = 0
for point in chunk.point_cloud.points:
pPos = chunk.crs.project(chunk.transform.matrix.mulp(point.coord[:3]))
GroundAlt += pPos.z
GroundAlt = GroundAlt / len(chunk.point_cloud.points)
# Get the average flight height above ground
FlightHeight = CamAlt - GroundAlt
# Get the average image width and height
ImageW = 0
ImageH = 0
for idx, camera in enumerate(chunk.cameras):
ImageW += chunk.cameras[idx].photo.image().width
ImageH += chunk.cameras[idx].photo.image().height
ImageW /= len(chunk.cameras)
ImageH /= len(chunk.cameras)
# Default sensor size (mm) for Phantom 4
# CMOS Sensor Sony EXMOR RS IMX377CQT Type 1/2.3" 12.35MP (4:3 or 1.33)
# FOV = 94.0
# Pixel size: 1.55 x 1.55 um
# Two recording modes:
# 1/2.3" = 12.32MP
# 1/2.5" = 9.05MP
# Total number of pixels 4152 x 3062 = 12.71MP
# Effective pixels (1/2.3") 4056 x 3046 = 12.35MP
# Active pixels (1/2.3") 4024 x 3036 = 12.22MP
# Sensor width = 6.24 mm - 4024 pixels
# Sensor height = 4.71 mm - 3036 pixels
# Diagonal size = 7.81 mm - 5040.82 pixels
# Recommended recording (1/2.3") 4000 x 3000 = 12.00MP (4:3 ratio)
# Pixel size in micrometers
PixelSize = 1.55 # square pixel (H = V)
# Calculated sensor dimensions in milimeters
SensorW = ImageW * PixelSize / 1000 # 4000 pixels
SensorH = ImageH * PixelSize / 1000 # 3000 pixels
SensorD = math.sqrt(SensorW ** 2 + SensorH ** 2) # 5000 pixels
# Calculate average GSD (cm/pix)
avgGSD = (SensorW * FlightHeight * 100) / (FocalLength * ImageW)
# Calculate average GSD (from focal lenght in pixels after sensor calib)
avgGSD2 = FlightHeight * 100 / FocalPix
# Calculate footprint width and height
FootprintW = (ImageW * avgGSD) / 100
FootprintH = (ImageH * avgGSD) / 100
print(55 * "-")
print(" GSD Report")
print(55 * "-")
print("Sensor Dimensions (mm): {} x {}, {}". format(SensorW, SensorH, SensorD))
print("Focal Length (mm): {}". format(FocalLength))
print("Focal Length - adjusted (pix): {:.2f}". format(FocalPix))
print("Average Camera Altitude (m): {:.2f}". format(CamAlt))
print("Average Ground Altitude (m): {:.2f}". format(GroundAlt))
print("Average Flight Height (m): {:.2f}". format(FlightHeight))
print("Image Dimensions (pix): {:.0f}x{:.0f}". format(ImageW, ImageH))
print("Ground Sampling Distance (cm/pix): {:.2f}". format(avgGSD))
print("Ground Sampling Distance - adjusted (cm/pix): {:.2f}". format(avgGSD2))
print("Footprint (m): {:.2f}x{:.2f}". format(FootprintW, FootprintH))
print(55 * "-")
getGSD(chunk)
def RMS_MAX_reprojection_error(chunk):
cameras = chunk.cameras
point_cloud = chunk.point_cloud
points = point_cloud.points
projections_per_camera = point_cloud.projections
tracks = point_cloud.tracks
point_squared_errors = [[] for i in range(len(points))]
point_key_point_size = [[] for i in range(len(points))]
track_cameras = [[] for i in range(len(tracks))]
track_projections = [[] for i in range(len(tracks))]
for camera_id, camera in enumerate(cameras):
if camera not in projections_per_camera:
continue
projections = projections_per_camera[camera]
for projection_id, projection in enumerate(projections):
track_id = projection.track_id
track_cameras[track_id].append(camera_id)
track_projections[track_id].append(projection_id)
for i, point in enumerate(points):
if point.valid is False: # se válido faz a estatística abaixo
continue
track_id = point.track_id
for idx in range(len(track_cameras[track_id])):
camera_id = track_cameras[track_id][idx]
projection_id = track_projections[track_id][idx]
camera = cameras[camera_id]
projections = projections_per_camera[camera]
projection = projections[projection_id]
key_point_size = projection.size
error = camera.error(point.coord, projection.coord) / key_point_size
point_squared_errors[i].append(error.norm() ** 2)
total_squared_error = sum([sum(el) for el in point_squared_errors])
# nº projeções
total_errors = sum([len(el) for el in point_squared_errors])
max_squared_error = max([max(el+[0])
for i, el in enumerate(point_squared_errors)])
rms_reprojection_error = math.sqrt(total_squared_error/total_errors)
max_reprojection_error = math.sqrt(max_squared_error)
return rms_reprojection_error, \
max_reprojection_error
Hello Luiz,
You can try to use this script sample:Code: [Select]file = open(path, "wt")
chunk = Metashape.app.document.chunk
photos_total = len(chunk.cameras) #number of photos in chunk
markers_total = len(chunk.markers) #number of markers in chunk
if (markers_total):
file.write(chunk.label + "\n")
for i in range (0, markers_total): #processing every marker in chunk
cur_marker = chunk.markers[i]
cur_projections = cur_marker.projections #list of marker projections
marker_name = cur_marker.label
for camera in cur_marker.projections.keys(): #
#for j in range (0, len(photos_list)): #processing every projection of the current marker
x, y = cur_projections[camera].coord #x and y coordinates of the current projection in pixels
label = camera.label
file.write(marker_name + "," + label + "," + "{:.4f}".format(x) + "," + "{:.4f}".format(y) + "\n") #writing output
file.write("#\n")
file.close()