text:
/****************************************************************************
 * Copyright (c) 2017-2021 by the ArborX authors                            *
 * All rights reserved.                                                     *
 *                                                                          *
 * This file is part of the ArborX library. ArborX is                       *
 * distributed under a BSD 3-clause license. For the licensing terms see    *
 * the LICENSE file in the top-level directory.                             *
 *                                                                          *
 * SPDX-License-Identifier: BSD-3-Clause                                    *
 ****************************************************************************/

#ifndef ARBORX_BOOST_RTREE_HELPERS_HPP
#define ARBORX_BOOST_RTREE_HELPERS_HPP

#include <ArborX_Config.hpp>

#include "ArborX_BoostGeometryAdapters.hpp"
#include "ArborX_BoostRangeAdapters.hpp"
#include <ArborX_Box.hpp>
#include <ArborX_DetailsKokkosExtAccessibilityTraits.hpp> // is_accessible_from_host
#include <ArborX_DetailsUtils.hpp> // exclusivePrefixSum, lastElement
#include <ArborX_Point.hpp>
#include <ArborX_Predicates.hpp>

#include <boost/range/adaptors.hpp>
#include <boost/range/algorithm/copy.hpp>
#include <boost/range/algorithm/transform.hpp>
#include <boost/range/combine.hpp>

#ifdef ARBORX_ENABLE_MPI
#include <mpi.h>
#endif

namespace BoostRTreeHelpers
{
// NOTE: The balancing algorithm does not really matter since the tree is
// created using packing algorithm and there is no later insertion of new
// objects into the tree.  The values of the maximum and the minimum number
// of elements could be adjusted in principle but I didn't bother setting
// this type as a template argument.
using Parameter = boost::geometry::index::linear<16>;

template <typename Value, typename Index>
struct PairMaker
{
  using result_type = std::pair<Value, Index>;
  template <typename T>
  inline result_type operator()(T const &v) const
  {
    return result_type(v.value(), v.index());
  }
};

template <typename Indexable>
using RTree =
    boost::geometry::index::rtree<std::pair<Indexable, int>, Parameter>;

template <typename View>
static RTree<typename View::value_type> makeRTree(View const &objects)
{
  using Indexable = typename View::value_type;
  return RTree<Indexable>(
      objects | boost::adaptors::indexed() |
      boost::adaptors::transformed(PairMaker<Indexable, int>()));
}

// NOTE: Boost.Config defines BOOST_NO_CXX11_VARIADIC_TEMPLATES for nvcc
// with the current version of CUDA we are using.  In consequence we are
// not able to use std::tuple<Indexable, ...> which is unfortunate :(
class AppendRankToPairObjectIndex
{
public:
  AppendRankToPairObjectIndex(int rank)
      : _rank(rank)
  {
  }

  template <typename T1, typename T2>
  inline boost::tuple<T1, T2, int> operator()(std::pair<T1, T2> const &p) const
  {
    return boost::make_tuple(std::get<0>(p), std::get<1>(p), _rank);
  }

private:
  int _rank = -1;
};

#ifdef ARBORX_ENABLE_MPI
template <typename Indexable>
using ParallelRTree =
    boost::geometry::index::rtree<boost::tuple<Indexable, int, int>, Parameter>;

template <typename View>
static ParallelRTree<typename View::value_type> makeRTree(MPI_Comm comm,
                                                          View const &objects)
{
  using Indexable = typename View::value_type;

  // Fill buffer with pair (object, index)
  std::vector<std::pair<Indexable, int>> buffer;
  boost::copy(objects | boost::adaptors::indexed() |
                  boost::adaptors::transformed(PairMaker<Indexable, int>()),
              std::back_inserter(buffer));

  // Gather all buffers
  int comm_size;
  MPI_Comm_size(comm, &comm_size);
  int comm_rank;
  MPI_Comm_rank(comm, &comm_rank);
  std::vector<int> counts(comm_size);
  counts[comm_rank] = buffer.size();
  MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, counts.data(), 1, MPI_INT,
                comm);
  std::vector<int> offsets(comm_size + 1);
  offsets[0] = 0;
  for (int i = 0; i < comm_size; ++i)
    offsets[i + 1] = counts[i] + offsets[i];
  decltype(buffer) all_buffers(offsets.back());
  auto const bytes_per_element = sizeof(typename decltype(buffer)::value_type);
  for (auto *pv : {&counts, &offsets})
    for (auto &x : *pv)
      x *= bytes_per_element;
  MPI_Allgatherv(buffer.data(), counts[comm_rank], MPI_BYTE, all_buffers.data(),
                 counts.data(), offsets.data(), MPI_BYTE, comm);
  for (auto &x : offsets)
    x /= bytes_per_element;

  std::vector<boost::tuple<Indexable, int, int>> all_objects;
  for (int i = 0; i < comm_size; ++i)
    boost::transform(std::make_pair(all_buffers.data() + offsets[i],
                                    all_buffers.data() + offsets[i + 1]),
                     std::back_inserter(all_objects),
                     AppendRankToPairObjectIndex(i));

  return ParallelRTree<Indexable>(all_objects);
}
#endif

template <typename Value>
struct UnaryPredicate
{
  using Function = std::function<bool(Value const &)>;
  UnaryPredicate(Function pred)
      : _pred(pred)
  {
  }
  inline bool operator()(Value const &val) const { return _pred(val); }
  Function _pred;
};

template <typename Value>
static auto translate(ArborX::Intersects<ArborX::Sphere> const &query)
{
  auto const sphere = getGeometry(query);
  auto const radius = sphere.radius();
  auto const centroid = sphere.centroid();
  ArborX::Box box;
  ArborX::Details::expand(box, sphere);
  return boost::geometry::index::intersects(box) &&
         boost::geometry::index::satisfies(
             UnaryPredicate<Value>([centroid, radius](Value const &val) {
               boost::geometry::index::indexable<Value> indexableGetter;
               auto const &geometry = indexableGetter(val);
               return boost::geometry::distance(centroid, geometry) <= radius;
             }));
}

template <typename Value, typename Geometry>
static auto translate(ArborX::Nearest<Geometry> const &query)
{
  auto const geometry = getGeometry(query);
  auto const k = getK(query);
  return boost::geometry::index::nearest(geometry, k);
}

template <typename Indexable, typename InputView,
          typename OutputView = Kokkos::View<int *, Kokkos::HostSpace>>
static std::tuple<OutputView, OutputView>
performQueries(RTree<Indexable> const &rtree, InputView const &queries)
{
  static_assert(KokkosExt::is_accessible_from_host<InputView>::value, "");

  using Value = typename RTree<Indexable>::value_type;
  auto const n_queries = queries.extent_int(0);
  OutputView offset("offset", n_queries + 1);
  std::vector<Value> returned_values;
  for (int i = 0; i < n_queries; ++i)
    offset(i) = rtree.query(translate<Value>(queries(i)),
                            std::back_inserter(returned_values));
  using ExecutionSpace = typename InputView::execution_space;
  ArborX::exclusivePrefixSum(ExecutionSpace{}, offset);
  auto const n_results = ArborX::lastElement(offset);
  OutputView indices("indices", n_results);
  for (int i = 0; i < n_queries; ++i)
    for (int j = offset(i); j < offset(i + 1); ++j)
      indices(j) = returned_values[j].second;
  return std::make_tuple(offset, indices);
}

#ifdef ARBORX_ENABLE_MPI
template <typename Indexable, typename InputView,
          typename OutputView1 =
              Kokkos::View<Kokkos::pair<int, int> *, Kokkos::HostSpace>,
          typename OutputView2 = Kokkos::View<int *, Kokkos::HostSpace>>
static std::tuple<OutputView2, OutputView1>
performQueries(ParallelRTree<Indexable> const &rtree, InputView const &queries)
{
  static_assert(KokkosExt::is_accessible_from_host<InputView>::value, "");
  using Value = typename ParallelRTree<Indexable>::value_type;
  auto const n_queries = queries.extent_int(0);
  OutputView2 offset("offset", n_queries + 1);
  std::vector<Value> returned_values;
  for (int i = 0; i < n_queries; ++i)
    offset(i) = rtree.query(translate<Value>(queries(i)),
                            std::back_inserter(returned_values));
  using ExecutionSpace = typename InputView::execution_space;
  ArborX::exclusivePrefixSum(ExecutionSpace{}, offset);
  auto const n_results = ArborX::lastElement(offset);
  OutputView1 values("values", n_results);
  for (int i = 0; i < n_queries; ++i)
    for (int j = offset(i); j < offset(i + 1); ++j)
    {
      int index;
      int rank;
      boost::tie(boost::tuples::ignore, index, rank) = returned_values[j];
      values(j) = {index, rank};
    }
  return std::make_tuple(offset, values);
}
#endif
} // end namespace BoostRTreeHelpers

namespace BoostExt
{

// FIXME Goal is to match the BVH interface
template <typename Indexable>
class RTree
{
public:
  using DeviceType = Kokkos::DefaultHostExecutionSpace::device_type;
  using memory_space = typename DeviceType::memory_space;

  template <typename ExecutionSpace>
  RTree(ExecutionSpace, Kokkos::View<Indexable *, DeviceType> const &values)
  {
    static_assert(Kokkos::is_execution_space<ExecutionSpace>::value, "");

    _tree = BoostRTreeHelpers::makeRTree(values);
  }

  // WARNING trailing pack will match anything :/
  template <typename ExecutionSpace, typename Predicates, typename InputView,
            typename... TrailingArgs>
  void query(ExecutionSpace const &, Predicates const &predicates,
             InputView &indices, InputView &offset, TrailingArgs &&...) const
  {
    static_assert(Kokkos::is_execution_space<ExecutionSpace>::value, "");

    std::tie(offset, indices) =
        BoostRTreeHelpers::performQueries(_tree, predicates);
  }

  template <typename ExecutionSpace, typename Predicates, typename Callback,
            typename... TrailingArgs>
  void query(ExecutionSpace const &, Predicates const &, Callback const &,
             TrailingArgs &&...) const
  {
    static_assert(Kokkos::is_execution_space<ExecutionSpace>::value, "");

    throw std::runtime_error(
        "Boost RTree does not support callback only query overload.");
  }

private:
  BoostRTreeHelpers::RTree<Indexable> _tree;
};

#ifdef ARBORX_ENABLE_MPI
template <typename Indexable>
class ParallelRTree
{
public:
  using DeviceType = Kokkos::DefaultHostExecutionSpace::device_type;
  using memory_space = typename DeviceType::memory_space;

  template <typename ExecutionSpace>
  ParallelRTree(MPI_Comm comm, ExecutionSpace const &,
                Kokkos::View<Indexable *, DeviceType> const &values)
  {
    static_assert(Kokkos::is_execution_space<ExecutionSpace>::value, "");

    _tree = BoostRTreeHelpers::makeRTree(comm, values);
  }

  // WARNING trailing pack will match anything :/
  template <typename ExecutionSpace, typename Predicates, typename InputView1,
            typename InputView2, typename... TrailingArgs>
  void query(ExecutionSpace const &, Predicates const &predicates,
             InputView1 &indices, InputView2 &offset, TrailingArgs &&...) const
  {
    static_assert(Kokkos::is_execution_space<ExecutionSpace>::value, "");

    std::tie(offset, indices) =
        BoostRTreeHelpers::performQueries(_tree, predicates);
  }

private:
  BoostRTreeHelpers::ParallelRTree<Indexable> _tree;
};
#endif

} // namespace BoostExt

namespace ArborX
{
// Specialization of ArborX::query
template <typename Indexable, typename ExecutionSpace, typename Predicates,
          typename InputView, typename... TrailingArgs>
inline void query(BoostExt::RTree<Indexable> const &rtree,
                  ExecutionSpace const &space, Predicates const &predicates,
                  InputView &indices, InputView &offset,
                  TrailingArgs &&... args)
{
  rtree.query(space, predicates, indices, offset,
              std::forward<TrailingArgs>(args)...);
}
} // namespace ArborX

#endif


text:
// Copyright (c) 2021 Graphcore Ltd. All rights reserved.
#include "CTCInferenceCodeletTestConnection.hpp"

#include <poplar/Engine.hpp>
#include <poplar/Graph.hpp>
#include <poplar/IPUModel.hpp>
#include <poplar/Type.hpp>

#include <poplibs_support/TestDevice.hpp>
#include <poplibs_test/CTCInference.hpp>
#include <poplibs_test/Util.hpp>
#include <poputil/VertexTemplates.hpp>

#include <boost/multi_array.hpp>

#include <limits>

using namespace poplar;
using namespace poplar::program;
using namespace poplibs_test::ctc;
using namespace poplibs_test;
using namespace poplibs_test::util;
using namespace poplibs_support;
using namespace poputil;

namespace poplibs_test {
namespace ctc {

template <typename PartialsType>
std::vector<Candidate<PartialsType>> runMergeCandidatesCodelet(
    Graph &graph, TestDevice &device, DeviceType deviceType, Type inType,
    Type partialsType,
    const std::vector<Candidate<PartialsType>> &extendCandidates,
    const Candidate<PartialsType> &copyCandidate, unsigned timestep,
    unsigned blankClass, const BeamHistory &beamHistory,
    const ArrayRef<unsigned> &outputLengths, unsigned lastBeamOutputSym,
    unsigned numClasses, bool profile) {

  const auto target = graph.getTarget();
  const auto beamwidth = beamHistory.symbols.size();
  const auto numExtendCandidates = extendCandidates.size();
  const auto maxT = beamHistory.symbols[0].size() - 1;

  auto beamAddend =
      graph.addVariable(UNSIGNED_INT, {maxT, beamwidth}, "beamAddend");
  auto beamParent =
      graph.addVariable(UNSIGNED_INT, {maxT, beamwidth}, "beamParent");
  auto lastBeamOutput =
      graph.addConstant(UNSIGNED_INT, {}, lastBeamOutputSym, "lastBeamOutput");

  auto beamLength = graph.addConstant(UNSIGNED_INT, {2 * beamwidth},
                                      outputLengths, "beamLength");

  auto currentTimestep = graph.addConstant(UNSIGNED_INT, {}, timestep);
  auto complete = graph.addConstant(UNSIGNED_INT, {}, 0u);

  graph.setTileMapping(beamAddend, 0);
  graph.setTileMapping(beamParent, 0);
  graph.setTileMapping(lastBeamOutput, 0);
  graph.setTileMapping(beamLength, 0);

  graph.setTileMapping(currentTimestep, 0);
  graph.setTileMapping(complete, 0);

  auto cs = graph.addComputeSet("cs");
  auto vertex = graph.addVertex(cs, templateVertex("popnn::CTCMergeCandidates",
                                                   partialsType, UNSIGNED_INT));
  graph.setTileMapping(vertex, 0);

  graph.connect(vertex["beamAddend"], beamAddend.flatten());
  graph.connect(vertex["beamParent"], beamParent.flatten());
  graph.connect(vertex["lastBeamOutput"], lastBeamOutput);
  graph.connect(vertex["beamLength"], beamLength);

  graph.connect(vertex["currentTimestep"], currentTimestep);
  graph.connect(vertex["complete"], complete);

  graph.setInitialValue(vertex["beamwidth"], beamwidth);
  graph.setInitialValue(vertex["blankClass"], blankClass);

  Sequence uploadProg, downloadProg;
  std::vector<std::pair<std::string, char *>> tmap;

  // InOut
  auto rawCopyCandidates =
      createAndConnectCandidates(graph, vertex, "copyCandidate", partialsType,
                                 {}, uploadProg, downloadProg, tmap);

  // Inputs
  auto rawExtendCandidates = createAndConnectCandidates(
      graph, vertex, "extendCandidate", partialsType, {numExtendCandidates},
      uploadProg, downloadProg, tmap, false);

  // Extend candidate inputs
  std::vector<unsigned> extendCandidateParentIn{};
  std::vector<unsigned> extendCandidateAddendIn{};
  std::vector<float> extendCandidateBeamProbNonBlankIn{};
  std::vector<float> extendCandidateBeamProbBlankIn{};

  for (unsigned c = 0; c < numExtendCandidates; c++) {
    extendCandidateParentIn.push_back(extendCandidates[c].beam);
    extendCandidateAddendIn.push_back(extendCandidates[c].addend);
    extendCandidateBeamProbNonBlankIn.push_back(extendCandidates[c].pnb);
    extendCandidateBeamProbBlankIn.push_back(extendCandidates[c].pb);
  }

  copy(target, extendCandidateParentIn, UNSIGNED_INT,
       rawExtendCandidates.parent.get());
  copy(target, extendCandidateAddendIn, UNSIGNED_INT,
       rawExtendCandidates.addend.get());
  copy(target, extendCandidateBeamProbNonBlankIn, partialsType,
       rawExtendCandidates.probNonBlank.get());

  // Copy candidate inputs (A single candidate)
  std::vector<unsigned> copyCandidateParentIn = {copyCandidate.beam};
  std::vector<unsigned> copyCandidateAddendIn = {copyCandidate.addend};
  std::vector<float> copyCandidateBeamProbNonBlankIn = {copyCandidate.pnb};
  std::vector<float> copyCandidateBeamProbBlankIn = {copyCandidate.pb};
  std::vector<float> copyCandidateBeamProbTotalIn = {copyCandidate.pTotal};

  copy(target, copyCandidateParentIn, UNSIGNED_INT,
       rawCopyCandidates.parent.get());
  copy(target, copyCandidateAddendIn, UNSIGNED_INT,
       rawCopyCandidates.addend.get());
  copy(target, copyCandidateBeamProbNonBlankIn, partialsType,
       rawCopyCandidates.probNonBlank.get());
  copy(target, copyCandidateBeamProbBlankIn, partialsType,
       rawCopyCandidates.probBlank.get());
  copy(target, copyCandidateBeamProbTotalIn, partialsType,
       rawCopyCandidates.probTotal.get().get());

  std::unique_ptr<char[]> rawBeamAddend, rawBeamParent;

  rawBeamAddend = allocateHostMemoryForTensor(beamAddend, "beamAddend", graph,
                                              uploadProg, downloadProg, tmap);
  rawBeamParent = allocateHostMemoryForTensor(beamParent, "beamParent", graph,
                                              uploadProg, downloadProg, tmap);

  std::vector<unsigned> beamAddendIn{};
  std::vector<unsigned> beamParentIn{};
  for (unsigned t = 0; t < maxT; t++) {
    for (unsigned b = 0; b < beamwidth; b++) {
      beamAddendIn.push_back(beamHistory.symbols[b][t]);
      if (beamHistory.parents[b][t]) {
        beamParentIn.push_back(*beamHistory.parents[b][t]);
      } else {
        // TODO need this column in real impl?
        beamParentIn.push_back(std::numeric_limits<unsigned>::max());
      }
    }
  }

  copy(target, beamAddendIn, UNSIGNED_INT, rawBeamAddend.get());
  copy(target, beamParentIn, UNSIGNED_INT, rawBeamParent.get());

  OptionFlags engineOptions;
  if (profile) {
    engineOptions.set("debug.instrumentCompute", "true");
  }
  Sequence prog;
  prog.add(Execute(cs));
  Engine engine(graph, Sequence{uploadProg, prog, downloadProg}, engineOptions);
  attachStreams(engine, tmap);
  device.bind([&](const Device &d) {
    engine.load(d);
    engine.run();
  });

  // TODO partialsType == float
  std::vector<float> copyCandidateBeamProbBlankOut(1);
  std::vector<float> copyCandidateBeamProbNonBlankOut(1);
  std::vector<float> copyCandidateBeamProbTotalOut(1);
  std::vector<unsigned> copyCandidateAddendOut(1);

  // This is all that should have changed - either a merge happened or it didn't
  // and if so the copy candidate is updated with probabilities.
  // The vector of extend beams can be shared between multiple vertices so is
  // not changed and the copy candidate addend indicates which was merged into
  // the copy.
  copy(target, partialsType, rawCopyCandidates.probNonBlank.get(),
       copyCandidateBeamProbNonBlankOut);
  copy(target, partialsType, rawCopyCandidates.probBlank.get(),
       copyCandidateBeamProbBlankOut);
  copy(target, partialsType, rawCopyCandidates.probTotal.get().get(),
       copyCandidateBeamProbTotalOut);
  copy(target, UNSIGNED_INT, rawCopyCandidates.addend.get(),
       copyCandidateAddendOut);

  if (profile && deviceType != DeviceType::Cpu) {
    engine.printProfileSummary(std::cout,
                               OptionFlags{{"showExecutionSteps", "true"}});
  }
  // Return a vector of candidates for comparison: The copy candidate first
  // and the the extend candidates, omitting the one that was merged if a merge
  // happened

  std::vector<Candidate<float>> mergedCandidates;
  mergedCandidates.push_back(
      {copyCandidateParentIn[0], copyCandidateAddendIn[0],
       copyCandidateBeamProbNonBlankOut[0], copyCandidateBeamProbBlankOut[0],
       copyCandidateBeamProbTotalOut[0]});
  for (unsigned i = 0; i < numExtendCandidates; i++) {
    if (extendCandidateAddendIn[i] != copyCandidateAddendOut[0]) {
      mergedCandidates.push_back(
          {extendCandidateParentIn[i], extendCandidateAddendIn[i],
           extendCandidateBeamProbNonBlankIn[i],
           extendCandidateBeamProbBlankIn[i], extendCandidates[i].pTotal});
    }
  }
  return mergedCandidates;
}

template std::vector<Candidate<float>> runMergeCandidatesCodelet(
    Graph &graph, TestDevice &device, DeviceType deviceType, Type inType,
    Type partialsType, const std::vector<Candidate<float>> &extendCandidates,
    const Candidate<float> &copyCandidate, unsigned timestep,
    unsigned blankClass, const BeamHistory &beamHistory,
    const ArrayRef<unsigned> &outputLengths, unsigned lastBeamOutputSym,
    unsigned numClasses, bool profile);

} // namespace ctc
} // namespace poplibs_test

ext: cpp

max_stars_repo_name: graphcore/poplibs

max_issues_repo_name: graphcore/poplibs

max_forks_repo_name: graphcore/poplibs


text:
#include "tutorial_shared_path.h"
#include <igl/fast_winding_number.h>
#include <igl/read_triangle_mesh.h>
#include <igl/slice_mask.h>
#include <Eigen/Geometry>
#include <igl/octree.h>
#include <igl/barycenter.h>
#include <igl/knn.h>
#include <igl/random_points_on_mesh.h>
#include <igl/bounding_box_diagonal.h>
#include <igl/per_face_normals.h>
#include <igl/copyleft/cgal/point_areas.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/get_seconds.h>
#include <iostream>
#include <cstdlib>

int main(int argc, char *argv[])
{
  const auto time = [](std::function<void(void)> func)->double
  {
    const double t_before = igl::get_seconds();
    func();
    const double t_after = igl::get_seconds();
    return t_after-t_before;
  };

  Eigen::MatrixXd V;
  Eigen::MatrixXi F;
  igl::read_triangle_mesh(argc>1?argv[1]:TUTORIAL_SHARED_PATH "/bunny.off",V,F);
  // Sample mesh for point cloud
  Eigen::MatrixXd P,N;
  {
    Eigen::VectorXi I;
    Eigen::SparseMatrix<double> B;
    igl::random_points_on_mesh(10000,V,F,B,I);
    P = B*V;
    Eigen::MatrixXd FN;
    igl::per_face_normals(V,F,FN);
    N.resize(P.rows(),3);
    for(int p = 0;p<I.rows();p++)
    {
      N.row(p) = FN.row(I(p));
    }
  }
  // Build octree
  std::vector<std::vector<int > > O_PI;
  Eigen::MatrixXi O_CH;
  Eigen::MatrixXd O_CN;
  Eigen::VectorXd O_W;
  igl::octree(P,O_PI,O_CH,O_CN,O_W);
  Eigen::VectorXd A;
  {
    Eigen::MatrixXi I;
    igl::knn(P,20,O_PI,O_CH,O_CN,O_W,I);
    // CGAL is only used to help get point areas
    igl::copyleft::cgal::point_areas(P,I,N,A);
  }

  if(argc<=1)
  {
    // corrupt mesh
    Eigen::MatrixXd BC;
    igl::barycenter(V,F,BC);
    Eigen::MatrixXd OV = V;
    V.resize(F.rows()*3,3);
    for(int f = 0;f<F.rows();f++)
    {
      for(int c = 0;c<3;c++)
      {
        int v = f+c*F.rows();
        // random rotation about barycenter
        Eigen::AngleAxisd R(
          0.5*static_cast <double> (rand()) / static_cast <double> (RAND_MAX),
          Eigen::Vector3d::Random(3,1));
        V.row(v) = (OV.row(F(f,c))-BC.row(f))*R.matrix()+BC.row(f);
        F(f,c) = v;
      }
    }
  }

  // Generate a list of random query points in the bounding box
  Eigen::MatrixXd Q = Eigen::MatrixXd::Random(1000000,3);
  const Eigen::RowVector3d Vmin = V.colwise().minCoeff();
  const Eigen::RowVector3d Vmax = V.colwise().maxCoeff();
  const Eigen::RowVector3d Vdiag = Vmax-Vmin;
  for(int q = 0;q<Q.rows();q++)
  {
    Q.row(q) = (Q.row(q).array()*0.5+0.5)*Vdiag.array() + Vmin.array();
  }

  // Positions of points inside of point cloud P
  Eigen::MatrixXd QiP;
  {
    Eigen::MatrixXd O_CM;
    Eigen::VectorXd O_R;
    Eigen::MatrixXd O_EC;
    printf("  point cloud precomputation: %g secs\n",
      time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,2,O_CM,O_R,O_EC);}));
    Eigen::VectorXd WiP;
    printf("      point cloud evaluation: %g secs\n",
      time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,O_CM,O_R,O_EC,Q,2,WiP);}));
    igl::slice_mask(Q,WiP.array()>0.5,1,QiP);
  }

  // Positions of points inside of triangle soup (V,F)
  Eigen::MatrixXd QiV;
  {
    igl::FastWindingNumberBVH fwn_bvh;
    printf("triangle soup precomputation: %g secs\n",
      time([&](){igl::fast_winding_number(V.cast<float>(),F,2,fwn_bvh);}));
    Eigen::VectorXf WiV;
    printf("    triangle soup evaluation: %g secs\n",
      time([&](){igl::fast_winding_number(fwn_bvh,2,Q.cast<float>(),WiV);}));
    igl::slice_mask(Q,WiV.array()>0.5,1,QiV);
  }

  // Visualization
  igl::opengl::glfw::Viewer viewer;
  // For dislpaying normals as little line segments
  Eigen::MatrixXd PN(2*P.rows(),3);
  Eigen::MatrixXi E(P.rows(),2);
  const double bbd = igl::bounding_box_diagonal(V);
  for(int p = 0;p<P.rows();p++)
  {
    E(p,0) = 2*p;
    E(p,1) = 2*p+1;
    PN.row(E(p,0)) = P.row(p);
    PN.row(E(p,1)) = P.row(p)+bbd*0.01*N.row(p);
  }

  bool show_P = false;
  int show_Q = 0;

  int query_data = 0;
  viewer.data_list[query_data].set_mesh(V,F);
  viewer.data_list[query_data].clear();
  viewer.data_list[query_data].point_size = 2;
  viewer.append_mesh();
  int object_data = 1;
  viewer.data_list[object_data].set_mesh(V,F);
  viewer.data_list[object_data].point_size = 5;

  const auto update = [&]()
  {
    viewer.data_list[query_data].clear();
    switch(show_Q)
    {
      case 1:
        // show all Q
        viewer.data_list[query_data].set_points(Q,Eigen::RowVector3d(0.996078,0.760784,0.760784));
        break;
      case 2:
        // show all Q inside
        if(show_P)
        {
          viewer.data_list[query_data].set_points(QiP,Eigen::RowVector3d(0.564706,0.847059,0.768627));
        }else
        {
          viewer.data_list[query_data].set_points(QiV,Eigen::RowVector3d(0.564706,0.847059,0.768627));
        }
        break;
    }

    viewer.data_list[object_data].clear();
    if(show_P)
    {
      viewer.data_list[object_data].set_points(P,Eigen::RowVector3d(1,1,1));
      viewer.data_list[object_data].set_edges(PN,E,Eigen::RowVector3d(0.8,0.8,0.8));
    }else
    {
      viewer.data_list[object_data].set_mesh(V,F);
    }
  };

  viewer.callback_key_pressed =
    [&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)
  {
    switch(key)
    {
      default:
        return false;
      case '1':
        show_P = !show_P;
        break;
      case '2':
        show_Q = (show_Q+1) % 3;
        break;
    }
    update();
    return true;
  };

  update();
  viewer.launch();

}