Forum

Author Topic: Reprojection error  (Read 6210 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.

`import Metashape as PSimport mathdoc = PS.app.documentchunk = doc.chunkdef 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)`