Program Listing for File triangulation.hpp

↰ Return to documentation for file (src/triangulation.hpp)

/*
// Copyright © 2019,2020 Greg Tucker <greg.tucker@stfc.ac.uk>
//
// This file is part of brille.
//
// brille is free software: you can redistribute it and/or modify it under the
// terms of the GNU Affero General Public License as published by the Free
// Software Foundation, either version 3 of the License, or (at your option)
// any later version.
//
// brille is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Affero General Public License for more details.
//
// You should have received a copy of the GNU Affero General Public License
// along with brille. If not, see <https://www.gnu.org/licenses/>.
*/

#ifndef BRILLE_TRIANGULATION_HPP_
#define BRILLE_TRIANGULATION_HPP_
#include <vector>
#include <array>
#include <cassert>
#include <algorithm>
#include "array_latvec.hpp" // defines bArray
#include "tetgen.h"
#include "debug.hpp"
#include "balltrellis.hpp"
#include "approx.hpp"
namespace brille {

template<class T, size_t N> static size_t find_first(const std::array<T,N>& x, const T val){
  auto at = std::find(x.begin(), x.end(), val);
  if (at == x.end()) throw std::logic_error("Value not found?!");
  return std::distance(x.begin(), at);
}

// // In case we need to store the input polyhedron and mesh-refining parameters:
// class TetgenInput{
//   bArray<double> _vertices;
//   std::vector<std::vector<int>> _vertices_per_facet;
//   double _max_cell_size;
//   double _min_dihedral_angle;
//   double _max_dihedral_angle;
//   double _radius_edge_ratio;
//   int _max_mesh_points;
// public:
//   TetgenInput(
//     const bArray<double>& v, const std::vector<std::vector<int>>& f,
//     const double mcs=-1.0, const double nda=-1.0, const double xda=-1.0,
//     const double rer=-1.0, const int mmp=-1):
//       _vertices(v), _vertices_per_facet(f), _max_cell_size(mcs),
//       _min_dihedral_angle(nda), _max_dihedral_angle(xda),
//       _radius_edge_ratio(rer), _max_mesh_points(mmp) {}
//   const bArray<double>& vertices() const { return _vertices; }
//   const std::vector<std::vector<int>>& vertices_per_facet() const {return _vertices_per_facet; }
//   double max_cell_size() const { return _max_cell_size; }
//   double max_dihedral_angle() const { return _max_dihedral_angle; }
//   double min_dihedral_angle() const { return _min_dihedral_angle; }
//   double radius_edge_ratio() const { return _radius_edge_ratio; }
//   int max_mesh_points() const { return _max_mesh_points; }
// };

// template<class T, template<class> class L, typename=typename std::enable_if<std::is_base_of<bArray<T>,L<T>>::value>::type>
class TetTri{
  // L<T> vertex_positions; // (nVertices, 3), bArray<T>, LQVec<T>, or LDVec<T>
  size_t nVertices;
  size_t nTetrahedra;
  bArray<double> vertex_positions; // (nVertices, 3)
  bArray<size_t> vertices_per_tetrahedron; // (nTetrahedra, 4)
  std::vector<std::vector<size_t>> tetrahedra_per_vertex; // (nVertices,)(1+,)
  std::vector<std::vector<size_t>> neighbours_per_tetrahedron; // (nTetrahedra,)(1+,)
  //tetgenio tgsource; // we need to store the output of tetgen so that we can refine the mesh
  // BallNode tetrahedraTree;
  // std::vector<BallLeaf> leaves;
  Trellis tetrahedraTrellis;
  std::vector<TrellisLeaf> leaves;
public:
  size_t number_of_vertices(void) const {return nVertices;}
  size_t number_of_tetrahedra(void) const {return nTetrahedra;}
  const bArray<double>& get_vertex_positions(void) const {return vertex_positions;}
  const bArray<size_t>& get_vertices_per_tetrahedron(void) const {return vertices_per_tetrahedron;}

  TetTri(void): nVertices(0), nTetrahedra(0), vertex_positions({0u,3u}), vertices_per_tetrahedron({0u,4u}){}
  TetTri(const tetgenio& tgio, const double fraction): vertex_positions({0u,3u}), vertices_per_tetrahedron({0u,4u}){ //, tgsource(tgio){
    nVertices = static_cast<size_t>(tgio.numberofpoints);
    nTetrahedra = static_cast<size_t>(tgio.numberoftetrahedra);
    // copy-over all vertex positions:
    vertex_positions.resize(nVertices);
    for (size_t i=0; i<nVertices; ++i)
    for (size_t j=0; j<3u; ++j)
    vertex_positions.val(i,j) = tgio.pointlist[3*i+j];
    // copy-over all tetrahedron vertex indices
    vertices_per_tetrahedron.resize(nTetrahedra);
    for (size_t i=0; i<nTetrahedra; ++i)
    for (size_t j=0; j<4u; ++j)
    vertices_per_tetrahedron.val(i,j) = static_cast<size_t>(tgio.tetrahedronlist[i*tgio.numberofcorners+j]);
    // for (size_t i=0; i<nTetrahedra; ++i)
    // Construct the tetrahedra per vertex vector of vectors
    tetrahedra_per_vertex.resize(nVertices);
    for (size_t i=0; i<nVertices; ++i)
    for (size_t j=0; j<nTetrahedra; ++j)
    for (size_t k=0; k<4u; ++k)
    if (vertices_per_tetrahedron.val(j,k)==i)
      tetrahedra_per_vertex[i].push_back(j);
    // Construct the neighbours per tetrahedron vector of vectors
    neighbours_per_tetrahedron.resize(nTetrahedra);
    for (size_t i=0; i<nTetrahedra; ++i)
    for (size_t j=0; j<4u; ++j)
    if (tgio.neighborlist[i*4u+j] >= 0)
      neighbours_per_tetrahedron[i].push_back(static_cast<size_t>(tgio.neighborlist[i*4u+j]));
    // ensure that all tetrahedra have positive (orient3d) volume
    this->correct_tetrahedra_vertex_ordering();
    // construct the tree for faster locating:
    this->make_balltree(fraction);
    // Create a string full of object information:
  }
  std::string to_string(void) const {
    std::string str;
    str  = std::to_string(nVertices) + " vertices";
    str += " in " + std::to_string(nTetrahedra) + " tetrahedra";
    str += " with a Trellis[" + tetrahedraTrellis.to_string() + "]";
    return str;
  }
  // emulate the locate functionality of CGAL -- for which we need an enumeration
  enum Locate_Type{ VERTEX, EDGE, FACET, CELL, OUTSIDE_CONVEX_HULL };
  size_t old_locate(const bArray<double>& x, Locate_Type& type, size_t& v0, size_t& v1) const{
    std::vector<size_t> v;
    std::vector<double> w;
    size_t idx = this->locate(x, v, w);
    switch (v.size()){
      case 0:
        type = OUTSIDE_CONVEX_HULL;
        return 0;
      case 1:
        type = VERTEX;
        v0 = v[0];
        break;
      case 2:
        type = EDGE;
        v0 = v[0];
        v1 = v[1];
        break;
      case 3:
        type = FACET;
        for (size_t i=0; i<4u; ++i){
          v0 = vertices_per_tetrahedron.val(idx, i);
          if (std::find(v.begin(), v.end(), v0) == v.end()) break;
        }
        break;
      case 4:
        type = CELL;
        break;
      default:
        throw std::logic_error("Intentionally unreachable switch statement..");
    }
    return idx;
  }
  // Make a new locator which slots into the interpolation routine more easily
  // size_t locate(const bArray<double>& x, std::vector<size_t>& v, std::vector<double>& w) const {
  //   if (x.numel() != 3u || x.size() != 1u)
  //     throw std::runtime_error("locate requires a single 3-element vector.");
  //   std::array<double,4> ws;
  //   v.clear();
  //   w.clear(); // make sure w is back to zero-elements
  //
  //   size_t tet_idx;
  //   for (tet_idx=0; tet_idx < nTetrahedra; ++tet_idx) if (this->might_contain(tet_idx, x)){
  //     this->weights(tet_idx, x, ws);
  //     // if all weights are greater or equal to ~zero, we can use this tetrahedron
  //     if (std::all_of(ws.begin(), ws.end(), [](double z){return z>0. || brille::approx::scalar(z,0.);})){
  //       for (size_t i=0; i<4u; ++i) if (!brille::approx::scalar(ws[i], 0.)){
  //         v.push_back(vertices_per_tetrahedron.getvalue(tet_idx, i));
  //         w.push_back(ws[i]);
  //       }
  //       break;
  //     }
  //   }
  //   return tet_idx;
  // }
  size_t locate(const bArray<double>& x, std::vector<size_t>& v, std::vector<double>& w) const {
    if (x.ndim()!=2u || x.size(0)!=1u || x.size(1)!=3u)
      throw std::runtime_error("locate requires a single 3-element vector.");
    std::array<double,4> ws;
    v.clear();
    w.clear(); // make sure w is back to zero-elements

    // verbose_update("Find which tetrahedra might contain the point ",x.to_string(0));
    // std::vector<size_t> tets_to_check = tetrahedraTree.all_containing_leaf_indexes(x);
    // verbose_update("The tree would have us search ",tets_to_check);
    // for (size_t tet_idx: tets_to_check) if (this->unsafe_might_contain(tet_idx, x)){

    for (size_t node: tetrahedraTrellis.nodes_to_search(x)) // returns a distance-sorted list of nodes which might have a containing leaf
    for (auto leaf: tetrahedraTrellis.node_leaves(node))
    if (this->unsafe_might_contain(leaf.index(), x)){
      this->weights(leaf.index(), x, ws);
      // if all weights are greater or equal to ~zero, we can use this tetrahedron
      if (std::all_of(ws.begin(), ws.end(), [](double z){return z>0. || brille::approx::scalar(z,0.);})){
        for (size_t i=0; i<4u; ++i) if (!brille::approx::scalar(ws[i], 0.)){
          v.push_back(vertices_per_tetrahedron[leaf.index(), i]);
          w.push_back(ws[i]);
        }
        return leaf.index();
      }
    }
    throw std::runtime_error("The containing tetrahedra was not found?!");
    return nTetrahedra;
  }
  // and a special version which doesn't return the weights
  size_t locate(const bArray<double>& x, std::vector<size_t>& v) const{
    std::vector<double> w;
    return this->locate(x, v, w);
  }

  // /* If the following function is more useful than the preceeding, it could be
  //    advantageous to replicate the above code in this function's for loop.
  //    Either way, this should probably be parallelised with OpenMP.
  // */
  // std::vector<std::vector<size_t>> locate_all_for_interpolation(const bArray<double>& x) const {
  //   if (x.numel()!=3u){
  //     std::string msg = "locate_all requires 3-element vector(s)";
  //     throw std::runtime_error(msg);
  //   }
  //   std::vector<std::vector<size_t>> all_idx(x.size());
  //   for (size_t i=0; i<x.size(); ++i)
  //     all_idx[i] = this->locate_for_interpolation(x.extract(i));
  //   return all_idx;
  // }

  /* Given a vertex in the mesh, return a vector of all of the other vertices to
     which it is connected.
  */
  std::vector<size_t> neighbours(const bArray<double>& x) const {
    std::vector<size_t> v;
    this->locate(x, v);
    if (v.size() != 1u){
      std::string msg = "The provided point is not a mesh vertex.";
      throw std::runtime_error(msg);
    }
    return this->neighbours(v[0]);
  }
  std::vector<size_t> neighbours(const size_t vert) const {
    if (vert >= this->nVertices){
      std::string msg = "The provided vertex index is out of bounds";
      throw std::out_of_range(msg);
    }
    std::vector<size_t> n;
    size_t v;
    for (size_t t: this->tetrahedra_per_vertex[vert])
    // for (size_t v: this->vertices_per_tetrahedron[t]) // would work if vertices_per_tetrahedron was a vector<array<size_t,4>>
    for (size_t j=0; j<4u; ++j){
      v = this->vertices_per_tetrahedron.val(t, j);
      if ( v!= vert && std::find(n.begin(), n.end(), v) == n.end() ) n.push_back(v);
    }
    return n;
  }
  double volume(const size_t tet) const {
    double v;
    const ind_t* i = vertices_per_tetrahedron.ptr(tet);
    v = orient3d(
      vertex_positions.ptr(i[0]), vertex_positions.ptr(i[1]),
      vertex_positions.ptr(i[2]), vertex_positions.ptr(i[3]))/6.0;return v;
  }
  bool might_contain(const size_t tet, const bArray<double>& x) const {
    if (x.ndim()!=2u || x.size(0)!=1u || x.size(1)!=3u)
      throw std::runtime_error("x must be a single 3-vector");
    if (tet >= nTetrahedra) return false;
    return this->unsafe_might_contain(tet, x);
  }
protected:
  bool unsafe_might_contain(const size_t tet, const bArray<double>& x) const {
    return leaves[tet].fuzzy_contains(x);
  }
  void weights(const size_t tet, const bArray<double>& x, std::array<double,4>& w) const {
    double vol6 = 6.0*this->volume(tet);
    const ind_t* i = vertices_per_tetrahedron.ptr(tet);
    const auto& p{vertex_positions};
    w[0] = orient3d(x.ptr(0),    p.ptr(i[1]), p.ptr(i[2]), p.ptr(i[3]) )/vol6;
    w[1] = orient3d(p.ptr(i[0]), x.ptr(0),    p.ptr(i[2]), p.ptr(i[3]) )/vol6;
    w[2] = orient3d(p.ptr(i[0]), p.ptr(i[1]), x.ptr(0),    p.ptr(i[3]) )/vol6;
    w[3] = orient3d(p.ptr(i[0]), p.ptr(i[1]), p.ptr(i[2]), x.ptr(0)    )/vol6;
  }
  void correct_tetrahedra_vertex_ordering(void){
    for (size_t i=0; i<nTetrahedra; ++i)
    if (std::signbit(this->volume(i))) // the volume of tetrahedra i is negative
    vertices_per_tetrahedron.swap(i, 0,1); // swap two vertices to switch sign
  }
  void make_balltree(const double fraction){
    // Construct a vector of BallLeaf objects for each tetrahedra:
    // std::vector<BallLeaf> leaves;
    leaves.clear();
    std::array<double,3> centre;
    double radius;
    verbose_update("Pull together the circumsphere information for all tetrahedra");
    tetgenmesh tgm; // to get access to circumsphere
    const auto& p{vertex_positions};
    for (size_t i=0; i<nTetrahedra; ++i){
      const ind_t* v = vertices_pertetrahedron.ptr(i);
      // use tetgen's circumsphere to find the centre and radius for each tetrahedra
      tgm.circumsphere(p.ptr(v[0]), p.ptr(v[1]), p.ptr(v[2]), p.ptr(v[3]), centre.data(), &radius);
      // leaves.push_back(BallLeaf(centre, radius, i));
      leaves.push_back(TrellisLeaf(centre, radius, i));
    }
    // // construct the full tree structure at once using an algorithm similar to
    // // Omohundro's Kd algorithm in 'Five Balltree Construction Algorithms'
    // if (nTetrahedra > 0){
    //   // we will construct a binary tree. If we have N leaves to place at the
    //   // end of the branches, they can occupy as many as N branches.
    //   // A tree with N final branches has log₂(N) levels, so we will never need
    //   // to exceed this.
    //   // If we reduce the number of branchings we increase the number of leaves
    //   // per terminal branch. There must be some trade-off between tree-granularity
    //   // and tree-size for how quickly a containing tetrahedra can be located.
    //   size_t maximum_branchings = static_cast<size_t>(std::log2(nTetrahedra))-2;
    //   maximum_branchings = 3;
    //   verbose_update("Construct a tree for tetrahedra location with up to ",maximum_branchings," branchings.");
    //   tetrahedraTree = construct_balltree(leaves, maximum_branchings);
    //   verbose_update("Tree for locating tetrahedra now exists");
    //   // info_update("Tetrahedra tree:\n",tetrahedraTree.to_string());
    // }
    tetrahedraTrellis = construct_trellis(leaves, fraction);
  }
};

template <class T>
TetTri
triangulate(
  const bArray<T>& verts,
  const std::vector<std::vector<int>>& vpf,
  const double max_cell_size=-1.0,
  const double min_dihedral=-1.0,
  const double max_dihedral=-1.0,
  const double radius_edge_ratio=-1.0,
  const int max_mesh_points=-1,
  const double fraction=1.0
) {
  assert(verts.ndim()==2 && verts.size(1)==3);// otherwise we can't make a 3-D triangulation
  // create the tetgenbehavior object which contains all options/switches for tetrahedralize
  verbose_update("Creating `tetgenbehavior` object");
  tetgenbehavior tgb;
  tgb.plc = 1; // we will always tetrahedralize a piecewise linear complex
  tgb.quality = 1; // we will (almost) always improve the tetrahedral mesh
  tgb.neighout = 1; // we *need* the neighbour information to be stored into tgo.
  if (max_cell_size > 0){
    // volume constraint with a specified size
    tgb.fixedvolume = 1;
    tgb.maxvolume = max_cell_size;
  } else{
    //volume constraint without a specified size?
    tgb.varvolume = 1;
  }
  if (max_mesh_points>0) tgb.steinerleft = max_mesh_points;
  if (radius_edge_ratio>0) tgb.minratio = radius_edge_ratio;
  if (min_dihedral>0) tgb.mindihedral = min_dihedral;
  if (max_dihedral>0) tgb.optmaxdihedral = max_dihedral;
  #ifndef VERBOSE_MESHING
  tgb.quiet = 1;
  #endif
  #ifdef VERBOSE_MESHING
  tgb.verbose = 10000;
  #endif

  // make the input and output tetgenio objects and fill the input with our polyhedron
  verbose_update("Creating input and output `tetgenio` objects");
  tetgenio tgi, tgo;
  // we have to handle initializing points/facets, but tetgenio has a destructor
  // which handles deleting all non-NULL fields.
  verbose_update("Initialize and fill the input object's pointlist parameter");
  tgi.numberofpoints = static_cast<int>(verts.size(0));
  tgi.pointlist = new double[3*tgi.numberofpoints];
  tgi.pointmarkerlist = new int[tgi.numberofpoints];
  //tgi.point2tetlist = new int[tgi.numberofpoints];
  int idx=0;
  for (size_t i=0; i<verts.size(0); ++i){
    tgi.pointmarkerlist[i] = static_cast<int>(i);
    for (size_t j=0; j<verts.size(1); ++j)
      tgi.pointlist[idx++] = verts.val(i,j);
  }
  verbose_update("Initialize and fill the input object's facetlist parameter");
  tgi.numberoffacets = static_cast<int>(vpf.size());
  tgi.facetlist = new tetgenio::facet[tgi.numberoffacets];
  tgi.facetmarkerlist = new int[tgi.numberoffacets];
  for (size_t i=0; i<vpf.size(); ++i){
    tgi.facetmarkerlist[i] = static_cast<int>(i);
    tgi.facetlist[i].numberofpolygons = 1;
    tgi.facetlist[i].polygonlist = new tetgenio::polygon[1];
    tgi.facetlist[i].polygonlist[0].numberofvertices = static_cast<int>(vpf[i].size());
    tgi.facetlist[i].polygonlist[0].vertexlist = new int[tgi.facetlist[i].polygonlist[0].numberofvertices];
    for (size_t j=0; j<vpf[i].size(); ++j)
      tgi.facetlist[i].polygonlist[0].vertexlist[j] = vpf[i][j];
  }
  // The input is now filled with the piecewise linear complex information.
  // so we can call tetrahedralize:
  verbose_update("Calling tetgen::tetrahedralize");
  try {
      tetrahedralize(&tgb, &tgi, &tgo);
  } catch (const std::logic_error& e) {
    std::string msg = "tetgen threw a logic error with message\n" + std::string(e.what());
    throw std::runtime_error(msg);
  } catch (const std::runtime_error& e) {
    std::string msg = "tetgen threw a runtime_error with message\n" + std::string(e.what());
    throw std::runtime_error(msg);
  } catch (...) {
    std::string msg = "tetgen threw an undetermined error";
    throw std::runtime_error(msg);
  }
  verbose_update("Constructing TetTri object");
  return TetTri(tgo, fraction);
}

} // end namespace brille
#endif // _TRIANGULATION_H_