Try to fix georeferencing accuracy

[originlake]

This issue was discussed earlier here https://community.opendronemap.org/t/model-accuracy-issues-related-to-geocoord-exports/13054/1

The horizontal error can get up to meters when the GCPs are a few kilometer away from the dataset center. The change would definitely slow down the processing speed a bit, but I think such big errors shouldn't be acceptable. My proposed fix is

1. Skip exporting geocoord in opensfm stage, so the file generated will all be in topocentric coordinate system.
2. In the following stages, openmvs, filter_pointcloud, odm_meshing, odm_texturing will all generate assets in topocentric coordinate system, assets will include _topocentric in name. openmvs and odm_texturing do require the camera parameters remain unchanged to guarantee accuracy, so all these stages have to be kept in topocentric
3. In odm_georeferencing stage, convert all the topocentric assets that are needed in the final output to UTM coordinate system, the conversion should be done point by point instead of using a rigid affine transformation matrix.
4. The orthophoto and dsm will be generated with assets in UTM coordinate system.

The current implementation is done by using pdal-python to read and write the points, and use numpy to do bulk conversions. For obj file, I simply followed the implementation in the alignment method. For efficiency, if we can use newer version of PROJ lib, it supports topocentric conversion method, so we could define a proj pipeline and use pdal solely to do the conversion https://proj.org/en/stable/operations/conversions/topocentric.html#geocentric-to-topocentric-conversion. This should provide faster speed.

Currently I implemented converting odm_georeferenced_model.laz, point_cloud.ply and odm_textured_model_geo.obj, I need to check 25dmeshing related files as I don't normally use it. Let me know if any other file should be included.

You could find some error calculations numerically in the post above. Below is a comparison between orthophotos processed in two ways, the gcp is around 1km away from the dataset center

Rigid transfrom(current way)

Point to point transform(proposed fix)

[pierotofy]

Hey thanks for the PR @originlake ! This is some good work on addressing the accuracy issue.

I'm concerned about introducing a topocentric version of all assets, would it not be possible to just do the point-by-point transformation in export_geocoords (and perhaps run bundle adjustment with cameras, GCPs and points fixed after the transformation to address small changes in camera rotations/angles, assuming that we can only transform coordinates and not rotations)?

[originlake]

I think the 3D coordinate system with UTM zone is not an accurate 3D cartesian coordinate system. The horizontal part is a traverse mercator projection, the warping is different from point to point, 1 meter on map could not be accurate 1 meter in real world, and the vertical part is the ellipsoid height with respect to ellipsoidal earth surface, which is also not a simple offset from topocentric coordinate system. I think doing bundle adjustment in such coordinate system could lead to worse camera parameter accuracy and hence affecting the dense point cloud accuracy? After all the 3d model accuracy itself is good, it's just the coordinate system conversion is not done accurately.

[pierotofy]

I'm not sure to what degree this would affect accuracy, but I have a feeling that it might possibly be better, currently the transformation is done post-SFM/MVS, so while the absolute accuracy is higher, your relative accuracy is probably a bit lower (the code warps the models). Bundle adjustment might be able to find an optimal middle ground between the two. It's something that I would need to test.

[originlake]

IMHO, the problem is not the 3D model being warped hence inaccurate.
The 3d reconstruction is designed to be done in topocentric coordinate system at the beginning because it's being an isotropic space. It's not suitable in projected coordinate system, projected coordinate system is non-isotropic, an extreme example, the notoriously mercator effect, the area of Greenland is larger than USA in web Mercator coordinate system. Assuming we have a dataset with drone taking photos covering the whole earth, we definitely can't use web mercator coordinates to do the bundle adjustment.
https://en.wikipedia.org/wiki/Web_Mercator_projection#/media/File:Web_maps_Mercator_projection_SW.jpg

This doesn't mean a coordinate in the web-mercator CRS is wrong, it's just a reference that we can use to find the actual location on earth, when needed for accurate 3d space calculation, we can always convert it back to topocentric coordinate system. Usually in order to get good measurement accuracy, the surveyor will choose a coordinate system where center is defined closed to the survey region, some will even create a custom crs with origin center defined as their site center for best accuracy. They don't do topocentric, as far as I can tell, is because topocentric is inconvenient to do mapping on earth surface.

Same for the geo coordinate exporting, at this stage, it's not about making the model accurate, it's about do the right conversion to the target CRS so we can still use it to accurately find the exact spot on the earth, if we need to export UTM, we do the correct conversion to UTM, if we need to export to US state plane, we do the correct conversion to US state plane, we don't need to do bundle adjustment each time exporting to a different coordinate system. The point location is always unique on earth, but the coordinate representation is different in different coordinate systems.

Back to the problem, the 3D model itself is not warped after the conversion, its exact location on earth remains the same, it's just being represented in a warped coordinate system, doing measurement in such CRS is inaccurate, but it's the surveyor's responsibility to choose a better accurate CRS if needed.

Anyway, these are just in theory, the scale of a normal datasets would probably not have large errors. We can probably compare the reprojection errors to see the difference.

[smathermather]

Back to the problem, the 3D model itself is not warped after the conversion, its exact location on earth remains the same, it's just being represented in a warped coordinate system, doing measurement in such CRS is inaccurate, but it's the surveyor's responsibility to choose a better accurate CRS if needed.

Anyway, these are just in theory, the scale of a normal datasets would probably not have large errors. We can probably compare the reprojection errors to see the difference.

Yeah, bit of a classic problem. Actual reprojection errors are (probably) sufficiently small to be meaningless on most projects of most sizes. But, folks who care will want to know that the math is done correctly so they don't have to worry about unanticipated edge cases.

[Saijin-Naib]

From the Community, I mostly only see this issue raised by folks who have/use RTK GCPs, or have known survey markers they are comparing to and expect survey-grade accuracy with RTK-enabled input data.

So, a cross-section of industry, commercial, and academic folks (mostly), and the occasional very savvy home-gamer.

I think there is potential for a lot of benefit here.

[pierotofy]

Indeed; does anyone have a small dataset with GCPs that clearly triggers this problem?

[Saijin-Naib]

https://community.opendronemap.org/t/call-for-datasets-gcps-drifting-at-boundaries/23934?u=saijin_naib

[Saijin-Naib]

We have seen the data, so lets hope some of it is readily sharable.

[pierotofy]

@originlake I've looked more into this. I've implemented Kabsch's algorithm here (pierotofy/OpenSfM@7366eec) which instead computes a rigid transformation by sampling points from the area covering the model (either via GCP locations, or a bbox computed from the camera poses).

I tested this on existing datasets with and without GCPs and seems to work on smaller areas, but I don't have a large dataset to test with, so perhaps you can try to replace opensfm/actions/export_geocoords.py with my copy and report back if it improves the results for your dataset?

If this works, then there's no need to keep separate models in topocentric coordinates.

[originlake]

OK, I'll test your updates and report back. Unfortunately I can't share my dataset

[pierotofy]

Also do note I might be wrong here, but it's difficult to try things without a dataset to test.

[originlake]

It's hard to show it without zooming around the map, below are some screenshots around the markers, the error still exists. The error vectors, from the actual marker on the ground to the gcp coordinate point, is generally pointing to the center of the dataset. There are scaling involved as well, so rigid transform could not solve the problem. The conversion from topocentric coordinate system to projected coordinate system is non-linear, so affine transformation(linear transformation) would not able to solve it accurately, it can be approximated to minimize the target loss, but probably still can't be precise enough.

[originlake]

I created a colab notebook can simulate the issue. The Idea is to create a plane with grids of coordinates in topocentric coordinate system, easting, northing are ranged from -2000m to 2000m, elevation is set to 1000. Every points in the grid will be converted to UTM representing the most accurate conversion as a baseline. Then can apply different method(Kabsch algorithm is used) to compute the errors between the baseline.

At the bottom, I plotted the errors between two conversions, separated horizontal error and elevation error for better visualization. You can see the both the horizontal errors and vertical errors are non-uniform, and the error vector's direction is different.

You can check it here, feel free to modify it https://colab.research.google.com/drive/1J3VBcmLbZ4YqVTrNXSisbjDr41E75ek3#scrollTo=lEMSfsirnm8T

[pierotofy]

Thanks for testing and the colab. It makes sense.

Keeping everything in topocentric seems like the way to go, I do have some concerns around performance (in particular for the PDAL pipeline stuff) and the extra 3D models, especially for multispectral datasets (which are missing currently) and the 2.5D only workflows where a 3D model is not built. Also disk usage concerns for downstream applications like NodeODM, which would currently include both textured models (topo and UTM), which is no good. Plus careful review and tons of testing, which is needed for changes like this.

As a heads up, it might take me a while to take an active role in doing all of the above, due to other work priorities, but I do think we need to include this fix because it's critical for accurate results. 🙏

[originlake]

understandable, I was unsure about whether to keep a standalone topocentric file or simply overwrite it. But in order to make rerun work properly, this is the only way I could think off. Maybe you can come up with a better idea.

For texture model conversion, only the obj file needs to be converted. The implementation I followed in the alignment feature, simply read the vertex and do conversion line by line, which should be a lot less efficient than using numpy to convert all vertices at once. Should probably introduce an efficient parser for better performance.

For point cloud conversion, if we can have a newer version of PROJ lib, version 8.0 at least, with this new feature combined with pdal filters.projpipeline , should allow us to get rid of the python part of conversion to improve the speed.

Let me know what I can help with.

[originlake]

Updated the code to convert all geomodel objs, the alignment part provided good reference for me to use. Still need to do GLTF conversion. But I think we should probably move the GLTF exporting to somewhere after the georeferencing stage, I noticed that the alignment section doesn't align GLTF file either.

[pierotofy]

I would expect a lot of things to be needing re-arrangement (and that's not a trivial thing to do, because there's subtle things that require assets to be generated in a certain order, depending on flags). That's why I was hoping that transforming the coordinate system + running bundle adjustment would work. It would make the fix a lot easier.

When I have time (and a dataset) I think I will try this approach. I know it's not ideal in theory, but we're talking about differences of centimeters (over a large area) and I'm just plain curious to see what would happen in practice.

[pierotofy]

Note to self: all calls to opensfm_reconstruction_average_gsd would break with this PR because it assumes a Z-up axis, but that can only be guaranteed after georeferencing.

Also possible memory/disk issues with 2.5D mesh generation, as the DEM generation assumes georeferenced units in order to not explode.

Edit: same for the spacing estimation in odm_filterpoints.

[nbnielsen]

I have a "small" dataset that I can probably share, that has this error showing. It is just under 2000 images

[VPA-Dave]

I have a dataset where I captured 4 x Rugby fields two ways - once just flying the programmed route and once with RTK. Both sets have been processed with GCPs and results are pretty close. I can provide all or some of the data. just let me know what you want and how much of it.

[smathermather]

Thanks. Is this a dataset or set of datasets that exhibit a behavior relevant to the discussion? If do, the smallest possible subset that does so would be useful.

[jasonowsley]

I'b be happy to create a set of data while im out in the field. Just let me know how many images, GSD, area etc you would like.

[jasonowsley]

I'b be happy to create a set of data while im out in the field. Just let me know how many images, GSD, area etc you would like.

[pierotofy]

Ideally, the smallest dataset you can get (ideally < 50-100 images) spanning the largest area possible (thus highest altitude legally allowed), with a rectangular/square coverage (again ideally, no corridor patterns).

Also 5-8 GCPs.

[jasonowsley]

Here is a link to a dataset with the gcp file

https://u.pcloud.link/publink/show?code=kZUwVu5Z7McgcJKgvCzGvbxAlMoFfSRVdaE7

[pierotofy]

Nice! I noticed the GCP file doesn't show the association to the pictures.

EPSG:4326
GCP1 -104.007793951133 32.1224660000333 903.75902
GCP2 -104.0077919565 32.1219366718333 903.70224
...

It's missing the im_x im_y image_name fields? https://docs.opendronemap.org/gcp/#gcp-file-format

[jasonowsley]

Yeah, I just made a reference file to input into WebODM’s georeferencer. Figured you'd want to align the GCPs in your workflow. I’d be happy to map it completely if you like.

[pierotofy]

I'm just unsure which point would correspond to which marker in the images, but any help would be appreciated!

[jasonowsley]

Ok ive added a fully referenced GCP file. one in the root and one in thephoto directory.

[pierotofy]

Just finished processing it; I'd say this is an excellent dataset with ground control points, but the area covered is too small, and thus it's difficult to highlight the issue at hand, which is accuracy over large areas. The GCPs align quite well here.

https://cloud1.webodm.net/public/task/299c48c9-ed84-4e5d-b36b-081268c108e7/map/?t=orthophoto

[jasonowsley]

I will get you a better set. Is it ok if its a large rectangular area? I can fly a tract possibly

[smathermather]

Do we have a sense of what kind of area we could start to measure the error in? If it's meters error in km distances, do we need a km wide dataset?

[originlake]

It also depends on the location of the survey region. The further away from the UTM zone center, the distortion gets bigger, so the requirement for the area would be smaller, vice versa. Could probably estimate the error expectations by modifying the script I created earlier with the actual center of the survey and distance. For example, use gcp1 as the center and define a 1kmx1km square plane at height 1km, the expected max error (should be one of the corners, which is around 700m away from the center) is 0.47m
https://colab.research.google.com/drive/1J3VBcmLbZ4YqVTrNXSisbjDr41E75ek3#scrollTo=lEMSfsirnm8T

[smathermather]

Ah nice. Couldn't we then force the use of an adjacent zone to simulate the error over a much smaller area and then calculate the expected impact of improvements? It's a bit synthetic, but arbitrarily large datasets with good GCPs are hard to come by. And just because they're hard to come by doesn't mean we shouldn't support, but it does make it harder to test for, and I wonder if synthesizing the problem by increasing it's magnitude by using an adjacent, otherwise inappropriate UTM zone could help us validate an improved approach.

[jasonowsley]

I just want to make sure I’m understanding this correctly. At the height and coverage area needed, it would take about 600 photos to cover a square kilometer, but the requirement is for 50 to 100 photos. I’m not sure I’ll be able to meet that with a drone as requested

[pierotofy]

but the requirement is for 50 to 100 photos

The photo count is an ideal type of requirement; if necessary, a larger dataset can be captured too. It's just much more difficult (and slow) to work with larger datasets.

[jasonowsley]

ok ill see what i can come up with this week

[jasonowsley]

I apologize if this was already answered somewhere in the thread, but I wanted to clarify something. Is this issue specific to using GCPs on autonomous flights, or does it also apply when using GCPs on RTK-enabled flights? I usually rely on surveyed coordinates and compare elevations using identical XY references to my models, and I haven’t really verified horizontal accuracy beyond visually checking alignment with the orthophoto overlay. The reason I’m asking is that with such a large area, setting a lot of GCPs would be time-consuming, and I’m wondering if combining RTK with fewer GCPs would still yield good results. Or would the RTK positioning essentially skew the experiment if the goal is to assess accuracy from GCPs alone?

[pierotofy]

The number of GCPs should still be in the 5-8 range, regardless of area. This problem affects all types of reconstructions, but for development/testing purposes, the GCP-only case is the one of interest because the problem is more detectable there (so don't enable RTK).

[peter5111]

Are we still in need of small datasets with GCP's ?

As I understand it you need a (roughly) square area with 5 to 8 GCP's.
Missions flown without RTK corrections to the imagery.
RGB only at this stage (MS not required yet)

What sort of area do you need ?

I can capture M3M and P4M imagery if it would help to have different datasets of the same GCP's ?

[pierotofy]

What sort of area do you need ?

Best would be a flat area, near the edge of a UTM zone, covering a large area.

[peter5111]

What sort of area do you need ?

Best would be a flat area, near the edge of a UTM zone, covering a large area.

Define 'large' ?
1 Ha, 10 Ha, 100 Ha?

Define 'Edge of a UTM zone' ?
Straddling, 1Km, 10Km, 100Km ?
I'm in 55G if that helps.

[pierotofy]

Edge of UTM zone => Further away from the UTM zone center.

I wish I could provide exact definitions for area, but at this point I'm not really sure. Perhaps @originlake can share more details about the dataset that triggered this issue.

[pierotofy]

I've found a dataset I can work with, so I will be testing the alternative approach in the upcoming days.

[smathermather]

Ah. Fantastic.

[YanNoun]

Hey,

Happy to have a look at this once we have a dataset. How I would solve it, is with bilinear interpolation that approximate the projection :

• Compute the extent on the data to get the bounds of the reconstruction topo frame.
• Sampling a grid of these bounds (grid of topo frames points)
• Transforms each grid's point using pyproj : this is the projected grid we'll use for bilinear interpolation.

Now, for transforming any point, we get its 4 neighbors in the projected grid and do bilinear-interpolate to get its projected value.

Or, for some points like the GCPs, we transform them individually, instead of using the linear (or the grid above). This could a first easy thing to try.

Cheers,

Yann

[alanterra]

Hi

I can add a dataset (actually 3 datasets) if helpful. I surveyed an area about about 500 m x 700 m with a GSD of about 2.8 cm (altitude 100 m) with a DJI Air 2S. Flight plans were 80%/80% overlap, azimuths 90°/270° and gimbal angle 5°, and 80%/65% overlap and gimbal angle 25°, giving ~1130 photos. I actually did these flights 3 times over various days with the same GCP markers. The GCPs were mapped with a BadElf Flex with RTK against a CORS site about 30 km distant, with a displayed error of ≤ 2 cm (plus any user error, which I tried to minimize).
Here are some images giving an idea of the scale of mismatch of GCPs. The GCP markers are 50 cm on a side.

[YanNoun]

Hey there,

I'm trying to improves the GCP handling these days and having access to a lot of datasets with RTK GCP (and GPS in image) is absolutely crucial to guide the improvments. And as of now, the data supply is quite scarce on that front.

So, yes, we definitely would love to have access to these datasets.

Cheers,

Yann

[alanterra]

Hi Yann
Please send an email to me at [email protected] and I'll reply with a Dropbox link to the photos (three surveys of the same place, each about 1150 photos). The photos are 43 GB and I am also including the WebODM results for each survey, which add another 59 GB. If this is too much data, I can trim it. Just let me know.