<commit_before>// __BEGIN_LICENSE__
//  Copyright (c) 2009-2013, United States Government as represented by the
//  Administrator of the National Aeronautics and Space Administration. All
//  rights reserved.
//
//  The NGT platform is licensed under the Apache License, Version 2.0 (the
//  "License"); you may not use this file except in compliance with the
//  License. You may obtain a copy of the License at
//  http://www.apache.org/licenses/LICENSE-2.0
//
//  Unless required by applicable law or agreed to in writing, software
//  distributed under the License is distributed on an "AS IS" BASIS,
//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//  See the License for the specific language governing permissions and
//  limitations under the License.
// __END_LICENSE__

// This tool uses libpointmatcher for alignment,
// https://github.com/ethz-asl/libpointmatcher
// Copyright (c) 2010--2012,
// <NAME> and <NAME>, ASL, ETHZ, Switzerland
// You can contact the authors at <f dot pomerleau at gmail dot com> and
// <stephane at magnenat dot net>

// All rights reserved.

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above copyright
//       notice, this list of conditions and the following disclaimer in the
//       documentation and/or other materials provided with the distribution.
//     * Neither the name of the <organization> nor the
//       names of its contributors may be used to endorse or promote products
//       derived from this software without specific prior written permission.

// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL ETH-ASL BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#include <vw/Core/Stopwatch.h>
#include <vw/FileIO/DiskImageView.h>
#include <vw/Math.h>
#include <vw/Image.h>
#include <vw/Cartography/Datum.h>
#include <vw/Cartography/GeoReference.h>
#include <vw/Cartography/PointImageManipulation.h>
#include <vw/FileIO/DiskImageUtils.h>
#include <vw/Core/CmdUtils.h>
#include <asp/Core/Common.h>
#include <asp/Core/Macros.h>
#include <asp/Core/PointUtils.h>
#include <asp/Core/InterestPointMatching.h>
#include <liblas/liblas.hpp>

#include <limits>
#include <cstring>

#include <pointmatcher/PointMatcher.h>
#include <ceres/ceres.h>
#include <ceres/loss_function.h>

#include <asp/Tools/pc_align_utils.h>

#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
namespace fs = boost::filesystem;
namespace po = boost::program_options;

using namespace vw;
using namespace std;
using namespace vw::cartography;
using namespace asp;

typedef double RealT; // We will use doubles in libpointmatcher.
typedef PointMatcher<RealT> PM;
typedef PM::DataPoints DP;
using namespace PointMatcherSupport;

const double BIG_NUMBER = 1e+300; // libpointmatcher does not like here the largest double

/// Options container for the pc_align tool
struct Options : public vw::cartography::GdalWriteOptions {
  // Input
  string reference, source, init_transform_file, alignment_method, config_file,
    datum, csv_format_str, csv_proj4_str, match_file, hillshade_options,
    ipfind_options, ipmatch_options;
  PointMatcher<RealT>::Matrix init_transform;
  int    num_iter,
         max_num_reference_points,
         max_num_source_points;
  double diff_translation_err,
         diff_rotation_err,
         max_disp,
         outlier_ratio,
         semi_major,
         semi_minor;
  bool   compute_translation_only, use_hillshading,
         dont_use_dem_distances,
         save_trans_source,
         save_trans_ref,
         highest_accuracy,
         verbose;
  std::string initial_ned_translation;
  
  // Output
  string out_prefix;

  Options() : max_disp(-1.0), verbose(true){}

  /// Return true if the reference file is a DEM file and this option is not disabled
  bool use_dem_distances() const { return ( (asp::get_cloud_type(this->reference) == "DEM") && !dont_use_dem_distances); }
};


void handle_arguments( int argc, char *argv[], Options& opt ) {
  po::options_description general_options("");
  general_options.add_options()
    ("initial-transform",        po::value(&opt.init_transform_file)->default_value(""),
                                 "The file containing the transform to be used as an initial guess. It can come from a previous run of the tool.")
    ("num-iterations",           po::value(&opt.num_iter)->default_value(1000),
                                 "Maximum number of iterations.")
    ("diff-rotation-error",      po::value(&opt.diff_rotation_err)->default_value(1e-8),
                                 "Change in rotation amount below which the algorithm will stop (if translation error is also below bound), in degrees.")
    ("diff-translation-error",   po::value(&opt.diff_translation_err)->default_value(1e-3),
                                 "Change in translation amount below which the algorithm will stop (if rotation error is also below bound), in meters.")
    ("max-displacement",         po::value(&opt.max_disp)->default_value(0.0),
                                 "Maximum expected displacement of source points as result of alignment, in meters (after the initial guess transform is applied to the source points). Used for removing gross outliers in the source point cloud.")
    ("outlier-ratio",            po::value(&opt.outlier_ratio)->default_value(0.75),
                                 "Fraction of source (movable) points considered inliers (after gross outliers further than max-displacement from reference points are removed).")
    ("max-num-reference-points", po::value(&opt.max_num_reference_points)->default_value(100000000),
                                 "Maximum number of (randomly picked) reference points to use.")
    ("max-num-source-points",    po::value(&opt.max_num_source_points)->default_value(100000),
                                 "Maximum number of (randomly picked) source points to use (after discarding gross outliers).")
    ("alignment-method",         po::value(&opt.alignment_method)->default_value("point-to-plane"),
                                 "The type of iterative closest point method to use. [point-to-plane, point-to-point, similarity-point-to-point, least-squares, similarity-least-squares]")
    ("highest-accuracy",         po::bool_switch(&opt.highest_accuracy)->default_value(false)->implicit_value(true),
                                 "Compute with highest accuracy for point-to-plane (can be much slower).")
    ("csv-format",               po::value(&opt.csv_format_str)->default_value(""), asp::csv_opt_caption().c_str())
    ("csv-proj4",                po::value(&opt.csv_proj4_str)->default_value(""),
                                 "The PROJ.4 string to use to interpret the entries in input CSV files.")
    ("datum",                    po::value(&opt.datum)->default_value(""),
                                 "Use this datum for CSV files instead of auto-detecting it. Options: WGS_1984, D_MOON (1,737,400 meters), D_MARS (3,396,190 meters), MOLA (3,396,000 meters), NAD83, WGS72, and NAD27. Also accepted: Earth (=WGS_1984), Mars (=D_MARS), Moon (=D_MOON).")
    ("semi-major-axis",          po::value(&opt.semi_major)->default_value(0),
                                 "Explicitly set the datum semi-major axis in meters.")
    ("semi-minor-axis",          po::value(&opt.semi_minor)->default_value(0),
                                 "Explicitly set the datum semi-minor axis in meters.")
    ("output-prefix,o",          po::value(&opt.out_prefix)->default_value("run/run"),
                                 "Specify the output prefix.")
    ("compute-translation-only", po::bool_switch(&opt.compute_translation_only)->default_value(false)->implicit_value(true),
                                 "Compute the transform from source to reference point cloud as a translation only (no rotation).")
    ("save-transformed-source-points", po::bool_switch(&opt.save_trans_source)->default_value(false)->implicit_value(true),
                                  "Apply the obtained transform to the source points so they match the reference points and save them.")
    ("save-inv-transformed-reference-points", po::bool_switch(&opt.save_trans_ref)->default_value(false)->implicit_value(true),
     "Apply the inverse of the obtained transform to the reference points so they match the source points and save them.")

    ("initial-ned-translation", po::value(&opt.initial_ned_translation)->default_value(""),
                                 "Initialize the alignment transform based on a translation with this vector in the local North-East-Down coordinate system. Specify it in quotes, separated by spaces or commas.")

    ("initial-transform-from-hillshading", po::bool_switch(&opt.use_hillshading)->default_value(false)->implicit_value(true), "If both input clouds are DEMs, find interest point matches among their hillshaded versions, and use them to compute an initial rotation + translation + scale transform to apply to the source cloud before proceeding with alignment.")
    ("hillshade-options", po::value(&opt.hillshade_options)->default_value("--azimuth 300 --elevation 20"), "Options to pass to the hillshade program when computing the transform from hillshading.")
    ("ipfind-options", po::value(&opt.ipfind_options)->default_value("--ip-per-image 500000 --interest-operator sift --descriptor-generator sift"), "Options to pass to the ipfind program when computing the transform from hillshading.")
    ("ipmatch-options", po::value(&opt.ipmatch_options)->default_value("--inlier-threshold 100 --ransac-constraint similarity"), "Options to pass to the ipmatch program when computing the transform from hillshading.")
    ("match-file", po::value(&opt.match_file)->default_value(""), "Compute a translation + rotation + scale transform from the source to the reference point cloud using manually selected point correspondences from the reference to the source (obtained for example using stereo_gui).")

    ("no-dem-distances",         po::bool_switch(&opt.dont_use_dem_distances)->default_value(false)->implicit_value(true),
                                 "For reference point clouds that are DEMs, don't take advantage of the fact that it is possible to interpolate into this DEM when finding the closest distance to it from a point in the source cloud and hence the error metrics.")

    ("config-file",              po::value(&opt.config_file)->default_value(""),
     "This is an advanced option. Read the alignment parameters from a configuration file, in the format expected by libpointmatcher, over-riding the command-line options.");

  //("verbose", po::bool_switch(&opt.verbose)->default_value(false)->implicit_value(true),
  // "Print debug information");

  general_options.add( vw::cartography::GdalWriteOptionsDescription(opt) );

  po::options_description positional("");
  positional.add_options()
    ("reference", po::value(&opt.reference), "The reference (fixed) point cloud/DEM.")
    ("source",    po::value(&opt.source),    "The source (movable) point cloud/DEM.");

  po::positional_options_description positional_desc;
  positional_desc.add("reference", 1);
  positional_desc.add("source",    1);

  string usage("--max-displacement arg [other options] <reference cloud> <source cloud> -o <output prefix>");<commit_msg>pc_align: A small tweak<commit_after>// __BEGIN_LICENSE__
//  Copyright (c) 2009-2013, United States Government as represented by the
//  Administrator of the National Aeronautics and Space Administration. All
//  rights reserved.
//
//  The NGT platform is licensed under the Apache License, Version 2.0 (the
//  "License"); you may not use this file except in compliance with the
//  License. You may obtain a copy of the License at
//  http://www.apache.org/licenses/LICENSE-2.0
//
//  Unless required by applicable law or agreed to in writing, software
//  distributed under the License is distributed on an "AS IS" BASIS,
//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//  See the License for the specific language governing permissions and
//  limitations under the License.
// __END_LICENSE__

// This tool uses libpointmatcher for alignment,
// https://github.com/ethz-asl/libpointmatcher
// Copyright (c) 2010--2012,
// <NAME> and <NAME>, ASL, ETHZ, Switzerland
// You can contact the authors at <f dot pomerleau at gmail dot com> and
// <stephane at magnenat dot net>

// All rights reserved.

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above copyright
//       notice, this list of conditions and the following disclaimer in the
//       documentation and/or other materials provided with the distribution.
//     * Neither the name of the <organization> nor the
//       names of its contributors may be used to endorse or promote products
//       derived from this software without specific prior written permission.

// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL ETH-ASL BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#include <vw/Core/Stopwatch.h>
#include <vw/FileIO/DiskImageView.h>
#include <vw/Math.h>
#include <vw/Image.h>
#include <vw/Cartography/Datum.h>
#include <vw/Cartography/GeoReference.h>
#include <vw/Cartography/PointImageManipulation.h>
#include <vw/FileIO/DiskImageUtils.h>
#include <vw/Core/CmdUtils.h>
#include <asp/Core/Common.h>
#include <asp/Core/Macros.h>
#include <asp/Core/PointUtils.h>
#include <asp/Core/InterestPointMatching.h>
#include <liblas/liblas.hpp>

#include <limits>
#include <cstring>

#include <pointmatcher/PointMatcher.h>
#include <ceres/ceres.h>
#include <ceres/loss_function.h>

#include <asp/Tools/pc_align_utils.h>

#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
namespace fs = boost::filesystem;
namespace po = boost::program_options;

using namespace vw;
using namespace std;
using namespace vw::cartography;
using namespace asp;

typedef double RealT; // We will use doubles in libpointmatcher.
typedef PointMatcher<RealT> PM;
typedef PM::DataPoints DP;
using namespace PointMatcherSupport;

const double BIG_NUMBER = 1e+300; // libpointmatcher does not like here the largest double

/// Options container for the pc_align tool
struct Options : public vw::cartography::GdalWriteOptions {
  // Input
  string reference, source, init_transform_file, alignment_method, config_file,
    datum, csv_format_str, csv_proj4_str, match_file, hillshade_options,
    ipfind_options, ipmatch_options;
  PointMatcher<RealT>::Matrix init_transform;
  int    num_iter,
         max_num_reference_points,
         max_num_source_points;
  double diff_translation_err,
         diff_rotation_err,
         max_disp,
         outlier_ratio,
         semi_major,
         semi_minor;
  bool   compute_translation_only, use_hillshading,
         dont_use_dem_distances,
         save_trans_source,
         save_trans_ref,
         highest_accuracy,
         verbose;
  std::string initial_ned_translation;
  
  // Output
  string out_prefix;

  Options() : max_disp(-1.0), verbose(true){}

  /// Return true if the reference file is a DEM file and this option is not disabled
  bool use_dem_distances() const { return ( (asp::get_cloud_type(this->reference) == "DEM") && !dont_use_dem_distances); }
};


void handle_arguments( int argc, char *argv[], Options& opt ) {
  po::options_description general_options("");
  general_options.add_options()
    ("initial-transform",        po::value(&opt.init_transform_file)->default_value(""),
                                 "The file containing the transform to be used as an initial guess. It can come from a previous run of the tool.")
    ("num-iterations",           po::value(&opt.num_iter)->default_value(1000),
                                 "Maximum number of iterations.")
    ("diff-rotation-error",      po::value(&opt.diff_rotation_err)->default_value(1e-8),
                                 "Change in rotation amount below which the algorithm will stop (if translation error is also below bound), in degrees.")
    ("diff-translation-error",   po::value(&opt.diff_translation_err)->default_value(1e-3),
                                 "Change in translation amount below which the algorithm will stop (if rotation error is also below bound), in meters.")
    ("max-displacement",         po::value(&opt.max_disp)->default_value(0.0),
                                 "Maximum expected displacement of source points as result of alignment, in meters (after the initial guess transform is applied to the source points). Used for removing gross outliers in the source point cloud.")
    ("outlier-ratio",            po::value(&opt.outlier_ratio)->default_value(0.75),
                                 "Fraction of source (movable) points considered inliers (after gross outliers further than max-displacement from reference points are removed).")
    ("max-num-reference-points", po::value(&opt.max_num_reference_points)->default_value(100000000),
                                 "Maximum number of (randomly picked) reference points to use.")
    ("max-num-source-points",    po::value(&opt.max_num_source_points)->default_value(100000),
                                 "Maximum number of (randomly picked) source points to use (after discarding gross outliers).")
    ("alignment-method",         po::value(&opt.alignment_method)->default_value("point-to-plane"),
                                 "The type of iterative closest point method to use. [point-to-plane, point-to-point, similarity-point-to-point, least-squares, similarity-least-squares]")
    ("highest-accuracy",         po::bool_switch(&opt.highest_accuracy)->default_value(false)->implicit_value(true),
                                 "Compute with highest accuracy for point-to-plane (can be much slower).")
    ("csv-format",               po::value(&opt.csv_format_str)->default_value(""), asp::csv_opt_caption().c_str())
    ("csv-proj4",                po::value(&opt.csv_proj4_str)->default_value(""),
                                 "The PROJ.4 string to use to interpret the entries in input CSV files.")
    ("datum",                    po::value(&opt.datum)->default_value(""),
                                 "Use this datum for CSV files instead of auto-detecting it. Options: WGS_1984, D_MOON (1,737,400 meters), D_MARS (3,396,190 meters), MOLA (3,396,000 meters), NAD83, WGS72, and NAD27. Also accepted: Earth (=WGS_1984), Mars (=D_MARS), Moon (=D_MOON).")
    ("semi-major-axis",          po::value(&opt.semi_major)->default_value(0),
                                 "Explicitly set the datum semi-major axis in meters.")
    ("semi-minor-axis",          po::value(&opt.semi_minor)->default_value(0),
                                 "Explicitly set the datum semi-minor axis in meters.")
    ("output-prefix,o",          po::value(&opt.out_prefix)->default_value("run/run"),
                                 "Specify the output prefix.")
    ("compute-translation-only", po::bool_switch(&opt.compute_translation_only)->default_value(false)->implicit_value(true),
                                 "Compute the transform from source to reference point cloud as a translation only (no rotation).")
    ("save-transformed-source-points", po::bool_switch(&opt.save_trans_source)->default_value(false)->implicit_value(true),
                                  "Apply the obtained transform to the source points so they match the reference points and save them.")
    ("save-inv-transformed-reference-points", po::bool_switch(&opt.save_trans_ref)->default_value(false)->implicit_value(true),
     "Apply the inverse of the obtained transform to the reference points so they match the source points and save them.")

    ("initial-ned-translation", po::value(&opt.initial_ned_translation)->default_value(""),
                                 "Initialize the alignment transform based on a translation with this vector in the local North-East-Down coordinate system. Specify it in quotes, separated by spaces or commas.")

    ("initial-transform-from-hillshading", po::bool_switch(&opt.use_hillshading)->default_value(false)->implicit_value(true), "If both input clouds are DEMs, find interest point matches among their hillshaded versions, and use them to compute an initial rotation + translation + scale transform to apply to the source cloud before proceeding with alignment.")
    ("hillshade-options", po::value(&opt.hillshade_options)->default_value("--azimuth 300 --elevation 20"), "Options to pass to the hillshade program when computing the transform from hillshading.")
    ("ipfind-options", po::value(&opt.ipfind_options)->default_value("--ip-per-image 1000000 --interest-operator sift --descriptor-generator sift"), "Options to pass to the ipfind program when computing the transform from hillshading.")
    ("ipmatch-options", po::value(&opt.ipmatch_options)->default_value("--inlier-threshold 100 --ransac-constraint similarity"), "Options to pass to the ipmatch program when computing the transform from hillshading.")
    ("match-file", po::value(&opt.match_file)->default_value(""), "Compute a translation + rotation + scale transform from the source to the reference point cloud using manually selected point correspondences from the reference to the source (obtained for example using stereo_gui).")

    ("no-dem-distances",         po::bool_switch(&opt.dont_use_dem_distances)->default_value(false)->implicit_value(true),
                                 "For reference point clouds that are DEMs, don't take advantage of the fact that it is possible to interpolate into this DEM when finding the closest distance to it from a point in the source cloud and hence the error metrics.")

    ("config-file",              po::value(&opt.config_file)->default_value(""),
     "This is an advanced option. Read the alignment parameters from a configuration file, in the format expected by libpointmatcher, over-riding the command-line options.");

  //("verbose", po::bool_switch(&opt.verbose)->default_value(false)->implicit_value(true),
  // "Print debug information");

  general_options.add( vw::cartography::GdalWriteOptionsDescription(opt) );

  po::options_description positional("");
  positional.add_options()
    ("reference", po::value(&opt.reference), "The reference (fixed) point cloud/DEM.")
    ("source",    po::value(&opt.source),    "The source (movable) point cloud/DEM.");

  po::positional_options_description positional_desc;
  positional_desc.add("reference", 1);
  positional_desc.add("source",    1);

  string usage("--max-displacement arg [other options] <reference cloud> <source cloud> -o <output prefix>");
{'fraction_non_alphanumeric': 0.10598702199730645, 'fraction_numerical': 0.009100926417173407, 'mean_word_length': 3.4552727272727273, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 139, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6406090', 'n_tokens_mistral': 6859, 'n_tokens_neox': 6459, 'n_words': 2226}