Implement feature matching with cascade hashing

Using cascade hashing for feature matching trades a slight drop in
matching accuracy for a speed-up of an order of magnitude.
Speed-up over bruteforce matching on the citywall-20140923 dataset:
1x Intel Core i7-3930K   : ~20x
2x Intel Xeon E5-2687W v3: ~18x

apps/sfmrecon/sfmrecon.cc

@@ -58,6 +58,7 @@ struct AppSettings
     int initial_pair_1 = -1;
     int initial_pair_2 = -1;
     int min_views_per_track = 3;
+    bool cascade_hashing = false;
 };
 
 void
@@ -111,6 +112,9 @@ features_and_matching (mve::Scene::Ptr scene, AppSettings const& conf,
     matching_opts.ransac_opts.verbose_output = false;
     matching_opts.use_lowres_matching = conf.lowres_matching;
     matching_opts.match_num_previous_frames = conf.video_matching;
+    matching_opts.matcher_type = conf.cascade_hashing
+        ? sfm::bundler::Matching::MATCHER_CASCADE_HASHING
+        : sfm::bundler::Matching::MATCHER_EXHAUSTIVE;
 
     std::cout << "Performing feature matching..." << std::endl;
     {
@@ -491,6 +495,7 @@ main (int argc, char** argv)
     args.add_option('\0', "track-thres-factor", true, "Error threshold factor for tracks [25]");
     args.add_option('\0', "use-2cam-tracks", false, "Triangulate tracks from only two cameras");
     args.add_option('\0', "initial-pair", true, "Manually specify initial pair IDs [-1,-1]");
+    args.add_option('\0', "cascade-hashing", false, "Use cascade hashing for matching [false]");
     args.parse(argc, argv);
 
     /* Setup defaults. */
@@ -549,6 +554,8 @@ main (int argc, char** argv)
             std::cout << "Using initial pair (" << conf.initial_pair_1
                 << "," << conf.initial_pair_2 << ")." << std::endl;
         }
+        else if (i->opt->lopt == "cascade-hashing")
+            conf.cascade_hashing = true;
         else
         {
             std::cerr << "Error: Unexpected option: "

libs/math/matrix.h

@@ -53,6 +53,7 @@ template <typename T, int N, int M>
 class Matrix
 {
 public:
+    typedef T ValueType;
 
     /* ------------------------ Constructors ---------------------- */
 

libs/math/vector.h

@@ -80,6 +80,7 @@ template <typename T, int N>
 class Vector
 {
 public:
+    typedef T ValueType;
 
     /* ------------------------ Constructors ---------------------- */
 

libs/sfm/bundler_matching.cc

@@ -18,6 +18,7 @@
 #include "sfm/sift.h"
 #include "sfm/ransac.h"
 #include "sfm/bundler_matching.h"
+#include "sfm/cascade_hashing.h"
 #include "sfm/exhaustive_matching.h"
 
 SFM_NAMESPACE_BEGIN
@@ -32,6 +33,9 @@ Matching::Matching (Options const& options, Progress* progress)
         case MATCHER_EXHAUSTIVE:
             this->matcher.reset(new ExhaustiveMatching());
             break;
+        case MATCHER_CASCADE_HASHING:
+            this->matcher.reset(new CascadeHashing());
+            break;
         default:
             throw std::runtime_error("Unhandled matcher type");
     }

libs/sfm/bundler_matching.h

@@ -51,7 +51,8 @@ class Matching
 public:
     enum MatcherType
     {
-        MATCHER_EXHAUSTIVE
+        MATCHER_EXHAUSTIVE,
+        MATCHER_CASCADE_HASHING
     };
 
     /** Options for feature matching. */

libs/sfm/cascade_hashing.cc

@@ -0,0 +1,260 @@
+/*
+ * Copyright (C) 2016, Andre Schulz
+ * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
+ * All rights reserved.
+ *
+ * This software may be modified and distributed under the terms
+ * of the BSD 3-Clause license. See the LICENSE.txt file for details.
+ */
+
+#include <algorithm>
+#include <cstdint>
+#include <vector>
+
+#include "sfm/cascade_hashing.h"
+
+SFM_NAMESPACE_BEGIN
+
+void
+CascadeHashing::GlobalData::generate_proj_matrices (Options const& opts)
+{
+    generate_proj_matrices(
+        &this->sift.prim_proj_mat,
+        &this->sift.sec_proj_mats,
+        opts);
+    generate_proj_matrices(
+        &this->surf.prim_proj_mat,
+        &this->surf.sec_proj_mats,
+        opts);
+}
+
+/* ---------------------------------------------------------------- */
+
+void
+CascadeHashing::init (bundler::ViewportList* viewports)
+{
+    ExhaustiveMatching::init(viewports);
+
+    util::WallTimer timer;
+    this->local_data_sift.clear();
+    this->local_data_sift.resize(viewports->size());
+    this->local_data_surf.clear();
+    this->local_data_surf.resize(viewports->size());
+
+    this->global_data.generate_proj_matrices(this->cashash_opts);
+
+    /* Convert feature descriptors to zero mean. */
+    math::Vec128f sift_avg;
+    math::Vec64f surf_avg;
+    compute_avg_descriptors(this->processed_feature_sets, &sift_avg, &surf_avg);
+
+#pragma omp parallel for schedule(dynamic)
+    for (std::size_t i = 0; i < viewports->size(); i++)
+    {
+        LocalData* ld_sift = &this->local_data_sift[i];
+        LocalData* ld_surf = &this->local_data_surf[i];
+
+        ProcessedFeatureSet const& pfs = this->processed_feature_sets[i];
+        std::vector<math::Vec128f> sift_zero_mean_descs;
+        std::vector<math::Vec64f> surf_zero_mean_descs;
+        this->compute_zero_mean_descs(&sift_zero_mean_descs,
+            &surf_zero_mean_descs, pfs.sift_descr, pfs.surf_descr, sift_avg,
+            surf_avg);
+
+        this->compute(ld_sift, ld_surf, sift_zero_mean_descs,
+            surf_zero_mean_descs, this->global_data, this->cashash_opts);
+    }
+    std::cout << "Computing cascade hashes took " << timer.get_elapsed()
+        << " ms" << std::endl;
+}
+
+void
+CascadeHashing::pairwise_match (int view_1_id, int view_2_id,
+    Matching::Result* result) const
+{
+    /* SIFT matching. */
+    ProcessedFeatureSet const& pfs_1 = this->processed_feature_sets[view_1_id];
+    ProcessedFeatureSet const& pfs_2 = this->processed_feature_sets[view_2_id];
+    LocalData const& ld_sift_1 = this->local_data_sift[view_1_id];
+    LocalData const& ld_sift_2 = this->local_data_sift[view_2_id];
+    Matching::Result sift_result;
+    if (pfs_1.sift_descr.size() > 0)
+    {
+        this->twoway_match(this->opts.sift_matching_opts,
+            ld_sift_1, ld_sift_2,
+            pfs_1.sift_descr, pfs_2.sift_descr,
+            &sift_result, this->cashash_opts);
+        Matching::remove_inconsistent_matches(&sift_result);
+    }
+
+    /* SURF matching. */
+    LocalData const& ld_surf_1 = this->local_data_surf[view_1_id];
+    LocalData const& ld_surf_2 = this->local_data_surf[view_2_id];
+    Matching::Result surf_result;
+    if (pfs_1.surf_descr.size() > 0)
+    {
+        this->twoway_match(this->opts.surf_matching_opts,
+            ld_surf_1, ld_surf_2,
+            pfs_1.surf_descr, pfs_2.surf_descr,
+            &surf_result, this->cashash_opts);
+        Matching::remove_inconsistent_matches(&surf_result);
+    }
+
+    /* TODO: The following code is the same as in
+     *       ExhaustiveMatching::pairwise_match() and could be moved into a
+     *       separate function to avoid code duplication.
+     */
+    /* Fix offsets in the matching result. */
+    std::size_t other_surf_offset = pfs_2.sift_descr.size();
+    if (other_surf_offset > 0)
+        for (std::size_t i = 0; i < surf_result.matches_1_2.size(); ++i)
+            if (surf_result.matches_1_2[i] >= 0)
+                surf_result.matches_1_2[i] += other_surf_offset;
+
+    std::size_t this_surf_offset = pfs_1.sift_descr.size();
+    if (this_surf_offset > 0)
+        for (std::size_t i = 0; i < surf_result.matches_2_1.size(); ++i)
+            if (surf_result.matches_2_1[i] >= 0)
+                surf_result.matches_2_1[i] += this_surf_offset;
+
+    /* Create a combined matching result. */
+    std::size_t this_num_descriptors = pfs_1.sift_descr.size()
+        + pfs_1.surf_descr.size();
+    std::size_t other_num_descriptors = pfs_2.sift_descr.size()
+        + pfs_2.surf_descr.size();
+
+    result->matches_1_2.clear();
+    result->matches_1_2.reserve(this_num_descriptors);
+    result->matches_1_2.insert(result->matches_1_2.end(),
+        sift_result.matches_1_2.begin(), sift_result.matches_1_2.end());
+    result->matches_1_2.insert(result->matches_1_2.end(),
+        surf_result.matches_1_2.begin(), surf_result.matches_1_2.end());
+
+    result->matches_2_1.clear();
+    result->matches_2_1.reserve(other_num_descriptors);
+    result->matches_2_1.insert(result->matches_2_1.end(),
+        sift_result.matches_2_1.begin(), sift_result.matches_2_1.end());
+    result->matches_2_1.insert(result->matches_2_1.end(),
+        surf_result.matches_2_1.begin(), surf_result.matches_2_1.end());
+}
+
+void
+CascadeHashing::compute(
+    LocalData* ld_sift, LocalData* ld_surf,
+    std::vector<math::Vec128f> const& sift_zero_mean_descs,
+    std::vector<math::Vec64f> const& surf_zero_mean_descs,
+    GlobalData const& cashash_global_data,
+    Options const& cashash_opts)
+{
+    /* Compute cascade hashes of zero mean SIFT/SURF descriptors. */
+    compute_cascade_hashes(
+        sift_zero_mean_descs,
+        &ld_sift->comp_hash_data,
+        &ld_sift->bucket_grps_bucket_ids,
+        cashash_global_data.sift.prim_proj_mat,
+        cashash_global_data.sift.sec_proj_mats,
+        cashash_opts);
+    compute_cascade_hashes(
+        surf_zero_mean_descs,
+        &ld_surf->comp_hash_data,
+        &ld_surf->bucket_grps_bucket_ids,
+        cashash_global_data.surf.prim_proj_mat,
+        cashash_global_data.surf.sec_proj_mats,
+        cashash_opts);
+
+    /* Build buckets. */
+    build_buckets(
+        &ld_sift->bucket_grps_feature_ids,
+        ld_sift->bucket_grps_bucket_ids,
+        sift_zero_mean_descs.size(),
+        cashash_opts);
+    build_buckets(
+        &ld_surf->bucket_grps_feature_ids,
+        ld_surf->bucket_grps_bucket_ids,
+        surf_zero_mean_descs.size(),
+        cashash_opts);
+}
+
+/* ---------------------------------------------------------------- */
+
+void
+CascadeHashing::compute_avg_descriptors (ProcessedFeatureSets const& pfs,
+    math::Vec128f* sift_avg, math::Vec64f* surf_avg)
+{
+    math::Vec128f sift_vec_sum(0.0f);
+    math::Vec64f surf_vec_sum(0.0f);
+    std::size_t num_sift_descs_total = 0;
+    std::size_t num_surf_descs_total = 0;
+
+    /* Compute sum of all SIFT/SURF descriptors. */
+    for (std::size_t i = 0; i < pfs.size(); i++)
+    {
+        SiftDescriptors const& sift_descr = pfs[i].sift_descr;
+        SurfDescriptors const& surf_descr = pfs[i].surf_descr;
+
+        std::size_t num_sift_descriptors = sift_descr.size();
+        std::size_t num_surf_descriptors = surf_descr.size();
+        num_sift_descs_total += num_sift_descriptors;
+        num_surf_descs_total += num_surf_descriptors;
+
+        for (std::size_t j = 0; j < num_sift_descriptors; j++)
+            for (int k = 0; k < 128; k++)
+                sift_vec_sum[k] += sift_descr[j][k] / 255.0f;
+
+        for (std::size_t j = 0; j < num_surf_descriptors; j++)
+            for (int k = 0; k < 64; k++)
+                surf_vec_sum[k] += surf_descr[j][k] / 127.0f;
+    }
+
+    /* Compute average vectors for SIFT/SURF. */
+    *sift_avg = sift_vec_sum / num_sift_descs_total;
+    *surf_avg = surf_vec_sum / num_surf_descs_total;
+}
+
+void
+CascadeHashing::compute_zero_mean_descs(
+    std::vector<math::Vec128f>* sift_zero_mean_descs,
+    std::vector<math::Vec64f>* surf_zero_mean_descs,
+    SiftDescriptors const& sift_descs, SurfDescriptors const& surf_descs,
+    math::Vec128f const& sift_avg, math::Vec64f const& surf_avg)
+{
+    /* Compute zero mean descriptors. */
+    sift_zero_mean_descs->resize(sift_descs.size());
+    for (std::size_t i = 0; i < sift_descs.size(); i++)
+        for (int j = 0; j < 128; j++)
+            (*sift_zero_mean_descs)[i][j] = sift_descs[i][j] / 255.0f - sift_avg[j];
+
+    surf_zero_mean_descs->resize(surf_descs.size());
+    for (std::size_t i = 0; i < surf_descs.size(); i++)
+        for (int j = 0; j < 64; j++)
+            (*surf_zero_mean_descs)[i][j] = surf_descs[i][j] / 127.0f - surf_avg[j];
+}
+
+/* ---------------------------------------------------------------- */
+
+void
+CascadeHashing::build_buckets(
+    BucketGroupsFeatures *bucket_grps_feature_ids,
+    BucketGroupsBuckets const& bucket_grps_bucket_ids,
+    size_t num_descs,
+    Options const& opts)
+{
+    uint8_t const num_bucket_grps = opts.num_bucket_groups;
+    uint8_t const num_bucket_bits = opts.num_bucket_bits;
+    uint32_t const num_buckets_per_group = 1 << num_bucket_bits;
+
+    bucket_grps_feature_ids->resize(num_bucket_grps);
+
+    for (uint8_t grp_idx = 0; grp_idx < num_bucket_grps; grp_idx++)
+    {
+        BucketGroupFeatures &bucket_grp_features = (*bucket_grps_feature_ids)[grp_idx];
+        bucket_grp_features.resize(num_buckets_per_group);
+        for (size_t i = 0; i < num_descs; i++)
+        {
+            uint16_t bucket_id = bucket_grps_bucket_ids[grp_idx][i];
+            bucket_grp_features[bucket_id].emplace_back(i);
+        }
+    }
+}
+
+SFM_NAMESPACE_END