Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Align by DGPS positions in images in SfmTransform - test with DJI drone images for example #1523

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from

Conversation

belveder79
Copy link

@belveder79 belveder79 commented Aug 30, 2023

Description

The changes were made to the SfmTransform executable mainly, plus building the dependency GeographicLib from public github repo, and having the feature ON/OFF in the main cmake configuration depending on GeographicLib being available (either built from the dependencies) or the system. The mode aligns the reconstruction using a Ceres implementation and a Huber Loss Function to the x/y/z coordinates given by the GPS coordinates (in metric UTM). The reconstruction ends up in a local coordinate system. The global center coordinates are stored in a separate file.

Features list

  • Adds a switchable dependency on GeographicLib which determines if the feature is available or not.
  • SfmTransform can now use an additional method gps2utm to align to good GPS positions in images. Good means for example from DJI drones, which are +/- 20cm off, or DGPS measurements work as well. Having all 3, latitude, longitude and altitude in the Exiv-Tags of the images is mandatory.
  • The alignment is based on a transform to an UTM coordinate system. The GeographicLib automatically returns the correct system when the latitude/longitude are provided. The mean of right/up/sky is calculated and removed to put the reconstruction into a local metric coordinate system again. The mean and UTM EPSG zone is stored into a new file
    localcenter.json in JSON format, which can be used later for global alignment.
  • The SfmTransform is supposed to plug in between the StructureFromMotion and DepthMapPrepare Node in Meshroom. It estimates scale, translation and rotation. Matlab code for verification is available too if needed.
  • For whatever reason the alignment is correct, but the "final" alignment of a textured reconstruction (i.e. at the end of the full Photogrammetry pipeline) has flipped axes. It is still right-handed, but Y/Z are inverted. Therefore the SfmTransform node applies a second transform right after the alignment to invert Y/Z. That gives the correctly textured 3D model which is in a local coordinate system, but globally aligned with correct right/up/sky coordinates. This can be verified in Meshlab (+z is sky, +y is up, i.e. towards north pole on the northern hemisphere, and +x is right).

Implementation remarks

Depending on the real use case, one might want to make additional modifications like (not implemented):

  • [] the final transform (inversion of Y/Z) optional with a flag (not implemented right now). If the goal is to do some localization, because this might be broken by inverting the Y/Z coordinates in the sense that the localization result is not the real UTM coordinate. The goal was to have the textured model aligned correctly, not the sfm.abc file at the end of the SfmTransform node. The best might be to have 2 files stored, one for making the textured model "correct" and one to have the correct alignment for localization.
  • [] passing a variable to the executable to store the local UTM center in a file named by the user
  • [] make the algorithm robust wrt. data spanning 2 zones, which might be rarely the case, however.
  • [] alignment chooses poses only which have valid GPS Exiv tags. Probably testing with missing tags would be an option, but this is really artificial if one uses DJI drone footage.
  • [] adding initial parameters to the algorithm could help in alignment (i.e. if an initial camera should be chosen as for the coordinate system instead of the mean of all camera coordinates). Again, this is merely optimization. One could think of more stuff.

The alignment was tested on 2 sets of images from DJI Mavic drone footage and shown to work very well.
kapelle00
kapelle01

Copy link
Member

@simogasp simogasp left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks for your contribution!
I left some comments and we really should clarify this axis inversion thing

set(ALICEVISION_HAVE_GEOGRAPHIC 0)

if(NOT ALICEVISION_USE_GEOGRAPHIC STREQUAL "OFF")
find_library(GeographicLib_LIBRARY NAMES Geographic GeographicLib DOC "GeographicLib library")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
find_library(GeographicLib_LIBRARY NAMES Geographic GeographicLib DOC "GeographicLib library")
find_package(GeographicLib CONFIG)

then the target GeographicLib::GeographicLib_SHARED can be used to express the dependency in the remaining cmake parts

@@ -194,6 +193,7 @@ alicevision_add_software(aliceVision_sfmTransform
aliceVision_sfmData
aliceVision_sfmDataIO
Boost::program_options
${GeographicLib_LIBRARY}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

GeographicLib_LIBRARY can be set to GeographicLib::GeographicLib_SHAREDif found or to empty string.

using ceres::Solve;
using ceres::Solver;

typedef Eigen::Matrix<double,3,4> Matrix34d;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
typedef Eigen::Matrix<double,3,4> Matrix34d;
using Matrix34d = Eigen::Matrix<double,3,4>;

GeographicLib::UTMUPS::Forward(lat, lon, zone, northp, x, y, gamma, k);
mean += Eigen::Vector3d(x, y, alt);

coords.push_back(Eigen::Vector3d(x, y, alt));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
coords.push_back(Eigen::Vector3d(x, y, alt));
coords.emplace_back(x, y, alt);

Comment on lines +284 to +293
{
FILE* localcoordfile = fopen("localcenter.json", "wt");
fprintf(localcoordfile,"{ \"Coords\": {\n");
fprintf(localcoordfile,"\"Center\": [ %.16f, %.16f, %.16f ],\n", mean(0), mean(1), mean(2));
fprintf(localcoordfile,"\"Zone\": %d,\n",zone);
fprintf(localcoordfile,"\"North\": %s,\n",northp ? "true" : "false");
fprintf(localcoordfile,"\"EPSGCode\": %d\n",northp ? 32600+zone : 32700+zone);
fprintf(localcoordfile,"}\n}");
fclose(localcoordfile);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it would be better to write a function that does that using boost.json

t = Tnew;
S = scale;

delete[] pose_for_alignment;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there are some other stuff to delete

delete[] Rn;
delete loss_function;

{
bool debug = false;

std::vector<Eigen::Vector3d> coords;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this seems to be filled but never used in any other part of the code, maybe remove it?

<< "No pose for view " << v.second->getViewId() << std::endl);
}
}
mean /= numValidPoses;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should check whether numValidPoses > 0 before otherwise return with a message/error as there are no valid poses

Comment on lines +254 to +282
std::vector<Eigen::Vector3d> coords;

Eigen::Vector3d mean = Eigen::Vector3d::Zero();
int numValidPoses = 0;
const sfmData::Poses poses = sfmData.getPoses();
int zone; bool northp;
for (auto v : sfmData.getViews()) {

IndexT poseId = v.second->getPoseId();
if(poses.find(poseId) != poses.end())
{
double lat; double lon; double alt;
v.second->getGpsPositionWGS84FromMetadata(lat, lon, alt);
// zone and northp should be returned!!!
double x, y, gamma, k;
GeographicLib::UTMUPS::Forward(lat, lon, zone, northp, x, y, gamma, k);
mean += Eigen::Vector3d(x, y, alt);

coords.push_back(Eigen::Vector3d(x, y, alt));

numValidPoses++;
}
else
{
ALICEVISION_LOG_INFO(std::setprecision(17)
<< "No pose for view " << v.second->getViewId() << std::endl);
}
}
mean /= numValidPoses;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would refactor this whole section into a separate function that computes the mean position and the number of valid poses, it would make the code more readable

Comment on lines +767 to +779

#if ALICEVISION_IS_DEFINED(ALICEVISION_HAVE_GEOGRAPHIC)
if (alignmentMethod == EAlignmentMethod::FROM_GPS2UTM)
{
ALICEVISION_LOG_INFO("Applying additional inversion to Y and Z to make textured mesh end up aligned correctly!");
Eigen::Matrix3d invYZ = Eigen::Matrix3d::Identity();
invYZ(1,1) = -1; invYZ(2,2) = -1;
Eigen::Vector3d zeroT = Eigen::Vector3d::Zero();

sfm::applyTransform(sfmData, 1.0, invYZ, zeroT);
}
#endif

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know why there is such inversion. This smells a lot like the infamous problem of conversion between the opengl reference frame and the computer vision reference frame.
@servantftechnicolor Maybe I'm not recalling well, but wasn't there a change in the way the obj are saved that has been introduced lately? Like to have the correct visualization there was this kind of inversion applied to all saved objs?

double * cPtr = centers;
for (auto v : sfmData.getViews())
{
IndexT poseId = v.second->getPoseId();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
IndexT poseId = v.second->getPoseId();
const IndexT poseId = v.second->getPoseId();

@belveder79
Copy link
Author

thx for the comments. I'm going to review your suggestions and update the request accordingly.

Yes, it has certainly something to do with the coordinate system convention used. It does not convert left-right handed or something, it is a pure 180° rotation around the x-axis. The inversion essentially makes z point to the sky and x/y correspond to right/up in the UTM coordinate system (which is still right-handed). From that point on it is also correct metric scale if it was not before.

Admittedly, the rotation actually makes the visualization in Meshroom look a bit odd, which (I believe) assumes z to be forward and y point down. In a MeshRoom pipeline, the SfmTransform would plug in right after the StructureFromMotion node and before any (I don't recall the exact names) PrepareDense node. That makes the entire textured reconstruction at the very end of the pipeline correct.

@fabiencastan
Copy link
Member

Hi @belveder79,
Thanks for your contribution. It's really interesting.
Would you have time to finalize it?
The best would be to rebase on develop to fix the conflicts (to avoid multiple merges with develop that complicate the git history) and do the updates from the review.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants