Converting coordinates

[arjunkg96]

Hello,

Given a 2D pixel coordinate on the orthomosaic, I want to do the following tasks:
- Convert it to world 2D coordinate
- project on DEM/model to gather Z coordinate
- reproject 3D point to source image
- get reprojection pixel coordinate on the source image (2D coordinates)

Could anyone please tell me how I can do these things with the Python API? My images are NOT geo referenced. I am doing this because I need a Python method to get the corresponding pixel coordinate on the source image given the coordinate on the ortho. I believe the above mentioned workflow is the way to do it.

Any help will be greatly appreciated.

Thank you.

[Paulo]

Hi arjunkg96,

Given a point p with (x, y) pixel coordinates on orthomosaic, then following code will print out its camera pixel coordinates (u, v) for each camera where it appears:

Code: [Select]

chunk = Metashape.app.document.chunk
ortho = chunk.orthomosaic
T = chunk.transform.matrix
x, y = (587.2549617344768, 702.7152715589801) # pixel coordinates of point p on orthomosaic
X, Y  = (ortho.left + ortho.resolution * x, ortho.top - ortho.resolution * y) # point p 2D coordimates  in ortho CS
if chunk.elevation:   # case DEM
    dem =  chunk.elevation
    Z = dem.altitude(Metashape.Vector((X,Y))) # Altitude in dem.crs (supposing dem.crs  = ortho.crs)
    # X, Y, Z  point p 3D coordimates  in ortho CS
    if ortho.crs.name[0:17] == 'Local Coordinates':
    p = T.inv().mulp(ortho.projection.matrix.inv().mulp(Metashape.Vector((X,Y,Z)))) # point p in internal coordinate system for case of Local CS
    else:
    p = T.inv().mulp(ortho.crs.unproject(Metashape.Vector((X,Y,Z)))) # point p in internal coordinate system
elif chunk.model:  # case 3D model
    surface = chunk.model
    if ortho.crs.name[0:17] == 'Local Coordinates':
    p1 = T.inv().mulp(ortho.projection.matrix.inv().mulp(Metashape.Vector((X,Y,0)))) # point p1 in internal coordinate system
    p2 = T.inv().mulp(ortho.projection.matrix.inv().mulp(Metashape.Vector((X,Y,10)))) # point p2 in internal coordinate system
    else:
    p1 = T.inv().mulp(ortho.crs.unproject(Metashape.Vector((X,Y,0)))) # point p1 in internal coordinate system
    p2 = T.inv().mulp(ortho.crs.unproject(Metashape.Vector((X,Y,10)))) # point p2 in internal coordinate system
    p = surface.pickPoint(p1,p2) # point p in internal coordinate system check order of p1, p2 in case of multi folded surface
elif chunk.dense_cloud:  # case dense cloud
    surface = chunk.dense_cloud
    if ortho.crs.name[0:17] == 'Local Coordinates':
    p1 = T.inv().mulp(ortho.projection.matrix.inv().mulp(Metashape.Vector((X,Y,0)))) # point p1 in internal coordinate system
    p2 = T.inv().mulp(ortho.projection.matrix.inv().mulp(Metashape.Vector((X,Y,10)))) # point p2 in internal coordinate system
    else:
    p1 = T.inv().mulp(ortho.crs.unproject(Metashape.Vector((X,Y,0)))) # point p1 in internal coordinate system
    p2 = T.inv().mulp(ortho.crs.unproject(Metashape.Vector((X,Y,10)))) # point p2 in internal coordinate system
    p = surface.pickPoint(p1,p2) # point p in internal coordinate system check order of p1, p2 in case of multi folded surface
else:
    raise Exception("No DEM, 3D model or Dense Cloud")
print("Camera pixel coordinates of point p with following ortho pixel coordinates: (",x,",",y,")")
for camera in chunk.cameras:
    if not camera.project(p):
    continue
    u = camera.project(p).x   # u pixel coordinates in camera
    v = camera.project(p).y # v pixel coordinates in camera
    if (u < 0 or u > camera.sensor.width or v < 0 or v > camera.sensor.height):
    continue
    print(" Camera: ",camera.label,"u: ",u,"\tv: ",v)
as seen in 1st screen capture where Ortho is in projected CS (WGS84 / UTM14N)...
and in 2nd capture where Ortho was done in Local CS (Front XZ) from a project with spherical images taken  inside a church...
Hope this can be helpful,
PS. I modified the above code to include case where there is model or dense_cloud surface but no DEM. In these cases, the Z component of point 2D (X,Y) ortho CS coordinates is found using surface.pickPoint along vertical line from p1 (Z=0) to p2 (Z=10). However the result depends on order of (p1,p2) start-end pick points and for some complex surfaces, it can give erroneous results.
As seen in 3rd attachment, pickPoint from p1 to p2 (shown by green line) will result in finding a point behind Christ's INRI tablet while the point p from ortho was clearly in front....

[arjunkg96]

This is what I was looking for. Thank you Paulo, thank you very much!