I want use the detect informations on 2D image , how can I mapping the 2D point to 3D point though OpenSfm?

[baihaozi12]

I can running this project, now I want to use the camera matrix and depth map to calculate the 3D(X,Y,Z) point postition on the point cloud? Or a approximation value is OK. Thank you

[mattip]

I have a start for this, but am also looking for an authoritative answer:

from opensfm import dataset
data = dataset.DataSet(path-to-data)
recs = data.load_reconstruction()[0]
shot = rec.shots[image-name]
cam = shot.camera

For simplicity, assume I want to get the coordinates of the 4 image corners 100 meters in front of the camera, what would be the next steps? I expected to see back_project() on the camera, but it seems that function is not exported to python?

[fabianschenk]

Hi @baihaozi12 and @mattip ,

In OpenSfM you can compute a bearing that gives you a unit vector in the direction of the 3D point. You then have to scale it by the depth and you're done. For one corner pt2D = [0.3, 0.4] it would work something like this (untested):

from opensfm import dataset
data = dataset.DataSet(path-to-data)
recs = data.load_reconstruction()[0]
shot = rec.shots[image-name]
cam = shot.camera
pt2D = [0.3, 0.4] # must be in normalized coordinates
# if you only have pixels, use cam.PixelToNormalizedCoordinates(..) to convert
bearing = cam.bearing(pt2D) # unit vector [x,y,z]
bearing /= bearing[2] # Now [x/z, y/z, 1] to step exactly 1 unit, e.g. meters, in the z direction
pt3D_100 = bearing * 100 # this point is in the camera coordinate system
pt3D_100_world = shot.pose.inverse().transform(pt3D_100)

For one particular point in berlin (tested):

     pose = pygeometry.Pose(
        np.array([1.2549079575524087, -0.40135469242527716, 0.5410499805923656]),
        np.array([-0.48666580475695764, 6.902214741262238, 7.5422937224223965]),
    )
    cam = pygeometry.Camera.create_perspective(
        0.9619961153009631, 0.023214677968683928, -0.00041315211971044893
    )
    cam.width = 3264
    cam.height = 2448
    pt3D_gt = np.array([5.13047795171054, 36.185508557877895, 22.927045680858683])
    pt2D = np.array([-0.376388, -0.141775])
    pt3D_cam = pose.transform(pt3D_gt)
    bear = cam.pixel_bearing(pt2D)
    bear = bear/bear[2] # now 1 m in z direction
    bear *= pt3D_cam[2]
    pt3D_world = pose.inverse().transform(bear)
    print(pt3D_world, "==", pt3D_gt) # very small numeric difference

[mattip]

Thanks!

[mattip]

Sorry to bother you. I made some progress, but something is off. I am using OpenSFM in the framework of OpenDroneMap, experimenting with their Brighton Beach dataset. I can do something like this to project the picnic table from pixel 3707, 718 in the original image DJI_0028.JPG (way over on the right edge):

import gdal, affine, re, pyproj
from opensfm import dataset
base_dir =<path-to-processed-data>
data = dataset.DataSet(base_dir + '/opensfm')
rec = data.load_reconstruction()[0]

# Get the offset from the openfsm coordinate system to the geotiff's
dt = gdal.Open(base_dir + '/odm_orthophoto/odm_orthophoto.tif')
gdal_proj = dt.GetProjection()
match = re.search('zone \d*[NS]', gdal_proj)
zone = match[0].replace(' ', '=')
proj = pyproj.Proj(f'+proj=utm +ellps=WGS84 +datum=WGS84 +{zone}')

delta = proj(rec.reference.lon, rec.reference.lat)

# Choose pixel. First get the camera bearing and height
pose = rec.shots['DJI_0028.JPG'].pose
camera = rec.shots['DJI_0028.JPG'].camera
rel_pix = camera.pixel_to_normalized_coordinates((3707, 718))
bearing = camera.pixel_bearing(rel_pix)
bearing /= bearing[2]
# Assume the geotiff is at sea-level since this is a beach
# so cam_pos[2] is the height
cam_pos = pose.transform_inverse((0, 0, 1))
rel_pos = pose.transform_inverse(bearing * cam_pos[2])
# Now rel_pos is in the opensfm coordinate system. Move it to WGS84
utm_pos = rel_pos[:2] + delta
# Get the GPS coordinate
picnic_gps = proj(utm_pos[0], utm_pos[1], inverse=True)
# Get the pixel on the geotiff
m = affine.Affine.from_gdal(*dt.GetGeoTransform())
m_inv = ~m
geotiff_pix = m_inv * utm_pos

picnic_gps is (-91.99252174422621, 46.84167993640242), which is in the sea on google maps. Maybe that is OK because all I care about is mapping the image pixel to the geotiff? Nope, the geotiff_pixel is off too: it is (2937.857387267053, 2802.11845625937) which is also off in the sea.

Am I using the rec.reference wrong? The geotiff projection is "WGS 84 / UTM zone 15N"

The relevant part of the reconstruction.json is below

Thanks again for any thoughts.

[
    {
        "cameras": {
            "v2 dji fc300s 4000 2250 brown 0.5555 rgb": {
                "projection_type": "brown",
                "width": 4000,
                "height": 2250,
                "focal_x": 0.5574477107783634,
                "focal_y": 0.5574477107783634,
                "c_x": -0.005111158539149004,
                "c_y": 0.006883416343247216,
                "k1": 0.0073820694983555875,
                "k2": 0.007197529667484939,
                "p1": 0.000886923323385389,
                "p2": 0.0014776283138854302,
                "k3": 0.00567955978276901
            }
        },
        "shots": {
            "DJI_0028.JPG": {
                "rotation": [
                    -2.8531926645589927,
                    1.2828698465726849,
                    -0.01623135819146115
                ],
                "translation": [
                    -4.813656862428607,
                    -23.46603285081447,
                    198.53652586066767
                ],
                "camera": "v2 dji fc300s 4000 2250 brown 0.5555 rgb",
                "orientation": 1,
                "capture_time": 1466699591.0,
                "gps_dop": 10.0,
                "gps_position": [
                    -15.377324314795588,
                    -15.05441859612256,
                    198.60896370653063
                ],
                "vertices": [],
                "faces": [],
                "scale": 1.0,
                "covariance": [],
                "merge_cc": 0
            }
        },
        "reference_lla": {
            "latitude": 46.842656666666656,
            "longitude": -91.99400802469134,
            "altitude": 0.0
        }
    }
]

[mattip]

Once this is working I would be happy to add it to the docs or a wiki if you think it is appropriate.

[baihaozi12]

Thank you!

[mattip]

Digging around in OpenDroneMap, they create a odm_georeferencing/coords.txt. The top lines there confirm my zone (line 0) and delta (line 1) calculation

WGS84 UTM 15N
576705 5188170

[fabianschenk]

Hi @mattip , OpenSfM uses a local reference instead of absolute WGS84 coordinates. The reference is the TopocentricConverter stored in the map. To get the absolute WGS84 coordinates from local gps coordinates gps_local, you would have to run something like (untested!):

gps_world = reconstructio.get_reference().to_lla(gps_local)

[mattip]

Thanks. That means I don't need the delta calculations, reducing the script to

pose = rec.shots['DJI_0028.JPG'].pose
camera = rec.shots['DJI_0028.JPG'].camera
rel_pix = camera.pixel_to_normalized_coordinates((3707, 718))
bearing = camera.pixel_bearing(rel_pix)
bearing /= bearing[2]
# From the DTM, ground height is ~160m
cam_pos = pose.transform_inverse((0, 0, 0))
rel_pos = pose.transform_inverse(bearing * (cam_pos[2] - 160))
picnic_gps= rec.get_reference().to_lla(*rel_pos)

That seems to align pretty well with the data. My mistake was thinking that the beach was at sea level, this is a lake at 160m above sea level, not the ocean.

[mattip]

My offer to make a PR with this snippet still stands: is there a need and a place for this in the docs?

[fabianschenk]

Hi @mattip ,

Thanks for you offer. I'm not a 100% sure where to put it but https://www.opensfm.org/docs/cam_coord_system.html might be a good place. The files for the webpage are in: https://github.com/mapillary/OpenSfM/tree/master/doc/source. Might be good to have a self-contained example, e.g. with the berlin dataset. Create a PR and then we will review.

Best, Fabian

[Oringa]

Hi @baihaozi12 and @mattip , For one particular point in berlin (tested):

     pose = pygeometry.Pose(
        np.array([1.2549079575524087, -0.40135469242527716, 0.5410499805923656]),
        np.array([-0.48666580475695764, 6.902214741262238, 7.5422937224223965]),
    )
    cam = pygeometry.Camera.create_perspective(
        0.9619961153009631, 0.023214677968683928, -0.00041315211971044893
    )
    cam.width = 3264
    cam.height = 2448
    pt3D_gt = np.array([5.13047795171054, 36.185508557877895, 22.927045680858683])
    pt2D = np.array([-0.376388, -0.141775])
    pt3D_cam = pose.transform(pt3D_gt)
    bear = cam.pixel_bearing(pt2D)
    bear = bear/bear[2] # now 1 m in z direction
    bear *= pt3D_cam[2]
    pt3D_world = pose.inverse().transform(bear)
    print(pt3D_world, "==", pt3D_gt) # very small numeric difference

@fabianschenk A small issue with this solution - in order to obtain the scale for the bearing, pt3D_cam[2] should be obtained. However, in order to obtain pt3D_cam[2] we must first have the 3D point in the point cloud and the corresponding pixel in the image. Is it possible to get the 3D point from the 2D pixel without having the 3D point beforehand?

[wilsonmichaelc]

Hi @fabianschenk

I am also trying to find the 3D point from 2D pixel. It seems that "load_clean_depthmap" might be useful for making this z adjustment but I am having trouble with the depthmap which seems to have a different size than my image. If I determine relative pixel coordinates (rel_pix = camera.pixel_to_normalized_coordinates((2115, 2577))) how does this map to the depthmap?

Cheers,

-Michael