Openmvg: How to Supply Already Known GCPs

[mchildsCO76]

How can I provide my own list of GCPs for a set of images? In this case a GCP will be a list of control points with real-world EXIF XYZ location and associated pixel locations on any images that control point was marked on.

I know there is an existing UI for this (http://openmvg.readthedocs.io/en/latest/software/ui/SfM/control_points_registration/GCP/), but I need to provide these myself somewhere in the pipeline (they will all be selected before any part of the pipeline runs). From what I can discern it seems like they are added as 'Landmarks' in the sfm_data.bin file, but I wasn't sure if there is an existing .exe for providing them as text or some other way to add to the .bin file.

[pmoulon]

There is no way to provide them yet by a file but feel free to contribute it so anyone can use it ;-)

It should be trivial to implement (I can provide guidance):

• Load a SfM_data file.
• Add the GCP Landmarks and observation from your file
  
  
  
  • https://github.com/openMVG/openMVG/blob/master/src/openMVG/sfm/sfm_data.hpp#L43
  
  
• Save the updated sfm_data

[mchildsCO76]

Thanks I will see if I can add something like this. Where in the pipeline would I provide these control points? I would expect to have the control points locations on the images at the start of the entire process. I don't currently have a sfm_data.bin file until after running Incremental SfM.

Would the Landmarks be used during SfM or are they something the openMVG_main_geodesy_registration_to_gps_position.exe would use?

[pmoulon]

You can use or configure the control points at any point in the pipeline.

For example you can use the sfm_data.json with the control point GUI app once the project is initialized (You don't need any 3D data to add them).

For the moment the control points constraints are only enabled for the BA for registration purpose:
see here: https://github.com/openMVG/openMVG/blob/6ce820b2844a98941a5d2543050f4e278a323800/src/software/ui/SfM/control_points_registration/mainLayout.cpp#L341

In the other places by default the default configuration will not use the control points constraint (but having the control points in the sfm_data will not be a problem) https://github.com/openMVG/openMVG/search?utf8=%E2%9C%93&q=Control_Point_Parameter&type=

If you want activate the control_points and the GPS data in openMVG_main_geodesy_registration_to_gps_position it's doable, just set the boolean that activate the constraint to true.

[mchildsCO76]

Ok it sounds like I want the provide the control points before IncrementalSfM (or GlobalSfM) is called. Does the JSON format for the SfM data support the Landmarks or must it be in the .bin format?

Do you know if there is a pre-build Windows binary for the control points UI app anywhere? It's sometimes a bit of a chore getting things building on VS2015 so if someone already has that piece ready to run it might be easier for testing.

[pmoulon]

OpenMVG support various IO serialization for sfm_data (JSON, XML, BIN). So the control point can be used in any of this format. Just use the OpenMVG API and choose which output you want.

There is not pre-build binary yet. But note that to develop with Qt you just have to download their precompiled package and directly use it.

Else you can build it in a Linux VM quite easily and use the GUI there.

[mchildsCO76]

I'll just toss some in to a test app in the OpenMVG .exe's that I'm already building without QT and write out the JSON there.

Thanks for the always quick responses!

[mchildsCO76]

I'm working on creating the Landmark (control points) in the sfm_data from my user-provided control points. The Observation class appears to be what holds the associated pixel location on an input image as the Landmark class holds a collection of these, as you would expect.

The Landmarks type stored by the Sfm_data structure is declared as follows:

/// Define a collection of landmarks are indexed by their TrackId
using Landmarks = Hash_Map;

I'm not sure what the TrackId that the Hash_Map is indexed by refers to. I can read in all of my control points with a defined ECEF XYZ location and a collection of the associated pixel locations on the input images, but I'm not sure how to create the map with the TrackId. Is that just an index that I increment so that each Landmark is uniquely identified?

[pmoulon]

You are right TrackId is just an ID in order to give a unique label to your GCP.
You can start from 0 and increment the Id for any new GCP you are adding.

[mchildsCO76]

Thanks I've got them added in to the SfM data structure in my branch now for testing.

Is it the case now that the control points are only used in the control points UI, so adding them to the sfm_data at the start of the pipeline wouldn't actually do anything? Earlier you said 'If you want activate the control_points and the GPS data in openMVG_main_geodesy_registration_to_gps_position it's doable, just set the boolean that activate the constraint to true.' I didn't see an obvious place to add that flag.

For the other places, I'm assuming where ControlPointParameter_t() is currently passed in I should change it to use something like ControlPointParameter_t( 20.0, true ) to enable their use?

[mchildsCO76]

Or I guess just change the default for the 2nd parameter to Control_Point_Parameter_t to true. I would think if control points are provided the user would want them used.

[mchildsCO76]

I have locally updated the OpenMVG apps to read in the CONTROL_POINTS from the SFM data file so they aren't lost as the data is passed through the pipeline if added at the beginning.

Once the Control_Point_Parameter_t flag is set to true the SfM is using the points, but was initially crashing in sfm_data_BA_ceres.cpp in the 'if (options.control_point_opt.bUse_control_points )' block of the Adjust function as the calls were assuming there were always non-empty mapped entries for every control point. The parameter constants also were not added.

I have updated the loop to look as follow:

if (options.control_point_opt.bUse_control_points)
  {
    // Use Ground Control Point:
    // - fixed 3D points with weighted observations
    for (auto & gcp_landmark_it : sfm_data.control_points)
    {
      const Observations & obs = gcp_landmark_it.second.obs;

      for (const auto & obs_it : obs)
      {
        // Build the residual block corresponding to the track observation:
        const View * view = sfm_data.views.at(obs_it.first).get();

        // Each Residual block takes a point and a camera as input and outputs a 2
        // dimensional residual. Internally, the cost function stores the observed
        // image location and compares the reprojection against the observation.
        ceres::CostFunction* cost_function =
          IntrinsicsToCostFunction(
            sfm_data.intrinsics[view->id_intrinsic].get(),
            obs_it.second.x,
            options.control_point_opt.weight);

        if (cost_function) {
            auto map_intrinsic_for_view = map_intrinsics.find( view->id_intrinsic );
            auto map_pose_for_view = map_poses.find(view->id_pose);
            if ((map_intrinsic_for_view == map_intrinsics.end()) || map_intrinsic_for_view->second.empty()) {
                std::cerr << "Missing map intrinsic for view intrinsic ID " << view->id_intrinsic << " for GCP id: " << gcp_landmark_it.first << std::endl;
            }
            else if ((map_pose_for_view == map_poses.end()) || map_pose_for_view->second.empty()) {
                std::cerr << "Missing map pose for view pose ID " << view->id_pose << " for GCP id: " << gcp_landmark_it.first << std::endl;
            }
            else {
                problem.AddResidualBlock(
                    cost_function,
                    nullptr,
                    &map_intrinsic_for_view->second[0],
                    &map_pose_for_view->second[0],
                    gcp_landmark_it.second.X.data());
            }
        }
      }
      if (obs.empty()) {
        std::cerr
          << "Cannot use this GCP id: " << gcp_landmark_it.first
          << ". There is not linked image observation." << std::endl;
      }
      else {
         // Set the 3D point as FIXED (it's a valid GCP)
         double* landmark_coords = gcp_landmark_it.second.X.data();
         problem.AddParameterBlock(landmark_coords, 3);
         problem.SetParameterBlockConstant(landmark_coords);
      }
    }
  }

The code now runs without crashing and when stepping through the values look ok, but the results are an empty point cloud. Should I even be enabling control points during the SfM process or should they only be used after the entire pipeline is run using the code chunk from the control point GUI at https://github.com/openMVG/openMVG/blob/6ce820b2844a98941a5d2543050f4e278a323800/src/software/ui/SfM/control_points_registration/mainLayout.cpp#L341 ?

If I just want to replicate what the control point UI is doing I'm not entirely sure where the right place to do that is? Should I perhaps just ignore the control points all the way until the openMVG_main_geodesy_registration_to_gps_position after the SfM has been performed?

Also is there a way to back out the XYZ location for an image pixel coordinate after the SfM process is complete? If so I could pull the SfM-calculated XYZ location for each control point and generate some info about how far off the control point location and the SfM-calculated location are.

[pmoulon]

The safer change consists in introducing your GCP once the SfM is performed (and perform the registration at the end).

I think the code was working. You hit a special case when you have GCP with no image measure.

I mean the code you added was required only in the case the cost_function is null else the parameter block will be added by the ceres API with the AddResidualBlock call.

double* landmark_coords = gcp_landmark_it.second.X.data();
problem.AddParameterBlock(landmark_coords, 3);

Note: The registration control point GUI is already computing the distance between the 3D points (SFM computation and control points).

[mchildsCO76]

I see now where the AddParameterBlock should really only be called if AddResidualBlock was called for any control point for the landmark and have adjusted that. In any case I don't think that code gets called if I don't use the control points during SfM.

Should I run that chunk of code from mainLayout.cpp on the sfm_data before or after main_geodesy_registration_to_gps_position (which I call right after Incremental or Global SfM)?

[mchildsCO76]

I'm getting good results running after the main_geodesy_registration_to_gps_position, so that seems to be the proper place.

[ywjia]

@mchildsCO76 Could you please summarize your current findings about how to add / where to add and how to use GCPs so that they can be applied in BA? Thanks.

[mchildsCO76]

I found good results applying GCPs (landmarks) after the main_geodesy_registration_to_gps_position. In my branch I just added the GCP stuff to that .exe after the adjustment using EXIF data since it seems to fit what that .exe is doing (and I was too lazy to create a whole new one).

I will also experiment with applying the GCPs before running SfM, but haven't tried that. I have local changes that allow passing a parameter to control whether GCPs and EXIF adjustment is done or just one or the other to allow easily trying different ordering.

[ywjia]

@mchildsCO76 Thanks for the summary.
I saw the the following code in "sfm_data_BA_ceres.cpp" (in the function "bool Bundle_Adjustment_Ceres::Adjust"). I thought that was for using GCPs for BA. Have you tested this part of the code, please? Thanks again.

if (options.control_point_opt.bUse_control_points)
{
...
...
}

[mchildsCO76]

You can see my version at https://github.com/mchildsCO76/openMVG/blob/develop_local_origin/src/software/Geodesy/registration_to_exif_gps_position.cpp . My registerProjectUsingControlPoints function has a parameter to control whether or not BA adjustment is done as well. I'm always passing true for now.

[pmoulon]

Have you tested this part of the code, please? Thanks again.

This part is tested in the unit test for BundleAdjusment https://github.com/openMVG/openMVG/blob/master/src/openMVG/sfm/sfm_data_BA_test.cpp#L188 and it is working.

[mchildsCO76]

The code I added is at a higher level, but I believe it gets down that far. It is working well in any case. I have a new function in registration_to_exif_gps_position.cpp that basically copies your code from the control point UI application and does the same thing. It looks like:

static bool
registerProjectUsingControlPoints
    (
    SfM_Data&   sfm_data,
    bool        bDoBundleAdjustment
    )
{
    // Make sure we have enough data
    if (sfm_data.control_points.size() < 3)
    {
        if (!sfm_data.control_points.empty()) 
        {
            std::cerr << "At least 3 control points are required to apply the adjustment." << std::endl;
        }
        return (false);
    }

    // Assert that control points can be triangulated
    for ( const auto& control_point : sfm_data.control_points )
    {
        if ( control_point.second.obs.size() < 2)
        {
            std::cerr << "Each control point must be defined in at least 2 pictures." << std::endl;
            return (false);
        }
    }

    //---
    // registration (coarse):
    // - compute the 3D points corresponding to the control point observation for the SfM scene
    // - compute a coarse registration between the controls points & the triangulated point
    // - transform the scene according the found transformation
    //---
    std::map<IndexT, Vec3> vec_control_points, vec_triangulated;
    std::map<IndexT, double> vec_triangulation_errors;
    for (auto iterCP = sfm_data.control_points.begin(); iterCP != sfm_data.control_points.end(); ++iterCP)
    {
        Landmark & landmark = iterCP->second;

        //Triangulate the point:
        Triangulation trianObj;
        const Observations & obs = landmark.obs;
        for (auto itObs = obs.begin(); itObs != obs.end(); ++itObs)
        {
            const View * view = sfm_data.views.at(itObs->first).get();
            if (!sfm_data.IsPoseAndIntrinsicDefined(view))
                continue;
            const openMVG::cameras::IntrinsicBase * cam = sfm_data.GetIntrinsics().at(view->id_intrinsic).get();
            const openMVG::geometry::Pose3 pose = sfm_data.GetPoseOrDie(view);
            trianObj.add(
                cam->get_projective_equivalent(pose),
                cam->get_ud_pixel(itObs->second.x));
        }

        // Compute the 3D point
        const Vec3 X = trianObj.compute();
        if (trianObj.minDepth() > 0) // Keep the point only if it have a positive depth
        {
            vec_triangulated[iterCP->first] = X;
            vec_control_points[iterCP->first] = landmark.X;
            vec_triangulation_errors[iterCP->first] = trianObj.error() / (double)trianObj.size();
        }
        else
        {
            std::cerr << "Invalid triangulation" << std::endl;
            return (false);
        }
    }

    if (vec_control_points.size() < 3)
    {
        std::cerr << "Insufficient number of triangulated control points." << std::endl;
        return (false);
    }

    // compute the similarity
    {
        // data conversion to appropriate container
        Mat x1(3, vec_control_points.size()),
            x2(3, vec_control_points.size());
        for (size_t i = 0; i < vec_control_points.size(); ++i)
        {
            x1.col(i) = vec_triangulated[i];
            x2.col(i) = vec_control_points[i];
        }

        std::cout
            << "Control points observation triangulations:\n"
            << x1 << std::endl << std::endl
            << "Control points coords:\n"
            << x2 << std::endl << std::endl;

        Vec3 t;
        Mat3 R;
        double S;
        if (openMVG::geometry::FindRTS(x1, x2, &S, &t, &R))
        {
            openMVG::geometry::Refine_RTS(x1, x2, &S, &t, &R);
            std::cout << "Found transform:\n"
                << " scale: " << S << "\n"
                << " rotation:\n" << R << "\n"
                << " translation: " << t.transpose() << std::endl;


            //--
            // Apply the found transformation as a 3D Similarity transformation matrix // S * R * X + t
            //--

            const openMVG::geometry::Similarity3 sim(geometry::Pose3(R, -R.transpose() * t / S), S);
            openMVG::sfm::ApplySimilarity(sim, sfm_data);

            // Display some statistics:
            std::stringstream os;
            for (auto iterL = sfm_data.control_points.begin(); iterL != sfm_data.control_points.end(); ++iterL)
            {
                const IndexT CPIndex = iterL->first;
                // If the control point has not been used, continue...
                if (vec_triangulation_errors.find(CPIndex) == vec_triangulation_errors.end())
                    continue;

                os
                    << "CP index: " << CPIndex << "\n"
                    << "CP triangulation error: " << vec_triangulation_errors[CPIndex] << " pixel(s)\n"
                    << "CP registration error: "
                    << (sim(vec_triangulated[CPIndex]) - vec_control_points[CPIndex]).norm() << " user unit(s)" << "\n\n";
            }
            std::cout << os.str();
        }
        else
        {
            std::cerr << "Registration failed. Please check your Control Points coordinates." << std::endl;
            return (false);
        }
    }

    // Bundle adjustment with GCP
    if ( bDoBundleAdjustment )
    {
        using namespace openMVG::sfm;
        std::cout << "Performing GCP bundle adjustment using " << sfm_data.control_points.size() << " control points." << std::endl;
        Bundle_Adjustment_Ceres::BA_Ceres_options options;
        Bundle_Adjustment_Ceres bundle_adjustment_obj(options);
        Control_Point_Parameter control_point_opt(20.0, true);
        bool bControlPointsUsed = bundle_adjustment_obj.Adjust(
            sfm_data,
            Optimize_Options
            (
                cameras::Intrinsic_Parameter_Type::NONE, // Keep intrinsic constant
                Extrinsic_Parameter_Type::ADJUST_ALL, // Adjust camera motion
                Structure_Parameter_Type::ADJUST_ALL, // Adjust structure
                control_point_opt // Use GCP and weight more their observation residuals
            )
        );
        if (!bControlPointsUsed)
        {
            std::cerr << "Error performing GCP bundle adjustment" << std::endl;
            return (false);
        }
    }

    return (true);
}

[pmoulon]

So if I understand well you would like to have a new binary that performs registration to control_points data...
So at the end we would have:
registration_to_exif_gps_position.cpp
and a new binary registration_to_control_points_position.cpp

We could also imagine to extend the registration binary to support those mode:

• GPS only,
• CONTROL_POINTS only,
• GPS + CONTROL_POINTS.

We can add RIGID and BA registration mode too.

[mchildsCO76]

I've actually already updated the existing registration_to_exif_gps_position binary in my branch to do all that and have options to control whether using GPS, control points, or both. I haven't pushed up my local changes to add the command-line options for the switches, but I'll do that.

My version is at https://github.com/mchildsCO76/openMVG/blob/develop_local_origin/src/software/Geodesy/registration_to_exif_gps_position.cpp .