Projecting from image coordinates to the orthomosaic

[vivekhsridhar]

Hi all,

I'm a biologist that studies animal movement and behaviour. As part of my work, I acquired UAV images which I used in Agisoft Metashape to reconstruct a georeferenced map of my animals' habitat. As a next step, I have videos I collected within this habitat. I have used computer vision algorithms to track my animals within my video, thus giving me their location in pixel coordinates. I would now like to convert trajectories of these animals from pixel space to global coordinates. To test if this worked, I first marked specific points on my orthomosaic. I then converted these coordinates to image space (from 3D to 2D) using the following code:

Code: [Select]

import numpy as np
import pandas as pd
import Metashape

def list_coordinates(file_path):
    data = []
    with open(file_path, 'r') as file:
        for line in file:
            values = line.strip().split(',')[1:-1]
            # Filter out empty strings before converting to float
            values = [float(val) if val else 0.0 for val in values]
            data.append(tuple(values))
    return data

def process_chunk(chunk, global_coordinates):
    T = chunk.transform.matrix

    data_list = []

    for point_idx, global_coord in enumerate(global_coordinates):
        p = T.inv().mulp(chunk.crs.unproject(Metashape.Vector(global_coord)))
        print(f"Image pixel coordinates of point {point_idx + 1} with global coordinates ({chunk.crs.name}): {global_coord}")

        for i, camera in enumerate(chunk.cameras):
            project_point = camera.project(p)

            if project_point:
                u = project_point.x  # u pixel coordinates in camera
                v = project_point.y  # v pixel coordinates in camera

                if 0 <= u <= camera.sensor.width and 0 <= v <= camera.sensor.height:
                    # Extract video and frame_seq from the camera label
                    camera_label_parts = camera.label.split('_')
                    video = '_'.join(camera_label_parts[:-1])
                    frame_seq = camera_label_parts[-1]

                    data_list.append([point_idx + 1, camera.label, video, frame_seq, u, v])

    columns = ['Point', 'Camera', 'video', 'frame_seq', 'u', 'v']
    df = pd.DataFrame(data_list, columns=columns)
    return df

# Define global coordinates for multiple points
points = list_coordinates('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/Solitude/output/points_xyz.txt')

# Use active Metashape document
doc = Metashape.app.document

# Iterate through all chunks in the document
df = pd.DataFrame()
for chunk in doc.chunks:
    if chunk.label != 'BaseMap':
        # Set the current chunk to the one being processed
        doc.chunk = chunk

        # Process the current chunk
        tmp = process_chunk(chunk, points)
        df = pd.concat([df, tmp], ignore_index=True)

# Save the results to a CSV file
df.to_csv('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/Solitude/output/points_uv.csv', index=False)

I know this code works because I plotted the 2D coordinates onto the respective images and they appear at the right location on the image. So next, I decided to convert these 2D coordinates back to 3D and check if I obtained coordinates of my original points. I used the following code to do this:

Code: [Select]

import Metashape
import datetime


doc = Metashape.app.document
chunk = doc.chunk
surface = chunk.point_cloud

chunk = doc.chunks[1]
print(chunk.label)

camera = chunk.cameras[0]
print(camera)

x=337.67729912777475
y=903.5153164982171

coords_2D = Metashape.Vector([x, y])
point_internal = surface.pickPoint(camera.center, camera.unproject(coords_2D))
point3D_world = chunk.crs.project(chunk.transform.matrix.mulp(point_internal))

print(point_internal)
print(point3D_world)

This however does not work and the newly obtained 3D coordinates do not match the location of the original point on the orthomosaic. Where have I gone wrong?

Thanks in advance for your help! I really appreciate any pointers in this direction.

Cheers!
Vivek

[Paulo]

Hello,

what occurs to me is to check if chunk.crs in both codes are the same?

[vivekhsridhar]

@Paul, thanks for the quick response. I just checked and the chunk.crs in both scripts are the same (CoordinateSystem 'WGS 84 (EPSG::4326)'). Let me know if you spot some other bug in my code / error in my thought process.

Cheers!
Vivek

[Paulo]

Hello vivekhsridar,

the code seems good to me... I would check in your first code the global coordinates how are they extracted? from orthomosaic? Then how did you get the height or Z? does the height reference correspond to point_cloud?

Many questions to investigate...

[vivekhsridhar]

The initial global coordinates were extracted from the orthomosaic. Here, I use the 'draw point' function to click on specific features on the orthomosaic. This created a shape layer with my marked points. I exported this shape layer as a .txt (file attached). These are the X,Y and Z coordinates that I start with. The subsequent processing is as described before. I use my first script to convert these X,Y,Z coordinates to U,V image coordinates and then try to obtain the X,Y,Z back by unprojecting the U,V coordinates.

[Alexey Pasumansky]

Hello Vivek,

What are the coordinate systems of shapes and orthomosaic (according to the Chunk Info dialog) and what systems are set for chunk, markers and cameras in the Reference pane settings dialog?

[vivekhsridhar]

All coordinate systems seem to match.

I wrote these two lines of code to output crs of the orthomosaic and the shapes

Code: [Select]

import Metashape

doc = Metashape.app.document
chunk = doc.chunk

print('crs of orthomosiac is: ', chunk.orthomosaic.crs)
print('crs of shapes is: ', chunk.shapes.crs)

This is the output on the Metashape console
2024-01-18 10:53:19 crs of orthomosiac is:  <CoordinateSystem 'WGS 84 (EPSG::4326)'>
2024-01-18 10:53:19 crs of shapes is:  <CoordinateSystem 'WGS 84 (EPSG::4326)'>

As suggested by you, I also checked this on the chunk info (images attached).

I also checked the settings on the reference pane for the crs of the photos and the markers and they all also seem to be WGS 84 (EPSG::4326) (Screenshot attached).

[Paulo]

Hello again,

wht version are u using of Metashape?

If version is 1.8 or less than surface chunk.point_cloud corresponds to sparse point cloud or tie points and a pickPoint on such surface would not be accurate compared with pickPoint on dense cloud surface... and could explain the discrepancy...

[vivekhsridhar]

Hi Paul, thanks for taking the time to try and sort this out. I'm using version 2.0.4. So this should still be the dense cloud then?

[Paulo]

Yes it is the dense cloud... for your orthomosaic what surface did you use? If DEM what type (terrain or surface mopdel)...?

[vivekhsridhar]

I used a DEM surface model. I haven't classified my ground points

[vivekhsridhar]

I might have found the issue. In the code where I convert the U,V image coordinates back to the X,Y,Z orthomosaic, the unprojection seems to work if I apply it to photos within the same chunk as the orthomosaic. However, this does not seem to work when I do the unprojection on images on a separate chunk. For example, in the code below, everything seems to work once I silence the line "chunk = doc.chunks[1]" (as done below).

Code: [Select]

import Metashape
import datetime


doc = Metashape.app.document
chunk = doc.chunk
surface = chunk.point_cloud

# chunk = doc.chunks[1]
print(chunk.label)

camera = chunk.cameras[0]
print(camera)

x=2066
y=544

coords_2D = Metashape.Vector([x, y])
point_internal = surface.pickPoint(camera.center, camera.unproject(coords_2D))
point3D_world = chunk.crs.project(chunk.transform.matrix.mulp(point_internal))

print(point_internal)
print(point3D_world)

I checked and the other chunk (which I'm trying to project onto the orthomosaic) is aligned to the base chunk. However, the photos have aren't calibrated (they have an NC next to them). Could this be the issue? Here is the code I use to align images on the new chunk with my 'BaseMap'.

Code: [Select]

import Metashape
import datetime
import glob
import os
import pathlib
import numpy as np
import helper_functions as hf

doc = Metashape.app.document
reference_chunk = doc.chunk.key

# The directory where videos are saved
VideoDirectory = "/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/Solitude/images/"

# The directory where frames must be imported
ImportDirectory = "/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/Solitude/images/"

def list_videos(directory):
    directory_path = pathlib.Path(directory)
    video_extensions = ['.mp4', '.avi', '.mkv', '.mov', '.wmv']
    video_files = [(str(file.resolve()), file.stem) for file in directory_path.iterdir() if file.suffix.lower() in video_extensions]
    return video_files
def import_frames(video_path, video_name, import_directory):
    # Create a folder for each video
    video_folder = os.path.join(import_directory, video_name)
    os.makedirs(video_folder, exist_ok=True)

    # Add a new chunk for each video
    chunk = doc.addChunk()
    chunk.label = video_name
    doc.chunks.append(chunk)

    # Define the image names pattern for frames
    image_names = os.path.join(video_folder, f'{video_name}_{{filenum:04}}.png')

    # Import frames with custom frame step
    chunk.importVideo(video_path, image_names, frame_step=Metashape.FrameStep.CustomFrameStep, custom_frame_step=20)


hf.log( "--- Starting workflow ---" )
hf.log( "Metashape version " + Metashape.Application().version )
hf.log_time()

# List videos in VideoDirectory
videos = list_videos(VideoDirectory)

# Import frames from listed videos
for video_path, video_name in videos:
    import_frames(video_path, video_name, ImportDirectory)


# Parameters for feature matching photos
match_photos_config = {
    'downscale': 1,
    'generic_preselection': True,
    'reference_preselection': True,
    'reference_preselection_mode': Metashape.ReferencePreselectionEstimated
}

chunk_dirs = hf.get_subdirectories(ImportDirectory, ignore_files=['.DS_Store'])

all_chunks = [reference_chunk]
# Match photos and align cameras
for chunk in doc.chunks:
    if chunk.label != 'BaseMap':
        hf.log( "Processing chunk" )

        chunk.matchPhotos(**match_photos_config)
        chunk.alignCameras()
        all_chunks.append(chunk.key)

doc.alignChunks(chunks=all_chunks, reference=reference_chunk)

hf.log_time()
hf.log( "--- Finished workflow ---")

doc.save('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/Solitude/Solitude.psx')

How would I go about fixing this? I've also attached a screenshot of my workspace if that helps.

Cheers!
Vivek

[vivekhsridhar]

Hi again,

I was wondering if someone could confirm that the lack of image calibration could in fact result in this issue? And if so, is there a way in metashape to import frames from a video while also keeping the intrinsic parameters for calibration? Alternatively, since it is the same camera / drone that was used to create the orthomosaic, would it be a possibility to use the same calibration parameters?

Thanks again for helping me with this.

Vivek

[Alexey Pasumansky]

Hello Vivek,

camera.project and camera.unproject methods are using calibration information, so you need to load the available calibration data (from another project, presumably) for the sensor corresponding to the cameras to be used.

Something like the following

Code: [Select]

sensor = camera.sensor
calib = Metashape.Calibration()
calib.load(path)
sensor.calibration = calib


[vivekhsridhar]

Hi all,

I'm still struggling with this. I'd really appreciate if someone could help me solve this problem and so, I'm recapping what I've done so far here.

I'm a biologist that studies animal movement and behaviour. As part of my work, I acquired UAV images which I used in Agisoft Metashape to reconstruct a georeferenced map (an orthomosaic) of my animals' habitat. As a next step, I have videos I collected within this habitat. I have used computer vision algorithms to track my animals within my video, thus giving me their location in pixel coordinates. I would now like to convert trajectories of these animals from pixel space to global coordinates. To test if this would work, I first marked specific points on my orthomosaic. I then converted these coordinates to image space (from 3D to 2D) using the following code (note the images are part of a separate chunk which I've aligned separately. I've also georeferenced this chunk using GCPs and the tiepoints align beautifully to the original chunk where I created the orthomosaic):

Code: [Select]

import Metashape
import numpy as np
import pandas as pd

def list_coordinates(file_path):
    data = []
    with open(file_path, 'r') as file:
        for line in file:
            values = line.strip().split(',')[1:-1]
            # Filter out empty strings before converting to float
            values = [float(val) if val else 0.0 for val in values]
            data.append(tuple(values))
    return data

def process_chunk(chunk, global_coordinates):
    T = chunk.transform.matrix

    data_list = []

    for point_idx, global_coord in enumerate(global_coordinates):
        p = T.inv().mulp(chunk.crs.unproject(Metashape.Vector(global_coord)))
        print(f"Image pixel coordinates of point {point_idx + 1} with global coordinates ({chunk.crs.name}): {global_coord}")

        for i, camera in enumerate(chunk.cameras):
            project_point = camera.project(p)

            if project_point:
                u = project_point.x  # u pixel coordinates in camera
                v = project_point.y  # v pixel coordinates in camera

                if 0 <= u <= camera.sensor.width and 0 <= v <= camera.sensor.height:
                    # Extract video and frame_seq from the camera label
                    camera_label_parts = camera.label.split('_')
                    video = '_'.join(camera_label_parts[:-1])
                    frame_seq = camera_label_parts[-1]

                    data_list.append([point_idx + 1, camera.label, video, frame_seq, u, v])

    columns = ['Point', 'Camera', 'video', 'frame_seq', 'u', 'v']
    df = pd.DataFrame(data_list, columns=columns)
    return df

# Define global coordinates for multiple points
points = list_coordinates('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/TalChhapar/output/ninety_territory_points_xyz.txt')

# Use active Metashape document
doc = Metashape.app.document

# Iterate through all chunks in the document
df = pd.DataFrame()
for chunk in doc.chunks:
    if chunk.label != 'Chunk 1' and chunk.label != 'Chunk 2':
        # Set the current chunk to the one being processed
        doc.chunk = chunk

        # Process the current chunk
        tmp = process_chunk(chunk, points)
        df = pd.concat([df, tmp], ignore_index=True)

# Save the results to a CSV file
df.to_csv('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/TalChhapar/output/p1_territory_points_uv.csv', index=False)

I know this code works because I plotted the 2D coordinates onto the respective images and they appear at the right location on the image (attached sample image). So next, I decided to convert these 2D coordinates back to 3D and check if I obtained coordinates of my original points. I used the following code to do this:

Code: [Select]

import Metashape
import datetime
import pandas as pd

doc = Metashape.app.document
chunk = doc.chunks[0]
surface = chunk.point_cloud

chunk = doc.chunks[2]
print(chunk.label)

camera = chunk.cameras[0]
print(camera)

df = pd.read_csv('/Users/vivekhsridhar/Library/Mobile Documents/com~apple~CloudDocs/Documents/Metashape/TalChhapar/output/p1_territory_points_uv.csv')
df = df[df['Camera'] == '20230310_SE_Lek1_P1D1_DJI_0190_0000']

for index, row in df.iterrows():
x = row['u']
y = row['v']

coords_2D = Metashape.Vector([x, y])
point_internal = surface.pickPoint(camera.center, camera.unproject(coords_2D))
point3D_world = chunk.crs.project(chunk.transform.matrix.mulp(point_internal))

print(point_internal)
print(point3D_world)

This however does not work and the newly obtained 3D coordinates do not match the location of the original point on the orthomosaic. Where have I gone wrong? Unlike my previous post, all my chunks now have an [R] flag next to them so all images are referenced as well. I'd appreciate any help in terms of how I can proceed with this.

Best,
Vivek