Manifold_tracing.h

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Siargey Kachanovich
 *
 *    Copyright (C) 2019 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef MANIFOLD_TRACING_H_
#define MANIFOLD_TRACING_H_
 
#include <gudhi/IO/output_debug_traces_to_html.h>  // for DEBUG_TRACES
#include <gudhi/Coxeter_triangulation/Query_result.h>
 
#include <boost/functional/hash.hpp>
 
#include <Eigen/Dense>
 
#include <queue>
#include <unordered_map>
 
namespace Gudhi {
 
namespace coxeter_triangulation {
 
template <class Triangulation_>
class Manifold_tracing {
 public:
  using Simplex_handle = typename Triangulation_::Simplex_handle;
 
  struct Simplex_hash {
    typedef Simplex_handle argument_type;
    typedef std::size_t result_type;
    result_type operator()(const argument_type& s) const noexcept {
      return boost::hash<typename Simplex_handle::Vertex>()(s.vertex());
    }
  };
 
 public:
  typedef std::unordered_map<Simplex_handle, Eigen::VectorXd, Simplex_hash> Out_simplex_map;
 
  template <class Point_range, class Intersection_oracle>
  void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation,
                                  const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) {
    std::size_t cod_d = oracle.cod_d();
    std::queue<Simplex_handle> queue;
 
    for (const auto& p : seed_points) {
      Simplex_handle full_simplex = triangulation.locate_point(p);
      for (Simplex_handle face : full_simplex.face_range(cod_d)) {
        Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
        if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) {
#ifdef DEBUG_TRACES
          mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
#endif
          queue.emplace(face);
          break;
        }
      }
    }
 
    while (!queue.empty()) {
      Simplex_handle s = queue.front();
      queue.pop();
      for (auto cof : s.coface_range(cod_d + 1)) {
        for (auto face : cof.face_range(cod_d)) {
          Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
          if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
        }
      }
    }
  }
 
  template <class Point_range, class Intersection_oracle>
  void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation,
                                  const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map,
                                  Out_simplex_map& boundary_simplex_map) {
    std::size_t cod_d = oracle.cod_d();
    std::queue<Simplex_handle> queue;
 
    for (const auto& p : seed_points) {
      Simplex_handle full_simplex = triangulation.locate_point(p);
      for (Simplex_handle face : full_simplex.face_range(cod_d)) {
        auto qr = oracle.intersects(face, triangulation);
#ifdef DEBUG_TRACES
        mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
#endif
        if (qr.success) {
          if (oracle.lies_in_domain(qr.intersection, triangulation)) {
            if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
          } else {
            for (Simplex_handle cof : face.coface_range(cod_d + 1)) {
              auto qrb = oracle.intersects_boundary(cof, triangulation);
#ifdef DEBUG_TRACES
              mt_seed_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
#endif
              if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection);
            }
          }
          // break;
        }
      }
    }
 
    while (!queue.empty()) {
      Simplex_handle s = queue.front();
      queue.pop();
      for (auto cof : s.coface_range(cod_d + 1)) {
        for (auto face : cof.face_range(cod_d)) {
          auto qr = oracle.intersects(face, triangulation);
#ifdef DEBUG_TRACES
          mt_inserted_list.push_back(MT_inserted_info(qr, face, false));
#endif
          if (qr.success) {
            if (oracle.lies_in_domain(qr.intersection, triangulation)) {
              if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
            } else {
              auto qrb = oracle.intersects_boundary(cof, triangulation);
#ifdef DEBUG_TRACES
              mt_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
#endif
              if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection);
            }
          }
        }
      }
    }
  }
 
  Manifold_tracing() {}
};
 
template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map>
void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation,
                                const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) {
  Manifold_tracing<Triangulation> mt;
  mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map);
}
 
template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map>
void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation,
                                const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map,
                                Out_simplex_map& boundary_simplex_map) {
  Manifold_tracing<Triangulation> mt;
  mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map);
}
 
}  // namespace coxeter_triangulation
 
}  // namespace Gudhi
 
#endif