Messages

1

General / Inversing brown-conrady equations (for camera distortion)

« on: February 25, 2021, 01:21:19 AM »

Hi,

I am working on an application where i need to project 2d image pixel coordinates into the 3D camera coordinate system. I understand that this can be done with the metashape .unproject method, for example:

Code: [Select]

chunk = Metashape.app.document.chunk
camera = chunk.cameras[0] # test using first camera
cam_calib = camera.sensor.calibration # camera 0 calibration values

test_pixel = Metashape.Vector([100, 100]) # random point in image pixel coordinates
X, Y, Z = cam_calib.unproject(test_pixel) # result using unproject function
However I would like to understand the equations behind this reverse projection. In particular i have radial and tangental distortion parameters which i need to include. I understand and can replicate the forward projection from camera crs to pixel, as below, using the equations for frame cameras at the end of the documentation https://www.agisoft.com/pdf/photoscan-pro_1_4_en.pdf. Using this i get the point (100, 100) as expected.

Code: [Select]

test_pixel = Metashape.Vector([100, 100]) # random point in image pixel coordinates
X, Y, Z = cam_calib.unproject(test_pixel) # result using unproject function

# Test projection from camera crs to test pixel...
print('Inputing ' + str(X) + ',' + str(Y) + ',' + str(Z))
print('Expecting ' + str(test_pixel))

x = X/Z
y = Y/Z

### distortion corrections
r = math.sqrt(x**2 + y**2)

x_rad = x*(1 + cam_calib.k1*r**2 + cam_calib.k2*r**4 + cam_calib.k3*r**6 + cam_calib.k4*r**8)
x_tang = (cam_calib.p1*(r**2 + 2*x**2) + 2*cam_calib.p2*x*y)*(1 + cam_calib.p3*r**2 + cam_calib.p4*r**4)
y_rad = y*(1 + cam_calib.k1*r**2 + cam_calib.k2*r**4 + cam_calib.k3*r**6 + cam_calib.k4*r**8)
y_tang = (cam_calib.p1*(r**2 + 2*y**2) + 2*cam_calib.p2*x*y)*(1 + cam_calib.p3*r**2 + cam_calib.p4*r**4)

x_dash = x_rad + x_tang
y_dash = y_rad + y_tang
###

# Equivalent to multiplying x_dash, y_dash by intrinsic matrix K
u = cam_calib.width*0.5 + cam_calib.cx + x_dash*cam_calib.f + x_dash*cam_calib.b1 + y_dash*cam_calib.b2
v = cam_calib.height *0.5 + cam_calib.cy + y_dash*cam_calib.f
# i.e K = [[f + b1, b2, cx + width*0.5], [0, fy, cy + height*0.5], [0,0,1]]
# so [u,v] = K.mulp(x_dash, y_dash)

print('Got: ' + str(u) + ',' + str(v))
I also understand that i can invert the intrinsic matrix K in the final stage, and have also separately saved depth maps so i can later use the z values. Really it is the section where you apply the brown-conrady distortion equations which i'm unsure how to perform in reverse. Would it be possible to give any insight into how this is performed in the .unproject method? I have parameters f, cx, cy, k1, k2, k3, p1 and p2 for this particular camera model.

Thanks in advance anyone who might be able to advise!

2

General / Re: Rendering model from camera viewpoints to get z-buffer/ depth map

« on: February 17, 2021, 02:10:15 PM »

Ah I see - thanks so much for your explanation!

3

General / Re: Rendering model from camera viewpoints to get z-buffer/ depth map

« on: February 17, 2021, 12:57:21 AM »

Thanks very much for your help Alexey, this works perfectly for exporting the depth maps stored in the project.

I was wondering if you could explain what the difference is between the project depth maps and the images you get from using the dense_cloud.renderDepth() or model.renderDepth() functions. I can see visually that they are different results (attached screen shot), and the depth maps have half the dimensions - (1500, 2000) compared to (3000, 4000). I am trying to use the depth values to get the z-values in order to project detections from the 2d camera images onto the 3d model, and i'm not sure which maps i need to use for this.

Sorry if this is a very obvious question i'm very new to working with these kinds of datasets.
Thanks for your help already

4

General / Rendering model from camera viewpoints to get z-buffer/ depth map

« on: February 15, 2021, 08:52:05 PM »

Hi

I was wondering if anyone could help - I'm a bit confused about different methods for exporting the depth maps. I need to get the depth value of the model rendered from the perspective of each camera in the dataset (i.e z in "real distance" rather than as a 0-255 greyscale image). I have tried directly exporting the depth maps as .exr files with this code (for example for the first camera)

Code: [Select]

chunk = Metashape.app.document.chunk
camera = chunk.cameras[0]
depth = chunk.depth_maps[camera].image()
depth.save(camera.label + ".exr")
However when i open this is in python the resulting array is not the same dimension as the image from the camera (it's (1500, 2000) rather than (3000, 4000)), so i'm not exactly sure what these depth maps show?

I have also tried using the renderDepth function as so

Code: [Select]

depth = chunk.model.renderDepth(camera.transform, camera.sensor.calibration)from the 3D model, and also

Code: [Select]

depth = chunk.dense_cloud.renderDepth(camera.transform, camera.sensor.calibration)from the dense cloud (which would be preferable as this means i don't have to generate the mesh)

When i do depth.save(camera.label + '.tif') though it seems like the images are blank, and i'm not able to open them in Python.

Just wondering if anyone could help with clarifying what i need to get the real-world depth values for every camera? Ideally straight from the existing depth maps or dense cloud, to avoid extra processing making the mesh. Sorry i think there are several posts already answering this but i wasn't able to understand finally which method i needed!

Thanks in advance!

 

5

General / Re: Projecting from raw images to orthomosaic using chunk xml transforms

« on: February 12, 2021, 04:05:46 PM »

Hi Paul,

Thanks again! My camera.center coordinates match with the last column of the camera.transform matrix exactly as you say. This would have taken me a long time to figure out!

Also thanks for explaining the geocentric -> geographic projection. I was wondering whether it also possible to get this information using just the information from the xml, i presume from this <reference> section? I've previously used the gdal/ogr library to convert results into long/lat. But i've never worked with geocentric coordinates, so i'm not sure what to put as the crs for the input source, which i need to project into wgs84 (or vice versa).

  <reference>GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4326"]]</reference>
  <region>
    <center>8.0295666577678944e+000 7.5971928366392634e+000 -6.9026469309386114e+000</center>
    <size>4.6588402261177400e+001 3.3074866737197382e+001 1.1623963295679715e+001</size>
    <R>-1.1289887338639719e-001 -9.9224287739757688e-001 5.2037646390469303e-002 9.9353856559191167e-001 -1.1212388138353775e-001 1.7588459428196531e-002 -1.1617360700660691e-002 5.3687125805611760e-002 9.9849022501629103e-001</R>

Many thanks again,
Ellie

6

General / Re: Projecting from raw images to orthomosaic using chunk xml transforms

« on: February 10, 2021, 09:48:57 PM »

Hi Paul,

Thank you so much for your reply! Sorry i have only just got back to you, i have been looking into this problem of getting the 3D surface which you pointed out to me.

I was wondering in your code where the camera.center vector is stored (or if it is!) in the xml?  Or is there another way i can export these values in a different format (e.g omega kappa phi) in agisoft? Also just to clarify does this camera.center vector locate the middle of the 2d camera image in the 3d model reference system?

Also I realised my data is exported in WGS 84, so will i need to perform another transformation to convert it to the correct coordinate system for applying these transformations stored in the xml? I think this may be the reason i'm not able to get the values to match up at the moment. I wonder if you have any advice on how to project to and from geographic coordinate systems (outside the metashape python library).

Thank you again for your help already,
Ellie

7

General / Re: Projecting from raw images to orthomosaic using chunk xml transforms

« on: January 29, 2021, 02:41:12 AM »

Hi, i'm wondering if anyone can help with this? I have a more specific example.

I would like to use the xml to calculate the projection from a pixel in a 2d image (u, v) to a location in the point cloud (X, Y, Z). In the files i have, a specific example is getting from pixel (2521.209961, 2522.959961) in image DJI_0217 to point (105.63217, -10.47468, 243.83644). The xml files corresponding to this camera are at the bottom of this post.

I have found a post which i think outlines what i need to do https://programmer.group/xml-interpretation-of-agisoft-photoscan-mateshape-camera-parameters.html . The equation they give is


I'm unsure what parameters below should be used for dx, dy. And also R - the rotation matrix. Should this be the transforms for the specific camera, or the chunk transform matrix? I've tried a number of different combinations using the example points, but can't get them to match. I've also tried in reverse to go from 3d to 2d, following this post https://www.agisoft.com/forum/index.php?topic=6437.0, but again am unclear on where i'm going wrong with picking values from the xml. Any help would be very appreciated! Thanks for your time, sorry for the long post

-------------------------------------------------------------------------------

for the calibration i have...

Code: [Select]

      <calibration class="adjusted" type="frame">
        <resolution height="3000" width="4000" />
        <f>3118.23872343487</f>
        <cx>21.6929492116479</cx>
        <cy>-15.0408107258771</cy>
        <k1>0.0303345331370718</k1>
        <k2>-0.0667270494782827</k2>
        <k3>0.0653506193462477</k3>
        <p1>-3.94887881926124e-005</p1>
        <p2>-0.000593847873789833</p2>
      </calibration>
for the specific camera mentioned in the example i have...

Code: [Select]

    </camera>
    <camera id="215" label="DJI_0217" sensor_id="0">
      <transform>-9.9845367866015977e-001 -4.7360369045444578e-002 2.9107507860923802e-002 2.0635994236918069e+000 -5.0010632995164908e-002 9.9389702762500165e-001 -9.8324132671528716e-002 -6.4615532779234619e-001 -2.4273198335146698e-002 -9.9627776859993980e-002 -9.9472866547642924e-001 -2.8160118158196525e-003 0 0 0 1</transform>
      <rotation_covariance>8.4460505284871782e-006 -1.4902639446243033e-006 -1.1296027835376122e-008 -1.4902639446243037e-006 2.1520441968010174e-005 -6.9017167960388795e-007 -1.1296027835376149e-008 -6.9017167960388763e-007 5.7679957063865655e-006</rotation_covariance>
      <location_covariance>1.2021170339207746e-003 -1.8157680112466801e-006 -1.3976334597317290e-005 -1.8157680112466801e-006 1.2033529966674376e-003 4.4792361160838019e-005 -1.3976334597317290e-005 4.4792361160838019e-005 2.0488157555429709e-003</location_covariance>
      <orientation>1</orientation>
      <reference enabled="true" x="105.632175972222" y="-10.4748751388889" z="313.197" />
    </camera>
the chunk transform is

Code: [Select]

<transform>
    <rotation locked="true">-8.1053840961487511e-001 -5.3146724142603119e-001 -2.4610984911212641e-001 -3.3158710760897486e-001 7.0029820158928557e-002 9.4082188237520392e-001 -4.8278098202653230e-001 8.4417912529244210e-001 -2.3298954443061240e-001</rotation>
    <translation locked="true">-1.6902721015857060e+006 6.0408286681674914e+006 -1.1519789192598241e+006</translation>
    <scale locked="true">1.0044411226433386e+001</scale>
  </transform>which i understand can be made into a 4x4 tranform matrix by doing rotation times scale, and combining with translation, i.e to get
 Matrix([[-8.141381100991113, -5.33827552626121, -2.4720285313576693, -1690272.101585706],
[-3.330597266208162, 0.7034083117894531, 9.45000187740369, 6040828.668167491],
[-4.8492507157758356, 8.479282283208121, -2.3402427957204432, -1151978.919259824],
[0.0, 0.0, 0.0, 1.0]])

8

General / Projecting from raw images to orthomosaic using chunk xml transforms

« on: January 20, 2021, 09:47:17 PM »

Hi!

I'm working on a detection method for a UAV dataset, which was collected and processed in Agisoft by a partner in the project. For our method we would like to run detection on the raw images captured by the drone (in x_pixel, y_pixel form), and then project these onto the orthomosaic so we can combine overlapping detections (e.g to long/lat). However to build this into the machine learning pipeline, we would like to calculate projections in python using affine transformations, outside the agisoft software.

I was wondering if it is possible to calculate the projections using the xml file stored with the project? These have been saved from when our partner proccessed the orthomosaics and annotated the data.  I understand from a few posts on the forum this is possible using the PhotoScan python library - but am i right this can only be imported/run in agisoft not in another python environment? I can see there are parameters for the entire chunk (e.g f, cx, cy, k1, k2, etc), as well as translation and rotation/location covariance matrices for each raw image (example below). I was wondering if it would be possible to explain a little how these parameters can be used to create the projections we need? Especially if the intrinsic matrix is stored in the xml, and under which tag?

Thanks in advance!

Code: [Select]

.............................
      <calibration class="adjusted" type="frame">
        <resolution height="3000" width="4000" />
        <f>3118.23872343487</f>
        <cx>21.6929492116479</cx>
        <cy>-15.0408107258771</cy>
        <k1>0.0303345331370718</k1>
        <k2>-0.0667270494782827</k2>
        <k3>0.0653506193462477</k3>
        <p1>-3.94887881926124e-005</p1>
        <p2>-0.000593847873789833</p2>
      </calibration>
      <covariance>
        <params>f cx cy k1 k2 k3 p1 p2</params>
        <coeffs>2.0468751869991163e+001 -5.3764593396939653e-001 -2.0849002785751725e+000 6.7189315228864807e-004 -1.6497207752060089e-003 2.4646804676160135e-003 -4.2935290026870297e-006 -3.9453328768781123e-005 -5.3764593396939653e-001 1.0322535496260234e-001 5.2120752064783622e-002 -9.0786128682199680e-006 3.9834494161823351e-005 -6.0501242778801419e-005 2.2683126222501854e-006 -6.7464192211726927e-007 -2.0849002785751725e+000 5.2120752064783622e-002 2.9360214203540624e-001 -6.8661956725641512e-005 1.7218488711043673e-004 -2.5479146665632337e-004 -2.6130556639053032e-007 1.3178895802044476e-005 6.7189315228864807e-004 -9.0786128682199680e-006 -6.8661956725641512e-005 1.6020123978103154e-007 -3.1356187592838213e-007 3.4307591807736369e-007 -7.7777465362305866e-010 -4.7568949608808401e-009 -1.6497207752060089e-003 3.9834494161823351e-005 1.7218488711043673e-004 -3.1356187592838213e-007 1.0816372717144955e-006 -1.2016747628465014e-006 3.3070610393767247e-010 2.9129344608067300e-009 2.4646804676160135e-003 -6.0501242778801419e-005 -2.5479146665632337e-004 3.4307591807736369e-007 -1.2016747628465014e-006 1.3942773284397268e-006 -3.3260914196314242e-010 -4.4361470605027930e-009 -4.2935290026870297e-006 2.2683126222501854e-006 -2.6130556639053032e-007 -7.7777465362305866e-010 3.3070610393767247e-010 -3.3260914196314242e-010 3.6717553086160803e-010 -4.6052236881848388e-011 -3.9453328768781123e-005 -6.7464192211726927e-007 1.3178895802044476e-005 -4.7568949608808401e-009 2.9129344608067300e-009 -4.4361470605027930e-009 -4.6052236881848388e-011 1.6945015932751167e-009</coeffs>
      </covariance>
    </sensor>
  </sensors>
  <cameras next_group_id="0" next_id="229">
    <camera id="0" label="DJI_0002" sensor_id="0">
      <transform>9.9367761520725317e-001 -1.1154450989110640e-001 1.2752229184717922e-002 1.6355136081158577e+001 -1.1174408706270952e-001 -9.9360520953338949e-001 1.6184764274795770e-002 -7.2377742564520382e+000 1.0865359752364273e-002 -1.7507424175531350e-002 -9.9978769449128035e-001 4.9736616107715687e-001 0 0 0 1</transform>
      <rotation_covariance>9.8217092719010480e-006 -4.2857770172187034e-007 3.2338520246152594e-008 -4.2857770172186997e-007 2.1911384793986368e-005 -1.1127430980809863e-006 3.2338520246152594e-008 -1.1127430980809863e-006 5.8461829356115879e-006</rotation_covariance>
      <location_covariance>2.0055242000696963e-003 1.8493823777608193e-006 -9.6731442372176409e-006 1.8493823777608193e-006 2.0218170118861392e-003 9.1746781498499200e-005 -9.6731442372176409e-006 9.1746781498499200e-005 3.4676077055425325e-003</location_covariance>
      <orientation>1</orientation>
      <reference enabled="true" x="105.633025333333" y="-10.4760358611111" z="313.197" />
.......................