alicevision / Meshroom

3D Reconstruction Software
http://alicevision.org
Other
11.04k stars 1.07k forks source link

[question] georeferencing #1539

Open andrea-bistacchi opened 2 years ago

andrea-bistacchi commented 2 years ago

Hello, to use the output of Meshroom (both point clouds and textured surfaces), georeferencing in some external coordinate system is needed. It is not clear to me if this can be performed in Meshroom.

When you import images collected with a drone, the GPS EXIF data are read:

image

Hover these are in latitude/longitude (degrees) and meters (elevation) so they are not directly suitable for georeferencing.

It would be much better to be able to (alternatively):

1) read a file with a table of camera file names and X Y Z coordinates (transformed in some suitable coordinate system with external software), and use these to set the output transformation.

2) manually set a transformation by picking targets in the scene and defining their coordinates.

These are the methods used in other photogrammetry applications (e.g. VSFM etc.) and they work very well.

If you need help with this, we can help (if it is possible to do this in Python), but we would need some guidance.

Note that I am assuming that the coordinate transformation is performed at the last step in the processing, as in VSFM, while all other steps are performed in "internal" "model" coordinates.

Thanks very much!

simogasp commented 2 years ago

Hi, some preliminary support for gps has been added recently (see #1477 and https://github.com/alicevision/AliceVision/pull/1069). There is no support for it during the sfm step (i.e. gps as a prior on the poses) but you can rescale the whole sfm scene using the SfMTransform node and selecting from_gps in the combo box. You can then choose if only apply scale or even the X Y Z world coordinates (in the latter case the model is likely to exceed the meshroom volume to be displayed correctly).

Please feel free to test it and give us feedback.

andrea-bistacchi commented 2 years ago

Hi, and thanks for super-prompt answer!

Doing the transformation at the end is the right way to do it, at least according to what all other SFM codes do. In fact (according for instance to VSFM's developers) doing the transformation earlier could decrease accuracy.

So I'll test the transformation in the SFMtransform node tomorrow!

However I already know that being able to use custom coordinates that you read from a table is extremely useful.

For instance in many projects I just drop the kilometer digits of UTM coordinates, or I shift them by some constant value, in order to obtain a local model that is easier to manage.

Is there a way to alter the coordinates that Meshroom reads from the EXIF?

And if there are no coordinates in the EXIF, e.g. in terrestrial projects with photos collected with a DLSR camera, is it possible to add them somehow?

I guess that the coordinates are hidden somewhere in some of the project files, and if it is an ASCII file altering the coordinates will be easy.

Thanks very much!

andrea-bistacchi commented 2 years ago

Hello, the problem is from_gpsis still missing.

image

I am using Meshroom 2021.1.0 for Windows downloaded yesterday.

Shall I download or compile a different release?

BTW, where should I look for documentation on the other transform methods? Just a couple are cited in the manual and not in a very informative way.

Thanks!

simogasp commented 2 years ago

For instance in many projects I just drop the kilometer digits of UTM coordinates, or I shift them by some constant value, in order to obtain a local model that is easier to manage.

Is there a way to alter the coordinates that Meshroom reads from the EXIF?

Yes, in the SfMTransform you can choose what to apply from the transformation computed from the GPS that brings everything in the world coordinates. It depends on what you select here image

If you only select the scale, the first camera stays in the origin and all the others are updated according to the scale. Otherwise with rotation and translation selected as well the cameras get the XYZ "world" coordinates.

And you can apply an additional scale to reduce the size of the model (if I remember correctly everything is in meters) so you can divide by 10^3 or more to get in a different unit.

simogasp commented 2 years ago

Yes it's in the develop branch, you need to wait for the upcoming release or rebuild everything on your own.

andrea-bistacchi commented 2 years ago

Uh... thanks. Any idea on a release date? Would it be worth rebuilding on my own?

simogasp commented 2 years ago

BTW, where should I look for documentation on the other transform methods? Just a couple are cited in the manual and not in a very informative way.

yes the documentation is a cry for help... This is what the node does from the documentation embedded in meshroom

This node allows changing the coordinate system of one SfM scene.
The transformation can be based on:
      transformation: Apply a given transformation
      auto_from_cameras: Fit all cameras into a box [-1,1]
      auto_from_landmarks: Fit all landmarks into a box [-1,1]
      from_single_camera: Use a specific camera as the origin of the coordinate system
      from_markers: Align specific markers to custom coordinates

Yes it's in the develop branch, you need to wait for the upcoming release or rebuild everything on your own.

@fabiencastan can have a better idea for the release As for rebuilding from scratch you need to be a little familiar with the developing environment (visual studio vcpkg and so on).

andrea-bistacchi commented 2 years ago

Yes, in the SfMTransform you can choose what to apply from the transformation computed from the GPS that brings everything in the world coordinates. It depends on what you select here image

If you only select the scale, the first camera stays in the origin and all the others are updated according to the scale. Otherwise with rotation and translation selected as well the cameras get the XYZ "world" coordinates.

And you can apply an additional scale to reduce the size of the model (if I remember correctly everything is in meters) so you can divide by 10^3 or more to get in a different unit.

This is all very good!

May I suggest also an option to apply a translation with "round numbers" more or less as in CloudCompare?

I explain with an example. Let's say that the first camera is at Lat/Long 46°14'54.89"N 8°14'30.16"E, that in UTM Z32T is X=441547.502 Y=5121946.513 (meters), with an elevation of 2345.087 m.

Dropping the lengthy UTM coordinates make sense because it makes your data lighter-weight and allows to increase the number of significant digits, so this is done quite often, particularly if you convert to mm.

However, if you set translation to False in SfMTransform, you will have to add the coordinates of the first camera to all your output/results in order to project it in real world. Something that can be important in applications.

An alternative is to drop just the kilometres, so that the first camera coordinates will be, in your "local" reference, X=547.502 Y=946.513 Z=345.087. Then converting to "real" UTM is easier because you just have to add deltaX=441000.0 deltaY=5121000.0 deltaZ=2000.0. This is how CloudCompare deals with this problem and I find it very simple and intuitive. And you can also combine this with a 10^3 conversion from m to mm and viceversa.

Let me know if I can contribute somehow on this. Knowing where to work on the code, it would be very easy.

Thanks!

andrea-bistacchi commented 2 years ago
This node allows changing the coordinate system of one SfM scene.
The transformation can be based on:
      transformation: Apply a given transformation
      auto_from_cameras: Fit all cameras into a box [-1,1]
      auto_from_landmarks: Fit all landmarks into a box [-1,1]
      from_single_camera: Use a specific camera as the origin of the coordinate system
      from_markers: Align specific markers to custom coordinates

BTW, how are the "markers" defined and collected?

Thanks again.

simogasp commented 2 years ago

This is all very good!

May I suggest also an option to apply a translation with "round numbers" more or less as in CloudCompare?

I explain with an example. Let's say that the first camera is at Lat/Long 46°14'54.89"N 8°14'30.16"E, that in UTM Z32T is X=441547.502 Y=5121946.513 (meters), with an elevation of 2345.087 m.

Dropping the lengthy UTM coordinates make sense because it makes your data lighter-weight and allows to increase the number of significant digits, so this is done quite often, particularly if you convert to mm.

However, if you set translation to False in SfMTransform, you will have to add the coordinates of the first camera to all your output/results in order to project it in real world. Something that can be important in applications.

An alternative is to drop just the kilometres, so that the first camera coordinates will be, in your "local" reference, X=547.502 Y=946.513 Z=345.087. Then converting to "real" UTM is easier because you just have to add deltaX=441000.0 deltaY=5121000.0 deltaZ=2000.0. This is how CloudCompare deals with this problem and I find it very simple and intuitive. And you can also combine this with a 10^3 conversion from m to mm and viceversa.

Let me know if I can contribute somehow on this. Knowing where to work on the code, it would be very easy.

Thanks!

That's really interesting as it would still make the model visible in the 3D viewer. We can maybe add the option to drop as a flag to the SfMTransform application (and a checkbox in meshroom). I don't know how and/or where to save the delta translation if it is needed later, though. Maybe we can embedded it into the SfM data file as a global data member to be used in case you want to get back to the real coordinates? What do you think or suggest? Is it something that could be needed later?

BTW the code for the transformations is all here https://github.com/alicevision/AliceVision/blob/develop/src/aliceVision/numeric/gps.hpp https://github.com/alicevision/AliceVision/blob/develop/src/aliceVision/numeric/gps.cpp and the SfMTransform that call those functions https://github.com/alicevision/AliceVision/blob/develop/src/software/utils/main_sfmTransform.cpp#L402

Please feel free to check the code as I tried my best to find the correct transformations to cartesian coordinates but I'm no expert on that ;)

simogasp commented 2 years ago

BTW, how are the "markers" defined and collected?

if the markers are detected and reconstructed in the scene you can pass the the IDs of the markers and their (metric) XYZ to rescale the reconstruction. In the markers field you can pass a list of ID:x,y,z

andrea-bistacchi commented 2 years ago

Please feel free to check the code as I tried my best to find the correct transformations to cartesian coordinates but I'm no expert on that ;)

From the code at:

https://github.com/alicevision/AliceVision/blob/2622a945737c65598e7cd40808fae0f3332c288b/src/aliceVision/numeric/gps.cpp#L15

it looks like you are converting Lat/Long/Elevation to earth-centered, earth-fixed geocentric XYZ, so coordinates in a system like this:

image

Is this just my impression or is it true?

FlachyJoe commented 2 years ago

if the markers are detected and reconstructed in the scene you can pass the the IDs of the markers and their (metric) XYZ to rescale the reconstruction. In the markers field you can pass a list of ID:x,y,z

Are ArUco markers usable in the same way / Does exist an ArUco detector in meshroom ?

simogasp commented 2 years ago

Is this just my impression or is it true?

yes that was the idea

simogasp commented 2 years ago

if the markers are detected and reconstructed in the scene you can pass the the IDs of the markers and their (metric) XYZ to rescale the reconstruction. In the markers field you can pass a list of ID:x,y,z

Are ArUco markers usable in the same way / Does exist an ArUco detector in meshroom ?

not at the moment, recently some support for aprilTag has been added. It would be nice to add aruco as well

andrea-bistacchi commented 2 years ago

Is this just my impression or is it true?

yes that was the idea

OK, I see now why you do not like the coordinate that you obtain. They will always be around the size of the Earth's radius.

In applications these coordinates are never used, also because the Z has nothing to do with the local "vertical".

What people do is to use geographic projections to convert from Lat/Long into meters.

This is a bit complicated because there are a multitude of geographic projections and the results of each projection (basically a coordinate transformation formula) are different if the datum changes. The datum in turn is given by an ellipsoid (semi-major and semi-minor axes) and the orientation and translation of that ellipsoid with respect to the standard WGS84 datum.

To learn almost all you need about this, there is a standard publication by the United States Geological Survey:

Snyder John P. , 1987, Map projections: A working manual, USGS Professional Paper 1395

But the good news are that you do not need to study all this! There is a library called PROJ that is the standard and manages all this for you. I really recommend to use it.

andrea-bistacchi commented 2 years ago

On markers, would it be possible to use "manual" targets that you simply pick in 2D from photos or in 3D (from dense cloud or textured mesh)? Quite often we use natural features or man-made landmarks that you can measure in the field with your GPS...

andrea-bistacchi commented 2 years ago

Alternatively to using PROJ, you can switch to the method used in VSFM, that in the end is very flexible.

Here the camera coordinates (as well as marker coordinates) are imported from a file (e.g. a CSV with camera, x, y, z).

In this case the problem of properly converting lat/long to projected coordinates is taken care of with some external software, such as a GIS, and the EXIF is not used directly. Then you can also apply scale, shift, etc very easily (e.g. in Excel).

In the end this is less automatic but more flexible and very easy to implement.

If I may, I would suggest to start with a solution like this, and then add PROJ when more time for this issue is available.

Thanks!!

andrea-bistacchi commented 2 years ago

Hello, this should be the relevant page in PROJ documentation:

[https://proj.org/development/reference/functions.html]()

Basically you have to define a transformation with proj_create_crs_to_crs, and then you can use it to transform your coordinates with proj_trans.

andrea-bistacchi commented 2 years ago

Or you can use the Python wrapper pyproj4.

In this case you just need to follow the example code here:

https://pyproj4.github.io/pyproj/stable/examples.html#converting-between-geographic-and-projection-coordinates-within-one-datum

I would do it myself, but I do not know in depth the Meshroom GUI and it will take forever for me to design the pull down lists to select the CRS (coordinate reference system) from those available in PROJ.

andrea-bistacchi commented 2 years ago

One possible source of errors with PROJ is that geographical coordinates in Latitude, Longitude are supposed to have this axis order, i.e. Y,X order. However it is possible to always use X, Y order for both projected and geographic coordinates - X, Y or Longitude, Latitude order - by using the always_xyoption when defining the transformation.

https://pyproj4.github.io/pyproj/stable/examples.html#step-2-create-transformer-to-convert-from-crs-to-crs

simogasp commented 2 years ago

Thanks for all this info!

When I started think about using the exif information I saw that there were way too many different projections that could be used. Since I'm not an expert and I didn't have time to spend on learning all the possible ways I went for the "easiest" one with the cartesian representation.

Is there a transformation that is used more than others or really each domain has its preferred way to represent the coordinates according to the application or to the scope?

Adding another dependency just for that is another step I didn't want to go through because we have already enough dependencies. I saw that boost also had a geometry module in which there are some projections implemented. https://www.boost.org/doc/libs/1_76_0/boost/geometry/srs/projections/proj/somerc.hpp https://www.boost.org/doc/libs/1_69_0/boost/geometry/srs/projections/proj/etmerc.hpp https://www.boost.org/doc/libs/1_67_0/boost/geometry/srs/projections/proj/tmerc.hpp https://www.boost.org/doc/libs/1_67_0/boost/geometry/srs/projections/proj/geos.hpp That would be nice because we already have this.

Anyway, maybe the idea of having the position passed with a csv file is more interesting so people can choose their own projection and just pass it to the file.

andrea-bistacchi commented 2 years ago

Hello, basically you can follow two strategies, that are followed also by other SFM codes and are well known by end users (as I am since around 10 years!).

(1) Use PROJ. As you said map projections is a terribly complicated thing, and finding bugs etc. is very difficult, so starting from scratch is not advisable. PROJ is THE standard that everybody uses (including all main GIS projects like QGIS) and there is no reason to use something different or (even worst) write new code. I understand your strategy not to add too many dependencies, but actually I would not worry about PROJ. It is a very mature and stable library, and many different projects use it (e.g. the huge QGIS project). - This is the strategy used by Agisoft -

(2) Do not use EXIF directly, and let the user add a simple table with just four columns: image_name, X, Y, Z. Then the user must find a way to obtain and convert meaningful coordinates, but this is easy. If you like I can write a short guide for you. Basically you just have to dump the EXIF coordinates in a table (there are command-line applications for this), do something to convert it to the projection you need (e.g. use PROJ as command line, or use PROJ within QGIS), then load back the coordinates table to Meshroom. The advantage here is simplicity for you who write the code, and flexibility for the end user (if I like to scale the coordinates, I can multiply them in Excel, if I like to shift them, I just subtract a constant in Excel, if I have some bad GPS point, I can drop them). Important to add even more flexibility: leave the possibility to have some not-georeferenced camera, that does not appear in the table. This is fundamental for instance when you combine both drone (with GPS) and terrestrial (no GPS) photos in your project. - This is the strategy used by VSFM -

Now you might wonder what I suggest.. for me, as an experienced user that knows something about map projections, option (2) is the best, because of flexibility. For an unexperienced user option (1) could be easier.

In the end I suggest to start from (2) - it will take minutes for you to implement it - and then maybe add also (1) when you have time.

In case you decide to use option (2) just let me know and I'll write a short tutorial on how to obtain a CSV table with coordinates.

andrea-bistacchi commented 2 years ago

Update: we always end with PROJ...

This is what I found in Boost:

// Boost.Geometry - extensions-gis-projections (based on PROJ4)

So you are already depending on PROJ.

Now you just have to see if it is easier to use it through Boost or directly.

But still the CSV table method could be useful and very simple to be implemented.

andrea-bistacchi commented 2 years ago

Hello, did you manage to test new implementations?

Thanks!

simogasp commented 2 years ago

Hi sorry I didn't have time to look into the boost/proj. The csv file, though, should be easy to implement and we can add it as a way to scale in the existing node/binary app. Do you happen to have an example of a (small) dataset that we can use to test the code? Only the .abc/.sfm with the cameras and the csv are needed, no images.

natowi commented 2 years ago

You can find a sample dataset with coordinates in EXIF and GCP here https://github.com/OpenDroneMap/odm_data_bellus/tree/master

I used exiftool to batch export the gps information in csv:

exiftool -csv -filename -GPSLatitude# -GPSLongitude# -GPSAltitude# -r D:\odm_data_bellus-master\odm_data_bellus-master\images > D:\odm_data_bellus-master\odm_data_bellus-master\images\coord1.csv

coord1.csv

Maybe we could adopt the format that is used by odm for gcp: https://docs.opendronemap.org/gcp/

Here are the .abc/.sfm files: abcsfm.zip