Reprojection Errors (script x chunk info)

[LFSantosgeo]

Hello!

I'm running the following script to get reprojection errors by point of the point cloud:

Code: [Select]

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.")
While it succeeds to calculate the reprojection error for each tie point, when I compare the max error value obtained with script (14,5876 pix) with the max error provided by chunk info (21.8388 pix) they don't match. Is that supposed to happen? Looking forward for the reply.

Regards!

[Alexey Pasumansky]

Hello Luiz,

Can you please check, if the following script gives you the correct values for max and RMS reprojection errors:

Code: [Select]

import Metashape
import time, math

def calc_reprojection(chunk):
point_cloud = chunk.point_cloud
points = point_cloud.points
npoints = len(points)
projections = chunk.point_cloud.projections
err_sum = 0
num = 0
maxe = 0

point_ids = [-1] * len(point_cloud.tracks)
point_errors = dict()
for point_id in range(0, npoints):
point_ids[points[point_id].track_id] = point_id

for camera in chunk.cameras:
if not camera.transform:
continue
for proj in projections[camera]:
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 = camera.error(point.coord, proj.coord).norm() ** 2
err_sum += error
num += 1
if point_id not in point_errors.keys():
point_errors[point_id] = [error]
else:
point_errors[point_id].append(error)
if error > maxe: maxe = error

sigma = math.sqrt(err_sum / num)
return (sigma, point_errors, maxe)

chunk = Metashape.app.document.chunk
t0 = time.time()
result = calc_reprojection(chunk)
print('RMS reprojection error (pix): ' + str(result[0]))
print("Max reprojection error (pix): " + str(math.sqrt(result[2])))
t0 = time.time() - t0
print("Script finished in " + "{:.2f}".format(float(t0)) + " seconds")

I think that the different comes from the confusion - max. error is calculated per projection and not per point.

[LFSantosgeo]

Hello Luiz,

Can you please check, if the following script gives you the correct values for max and RMS reprojection errors:

Yes the code you provided calculates precisely RMS and max reprojection for the chunk.

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?

Just to clear a doubt of mine: the gradual selection (API and GUI) for reprojection error  goes through each one of the projections and removes them by the defined threshold? How does it select points then?

Thank you

[LFSantosgeo]

So I've changed my code to get point errors...

Code: [Select]

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()
But when I average the point errors they differ a bit (0.55pix) from the chunk info (0.53pix).

[LFSantosgeo]

Anything on this? I guess my code is wrong but I've double checked it and can't figure out what is wrong with it.
I trying to get reprojection error by points to text file with: point_id, X, Y, Z, error

[Alexey Pasumansky]

Hello Luiz,

The gradual selection dialog should select the tie points that have at least one projection with the error above the selected threshold.

Meanwhile I'll try to check your code on some sample projects and see, if any difference can be observed.

[LFSantosgeo]

Thank you Alexey!

Looking forward for your test results.

[Dora]

Hi Luiz,

It looks like the rms reprojection error in chunk info is calculated from the errors of all projections, not from point errors. To get the global rms value divide the sum of all projection errors by the total number of projections.

Code: [Select]

# Get point errors
error = math.sqrt(sum(point_errors[point_id]) / len(point_errors[point_id]))

# Get rms
# at start set err_sum = 0, num = 0
err_sum += sum(point_errors[point_id])
num += len(point_errors[point_id])

rms_error = math.sqrt(err_sum / num)

From this you can also get max reprojection error from chunk info. Now this leads to another question:
I understand how to remove tie points with a maximum reprojection error above a given threshold. This makes sense for points with only 2 projections, but for those with many projections, wouldn't it make sense to remove the projections with high errors and keep the rest? There is a way to remove tie points, and to turn off cameras, but is it possible to disable some projections associated with a given camera/tie point?

Any ideas?

Cheers!

[Alexey Pasumansky]

Hello Luiz,

In order to calculate the same value, as shown in chunk info you should keep (sum) the following values per point:
total += sum(point_errors[point_id])
n_total += len(point_errors[point_id])
then use sqrt(total / n_total).
It's necessary, as you should calculate RMS for projections. So I confirm the solution provided by Dora above.