Reworked the MNIST example to include the SIFT dataset.

Added support for reading the SIFT dataset to pico_toolshed.
Renamed the mnist example to kd_forest.
Documentation updated.

Author: Jaybro

README.md

@@ -61,7 +61,7 @@ PicoTree can interface with different types of points and point sets through tra
 * Creating a [custom search visitor](./examples/kd_tree/kd_tree_custom_search_visitor.cpp).
 * [Saving and loading](./examples/kd_tree/kd_tree_save_and_load.cpp) a KdTree to and from a file.
 * Support for [Eigen](./examples/eigen/eigen.cpp) and [OpenCV](./examples/opencv/opencv.cpp) data types.
-* Running the KdTree on the [MNIST](./examples/mnist/mnist.cpp) [database](http://yann.lecun.com/exdb/mnist/).
+* [Running the KdTree and KdForest](./examples/kd_forest/kd_forest.cpp) on the [MNIST](http://yann.lecun.com/exdb/mnist/) and [SIFT](http://corpus-texmex.irisa.fr/) datasets.
 * How to use the [KdTree with Python](./examples/python/kd_tree.py).
 
 # Requirements
@@ -113,9 +113,11 @@ $ pip install ./pico_tree
 
 # References
 
-* [Computational Geometry - Algorithms and Applications.](https://www.springer.com/gp/book/9783540779735) Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars, Springer-Verlag, third edition, 2008.
-* S. Maneewongvatana and D. M. Mount. [It's okay to be skinny, if your friends are fat.](http://www.cs.umd.edu/~mount/Papers/cgc99-smpack.pdf) 4th Annual CGC Workshop on Computational Geometry, 1999.
-* S. Arya and H. Y. Fu. [Expected-case complexity of approximate nearest neighbor searching.](https://www.cse.ust.hk/faculty/arya/pub/exp.pdf) InProceedings of the 11th ACM-SIAM Symposium on Discrete Algorithms, 2000.
-* S. Arya and D. M. Mount. [Algorithms for fast vector quantization.](https://www.cs.umd.edu/~mount/Papers/DCC.pdf) In IEEE Data Compression Conference, pages 381–390, March 1993.
-* N. Sample, M. Haines, M. Arnold and T. Purcell. [Optimizing Search Strategies in k-d Trees.](http://infolab.stanford.edu/~nsample/pubs/samplehaines.pdf) In: 5th WSES/IEEE World Multiconference on Circuits, Systems, Communications & Computers (CSCC 2001), July 2001.
-* A. Yershova and S. M. LaValle, [Improving Motion-Planning Algorithms by Efficient Nearest-Neighbor Searching.](http://msl.cs.uiuc.edu/~lavalle/papers/YerLav06.pdf) In IEEE Transactions on Robotics, vol. 23, no. 1, pp. 151-157, Feb. 2007.
+* J. L. Bentley, [Multidimensional binary search trees used for associative searching](https://dl.acm.org/doi/pdf/10.1145/361002.361007), Communications of the ACM, vol. 18, no. 9, pp. 509–517, 1975.
+* S. Arya and D. M. Mount, [Algorithms for fast vector quantization](https://www.cs.umd.edu/~mount/Papers/DCC.pdf), In IEEE Data Compression Conference, pp. 381–390, March 1993.
+* S. Maneewongvatana and D. M. Mount, [It's okay to be skinny, if your friends are fat](http://www.cs.umd.edu/~mount/Papers/cgc99-smpack.pdf), 4th Annual CGC Workshop on Computational Geometry, 1999.
+* S. Arya and H. Y. Fu, [Expected-case complexity of approximate nearest neighbor searching](https://www.cse.ust.hk/faculty/arya/pub/exp.pdf), InProceedings of the 11th ACM-SIAM Symposium on Discrete Algorithms, 2000.
+* N. Sample, M. Haines, M. Arnold and T. Purcell, [Optimizing Search Strategies in k-d Trees](http://infolab.stanford.edu/~nsample/pubs/samplehaines.pdf), In: 5th WSES/IEEE World Multiconference on Circuits, Systems, Communications & Computers (CSCC 2001), July 2001.
+* A. Yershova and S. M. LaValle, [Improving Motion-Planning Algorithms by Efficient Nearest-Neighbor Searching](http://msl.cs.uiuc.edu/~lavalle/papers/YerLav06.pdf), In IEEE Transactions on Robotics, vol. 23, no. 1, pp. 151-157, Feb. 2007.
+* M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars, [Computational Geometry - Algorithms and Applications](https://www.springer.com/gp/book/9783540779735), Springer-Verlag, third edition, 2008.
+* C. Silpa-Anan and R. Hartley, [Optimised KD-trees for fast image descriptor matching](http://vigir.missouri.edu/~gdesouza/Research/Conference_CDs/IEEE_CVPR_2008/data/papers/298.pdf), In CVPR, 2008.

examples/CMakeLists.txt

@@ -8,7 +8,7 @@ add_subdirectory(pico_understory)
 
 add_subdirectory(kd_tree)
 
-add_subdirectory(mnist)
+add_subdirectory(kd_forest)
 
 find_package(Eigen3 QUIET)
 

examples/kd_forest/CMakeLists.txt

@@ -0,0 +1,3 @@
+add_executable(kd_forest kd_forest.cpp)
+set_default_target_properties(kd_forest)
+target_link_libraries(kd_forest PUBLIC pico_toolshed pico_understory)

examples/kd_forest/kd_forest.cpp

@@ -0,0 +1,88 @@
+#include <iostream>
+#include <pico_toolshed/format/format_bin.hpp>
+#include <pico_toolshed/scoped_timer.hpp>
+#include <pico_tree/array_traits.hpp>
+#include <pico_tree/kd_tree.hpp>
+#include <pico_tree/vector_traits.hpp>
+#include <pico_understory/kd_forest.hpp>
+
+#include "mnist.hpp"
+#include "sift.hpp"
+
+template <typename Dataset>
+void RunDataset(
+    std::size_t max_leaf_size_exact,
+    std::size_t max_leaf_size_apprx,
+    std::size_t forest_size) {
+  using Point = typename Dataset::PointType;
+  using Space = std::reference_wrapper<std::vector<Point>>;
+  using Scalar = typename Point::value_type;
+
+  auto train = Dataset::ReadTrain();
+  auto test = Dataset::ReadTest();
+  std::size_t count = test.size();
+  std::vector<pico_tree::Neighbor<int, Scalar>> nns(count);
+  std::string fn_nns_gt = Dataset::kDatasetName + "_nns_gt.bin";
+
+  if (!std::filesystem::exists(fn_nns_gt)) {
+    auto kd_tree = [&train, &max_leaf_size_exact]() {
+      ScopedTimer t0("kd_tree build");
+      return pico_tree::KdTree<Space>(train, max_leaf_size_exact);
+    }();
+
+    {
+      ScopedTimer t1("kd_tree query");
+      for (std::size_t i = 0; i < nns.size(); ++i) {
+        kd_tree.SearchNn(test[i], nns[i]);
+      }
+    }
+
+    std::cout << "Writing " << fn_nns_gt << "." << std::endl;
+    pico_tree::WriteBin(fn_nns_gt, nns);
+  } else {
+    std::cout << "Reading " << fn_nns_gt << "." << std::endl;
+    pico_tree::ReadBin(fn_nns_gt, nns);
+  }
+
+  std::size_t equal = 0;
+
+  // Building the KdForest takes roughly forest_size times longer compared to
+  // building the regular KdTree. However, it is usually a lot faster.
+  {
+    auto rkd_tree = [&train, &max_leaf_size_apprx, &forest_size]() {
+      ScopedTimer t0("kd_forest build");
+      return pico_tree::KdForest<Space>(
+          train, max_leaf_size_apprx, forest_size);
+    }();
+
+    ScopedTimer t1("kd_forest query");
+    pico_tree::Neighbor<int, Scalar> nn;
+    for (std::size_t i = 0; i < nns.size(); ++i) {
+      rkd_tree.SearchNn(test[i], nn);
+
+      if (nns[i].index == nn.index) {
+        ++equal;
+      }
+    }
+  }
+
+  std::cout << "Precision: "
+            << (static_cast<float>(equal) / static_cast<float>(count))
+            << std::endl;
+}
+
+int main() {
+  // max_leaf_size_apprx = 128:
+  // forest_size 8: a precision of around 79%.
+  // forest_size 16: a precision of around 93%.
+  // forest_size 32: a precision of around 98%.
+  RunDataset<Mnist>(16, 128, 8);
+  // max_leaf_size_apprx = 1024:
+  // forest_size 8: a precision of around 58%.
+  // forest_size 16: a precision of around 68%.
+  // forest_size 32: a precision of around 80%.
+  // forest_size 64: a precision of around 87%.
+  // forest_size 128: out of memory :'(
+  RunDataset<Sift>(16, 1024, 8);
+  return 0;
+}

examples/kd_forest/mnist.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <algorithm>
+#include <filesystem>
+#include <pico_toolshed/format/format_mnist.hpp>
+
+template <typename U, typename T, std::size_t N>
+std::array<U, N> Cast(std::array<T, N> const& i) {
+  std::array<U, N> c;
+  std::transform(i.begin(), i.end(), c.begin(), [](T a) -> U {
+    return static_cast<U>(a);
+  });
+  return c;
+}
+
+template <typename U, typename T, std::size_t N>
+std::vector<std::array<U, N>> Cast(std::vector<std::array<T, N>> const& i) {
+  std::vector<std::array<U, N>> c;
+  std::transform(
+      i.begin(),
+      i.end(),
+      std::back_inserter(c),
+      [](std::array<T, N> const& a) -> std::array<U, N> { return Cast<U>(a); });
+  return c;
+}
+
+class Mnist {
+ private:
+  using Scalar = float;
+  using ImageByte = std::array<std::byte, 28 * 28>;
+  using ImageFloat = std::array<Scalar, 28 * 28>;
+
+  static std::vector<ImageFloat> ReadImages(std::string const& filename) {
+    if (!std::filesystem::exists(filename)) {
+      throw std::runtime_error(filename + " doesn't exist.");
+    }
+
+    std::vector<ImageByte> images_u8;
+    pico_tree::ReadMnistImages(filename, images_u8);
+    return Cast<Scalar>(images_u8);
+  }
+
+ public:
+  using PointType = ImageFloat;
+
+  static std::string const kDatasetName;
+
+  static std::vector<PointType> ReadTrain() {
+    std::string fn_images_train = "train-images.idx3-ubyte";
+    return ReadImages(fn_images_train);
+  }
+
+  static std::vector<PointType> ReadTest() {
+    std::string fn_images_test = "t10k-images.idx3-ubyte";
+    return ReadImages(fn_images_test);
+  }
+};
+
+std::string const Mnist::kDatasetName = "mnist";

examples/kd_forest/sift.hpp

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <filesystem>
+#include <pico_toolshed/format/format_xvecs.hpp>
+
+class Sift {
+ private:
+  using VectorFloat = std::array<float, 128>;
+
+  static std::vector<VectorFloat> ReadVectors(std::string const& filename) {
+    if (!std::filesystem::exists(filename)) {
+      throw std::runtime_error(filename + " doesn't exist.");
+    }
+
+    std::vector<VectorFloat> vectors;
+    pico_tree::ReadXvecs(filename, vectors);
+    return vectors;
+  }
+
+ public:
+  using PointType = VectorFloat;
+
+  static std::string const kDatasetName;
+
+  static std::vector<PointType> ReadTrain() {
+    std::string fn_images_train = "sift_base.fvecs";
+    return ReadVectors(fn_images_train);
+  }
+
+  static std::vector<PointType> ReadTest() {
+    std::string fn_images_test = "sift_query.fvecs";
+    return ReadVectors(fn_images_test);
+  }
+};
+
+std::string const Sift::kDatasetName = "sift";

examples/mnist/CMakeLists.txt

@@ -1,6 +1,6 @@
 #pragma once
 
-#include <limits>
+#include <filesystem>
 #include <pico_tree/internal/stream.hpp>
 
 namespace pico_tree {
@@ -22,20 +22,10 @@ void ReadBin(std::string const& filename, std::vector<T>& v) {
   std::fstream stream =
       internal::OpenStream(filename, std::ios::in | std::ios::binary);
 
-  // The four lines below are used when determining the file size.
-  // C++17 is not used here in order to keep the benchmark C++11.
-  stream.ignore(std::numeric_limits<std::streamsize>::max());
-  std::streamsize byte_count = stream.gcount();
-  stream.clear();
-  stream.seekg(0, std::ios_base::beg);
-
-  if (byte_count == 0) {
-    return;
-  }
-
+  auto bytes = std::filesystem::file_size(filename);
   std::size_t const element_size = sizeof(T);
   std::size_t const element_count =
-      static_cast<std::size_t>(byte_count / element_size);
+      static_cast<std::size_t>(bytes) / element_size;
   v.resize(element_count);
   stream.read(reinterpret_cast<char*>(&v[0]), element_size * element_count);
 }