Finding 3D (X, Y, Z) world coordinates from multiple 2D image points

[neavrys]

Hi,

Through the Metashape API 1.6.2, I am actively looking to develop a script that will automatically process a list of 2D coordinates on all cameras/image in a list of chunks. The output result will be the list of 3D point positions in the project cartographic coordinate, mostly UTM.

The 2D image coordinate list would be constant and contain 14336 (X,Y) pairs. Yes, that a lot.
This script should be executed after of all procedures (Camera alignment, optimization, Camera alignment, dense cloud, mesh, DEM and ortho).

I'm having trouble finding my way around. There are many point projection functions in the Metashape API.

Below is the beginning of a first draft script:

Code: [Select]

import Metashape
import numpy as np

#Load coord list from external numpy file:
Img2DCoordList = np.load("***.npy")

chunk2Process = ["chunk name #1", "chunk name #2"]
chunks = Metashape.app.document.chunks

if not len(chunks):
    Metashape.app.messageBox("No chunk!")
    exit()

Metashape.app.messageBox("-Start-")
for chunk in chunks:
    if chunk.label in chunk2Process:
        Metashape.app.messageBox("Process with Chunk :{0}".format(chunk.label))
        cameras = chunk.cameras
        if not len(cameras):
            Metashape.app.messageBox("No camera in this chunk called {0}!".format(chunk.label))
            continue
        for cam in cameras:
            pnt = Img2DCoordList[0]  # 1 point for testing
            ### exploration ###
            # 1.
            # Parameter : Coordinates of the point to be projected (as Vector);
            # Output : Returns a 2D coordinates, as vector, of the point projection on the photo;
            cam.project(pnt)
            # 2.
            # Parameter : Projection coordinates (as Vector);
            # Output : returns a 3D coordinates, as vector, which will have specified projected coordinates;
            cam.unproject(pnt)
            # 3.
            # 3.1 Create a marker from pnt
            # 3.2 Project the marker from a way or another
            # 3.3 Retrieve the 3D coordinates
            # [code]
            # 4.
            # Any ideas?


Many thanks,

[neavrys]

From the GUI, creating a marker works very well. It provides the necessary information which is easy and simple to export in several formats.

Is there a way to automate the creation of markers from images coordinates?

Thanks,

[Alexey Pasumansky]

Hello neavrys,

Have you checked the solution related to the marker placement that I have recently published in the similar thread:
https://www.agisoft.com/forum/index.php?topic=12160.msg54374#msg54374

Does it solve your task?

[neavrys]

Hi Alexey,

Yes it does! I pass a huge 2D coords list related to each camera from a dict and it's returning a 3D coords list that I save into a dict.

Now the problem is to reduce the processing time. The iteration over the list is unfortunately an option for this task. With the entire code bellow, I process 50 marker in 3.3 seconds and 500 in around 35 seconds. So, my estimation for ­~43000 pixels coords over ~5500 camera will take about ~163 days!!! 

Please, tell me that I can either do vectorization, broadcasting, multiprocessing, process on GPU, etc. to reduce the processing time? A way to avoid this for loop which is the main performance bottleneck?

Main performance bottleneck:

Code: [Select]

...
coord3DList = []
for coord in Img2DCoordList[0:50]:  # 5 point only for testing
    pnt = chunk.point_cloud.pickPoint(cam.center, cam.transform.mulp(cam.sensor.calibration.unproject(vec(coord))))
    mk = chunk.addMarker(pnt)
    #print("The point {0} for the {1} is : {2}".format(coord, cam.label, mk.reference.location))
    coord3DList.append(np.asarray(mk.reference.location))
...



All code:

Code: [Select]

import Metashape as MS
import numpy as np
from Metashape import Vector as vec
import pickle
import time

start = time.time()
sourcePath = r"G:\..."
sourceName = r"....p"
outputPath = r"G:\Projection\..."
outputCRS = "EPSG::4326"

chunk2Process = ["Merged Chunk_Done"]

imgCoordProject = {}

## import Coord from dictio
with open(r'{}\{}'.format(sourcePath,sourceName), 'rb') as fp:
    img2DCoordDict = pickle.load(fp)

MS.app.messageBox("-Start-")
start2 = time.time()
#In case of running from outside


#if __name__ == '__main__':
doc = MS.app.document
chunks = doc.chunks

if not len(chunks):
    MS.app.messageBox("No chunk!")
    exit()

for chunk in chunks:  # chunk = chunks[-1]
    if chunk.label in chunk2Process:
        cameras = chunk.cameras
        processed = 0
        if not len(cameras):
            MS.app.messageBox("No camera in this chunk called {0}!".format(chunk.label))
            continue
        chunk.remove(chunk.markers)  # Clear all markers
        for cam in cameras[0:2]:  # 2 camera only for testing
            Img2DCoordList = img2DCoordDict[cam.label]
            coord3DList = []
            for coord in Img2DCoordList[0:50]:  # 50 point only for testing
                pnt = chunk.point_cloud.pickPoint(cam.center, cam.transform.mulp(cam.sensor.calibration.unproject(vec(coord))))
                mk = chunk.addMarker(pnt)
                print("The point {0} for the {1} is : {2}".format(coord, cam.label, mk.reference.location))
                coord3DList.append(np.asarray(mk.reference.location))
                #or
                #pnt2 = chunk.model.pickPoint(cam.center, cam.transform.mulp(cam.sensor.calibration.unproject(vec(coord))))
                #mk2 = chunk.addMarker(pnt2)
                #mk2.reference.location
                #coord3DList.append(np.asarray(mk2.reference.location))

            imgCoordProject[cam.label] = np.asarray(coord3DList)
            # Save markers
            #chunk.exportMarkers("{0}\{1}.xml".format(outputPath, cam.label), "xml", outputCRS)
            #chunk.exportMarkers("{0}\{1}.shp".format(outputPath, cam.label), "shp", outputCRS)
            # Remove all markers
            chunk.remove(chunk.markers)

with open(outputPath + r'\imgCoordProject.p', 'wb') as fp:
    pickle.dump(imgCoordProject, fp, protocol=pickle.HIGHEST_PROTOCOL)

end = time.time()
print("Process Time after packages import: {:.2f} second".format(end - start))
print("Process Time after datasource import: {:.2f} second".format(end - start2))


Updated:
The addition of a cameras filter restrict the process on image that cover a specific area. This reduce significantly the amount of image to process to ~1900 images.
Process time estimation would still be around  ~62 days!!! 

Also, I look into a Multiprocessing method similar to the footprints_to_shapes.py script in "agisoft-llc/metashape-scripts". Not sure how to handle the output, but this should reduce the processing time by a factor of 28 (CPU physical core).

Code: [Select]

with concurrent.futures.ThreadPoolExecutor(multiprocessing.cpu_count()) as executor:
    executor.map(lambda camera: process_camera(chunk, camera), chunk.cameras)

Any idea how implement this?

[neavrys]

The implementation of multiprocessing in the previous code cause opening of multiple empty Metashape instances, command console say "invalid parameter -e" or something similar and the process wait indefinitely for multiprocess output return which will never be produce.

The multiprocessing implementation:

Code: [Select]

def project_point(data):
    chunk, cam, coord = data[0], data[1], data[2]
    mk = chunk.addMarker(chunk.point_cloud.pickPoint(cam.center, cam.transform.mulp(cam.sensor.calibration.unproject(vec(coord)))))
    return mk.reference.location
...

if __name__ == '__main__':
...
pool = multiprocessing.Pool(processes=multiprocessing.cpu_count())
result = pool.map(project_point, np.concatenate((np.full((Img2DCoords.shape[0], 1), chunk), np.full((Img2DCoords.shape[0], 1), cam), Img2DCoords), axis=-1))
...

Questions:

Is a multiprocessing implementation work within a Metashape console?

Is that will work from a Python IDE console outside Metashape?

Is only multithreading is "compatible" with Metashape?

Many thanks,