Reprojection error

[LFSantosgeo]

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?

[rongilad]

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]

Thank you for the reply rongilad!

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.

[rongilad]

You should multiply it by the GSD.

[LFSantosgeo]

You should multiply it by the GSD.

Thinking of what rongilad wrote:

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?