Show Posts

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.

Messages - LFSantosgeo

Pages: [1] 2 3 ... 5
Python Scripting / Re: Reprojection error
« on: March 11, 2019, 05:47:34 AM »
You should multiply it by the GSD.

Thinking of what rongilad wrote:

Code: [Select]
import Metashape as PS
import math

doc =
chunk = doc.chunk

def getGSD(chunk):
    Based on:
    # 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 =
            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 =[: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 * "-")


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?

Python Scripting / Re: Reprojection error
« on: March 10, 2019, 02:25:08 PM »
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.

Python Scripting / Reprojection error
« on: March 10, 2019, 01:46:14 PM »

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:

        projections = projections_per_camera[camera]

        for projection_id, projection in enumerate(projections):

            track_id = projection.track_id

    for i, point in enumerate(points):
        if point.valid is False:  # se válido faz a estatística abaixo

        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, \

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?

Python Scripting / Re: Markers XML to Text File
« on: February 14, 2019, 06:26:41 PM »
Hello Alexey,

It works! Thank you very much!

Hello Luiz,

You can try to use this script sample:

Code: [Select]
file = open(path, "wt")
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


Python Scripting / Markers XML to Text File
« on: February 14, 2019, 06:08:40 AM »

Is there any script to export MetaShape XML markers into text file with marker labels, photo labels and 2D coordinates from markers projections in the photos? Output something like:

#marker_label  #photo_label  #img_x  # img_y

I've run into lots of exemples of the opposite procedure...
Thanks in advance!

Thank you Alexey!

Hello Luiz,

Please check:
Code: [Select]

Hey mks_gis,

What do you think of the following code structure, for instance, to range from a value to another of filter levels and try to remove points if a condition is True? It also can restart searching through the same range values if one is found.  I'm working on still, just thought of sharing this idea to implement filtering for the USGS workflow (thinking of reconstruction uncertainty running more than once till no points are selected based on the chosen range). This of course would rapidly decrease points in point cloud.

Code: [Select]
threshold = range(0, 11, 1)  # here would be the filter levels range
listx = []  # throwing a list just for the exemple
for i in threshold:

restart = 0
restartLoop = True
while restartLoop:
    restartLoop = False
    for idx, i in enumerate(listx):
        print("do something as printing i:", i)

        if i > 5:  # if this condition: remove points and restart loop to search again in range
            print("found value for condition: ", i)
            del listx[idx]  # simulates removing points for a certain i value

            #  optimization would happen here

            restartLoop = True
            print("RESTARTING LOOP\n")
            restart += 1
            break  # break inner while and restart for loop
            # continue if the inner loop wasn't broken

print("restart - outer while", restart)

Looking forward to see your new version of the code!

Any news about this feature in Metashape?

Just found out. Python tricks...
Code: [Select]
points = point_cloud.points
Gives the total valid points.
Then Gradual selection is done.
And then:
Code: [Select]
chunk = doc.chunk
Gives the valid tiepoints number after gradual selection and for that you can assign to a new variable. You could reassign the variable points but you will loose the total points value before gradual selection.

Hey mks_gis nice coding.

    len(chunk.point_cloud.tracks) - counts valid and invalid tiepoints right?
    len(chunk.point_cloud.points) - and this counts the initial number of valid tiepoints before gradual selection

Have you managed to count tiepoints with Python after gradual selection? Every time I use len(chunk.point_cloud.points) I get the same initial value

Currently needs to be run as one session, valid ties not accessible after
    gradual selection. To be fixed in next version using tracks.
    len(chunk.point_cloud.points) - valid ties
    len(chunk.point_cloud.tracks) - all ties


How do I know by python if the rolling shutter is enabled? Or if it was done?


Try this!
Code: [Select]
Hope it helps!

Python Scripting / Re: Height Above Ground
« on: May 16, 2018, 07:33:27 PM »
Thank you Alexey and Erik! I'll give a try on your thoughts... And I'll post it here.

Python Scripting / Filtering Sparse Point Cloud
« on: May 16, 2018, 07:32:34 PM »
Hello! Is there a consensus about the gradual selection and build points to filter reprojection error of the sparse point cloud?

I've checked the following topics from which I extracted some quotes bellow:

GUI: Gradual Selection=X is not the same as buildPoints(error=X)
So they are not the same. And they are not equivalent for the filtering task. What's the difference?

For buildPoints() -- Build Points.
Code: [Select]
buildPoints(error=10[, min_image ][, progress])
Rebuild point cloud for the chunk.
• error (float) – Reprojection error threshold.
• min_image (int) – Minimum number of point projections.
• progress (Callable[[float], None]) – Progress callback.

The perk of using this [buildPoints], apart from that it's shorter, is that statistical 'outliers' that were removed in previous stages could be reinstated if their errors are reduced to within the threshold.

As for the PointCloud.Filter()  -- Gradual Selection:
Code: [Select]
threshold = 0.5
>>> f = PhotoScan.PointCloud.Filter()
>>> f.init(chunk, criterion = PhotoScan.PointCloud.Filter.ReprojectionError)
>>> f.selectPoints(threshold)

I understand that all matches ,even filtered ones, are kept within the database which means that chunk.buildPoints(error=X) {'coordinates applied...'} will 'resuscitate' some of them and is by no means equivalent to Gradual Selection:Reprojection Error

Python Scripting / Height Above Ground
« on: May 10, 2018, 11:26:59 PM »

Is there a way to access the estimated(calculated) height above ground through Python (if already available after alignment)?

Luiz Fernando

Pages: [1] 2 3 ... 5