# Forum

### Author Topic: Reprojection error  (Read 4402 times)

#### LFSantosgeo

• Jr. Member
• Posts: 70
##### Reprojection error
« on: March 10, 2019, 01:46:14 PM »
Hello,

I got the following code to calculate max and rms reprojection errors:

Code: [Select]

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

Sometimes it doesn't work at all at getting the values and sometimes it does. I haven't manage to figure out why. Can anyone help out?
Luiz Fernando

• Newbie
• Posts: 8
##### Re: Reprojection error
« Reply #1 on: March 10, 2019, 02:04:07 PM »
Try this...

Code: [Select]
doc = Metashape.app.document

# Compute the re-projection error
chunk = doc.chunk
point_cloud = chunk.point_cloud
points = point_cloud.points
error, tie_points = [], []

for camera in [cam for cam in doc.chunk.cameras if cam.transform]:
point_index = 0
photo_num = 0
for proj in doc.chunk.point_cloud.projections[camera]:
track_id = proj.track_id
while point_index < len(points) and points[point_index].track_id < track_id:
point_index += 1
if point_index < len(points) and points[point_index].track_id == track_id:
if not points[point_index].valid:
continue

dist = camera.error(points[point_index].coord, proj.coord).norm() ** 2
error.append(dist)
photo_num += 1

tie_points.append(photo_num)

reprojection_rmse = round(math.sqrt(sum(error) / len(error)), 2)
reprojection_max = round(max(error) , 2)
reprojection_std = round(statistics.stdev(error), 2)
tie_points_per_image = round(sum(tie_points) / len(tie_points), 0)

print("Average tie point residual error: " + str(reprojection_rmse))
print("Maxtie point residual error: " + str(reprojection_max))
print("Standard deviation for tie point residual error: " + str(reprojection_std))
print("Average number of tie points per image: " + str(tie_points_per_image))

#### LFSantosgeo

• Jr. Member
• Posts: 70
##### Re: Reprojection error
« Reply #2 on: March 10, 2019, 02:25:08 PM »

You script works perfectly for pixel values for RMSE and Max reprojection errors, despite the max value is a bit off by 0.03 as shown on the "show info...". You should use math.sqrt(max(error)), then the values are exactly the same.

But I'm trying to get the non pixel values for the reprojection errors.
Luiz Fernando

• Newbie
• Posts: 8
##### Re: Reprojection error
« Reply #3 on: March 10, 2019, 04:03:02 PM »
You should multiply it by the GSD.

#### LFSantosgeo

• Jr. Member
• Posts: 70
##### Re: Reprojection error
« Reply #4 on: March 11, 2019, 05:47:34 AM »
You should multiply it by the GSD.

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

From above script the calculated GSD is different from what shows in PDF report. Calculated GSD is a higher value. I guess it is because my average flight height (depending on average cameras and points altitudes -- coordinates z) is greater than what it actually is. Any insight on calculating GSD?
Luiz Fernando