diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h new file mode 100644 index 00000000000..8f2279fb3b8 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h @@ -0,0 +1,111 @@ +#ifndef CGAL_CARTESIAN_GRID_3_H +#define CGAL_CARTESIAN_GRID_3_H + +#include +#include + +#include + +namespace CGAL { + +// TODO: not sure if anything other than float works +template +class Cartesian_grid_3 { +public: + typedef typename Traits::FT FT; + +public: + Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox) { + create_image(xdim, ydim, zdim, bbox); + } + + Cartesian_grid_3(const Image_3 &image) : grid(image) {} + + FT value(const std::size_t x, const std::size_t y, const std::size_t z) const { + return grid.value(x, y, z); + } + + FT &value(const std::size_t x, const std::size_t y, const std::size_t z) { + FT *data = (FT *)grid.image()->data; + return data[(z * ydim() + y) * xdim() + x]; + } + + std::size_t xdim() const { + return grid.xdim(); + } + std::size_t ydim() const { + return grid.ydim(); + } + std::size_t zdim() const { + return grid.zdim(); + } + + // TODO: better return types + double voxel_x() const { + return grid.vx(); + } + double voxel_y() const { + return grid.vy(); + } + double voxel_z() const { + return grid.vz(); + } + + float offset_x() const { + return grid.tx(); + } + float offset_y() const { + return grid.ty(); + } + float offset_z() const { + return grid.tz(); + } + +private: + void create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox); + +private: + Image_3 grid; +}; + +template +void Cartesian_grid_3::create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, + const Bbox_3 &bbox) { + + WORD_KIND wordkind; + if (std::is_floating_point::value) + wordkind = WK_FLOAT; + else + wordkind = WK_FIXED; + + SIGN sign; + if (std::is_signed::value) + sign = SGN_SIGNED; + else + sign = SGN_UNSIGNED; + + const double vx = bbox.x_span() / xdim; + const double vy = bbox.y_span() / ydim; + const double vz = bbox.z_span() / zdim; + + _image *im = _createImage(xdim, ydim, zdim, + 1, // vectorial dimension + vx, vy, vz, // voxel size + sizeof(FT), // image word size in bytes + wordkind, // image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN + sign); // image word sign + + if (im == nullptr || im->data == nullptr) { + throw std::bad_alloc(); // TODO: idk? + } + + im->tx = bbox.xmin(); + im->ty = bbox.ymin(); + im->tz = bbox.zmin(); + + grid = Image_3(im, Image_3::OWN_THE_DATA); +} + +} // end namespace CGAL + +#endif // CGAL_CARTESIAN_GRID_3_H \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h index e46325d0a15..74c09114b6f 100644 --- a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h @@ -1,31 +1,22 @@ #ifndef CGAL_CARTESIAN_GRID_DOMAIN_H #define CGAL_CARTESIAN_GRID_DOMAIN_H +#include +#include +#include #include -#include "Cartesian_grid_3.h" -#include "Isosurfacing_3/internal/Tables.h" - namespace CGAL { namespace Isosurfacing { template -class Cartesian_grid_domain { +class Cartesian_grid_domain : public Cartesian_topology_base { public: typedef GeomTraits Geom_traits; typedef typename Geom_traits::FT FT; typedef typename Geom_traits::Point_3 Point; typedef typename Geom_traits::Vector_3 Vector; - typedef std::array Vertex_handle; - typedef std::array Edge_handle; - typedef std::array Cell_handle; - - typedef std::array Edge_vertices; - typedef std::array Cells_incident_to_edge; - typedef std::array Cell_vertices; - typedef std::array Cell_edges; - public: Cartesian_grid_domain(const Cartesian_grid_3& grid) : grid(&grid) {} @@ -50,61 +41,19 @@ public: return grid->value(v[0], v[1], v[2]); } - Edge_vertices edge_vertices(const Edge_handle& e) const { - Edge_vertices ev; - ev[0] = {e[0], e[1], e[2]}; - ev[1] = {e[0], e[1], e[2]}; - ev[1][e[3]] += 1; - return ev; - } - - Cells_incident_to_edge cells_incident_to_edge(const Edge_handle& e) const { - const int local = internal::Cube_table::edge_store_index[e[3]]; - auto neighbors = internal::Cube_table::edge_to_voxel_neighbor[local]; - - Cells_incident_to_edge cite; - for (std::size_t i = 0; i < cite.size(); i++) { - for (std::size_t j = 0; j < cite[i].size(); j++) { - cite[i][j] = e[j] + neighbors[i][j]; - } - } - return cite; - } - - Cell_vertices cell_vertices(const Cell_handle& v) const { - Cell_vertices cv; - for (std::size_t i = 0; i < cv.size(); i++) { - for (std::size_t j = 0; j < v.size(); j++) { - cv[i][j] = v[j] + internal::Cube_table::local_vertex_position[i][j]; - } - } - return cv; - } - - Cell_edges cell_edges(const Cell_handle& v) const { - Cell_edges ce; - for (std::size_t i = 0; i < ce.size(); i++) { - for (std::size_t j = 0; j < v.size(); j++) { - ce[i][j] = v[j] + internal::Cube_table::global_edge_id[i][j]; - } - ce[i][3] = internal::Cube_table::global_edge_id[i][3]; - } - return ce; + template + void iterate_vertices(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_vertices_base(f, tag, grid->xdim(), grid->ydim(), grid->zdim()); } template - void iterate_vertices(Functor& f, Sequential_tag) const { - const std::size_t size_x = grid->xdim(); - const std::size_t size_y = grid->ydim(); - const std::size_t size_z = grid->zdim(); + void iterate_edges(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_edges_base(f, tag, grid->xdim(), grid->ydim(), grid->zdim()); + } - for (std::size_t x = 0; x < size_x - 1; x++) { - for (std::size_t y = 0; y < size_y - 1; y++) { - for (std::size_t z = 0; z < size_z - 1; z++) { - f({x, y, z}); - } - } - } + template + void iterate_cells(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_cells_base(f, tag, grid->xdim(), grid->ydim(), grid->zdim()); } #ifdef CGAL_LINKED_WITH_TBB @@ -126,26 +75,7 @@ public: tbb::parallel_for(tbb::blocked_range(0, size_x - 1), iterator); } -#endif // CGAL_LINKED_WITH_TBB - template - void iterate_edges(Functor& f, Sequential_tag) const { - const std::size_t size_x = grid->xdim(); - const std::size_t size_y = grid->ydim(); - const std::size_t size_z = grid->zdim(); - - for (std::size_t x = 0; x < size_x - 1; x++) { - for (std::size_t y = 0; y < size_y - 1; y++) { - for (std::size_t z = 0; z < size_z - 1; z++) { - f({x, y, z, 0}); - f({x, y, z, 1}); - f({x, y, z, 2}); - } - } - } - } - -#ifdef CGAL_LINKED_WITH_TBB template void iterate_edges(Functor& f, Parallel_tag) const { const std::size_t size_x = grid->xdim(); @@ -166,24 +96,7 @@ public: tbb::parallel_for(tbb::blocked_range(0, size_x - 1), iterator); } -#endif // CGAL_LINKED_WITH_TBB - template - void iterate_cells(Functor& f, Sequential_tag) const { - const std::size_t size_x = grid->xdim(); - const std::size_t size_y = grid->ydim(); - const std::size_t size_z = grid->zdim(); - - for (std::size_t x = 0; x < size_x - 1; x++) { - for (std::size_t y = 0; y < size_y - 1; y++) { - for (std::size_t z = 0; z < size_z - 1; z++) { - f({x, y, z}); - } - } - } - } - -#ifdef CGAL_LINKED_WITH_TBB template void iterate_cells(Functor& f, Parallel_tag) const { const std::size_t size_x = grid->xdim(); diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain_old.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain_old.h new file mode 100644 index 00000000000..a9c55ab41ad --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain_old.h @@ -0,0 +1,48 @@ +#ifndef CGAL_CARTESIAN_GRID_DOMAIN_OLD_H +#define CGAL_CARTESIAN_GRID_DOMAIN_OLD_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Cartesian_grid_domain_old { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point_3; + +public: + Cartesian_grid_domain_old(const Cartesian_grid_3& grid) : grid(&grid) {} + + std::size_t size_x() const { + return grid->xdim(); + } + std::size_t size_y() const { + return grid->ydim(); + } + std::size_t size_z() const { + return grid->zdim(); + } + + Point_3 position(const std::size_t x, const std::size_t y, const std::size_t z) const { + const FT vx = grid->voxel_x(); + const FT vy = grid->voxel_y(); + const FT vz = grid->voxel_z(); + + return Point_3(x * vx + grid->offset_x(), y * vy + grid->offset_y(), z * vz + grid->offset_z()); + } + + FT value(const std::size_t x, const std::size_t y, const std::size_t z) const { + return grid->value(x, y, z); + } + +private: + const Cartesian_grid_3* grid; +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_CARTESIAN_GRID_DOMAIN_OLD_H diff --git a/Isosurfacing_3/include/CGAL/Cartesian_topology_base.h b/Isosurfacing_3/include/CGAL/Cartesian_topology_base.h new file mode 100644 index 00000000000..42be536ed95 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Cartesian_topology_base.h @@ -0,0 +1,117 @@ +#ifndef CGAL_CARTESIAN_TOPOLOGY_BASE_H +#define CGAL_CARTESIAN_TOPOLOGY_BASE_H + +#include +#include +#include + +#include + +namespace CGAL { +namespace Isosurfacing { + +class Cartesian_topology_base { +public: + typedef std::array Vertex_handle; + typedef std::array Edge_handle; + typedef std::array Cell_handle; + + static constexpr Cell_type CELL_TYPE = CUBICAL_CELL; + static constexpr std::size_t VERTICES_PER_CELL = 8; + static constexpr std::size_t EDGES_PER_CELL = 12; + + typedef std::array Edge_vertices; + typedef std::array Cells_incident_to_edge; + typedef std::array Cell_vertices; + typedef std::array Cell_edges; + +public: + Edge_vertices edge_vertices(const Edge_handle& e) const { + Edge_vertices ev; + ev[0] = {e[0], e[1], e[2]}; + ev[1] = {e[0], e[1], e[2]}; + ev[1][e[3]] += 1; + return ev; + } + + Cells_incident_to_edge cells_incident_to_edge(const Edge_handle& e) const { + const int local = internal::Cube_table::edge_store_index[e[3]]; + auto neighbors = internal::Cube_table::edge_to_voxel_neighbor[local]; + + Cells_incident_to_edge cite; + for (std::size_t i = 0; i < cite.size(); i++) { + for (std::size_t j = 0; j < cite[i].size(); j++) { + cite[i][j] = e[j] + neighbors[i][j]; + } + } + return cite; + } + + Cell_vertices cell_vertices(const Cell_handle& v) const { + Cell_vertices cv; + for (std::size_t i = 0; i < cv.size(); i++) { + for (std::size_t j = 0; j < v.size(); j++) { + cv[i][j] = v[j] + internal::Cube_table::local_vertex_position[i][j]; + } + } + return cv; + } + + Cell_edges cell_edges(const Cell_handle& v) const { + Cell_edges ce; + for (std::size_t i = 0; i < ce.size(); i++) { + for (std::size_t j = 0; j < v.size(); j++) { + ce[i][j] = v[j] + internal::Cube_table::global_edge_id[i][j]; + } + ce[i][3] = internal::Cube_table::global_edge_id[i][3]; + } + return ce; + } + +protected: + template + void iterate_vertices_base(Functor& f, Sequential_tag, const std::size_t size_x, const std::size_t size_y, + const std::size_t size_z) const { + + for (std::size_t x = 0; x < size_x; x++) { + for (std::size_t y = 0; y < size_y; y++) { + for (std::size_t z = 0; z < size_z; z++) { + f({x, y, z}); + } + } + } + } + + template + void iterate_edges_base(Functor& f, Sequential_tag, const std::size_t size_x, const std::size_t size_y, + const std::size_t size_z) const { + + for (std::size_t x = 0; x < size_x - 1; x++) { + for (std::size_t y = 0; y < size_y - 1; y++) { + for (std::size_t z = 0; z < size_z - 1; z++) { + f({x, y, z, 0}); + f({x, y, z, 1}); + f({x, y, z, 2}); + } + } + } + } + + template + void iterate_cells_base(Functor& f, Sequential_tag, const std::size_t size_x, const std::size_t size_y, + const std::size_t size_z) const { + + for (std::size_t x = 0; x < size_x - 1; x++) { + for (std::size_t y = 0; y < size_y - 1; y++) { + for (std::size_t z = 0; z < size_z - 1; z++) { + f({x, y, z}); + } + } + } + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_CARTESIAN_TOPOLOGY_BASE_H diff --git a/Isosurfacing_3/include/CGAL/Cell_type.h b/Isosurfacing_3/include/CGAL/Cell_type.h new file mode 100644 index 00000000000..b897569a493 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Cell_type.h @@ -0,0 +1,20 @@ +#ifndef CGAL_DOMAIN_CELL_TYPE +#define CGAL_DOMAIN_CELL_TYPE + +#include + +namespace CGAL { +namespace Isosurfacing { + +typedef std::size_t Cell_type; + +static constexpr Cell_type ANY_CELL = std::numeric_limits::max(); + +static constexpr Cell_type POLYHERDAL_CELL = ((Cell_type)1) << 0; +static constexpr Cell_type TETRAHEDRAL_CELL = ((Cell_type)1) << 1; +static constexpr Cell_type CUBICAL_CELL = ((Cell_type)1) << 2; + +} +} + +#endif // CGAL_DOMAIN_CELL_TYPE \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Dual_contouring_3.h b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h new file mode 100644 index 00000000000..ee3ce016508 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h @@ -0,0 +1,44 @@ +#ifndef CGAL_DUAL_CONTOURING_3_H +#define CGAL_DUAL_CONTOURING_3_H + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template > +void make_quad_mesh_using_dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value, + PointRange& points, PolygonRange& polygons, + const Positioning& positioning = Positioning()) { + + // static_assert(Domain_::CELL_TYPE & ANY_CELL); + + internal::Dual_contouring_position_functor pos_func(domain, iso_value, positioning); + domain.iterate_cells(pos_func, Concurrency_tag()); + + internal::Dual_contouring_quads_functor quad_func(domain, iso_value); + domain.iterate_edges(quad_func, Concurrency_tag()); + + // write points and quads in ranges + points.resize(pos_func.points_counter); + for (const auto& vtop : pos_func.map_voxel_to_point) { + points[pos_func.map_voxel_to_point_id[vtop.first]] = vtop.second; + } + + polygons.reserve(quad_func.quads.size()); + for (const auto& q : quad_func.quads) { + std::vector vertex_ids; + for (const auto& v_id : q.second) { + vertex_ids.push_back(pos_func.map_voxel_to_point_id[v_id]); + } + polygons.push_back(vertex_ids); + } +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_DUAL_CONTOURING_3_H \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Implicit_domain.h b/Isosurfacing_3/include/CGAL/Implicit_domain.h new file mode 100644 index 00000000000..2e9124e721a --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Implicit_domain.h @@ -0,0 +1,132 @@ +#ifndef CGAL_IMPLICIT_DOMAIN_H +#define CGAL_IMPLICIT_DOMAIN_H + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Implicit_domain : public Cartesian_topology_base { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Grid_spacing; + +public: + Implicit_domain(const Function& func, const Bbox_3& domain, const Grid_spacing& spacing) + : func(&func), bbox(domain), spacing(spacing) { + + sizes[0] = domain.x_span() / spacing.x(); + sizes[1] = domain.y_span() / spacing.y(); + sizes[2] = domain.z_span() / spacing.z(); + } + + Point position(const Vertex_handle& v) const { + return Point(v[0] * spacing.x() + bbox.xmin(), v[1] * spacing.y() + bbox.ymin(), + v[2] * spacing.z() + bbox.zmin()); + } + + FT value(const Vertex_handle& v) const { + return func->operator()(position(v)); + } + + template + void iterate_vertices(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_vertices_base(f, tag, sizes[0], sizes[1], sizes[2]); + } + + template + void iterate_edges(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_edges_base(f, tag, sizes[0], sizes[1], sizes[2]); + } + + template + void iterate_cells(Functor& f, Sequential_tag tag = Sequential_tag()) const { + iterate_cells_base(f, tag, sizes[0], sizes[1], sizes[2]); + } + +#ifdef CGAL_LINKED_WITH_TBB + template + void iterate_vertices(Functor& f, Parallel_tag) const { + const std::size_t size_x = sizes[0]; + const std::size_t size_y = sizes[1]; + const std::size_t size_z = sizes[2]; + + auto iterator = [=, &f](const tbb::blocked_range& r) { + for (std::size_t x = r.begin(); x != r.end(); x++) { + for (std::size_t y = 0; y < size_y; y++) { + for (std::size_t z = 0; z < size_z; z++) { + f({x, y, z}); + } + } + } + }; + + tbb::parallel_for(tbb::blocked_range(0, size_x), iterator); + } + + template + void iterate_edges(Functor& f, Parallel_tag) const { + const std::size_t size_x = sizes[0]; + const std::size_t size_y = sizes[1]; + const std::size_t size_z = sizes[2]; + + auto iterator = [=, &f](const tbb::blocked_range& r) { + for (std::size_t x = r.begin(); x != r.end(); x++) { + for (std::size_t y = 0; y < size_y - 1; y++) { + for (std::size_t z = 0; z < size_z - 1; z++) { + f({x, y, z, 0}); + f({x, y, z, 1}); + f({x, y, z, 2}); + } + } + } + }; + + tbb::parallel_for(tbb::blocked_range(0, size_x - 1), iterator); + } + + template + void iterate_cells(Functor& f, Parallel_tag) const { + const std::size_t size_x = sizes[0]; + const std::size_t size_y = sizes[1]; + const std::size_t size_z = sizes[2]; + + auto iterator = [=, &f](const tbb::blocked_range& r) { + for (std::size_t x = r.begin(); x != r.end(); x++) { + for (std::size_t y = 0; y < size_y - 1; y++) { + for (std::size_t z = 0; z < size_z - 1; z++) { + f({x, y, z}); + } + } + } + }; + + tbb::parallel_for(tbb::blocked_range(0, size_x - 1), iterator); + } +#endif // CGAL_LINKED_WITH_TBB + +private: + const Function* func; + + Bbox_3 bbox; + Grid_spacing spacing; + + std::array sizes; +}; + + +template +Implicit_domain create_implicit_domain(const Function& func, const Bbox_3& domain, + const typename GeomTraits::Vector_3& spacing) { + return Implicit_domain(func, domain, spacing); +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_IMPLICIT_DOMAIN_H diff --git a/Isosurfacing_3/include/CGAL/Implicit_domain_old.h b/Isosurfacing_3/include/CGAL/Implicit_domain_old.h new file mode 100644 index 00000000000..13e59a0677d --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Implicit_domain_old.h @@ -0,0 +1,64 @@ +#ifndef CGAL_IMPLICIT_DOMAIN_OLD_H +#define CGAL_IMPLICIT_DOMAIN_OLD_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Implicit_domain_old { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point_3; + typedef typename Geom_traits::Vector_3 Vector_3; + +public: + Implicit_domain_old(const Function& func, const CGAL::Bbox_3& domain, const Vector_3& resolution) + : func(func), bbox(domain), resolution(resolution) { + + sizes[0] = domain.x_span() / resolution.x(); + sizes[1] = domain.y_span() / resolution.y(); + sizes[2] = domain.z_span() / resolution.z(); + } + + std::size_t size_x() const { + return sizes[0]; + } + std::size_t size_y() const { + return sizes[1]; + } + std::size_t size_z() const { + return sizes[2]; + } + + Point_3 position(const std::size_t x, const std::size_t y, const std::size_t z) const { + return Point_3(x * resolution.x() + bbox.xmin(), y * resolution.y() + bbox.ymin(), + z * resolution.z() + bbox.zmin()); + } + + FT value(const std::size_t x, const std::size_t y, const std::size_t z) const { + return func(position(x, y, z)); + } + +private: + Function func; + + CGAL::Bbox_3 bbox; + Vector_3 resolution; + + std::array sizes; +}; + + +template +Implicit_domain_old create_implicit_domain_old(const Function& func, const CGAL::Bbox_3& domain, + const typename GeomTraits::Vector_3& resolution) { + return Implicit_domain_old(func, domain, resolution); +} + +} // namespace Isosurfacing +} // end namespace CGAL + +#endif // CGAL_IMPLICIT_DOMAIN_OLD_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Dual_contouring_internal.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Dual_contouring_internal.h new file mode 100644 index 00000000000..678a247def4 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Dual_contouring_internal.h @@ -0,0 +1,347 @@ +#ifndef CGAL_DUAL_CONTOURING_3_INTERNAL_DUAL_CONTOURING_3_H +#define CGAL_DUAL_CONTOURING_3_INTERNAL_DUAL_CONTOURING_3_H + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +namespace Positioning { +template +class QEM_SVD { +public: + /// + /// Compute vertex position for Dual Contouring + /// + /// + /// + /// + /// + /// + /// + /// true, if there is a point in the cell + template + bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh, + typename Domain_::Point& point) const { + typedef typename Domain_::Point Point; + typedef typename Domain_::Geom_traits::Vector_3 Vector; + typedef typename Domain_::FT FT; + + typename Domain_::Cell_vertices vertices = domain.cell_vertices(vh); + + namespace Tables = internal::Cube_table; + + std::array s; + std::transform(vertices.begin(), vertices.end(), s.begin(), [&](const auto& v) { return domain.value(v); }); + + std::array b; + std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; }); + + unsigned int cubeindex = 0; + // set bit if corresponding corner is below iso + for (int i = 0; i < Tables::N_VERTICES; ++i) { + cubeindex |= b[i] << i; + } + + if (cubeindex == 0 || cubeindex == 255) { + return false; + } + + std::array pos; + std::transform(vertices.begin(), vertices.end(), pos.begin(), + [&](const auto& v) { return domain.position(v) - CGAL::ORIGIN; }); + + point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0])); // set point to voxel center + + std::array normals; + std::transform(vertices.begin(), vertices.end(), normals.begin(), + [&](const auto& v) { return domain.gradient(v); }); + + // compute edge intersections + std::vector edge_intersections; + std::vector edge_intersection_normals; + + for (int i = 0; i < Tables::N_EDGES; ++i) { + const auto& v0 = Tables::edge_to_vertex[i][0]; + const auto& v1 = Tables::edge_to_vertex[i][1]; + + if (b[v0] != b[v1]) { // e0 + const FT u = (s[v0] - iso_value) / (s[v0] - s[v1]); + const Point p_lerp = CGAL::ORIGIN + ((1 - u) * pos[v0] + u * pos[v1]); + edge_intersections.push_back(p_lerp); + const Vector n_lerp = (1 - u) * normals[v0] + u * normals[v1]; + edge_intersection_normals.push_back(n_lerp); + } + } + + // MC Polygon Center of Mass + if (false) { + Vector com_vec(0, 0, 0); + + for (int i = 0; i < edge_intersections.size(); ++i) { + com_vec += edge_intersections[i] - CGAL::ORIGIN; + } + + Point p = CGAL::ORIGIN + com_vec / edge_intersections.size(); + point = p; + } + + // SVD QEM + if (true) { + Eigen::Matrix3d A; + A.setZero(); + Eigen::Vector3d b; + b.setZero(); + for (std::size_t i = 0; i < edge_intersections.size(); ++i) { + Eigen::Vector3d n_k = {edge_intersection_normals[i].x(), edge_intersection_normals[i].y(), + edge_intersection_normals[i].z()}; + Eigen::Vector3d p_k = {edge_intersections[i].x(), edge_intersections[i].y(), edge_intersections[i].z()}; + double d_k = n_k.transpose() * p_k; + + Eigen::Matrix3d A_k = n_k * n_k.transpose(); + Eigen::Vector3d b_k = d_k * n_k; + A += A_k; + b += b_k; + } + + Eigen::JacobiSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); + // set threshold as in Peter Lindstrom's paper, "Out-of-Core + // Simplification of Large Polygonal Models" + svd.setThreshold(1e-3); + + // Init x hat + Eigen::Vector3d x_hat; + x_hat << point.x(), point.y(), point.z(); + + // Lindstrom formula for QEM new position for singular matrices + Eigen::VectorXd v_svd = x_hat + svd.solve(b - A * x_hat); + + point = Point(v_svd[0], v_svd[1], v_svd[2]); + } + + // bbox + if (use_bbox) { + CGAL::Bbox_3 bbox = (CGAL::ORIGIN + pos[0]).bbox() + (CGAL::ORIGIN + pos[7]).bbox(); + + FT x = std::min(std::max(point.x(), bbox.xmin()), bbox.xmax()); + FT y = std::min(std::max(point.y(), bbox.ymin()), bbox.ymax()); + FT z = std::min(std::max(point.z(), bbox.zmin()), bbox.zmax()); + point = Point(x, y, z); + } + + return true; + } +}; + +class Voxel_center { +public: + /// + /// Compute vertex position for Dual Contouring + /// + /// + /// + /// + /// + /// + /// + /// true, if there is a point in the cell + template + bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh, + typename Domain_::Point& point) const { + typedef typename Domain_::Point Point; + typedef typename Domain_::Vector_3 Vector; + + namespace Tables = internal::Cube_table; + + std::array s = domain.voxel_values(vh); + + std::array b; + std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; }); + + unsigned int cubeindex = 0; + // set bit if corresponding corner is below iso + for (int i = 0; i < Tables::N_VERTICES; ++i) { + cubeindex |= b[i] << i; + } + + if (cubeindex == 0 || cubeindex == 255) { + return false; + } + + std::array p = domain.voxel_vertex_positions(vh); + std::array pos; + std::transform(p.begin(), p.end(), pos.begin(), [](const auto& e) { return e - CGAL::ORIGIN; }); + + point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0])); // set point to voxel center + + return true; + } +}; + +class MC_polygon_center { +public: + /// + /// Compute vertex position for Dual Contouring + /// + /// + /// + /// + /// + /// + /// + /// true, if there is a point in the cell + template + bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh, + typename Domain_::Point& point) const { + typedef typename Domain_::Point Point; + typedef typename Domain_::Vector_3 Vector; + typedef typename Domain_::FT FT; + + namespace Tables = internal::Cube_table; + + std::array s = domain.voxel_values(vh); + + std::array b; + std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; }); + + unsigned int cubeindex = 0; + // set bit if corresponding corner is below iso + for (int i = 0; i < Tables::N_VERTICES; ++i) { + cubeindex |= b[i] << i; + } + + if (cubeindex == 0 || cubeindex == 255) { + return false; + } + + std::array p = domain.voxel_vertex_positions(vh); + std::array pos; + std::transform(p.begin(), p.end(), pos.begin(), [](const auto& e) { return e - CGAL::ORIGIN; }); + + point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0])); // set point to voxel center + + std::array normals = domain.gradient(vh); + + // compute edge intersections + std::vector edge_intersections; + std::vector edge_intersection_normals; + + for (int i = 0; i < Tables::N_EDGES; ++i) { + const auto& v0 = Tables::edge_to_vertex[i][0]; + const auto& v1 = Tables::edge_to_vertex[i][1]; + + if (b[v0] != b[v1]) { // e0 + const FT u = (s[v0] - iso_value) / (s[v0] - s[v1]); + const Point p_lerp = CGAL::ORIGIN + ((1 - u) * pos[v0] + u * pos[v1]); + edge_intersections.push_back(p_lerp); + const Vector n_lerp = (1 - u) * normals[v0] + u * normals[v1]; + edge_intersection_normals.push_back(n_lerp); + } + } + + // MC Polygon Center of Mass + Vector com_vec(0, 0, 0); + + for (int i = 0; i < edge_intersections.size(); ++i) { + com_vec += edge_intersections[i] - CGAL::ORIGIN; + } + + point = CGAL::ORIGIN + com_vec / edge_intersections.size(); + + return true; + } +}; +} // namespace Positioning + +template +class Dual_contouring_position_functor { +private: + typedef Domain_ Domain; + typedef Positioning_ Positioning; + + typedef typename Domain::FT FT; + typedef typename Domain::Point Point; + typedef typename Domain::Cell_handle Cell_handle; + +public: + Dual_contouring_position_functor(const Domain& domain, FT iso_value, const Positioning& positioning) + : domain(domain), iso_value(iso_value), positioning(positioning), points_counter(0) {} + + void operator()(const Cell_handle& v) { + // compute dc-vertices + Point p; + if (positioning.position(domain, iso_value, v, p)) { + + std::lock_guard lock(mutex); + map_voxel_to_point[v] = p; + map_voxel_to_point_id[v] = points_counter++; + } + } + + // private: + const Domain& domain; + FT iso_value; + const Positioning& positioning; + + std::map map_voxel_to_point_id; + std::map map_voxel_to_point; + std::size_t points_counter; + + std::mutex mutex; +}; + +template +class Dual_contouring_quads_functor { +private: + typedef Domain_ Domain; + + typedef typename Domain_::FT FT; + typedef typename Domain_::Edge_handle Edge_handle; + typedef typename Domain_::Cell_handle Cell_handle; + +public: + Dual_contouring_quads_functor(const Domain& domain, FT iso_value) : domain(domain), iso_value(iso_value) {} + + void operator()(const Edge_handle& e) { + // save all quads + const auto& vertices = domain.edge_vertices(e); + const FT s0 = domain.value(vertices[0]); + const FT s1 = domain.value(vertices[1]); + + if (s0 <= iso_value && s1 > iso_value) { + const auto voxels = domain.cells_incident_to_edge(e); + + std::lock_guard lock(mutex); + quads[e] = {voxels[0], voxels[1], voxels[2], voxels[3]}; + + } else if (s1 <= iso_value && s0 > iso_value) { + const auto voxels = domain.cells_incident_to_edge(e); + + std::lock_guard lock(mutex); + quads[e] = {voxels[0], voxels[3], voxels[2], voxels[1]}; + } + } + + // private: + std::map> quads; + + const Domain& domain; + FT iso_value; + + std::mutex mutex; +}; + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_DUAL_CONTOURING_3_INTERNAL_DUAL_CONTOURING_3_H \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h new file mode 100644 index 00000000000..9a7fa587cb0 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h @@ -0,0 +1,167 @@ +#ifndef CGAL_MARCHING_CUBES_3_INTERNAL_MARCHING_CUBES_3_H +#define CGAL_MARCHING_CUBES_3_INTERNAL_MARCHING_CUBES_3_H + +#include + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +template +Point_3 vertex_interpolation(const Point_3& p0, const Point_3& p1, const FT d0, const FT d1, const FT iso_value) { + + FT mu; + + // don't divide by 0 + if (abs(d1 - d0) < 0.000001) { + mu = 0.5; // if both points have the same value, assume isolevel is in the middle + } else { + mu = (iso_value - d0) / (d1 - d0); + } + + assert(mu >= 0.0 || mu <= 1.0); + + // linear interpolation + return Point_3(p1.x() * mu + p0.x() * (1 - mu), p1.y() * mu + p0.y() * (1 - mu), p1.z() * mu + p0.z() * (1 - mu)); +} + +template +std::size_t get_cell_corners(const Domain_& domain, const typename Domain_::Cell_handle& cell, + const typename Domain_::FT iso_value, Corners_& corners, Values_& values) { + + typedef typename Domain_::Vertex_handle Vertex_handle; + + // collect function values and build index + std::size_t v_id = 0; + std::bitset index = 0; // TODO: get size from domain + for (const Vertex_handle& v : domain.cell_vertices(cell)) { + // collect scalar values and computex index + corners[v_id] = domain.position(v); + values[v_id] = domain.value(v); + + if (values[v_id] >= iso_value) { + index.set(v_id); + } + // next cell vertex + v_id++; + } + + return static_cast(index.to_ullong()); +} + +template +void mc_construct_vertices(const CellEdges& cell_edges, const FT iso_value, const std::size_t i_case, + const Corners_& corners, const Values_& values, Vertices_& vertices) { + + // compute for this case the vertices + std::size_t flag = 1; + std::size_t e_id = 0; + + for (const auto& edge : cell_edges) { + (void)edge; // TODO + + if (flag & Cube_table::intersected_edges[i_case]) { + + // generate vertex here, do not care at this point if vertex already exist + + const int v0 = Cube_table::edge_to_vertex[e_id][0]; + const int v1 = Cube_table::edge_to_vertex[e_id][1]; + + vertices[e_id] = vertex_interpolation(corners[v0], corners[v1], values[v0], values[v1], iso_value); + } + flag <<= 1; + e_id++; + } +} + +template +void mc_construct_triangles(const int i_case, const Vertices_& vertices, PointRange& points, PolygonRange& polygons) { + // construct triangles + for (int t = 0; t < 16; t += 3) { + + const int t_index = i_case * 16 + t; + // if (e_tris_list[t_index] == 0x7f) + if (Cube_table::triangle_cases[t_index] == -1) break; + + const int eg0 = Cube_table::triangle_cases[t_index + 0]; // TODO: move more of this stuff into the table + const int eg1 = Cube_table::triangle_cases[t_index + 1]; + const int eg2 = Cube_table::triangle_cases[t_index + 2]; + + const std::size_t p0_idx = points.size(); // TODO: not allowed + + points.push_back(vertices[eg0]); + points.push_back(vertices[eg1]); + points.push_back(vertices[eg2]); + + // insert new triangle in list + polygons.push_back({}); + auto& triangle = polygons.back(); + + triangle.push_back(p0_idx + 2); + triangle.push_back(p0_idx + 1); + triangle.push_back(p0_idx + 0); + } +} + +template +class Marching_cubes_functor { +private: + typedef Domain_ Domain; + typedef PointRange Point_range; + typedef PolygonRange Polygon_range; + + typedef typename Domain::FT FT; + typedef typename Domain::Point Point; + typedef typename Domain::Vertex_handle Vertex_handle; + typedef typename Domain::Edge_handle Edge_handle; + typedef typename Domain::Cell_handle Cell_handle; + +public: + Marching_cubes_functor(const Domain& domain, const FT iso_value, Point_range& points, Polygon_range& polygons) + : domain(domain), iso_value(iso_value), points(points), polygons(polygons) {} + + + void operator()(const Cell_handle& cell) { + + assert(domain.cell_vertices(cell).size() == 8); + assert(domain.cell_edges(cell).size() == 12); + + FT values[8]; + Point corners[8]; + const int i_case = get_cell_corners(domain, cell, iso_value, corners, values); + + const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1 + if (i_case == 0 || i_case == all_bits_set) { + return; + } + + std::array vertices; + mc_construct_vertices(domain.cell_edges(cell), iso_value, i_case, corners, values, vertices); + + std::lock_guard lock(mutex); + mc_construct_triangles(i_case, vertices, points, polygons); + } + +private: + const Domain& domain; + FT iso_value; + + Point_range& points; + Polygon_range& polygons; + + // compute a unique global index for vertices + // use as key the unique edge number + std::map vertex_map; + + std::mutex mutex; +}; + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_MARCHING_CUBES_3_INTERNAL_MARCHING_CUBES_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables.h new file mode 100644 index 00000000000..9125ef48a8b --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables.h @@ -0,0 +1,640 @@ +#ifndef CGAL_MARCHING_CUBES_3_INTERNAL_TABLES_H +#define CGAL_MARCHING_CUBES_3_INTERNAL_TABLES_H + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +namespace Cube_table { +/* + * Naming convention from "A parallel dual marching cubes approach to quad only surface reconstruction - Grosso & + * Zint" + * + * ^ y + * | + * v2------e2------v3 + * /| /| + * e11| e10| + * / e3 / e1 + * v6------e6------v7 | + * | | | | + * | v0------e0--|---v1 --> x + * e7 / e5 / + * | e8 | e9 + * |/ |/ + * v4------e4------v5 + * / + * < z + */ + +constexpr int N_VERTICES = 8; +constexpr int N_EDGES = 12; + +// This table iterates around an edge of a voxel in positive direction, starting from the given voxel (0,0,0). The +// iteration is described in coordinates relative to the given voxel. The last number is the local edge index. +constexpr int edge_to_voxel_neighbor[N_EDGES][4][4] = { + {{0, 0, 0, 0}, {0, -1, 0, 2}, {0, -1, -1, 6}, {0, 0, -1, 4}}, // e0 + {{0, 0, 0, 1}, {1, 0, 0, 3}, {1, 0, -1, 7}, {0, 0, -1, 5}}, // e1 + {{0, 0, 0, 2}, {0, 0, -1, 6}, {0, 1, -1, 4}, {0, 1, 0, 0}}, // e2 + {{0, 0, 0, 3}, {0, 0, -1, 7}, {-1, 0, -1, 5}, {-1, 0, 0, 1}}, // e3 + {{0, 0, 0, 4}, {0, 0, 1, 0}, {0, -1, 1, 2}, {0, -1, 0, 6}}, // e4 + {{0, 0, 0, 5}, {0, 0, 1, 1}, {1, 0, 1, 3}, {1, 0, 0, 7}}, // e5 + {{0, 0, 0, 6}, {0, 1, 0, 4}, {0, 1, 1, 0}, {0, 0, 1, 2}}, // e6 + {{0, 0, 0, 7}, {-1, 0, 0, 5}, {-1, 0, 1, 1}, {0, 0, 1, 3}}, // e7 + {{0, 0, 0, 8}, {-1, 0, 0, 9}, {-1, -1, 0, 10}, {0, -1, 0, 11}}, // e8 + {{0, 0, 0, 9}, {0, -1, 0, 10}, {1, -1, 0, 11}, {1, 0, 0, 8}}, // e9 + {{0, 0, 0, 10}, {1, 0, 0, 11}, {1, 1, 0, 8}, {0, 1, 0, 9}}, // e10 + {{0, 0, 0, 11}, {0, 1, 0, 8}, {-1, 1, 0, 9}, {-1, 0, 0, 10}} // e11 +}; + +/* The global edge index consists of the lexicographical index of the v0 vertex of a voxel, and an index that + * represents the axis. This table maps from the axis index to the local edge index: 0 = x-axis --> 0 1 = y-axis --> + * 3 2 = z-axis --> 8 + */ +constexpr int edge_store_index[3] = {0, 3, 8}; + +// The local vertex indices of an edge. The indices are sorted by axis direction. +constexpr int edge_to_vertex[N_EDGES][2] = { + {0, 1}, // e0 + {1, 3}, // e1 + {2, 3}, // e2 + {0, 2}, // e3 + {4, 5}, // e4 + {5, 7}, // e5 + {6, 7}, // e6 + {4, 6}, // e7 + {0, 4}, // e8 + {1, 5}, // e9 + {3, 7}, // e10 + {2, 6} // e11 +}; + +// The local vertex coordinates within a voxel. +constexpr int local_vertex_position[N_VERTICES][3] = { + {0, 0, 0}, // v0 + {1, 0, 0}, // v1 + {0, 1, 0}, // v2 + {1, 1, 0}, // v3 + {0, 0, 1}, // v4 + {1, 0, 1}, // v5 + {0, 1, 1}, // v6 + {1, 1, 1} // v7 +}; + +// edges are uniquely characterized by the two end vertices, which have a unique vertex id +// the end vertices of the edge are computed in the cell by giving the indices (i,j,k). +// These indices are obtained from the cell index by adding 0 or 1 to i, j or k respectively +// Example: edge 0: (i,j,k) - (i+1,j,k) +// edge 1: (i+1,j,k) - (i+1,j+1,k) +// The first 3 indices are for the first vertex and the second 3 for the second vertex. +// there are 12 edges, assign to each vertex three edges, the global edge numbering +// consist of 3*global_vertex_id + edge_offset. +constexpr int global_edge_id[][4] = {{0, 0, 0, 0}, {1, 0, 0, 1}, {0, 1, 0, 0}, {0, 0, 0, 1}, + {0, 0, 1, 0}, {1, 0, 1, 1}, + {0, 1, 1, 0}, {0, 0, 1, 1}, {0, 0, 0, 2}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 2}}; + +// probably a list without errors +// indicates which edges has to be intersected +constexpr int intersected_edges[256] = { + 0, 265, 515, 778, 2060, 2309, 2575, 2822, 1030, 1295, 1541, 1804, 3082, 3331, 3593, 3840, 400, 153, 915, + 666, 2460, 2197, 2975, 2710, 1430, 1183, 1941, 1692, 3482, 3219, 3993, 3728, 560, 825, 51, 314, 2620, 2869, + 2111, 2358, 1590, 1855, 1077, 1340, 3642, 3891, 3129, 3376, 928, 681, 419, 170, 2988, 2725, 2479, 2214, 1958, + 1711, 1445, 1196, 4010, 3747, 3497, 3232, 2240, 2505, 2755, 3018, 204, 453, 719, 966, 3270, 3535, 3781, 4044, + 1226, 1475, 1737, 1984, 2384, 2137, 2899, 2650, 348, 85, 863, 598, 3414, 3167, 3925, 3676, 1370, 1107, 1881, + 1616, 2800, 3065, 2291, 2554, 764, 1013, 255, 502, 3830, 4095, 3317, 3580, 1786, 2035, 1273, 1520, 2912, 2665, + 2403, 2154, 876, 613, 367, 102, 3942, 3695, 3429, 3180, 1898, 1635, 1385, 1120, 1120, 1385, 1635, 1898, 3180, + 3429, 3695, 3942, 102, 367, 613, 876, 2154, 2403, 2665, 2912, 1520, 1273, 2035, 1786, 3580, 3317, 4095, 3830, + 502, 255, 1013, 764, 2554, 2291, 3065, 2800, 1616, 1881, 1107, 1370, 3676, 3925, 3167, 3414, 598, 863, 85, + 348, 2650, 2899, 2137, 2384, 1984, 1737, 1475, 1226, 4044, 3781, 3535, 3270, 966, 719, 453, 204, 3018, 2755, + 2505, 2240, 3232, 3497, 3747, 4010, 1196, 1445, 1711, 1958, 2214, 2479, 2725, 2988, 170, 419, 681, 928, 3376, + 3129, 3891, 3642, 1340, 1077, 1855, 1590, 2358, 2111, 2869, 2620, 314, 51, 825, 560, 3728, 3993, 3219, 3482, + 1692, 1941, 1183, 1430, 2710, 2975, 2197, 2460, 666, 915, 153, 400, 3840, 3593, 3331, 3082, 1804, 1541, 1295, + 1030, 2822, 2575, 2309, 2060, 778, 515, 265, 0}; + +// list of triangles for Marching Cubes case t, position at t*16 + tri +constexpr int triangle_cases[4096] = { + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 0 <-> mc: 0, class rep: 0 + 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 1 <-> mc: 1, class rep: 1 + 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 2 <-> mc: 2, class rep: 1 + 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 3 <-> mc: 3, class rep: 3 + 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 4 <-> mc: 8, class rep: 1 + 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 5 <-> mc: 9, class rep: 3 + 1, 0, 9, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 6 <-> mc: 10, class rep: 6 + 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 7 <-> mc: 11, class rep: 7 + 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 8 <-> mc: 4, class rep: 1 + 0, 3, 8, 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 9 <-> mc: 5, class rep: 6 + 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 10 <-> mc: 6, class rep: 3 + 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 11 <-> mc: 7, class rep: 7 + 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 12 <-> mc: 12, class rep: 3 + 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 13 <-> mc: 13, class rep: 7 + 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 14 <-> mc: 14, class rep: 7 + 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 15 <-> mc: 15, class rep: 15 + 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 16 <-> mc: 16, class rep: 1 + 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 17 <-> mc: 17, class rep: 3 + 0, 9, 1, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 18 <-> mc: 18, class rep: 6 + 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 19 <-> mc: 19, class rep: 7 + 8, 7, 4, 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 20 <-> mc: 24, class rep: 6 + 11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 21 <-> mc: 25, class rep: 7 + 9, 1, 0, 8, 7, 4, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 22 <-> mc: 26, class rep: 22 + 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1, // quitte: 23 <-> mc: 27, class rep: 23 + 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 24 <-> mc: 20, class rep: 24 + 3, 7, 4, 3, 4, 0, 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 25 <-> mc: 21, class rep: 25 + 9, 10, 2, 9, 2, 0, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 26 <-> mc: 22, class rep: 25 + 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1, // quitte: 27 <-> mc: 23, class rep: 27 + 3, 1, 10, 3, 10, 11, 7, 4, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 28 <-> mc: 28, class rep: 25 + 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1, // quitte: 29 <-> mc: 29, class rep: 29 + 4, 8, 7, 9, 11, 0, 9, 10, 11, 11, 3, 0, -1, -1, -1, -1, // quitte: 30 <-> mc: 30, class rep: 30 + 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 31 <-> mc: 31, class rep: 7 + 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 32 <-> mc: 32, class rep: 1 + 9, 4, 5, 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 33 <-> mc: 33, class rep: 6 + 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 34 <-> mc: 34, class rep: 3 + 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 35 <-> mc: 35, class rep: 7 + 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 36 <-> mc: 40, class rep: 24 + 0, 2, 11, 0, 11, 8, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 37 <-> mc: 41, class rep: 25 + 0, 4, 5, 0, 5, 1, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 38 <-> mc: 42, class rep: 25 + 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1, // quitte: 39 <-> mc: 43, class rep: 29 + 1, 10, 2, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 40 <-> mc: 36, class rep: 6 + 3, 8, 0, 1, 10, 2, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 41 <-> mc: 37, class rep: 22 + 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 42 <-> mc: 38, class rep: 7 + 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1, // quitte: 43 <-> mc: 39, class rep: 23 + 10, 11, 3, 10, 3, 1, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 44 <-> mc: 44, class rep: 25 + 4, 5, 9, 0, 1, 8, 8, 1, 10, 8, 10, 11, -1, -1, -1, -1, // quitte: 45 <-> mc: 45, class rep: 30 + 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1, // quitte: 46 <-> mc: 46, class rep: 27 + 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 47 <-> mc: 47, class rep: 7 + 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 48 <-> mc: 48, class rep: 3 + 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 49 <-> mc: 49, class rep: 7 + 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 50 <-> mc: 50, class rep: 7 + 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 51 <-> mc: 51, class rep: 15 + 7, 5, 9, 7, 9, 8, 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 52 <-> mc: 56, class rep: 25 + 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1, // quitte: 53 <-> mc: 57, class rep: 27 + 2, 11, 3, 0, 8, 1, 1, 8, 7, 1, 7, 5, -1, -1, -1, -1, // quitte: 54 <-> mc: 58, class rep: 30 + 11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 55 <-> mc: 59, class rep: 7 + 9, 8, 7, 9, 7, 5, 10, 2, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 56 <-> mc: 52, class rep: 25 + 10, 2, 1, 9, 0, 5, 5, 0, 3, 5, 3, 7, -1, -1, -1, -1, // quitte: 57 <-> mc: 53, class rep: 30 + 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1, // quitte: 58 <-> mc: 54, class rep: 29 + 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 59 <-> mc: 55, class rep: 7 + 9, 8, 5, 8, 7, 5, 10, 3, 1, 10, 11, 3, -1, -1, -1, -1, // quitte: 60 <-> mc: 60, class rep: 60 + 5, 0, 7, 5, 9, 0, 7, 0, 11, 1, 10, 0, 11, 0, 10, -1, // quitte: 61 <-> mc: 61, class rep: 25 + 11, 0, 10, 11, 3, 0, 10, 0, 5, 8, 7, 0, 5, 0, 7, -1, // quitte: 62 <-> mc: 62, class rep: 25 + 11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 63 <-> mc: 63, class rep: 3 + 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 64 <-> mc: 128, class rep: 1 + 3, 8, 0, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 65 <-> mc: 129, class rep: 6 + 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 66 <-> mc: 130, class rep: 24 + 8, 9, 1, 8, 1, 3, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 67 <-> mc: 131, class rep: 25 + 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 68 <-> mc: 136, class rep: 3 + 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 69 <-> mc: 137, class rep: 7 + 2, 6, 7, 2, 7, 3, 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 70 <-> mc: 138, class rep: 25 + 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1, // quitte: 71 <-> mc: 139, class rep: 27 + 10, 2, 1, 6, 7, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 72 <-> mc: 132, class rep: 6 + 1, 10, 2, 3, 8, 0, 6, 7, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 73 <-> mc: 133, class rep: 22 + 2, 0, 9, 2, 9, 10, 6, 7, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 74 <-> mc: 134, class rep: 25 + 6, 7, 11, 2, 3, 10, 10, 3, 8, 10, 8, 9, -1, -1, -1, -1, // quitte: 75 <-> mc: 135, class rep: 30 + 10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 76 <-> mc: 140, class rep: 7 + 10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1, // quitte: 77 <-> mc: 141, class rep: 23 + 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1, // quitte: 78 <-> mc: 142, class rep: 29 + 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 79 <-> mc: 143, class rep: 7 + 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 80 <-> mc: 144, class rep: 3 + 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 81 <-> mc: 145, class rep: 7 + 8, 11, 6, 8, 6, 4, 9, 1, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 82 <-> mc: 146, class rep: 25 + 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1, // quitte: 83 <-> mc: 147, class rep: 29 + 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 84 <-> mc: 152, class rep: 7 + 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 85 <-> mc: 153, class rep: 15 + 1, 0, 9, 2, 4, 3, 2, 6, 4, 4, 8, 3, -1, -1, -1, -1, // quitte: 86 <-> mc: 154, class rep: 30 + 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 87 <-> mc: 155, class rep: 7 + 6, 4, 8, 6, 8, 11, 2, 1, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 88 <-> mc: 148, class rep: 25 + 1, 10, 2, 3, 11, 0, 0, 11, 6, 0, 6, 4, -1, -1, -1, -1, // quitte: 89 <-> mc: 149, class rep: 30 + 4, 8, 11, 4, 11, 6, 0, 9, 2, 2, 9, 10, -1, -1, -1, -1, // quitte: 90 <-> mc: 150, class rep: 60 + 10, 3, 9, 10, 2, 3, 9, 3, 4, 11, 6, 3, 4, 3, 6, -1, // quitte: 91 <-> mc: 151, class rep: 25 + 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1, // quitte: 92 <-> mc: 156, class rep: 27 + 10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 93 <-> mc: 157, class rep: 7 + 4, 3, 6, 4, 8, 3, 6, 3, 10, 0, 9, 3, 10, 3, 9, -1, // quitte: 94 <-> mc: 158, class rep: 25 + 10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 95 <-> mc: 159, class rep: 3 + 4, 5, 9, 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 96 <-> mc: 160, class rep: 6 + 0, 3, 8, 4, 5, 9, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 97 <-> mc: 161, class rep: 22 + 5, 1, 0, 5, 0, 4, 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 98 <-> mc: 162, class rep: 25 + 11, 6, 7, 8, 4, 3, 3, 4, 5, 3, 5, 1, -1, -1, -1, -1, // quitte: 99 <-> mc: 163, class rep: 30 + 7, 3, 2, 7, 2, 6, 5, 9, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 100 <-> mc: 168, class rep: 25 + 9, 4, 5, 0, 6, 8, 0, 2, 6, 6, 7, 8, -1, -1, -1, -1, // quitte: 101 <-> mc: 169, class rep: 30 + 3, 2, 6, 3, 6, 7, 1, 0, 5, 5, 0, 4, -1, -1, -1, -1, // quitte: 102 <-> mc: 170, class rep: 60 + 6, 8, 2, 6, 7, 8, 2, 8, 1, 4, 5, 8, 1, 8, 5, -1, // quitte: 103 <-> mc: 171, class rep: 25 + 9, 4, 5, 10, 2, 1, 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 104 <-> mc: 164, class rep: 22 + 6, 7, 11, 1, 10, 2, 0, 3, 8, 4, 5, 9, -1, -1, -1, -1, // quitte: 105 <-> mc: 165, class rep: 105 + 7, 11, 6, 5, 10, 4, 4, 10, 2, 4, 2, 0, -1, -1, -1, -1, // quitte: 106 <-> mc: 166, class rep: 30 + 3, 8, 4, 3, 4, 5, 3, 5, 2, 10, 2, 5, 11, 6, 7, -1, // quitte: 107 <-> mc: 167, class rep: 22 + 9, 4, 5, 10, 6, 1, 1, 6, 7, 1, 7, 3, -1, -1, -1, -1, // quitte: 108 <-> mc: 172, class rep: 30 + 1, 10, 6, 1, 6, 7, 1, 7, 0, 8, 0, 7, 9, 4, 5, -1, // quitte: 109 <-> mc: 173, class rep: 22 + 4, 10, 0, 4, 5, 10, 0, 10, 3, 6, 7, 10, 3, 10, 7, -1, // quitte: 110 <-> mc: 174, class rep: 25 + 7, 10, 6, 7, 8, 10, 5, 10, 4, 4, 10, 8, -1, -1, -1, -1, // quitte: 111 <-> mc: 175, class rep: 6 + 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 112 <-> mc: 176, class rep: 7 + 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1, // quitte: 113 <-> mc: 177, class rep: 23 + 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1, // quitte: 114 <-> mc: 178, class rep: 27 + 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 115 <-> mc: 179, class rep: 7 + 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1, // quitte: 116 <-> mc: 184, class rep: 29 + 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 117 <-> mc: 185, class rep: 7 + 1, 8, 5, 1, 0, 8, 5, 8, 6, 3, 2, 8, 6, 8, 2, -1, // quitte: 118 <-> mc: 186, class rep: 25 + 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 119 <-> mc: 187, class rep: 3 + 1, 10, 2, 9, 11, 5, 9, 8, 11, 11, 6, 5, -1, -1, -1, -1, // quitte: 120 <-> mc: 180, class rep: 30 + 0, 3, 11, 0, 11, 6, 0, 6, 9, 5, 9, 6, 1, 10, 2, -1, // quitte: 121 <-> mc: 181, class rep: 22 + 11, 5, 8, 11, 6, 5, 8, 5, 0, 10, 2, 5, 0, 5, 2, -1, // quitte: 122 <-> mc: 182, class rep: 25 + 6, 3, 11, 6, 5, 3, 2, 3, 10, 10, 3, 5, -1, -1, -1, -1, // quitte: 123 <-> mc: 183, class rep: 6 + 1, 6, 3, 1, 10, 6, 3, 6, 8, 5, 9, 6, 8, 6, 9, -1, // quitte: 124 <-> mc: 188, class rep: 25 + 10, 0, 1, 10, 6, 0, 9, 0, 5, 5, 0, 6, -1, -1, -1, -1, // quitte: 125 <-> mc: 189, class rep: 6 + 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 126 <-> mc: 190, class rep: 24 + 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 127 <-> mc: 191, class rep: 1 + 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 128 <-> mc: 64, class rep: 1 + 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 129 <-> mc: 65, class rep: 24 + 9, 1, 0, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 130 <-> mc: 66, class rep: 6 + 1, 3, 8, 1, 8, 9, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 131 <-> mc: 67, class rep: 25 + 2, 11, 3, 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 132 <-> mc: 72, class rep: 6 + 11, 8, 0, 11, 0, 2, 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 133 <-> mc: 73, class rep: 25 + 0, 9, 1, 2, 11, 3, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 134 <-> mc: 74, class rep: 22 + 5, 6, 10, 1, 2, 9, 9, 2, 11, 9, 11, 8, -1, -1, -1, -1, // quitte: 135 <-> mc: 75, class rep: 30 + 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 136 <-> mc: 68, class rep: 3 + 1, 5, 6, 1, 6, 2, 3, 8, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 137 <-> mc: 69, class rep: 25 + 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 138 <-> mc: 70, class rep: 7 + 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1, // quitte: 139 <-> mc: 71, class rep: 29 + 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 140 <-> mc: 76, class rep: 7 + 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1, // quitte: 141 <-> mc: 77, class rep: 27 + 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1, // quitte: 142 <-> mc: 78, class rep: 23 + 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 143 <-> mc: 79, class rep: 7 + 5, 6, 10, 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 144 <-> mc: 80, class rep: 6 + 4, 0, 3, 4, 3, 7, 6, 10, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 145 <-> mc: 81, class rep: 25 + 1, 0, 9, 5, 6, 10, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 146 <-> mc: 82, class rep: 22 + 10, 5, 6, 1, 7, 9, 1, 3, 7, 7, 4, 9, -1, -1, -1, -1, // quitte: 147 <-> mc: 83, class rep: 30 + 3, 2, 11, 7, 4, 8, 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 148 <-> mc: 88, class rep: 22 + 5, 6, 10, 4, 2, 7, 4, 0, 2, 2, 11, 7, -1, -1, -1, -1, // quitte: 149 <-> mc: 89, class rep: 30 + 0, 9, 1, 4, 8, 7, 2, 11, 3, 5, 6, 10, -1, -1, -1, -1, // quitte: 150 <-> mc: 90, class rep: 105 + 9, 1, 2, 9, 2, 11, 9, 11, 4, 7, 4, 11, 5, 6, 10, -1, // quitte: 151 <-> mc: 91, class rep: 22 + 6, 2, 1, 6, 1, 5, 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 152 <-> mc: 84, class rep: 25 + 1, 5, 2, 5, 6, 2, 3, 4, 0, 3, 7, 4, -1, -1, -1, -1, // quitte: 153 <-> mc: 85, class rep: 60 + 8, 7, 4, 9, 5, 0, 0, 5, 6, 0, 6, 2, -1, -1, -1, -1, // quitte: 154 <-> mc: 86, class rep: 30 + 7, 9, 3, 7, 4, 9, 3, 9, 2, 5, 6, 9, 2, 9, 6, -1, // quitte: 155 <-> mc: 87, class rep: 25 + 8, 7, 4, 3, 5, 11, 3, 1, 5, 5, 6, 11, -1, -1, -1, -1, // quitte: 156 <-> mc: 92, class rep: 30 + 5, 11, 1, 5, 6, 11, 1, 11, 0, 7, 4, 11, 0, 11, 4, -1, // quitte: 157 <-> mc: 93, class rep: 25 + 0, 9, 5, 0, 5, 6, 0, 6, 3, 11, 3, 6, 8, 7, 4, -1, // quitte: 158 <-> mc: 94, class rep: 22 + 6, 9, 5, 6, 11, 9, 4, 9, 7, 7, 9, 11, -1, -1, -1, -1, // quitte: 159 <-> mc: 95, class rep: 6 + 10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 160 <-> mc: 96, class rep: 3 + 4, 6, 10, 4, 10, 9, 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 161 <-> mc: 97, class rep: 25 + 10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 162 <-> mc: 98, class rep: 7 + 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1, // quitte: 163 <-> mc: 99, class rep: 27 + 10, 9, 4, 10, 4, 6, 11, 3, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 164 <-> mc: 104, class rep: 25 + 0, 2, 8, 2, 11, 8, 4, 10, 9, 4, 6, 10, -1, -1, -1, -1, // quitte: 165 <-> mc: 105, class rep: 60 + 3, 2, 11, 0, 6, 1, 0, 4, 6, 6, 10, 1, -1, -1, -1, -1, // quitte: 166 <-> mc: 106, class rep: 30 + 6, 1, 4, 6, 10, 1, 4, 1, 8, 2, 11, 1, 8, 1, 11, -1, // quitte: 167 <-> mc: 107, class rep: 25 + 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 168 <-> mc: 100, class rep: 7 + 3, 8, 0, 1, 9, 2, 2, 9, 4, 2, 4, 6, -1, -1, -1, -1, // quitte: 169 <-> mc: 101, class rep: 30 + 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 170 <-> mc: 102, class rep: 15 + 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 171 <-> mc: 103, class rep: 7 + 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1, // quitte: 172 <-> mc: 108, class rep: 29 + 8, 1, 11, 8, 0, 1, 11, 1, 6, 9, 4, 1, 6, 1, 4, -1, // quitte: 173 <-> mc: 109, class rep: 25 + 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1, // quitte: 174 <-> mc: 110, class rep: 7 + 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 175 <-> mc: 111, class rep: 3 + 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 176 <-> mc: 112, class rep: 7 + 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1, // quitte: 177 <-> mc: 113, class rep: 29 + 10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1, // quitte: 178 <-> mc: 114, class rep: 23 + 10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 179 <-> mc: 115, class rep: 7 + 2, 11, 3, 10, 8, 6, 10, 9, 8, 8, 7, 6, -1, -1, -1, -1, // quitte: 180 <-> mc: 120, class rep: 30 + 2, 7, 0, 2, 11, 7, 0, 7, 9, 6, 10, 7, 9, 7, 10, -1, // quitte: 181 <-> mc: 121, class rep: 25 + 1, 0, 8, 1, 8, 7, 1, 7, 10, 6, 10, 7, 2, 11, 3, -1, // quitte: 182 <-> mc: 122, class rep: 22 + 11, 1, 2, 11, 7, 1, 10, 1, 6, 6, 1, 7, -1, -1, -1, -1, // quitte: 183 <-> mc: 123, class rep: 6 + 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1, // quitte: 184 <-> mc: 116, class rep: 27 + 2, 9, 6, 2, 1, 9, 6, 9, 7, 0, 3, 9, 7, 9, 3, -1, // quitte: 185 <-> mc: 117, class rep: 25 + 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1, // quitte: 186 <-> mc: 118, class rep: 7 + 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 187 <-> mc: 119, class rep: 3 + 8, 6, 9, 8, 7, 6, 9, 6, 1, 11, 3, 6, 1, 6, 3, -1, // quitte: 188 <-> mc: 124, class rep: 25 + 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 189 <-> mc: 125, class rep: 24 + 7, 0, 8, 7, 6, 0, 3, 0, 11, 11, 0, 6, -1, -1, -1, -1, // quitte: 190 <-> mc: 126, class rep: 6 + 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 191 <-> mc: 127, class rep: 1 + 11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 192 <-> mc: 192, class rep: 3 + 11, 10, 5, 11, 5, 7, 8, 0, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 193 <-> mc: 193, class rep: 25 + 5, 7, 11, 5, 11, 10, 1, 0, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 194 <-> mc: 194, class rep: 25 + 10, 5, 7, 10, 7, 11, 9, 1, 8, 8, 1, 3, -1, -1, -1, -1, // quitte: 195 <-> mc: 195, class rep: 60 + 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 196 <-> mc: 200, class rep: 7 + 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1, // quitte: 197 <-> mc: 201, class rep: 29 + 9, 1, 0, 5, 3, 10, 5, 7, 3, 3, 2, 10, -1, -1, -1, -1, // quitte: 198 <-> mc: 202, class rep: 30 + 9, 2, 8, 9, 1, 2, 8, 2, 7, 10, 5, 2, 7, 2, 5, -1, // quitte: 199 <-> mc: 203, class rep: 25 + 11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 200 <-> mc: 196, class rep: 7 + 0, 3, 8, 1, 7, 2, 1, 5, 7, 7, 11, 2, -1, -1, -1, -1, // quitte: 201 <-> mc: 197, class rep: 30 + 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1, // quitte: 202 <-> mc: 198, class rep: 27 + 7, 2, 5, 7, 11, 2, 5, 2, 9, 3, 8, 2, 9, 2, 8, -1, // quitte: 203 <-> mc: 199, class rep: 25 + 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 204 <-> mc: 204, class rep: 15 + 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1, // quitte: 205 <-> mc: 205, class rep: 7 + 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1, // quitte: 206 <-> mc: 206, class rep: 7 + 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 207 <-> mc: 207, class rep: 3 + 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 208 <-> mc: 208, class rep: 7 + 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1, // quitte: 209 <-> mc: 209, class rep: 27 + 0, 9, 1, 8, 10, 4, 8, 11, 10, 10, 5, 4, -1, -1, -1, -1, // quitte: 210 <-> mc: 210, class rep: 30 + 10, 4, 11, 10, 5, 4, 11, 4, 3, 9, 1, 4, 3, 4, 1, -1, // quitte: 211 <-> mc: 211, class rep: 25 + 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1, // quitte: 212 <-> mc: 216, class rep: 23 + 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1, // quitte: 213 <-> mc: 217, class rep: 7 + 3, 2, 10, 3, 10, 5, 3, 5, 8, 4, 8, 5, 0, 9, 1, -1, // quitte: 214 <-> mc: 218, class rep: 22 + 5, 2, 10, 5, 4, 2, 1, 2, 9, 9, 2, 4, -1, -1, -1, -1, // quitte: 215 <-> mc: 219, class rep: 6 + 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1, // quitte: 216 <-> mc: 212, class rep: 29 + 0, 11, 4, 0, 3, 11, 4, 11, 5, 2, 1, 11, 5, 11, 1, -1, // quitte: 217 <-> mc: 213, class rep: 25 + 0, 5, 2, 0, 9, 5, 2, 5, 11, 4, 8, 5, 11, 5, 8, -1, // quitte: 218 <-> mc: 214, class rep: 25 + 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 219 <-> mc: 215, class rep: 24 + 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1, // quitte: 220 <-> mc: 220, class rep: 7 + 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 221 <-> mc: 221, class rep: 3 + 8, 5, 4, 8, 3, 5, 9, 5, 0, 0, 5, 3, -1, -1, -1, -1, // quitte: 222 <-> mc: 222, class rep: 6 + 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 223 <-> mc: 223, class rep: 1 + 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 224 <-> mc: 224, class rep: 7 + 0, 3, 8, 4, 7, 9, 9, 7, 11, 9, 11, 10, -1, -1, -1, -1, // quitte: 225 <-> mc: 225, class rep: 30 + 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1, // quitte: 226 <-> mc: 226, class rep: 29 + 3, 4, 1, 3, 8, 4, 1, 4, 10, 7, 11, 4, 10, 4, 11, -1, // quitte: 227 <-> mc: 227, class rep: 25 + 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1, // quitte: 228 <-> mc: 232, class rep: 27 + 9, 7, 10, 9, 4, 7, 10, 7, 2, 8, 0, 7, 2, 7, 0, -1, // quitte: 229 <-> mc: 233, class rep: 25 + 3, 10, 7, 3, 2, 10, 7, 10, 4, 1, 0, 10, 4, 10, 0, -1, // quitte: 230 <-> mc: 234, class rep: 25 + 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 231 <-> mc: 235, class rep: 24 + 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1, // quitte: 232 <-> mc: 228, class rep: 23 + 9, 4, 7, 9, 7, 11, 9, 11, 1, 2, 1, 11, 0, 3, 8, -1, // quitte: 233 <-> mc: 229, class rep: 22 + 11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1, // quitte: 234 <-> mc: 230, class rep: 7 + 11, 4, 7, 11, 2, 4, 8, 4, 3, 3, 4, 2, -1, -1, -1, -1, // quitte: 235 <-> mc: 231, class rep: 6 + 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1, // quitte: 236 <-> mc: 236, class rep: 7 + 4, 1, 9, 4, 7, 1, 0, 1, 8, 8, 1, 7, -1, -1, -1, -1, // quitte: 237 <-> mc: 237, class rep: 6 + 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 238 <-> mc: 238, class rep: 3 + 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 239 <-> mc: 239, class rep: 1 + 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 240 <-> mc: 240, class rep: 15 + 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1, // quitte: 241 <-> mc: 241, class rep: 7 + 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1, // quitte: 242 <-> mc: 242, class rep: 7 + 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 243 <-> mc: 243, class rep: 3 + 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1, // quitte: 244 <-> mc: 248, class rep: 7 + 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 245 <-> mc: 249, class rep: 3 + 2, 8, 3, 2, 10, 8, 0, 8, 1, 1, 8, 10, -1, -1, -1, -1, // quitte: 246 <-> mc: 250, class rep: 6 + 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 247 <-> mc: 251, class rep: 1 + 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1, // quitte: 248 <-> mc: 244, class rep: 7 + 3, 9, 0, 3, 11, 9, 1, 9, 2, 2, 9, 11, -1, -1, -1, -1, // quitte: 249 <-> mc: 245, class rep: 6 + 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 250 <-> mc: 246, class rep: 3 + 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 251 <-> mc: 247, class rep: 1 + 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 252 <-> mc: 252, class rep: 3 + 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 253 <-> mc: 253, class rep: 1 + 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, // quitte: 254 <-> mc: 254, class rep: 1 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 // quitte: 255 <-> mc: 255, class rep: 0 +}; + +// list of ambiguous cases +constexpr int t_ambig[256] = { + 0, // quitte: 0 <-> mc: 0, class representative: 0 + 1, // quitte: 1 <-> mc: 1, class representative: 1 + 2, // quitte: 2 <-> mc: 2, class representative: 1 + 3, // quitte: 3 <-> mc: 3, class representative: 3 + 4, // quitte: 4 <-> mc: 8, class representative: 1 + 5, // quitte: 5 <-> mc: 9, class representative: 3 + 105, // quitte: 6 <-> mc: 10, class representative: 6 + 7, // quitte: 7 <-> mc: 11, class representative: 7 + 8, // quitte: 8 <-> mc: 4, class representative: 1 + 105, // quitte: 9 <-> mc: 5, class representative: 6 + 10, // quitte: 10 <-> mc: 6, class representative: 3 + 11, // quitte: 11 <-> mc: 7, class representative: 7 + 12, // quitte: 12 <-> mc: 12, class representative: 3 + 13, // quitte: 13 <-> mc: 13, class representative: 7 + 14, // quitte: 14 <-> mc: 14, class representative: 7 + 15, // quitte: 15 <-> mc: 15, class representative: 15 + 16, // quitte: 16 <-> mc: 16, class representative: 1 + 17, // quitte: 17 <-> mc: 17, class representative: 3 + 105, // quitte: 18 <-> mc: 18, class representative: 6 + 19, // quitte: 19 <-> mc: 19, class representative: 7 + 105, // quitte: 20 <-> mc: 24, class representative: 6 + 21, // quitte: 21 <-> mc: 25, class representative: 7 + 105, // quitte: 22 <-> mc: 26, class representative: 22 + 23, // quitte: 23 <-> mc: 27, class representative: 23 + 105, // quitte: 24 <-> mc: 20, class representative: 24 + 105, // quitte: 25 <-> mc: 21, class representative: 25 + 105, // quitte: 26 <-> mc: 22, class representative: 25 + 27, // quitte: 27 <-> mc: 23, class representative: 27 + 105, // quitte: 28 <-> mc: 28, class representative: 25 + 29, // quitte: 29 <-> mc: 29, class representative: 29 + 105, // quitte: 30 <-> mc: 30, class representative: 30 + 31, // quitte: 31 <-> mc: 31, class representative: 7 + 32, // quitte: 32 <-> mc: 32, class representative: 1 + 105, // quitte: 33 <-> mc: 33, class representative: 6 + 34, // quitte: 34 <-> mc: 34, class representative: 3 + 35, // quitte: 35 <-> mc: 35, class representative: 7 + 105, // quitte: 36 <-> mc: 40, class representative: 24 + 105, // quitte: 37 <-> mc: 41, class representative: 25 + 105, // quitte: 38 <-> mc: 42, class representative: 25 + 39, // quitte: 39 <-> mc: 43, class representative: 29 + 105, // quitte: 40 <-> mc: 36, class representative: 6 + 105, // quitte: 41 <-> mc: 37, class representative: 22 + 42, // quitte: 42 <-> mc: 38, class representative: 7 + 43, // quitte: 43 <-> mc: 39, class representative: 23 + 105, // quitte: 44 <-> mc: 44, class representative: 25 + 105, // quitte: 45 <-> mc: 45, class representative: 30 + 46, // quitte: 46 <-> mc: 46, class representative: 27 + 47, // quitte: 47 <-> mc: 47, class representative: 7 + 48, // quitte: 48 <-> mc: 48, class representative: 3 + 49, // quitte: 49 <-> mc: 49, class representative: 7 + 50, // quitte: 50 <-> mc: 50, class representative: 7 + 51, // quitte: 51 <-> mc: 51, class representative: 15 + 105, // quitte: 52 <-> mc: 56, class representative: 25 + 53, // quitte: 53 <-> mc: 57, class representative: 27 + 105, // quitte: 54 <-> mc: 58, class representative: 30 + 55, // quitte: 55 <-> mc: 59, class representative: 7 + 105, // quitte: 56 <-> mc: 52, class representative: 25 + 105, // quitte: 57 <-> mc: 53, class representative: 30 + 58, // quitte: 58 <-> mc: 54, class representative: 29 + 59, // quitte: 59 <-> mc: 55, class representative: 7 + 105, // quitte: 60 <-> mc: 60, class representative: 60 + 105, // quitte: 61 <-> mc: 61, class representative: 25 + 105, // quitte: 62 <-> mc: 62, class representative: 25 + 63, // quitte: 63 <-> mc: 63, class representative: 3 + 64, // quitte: 64 <-> mc: 128, class representative: 1 + 105, // quitte: 65 <-> mc: 129, class representative: 6 + 105, // quitte: 66 <-> mc: 130, class representative: 24 + 105, // quitte: 67 <-> mc: 131, class representative: 25 + 68, // quitte: 68 <-> mc: 136, class representative: 3 + 69, // quitte: 69 <-> mc: 137, class representative: 7 + 105, // quitte: 70 <-> mc: 138, class representative: 25 + 71, // quitte: 71 <-> mc: 139, class representative: 27 + 105, // quitte: 72 <-> mc: 132, class representative: 6 + 105, // quitte: 73 <-> mc: 133, class representative: 22 + 105, // quitte: 74 <-> mc: 134, class representative: 25 + 105, // quitte: 75 <-> mc: 135, class representative: 30 + 76, // quitte: 76 <-> mc: 140, class representative: 7 + 77, // quitte: 77 <-> mc: 141, class representative: 23 + 78, // quitte: 78 <-> mc: 142, class representative: 29 + 79, // quitte: 79 <-> mc: 143, class representative: 7 + 80, // quitte: 80 <-> mc: 144, class representative: 3 + 81, // quitte: 81 <-> mc: 145, class representative: 7 + 105, // quitte: 82 <-> mc: 146, class representative: 25 + 83, // quitte: 83 <-> mc: 147, class representative: 29 + 84, // quitte: 84 <-> mc: 152, class representative: 7 + 85, // quitte: 85 <-> mc: 153, class representative: 15 + 105, // quitte: 86 <-> mc: 154, class representative: 30 + 87, // quitte: 87 <-> mc: 155, class representative: 7 + 105, // quitte: 88 <-> mc: 148, class representative: 25 + 105, // quitte: 89 <-> mc: 149, class representative: 30 + 105, // quitte: 90 <-> mc: 150, class representative: 60 + 105, // quitte: 91 <-> mc: 151, class representative: 25 + 92, // quitte: 92 <-> mc: 156, class representative: 27 + 93, // quitte: 93 <-> mc: 157, class representative: 7 + 105, // quitte: 94 <-> mc: 158, class representative: 25 + 95, // quitte: 95 <-> mc: 159, class representative: 3 + 105, // quitte: 96 <-> mc: 160, class representative: 6 + 105, // quitte: 97 <-> mc: 161, class representative: 22 + 105, // quitte: 98 <-> mc: 162, class representative: 25 + 105, // quitte: 99 <-> mc: 163, class representative: 30 + 105, // quitte: 100 <-> mc: 168, class representative: 25 + 105, // quitte: 101 <-> mc: 169, class representative: 30 + 105, // quitte: 102 <-> mc: 170, class representative: 60 + 105, // quitte: 103 <-> mc: 171, class representative: 25 + 105, // quitte: 104 <-> mc: 164, class representative: 22 + 105, // quitte: 105 <-> mc: 165, class representative: 105 + 105, // quitte: 106 <-> mc: 166, class representative: 30 + 105, // quitte: 107 <-> mc: 167, class representative: 22 + 105, // quitte: 108 <-> mc: 172, class representative: 30 + 105, // quitte: 109 <-> mc: 173, class representative: 22 + 105, // quitte: 110 <-> mc: 174, class representative: 25 + 105, // quitte: 111 <-> mc: 175, class representative: 6 + 112, // quitte: 112 <-> mc: 176, class representative: 7 + 113, // quitte: 113 <-> mc: 177, class representative: 23 + 114, // quitte: 114 <-> mc: 178, class representative: 27 + 115, // quitte: 115 <-> mc: 179, class representative: 7 + 116, // quitte: 116 <-> mc: 184, class representative: 29 + 117, // quitte: 117 <-> mc: 185, class representative: 7 + 105, // quitte: 118 <-> mc: 186, class representative: 25 + 119, // quitte: 119 <-> mc: 187, class representative: 3 + 105, // quitte: 120 <-> mc: 180, class representative: 30 + 105, // quitte: 121 <-> mc: 181, class representative: 22 + 105, // quitte: 122 <-> mc: 182, class representative: 25 + 105, // quitte: 123 <-> mc: 183, class representative: 6 + 105, // quitte: 124 <-> mc: 188, class representative: 25 + 105, // quitte: 125 <-> mc: 189, class representative: 6 + 105, // quitte: 126 <-> mc: 190, class representative: 24 + 127, // quitte: 127 <-> mc: 191, class representative: 1 + 128, // quitte: 128 <-> mc: 64, class representative: 1 + 105, // quitte: 129 <-> mc: 65, class representative: 24 + 105, // quitte: 130 <-> mc: 66, class representative: 6 + 105, // quitte: 131 <-> mc: 67, class representative: 25 + 105, // quitte: 132 <-> mc: 72, class representative: 6 + 105, // quitte: 133 <-> mc: 73, class representative: 25 + 105, // quitte: 134 <-> mc: 74, class representative: 22 + 105, // quitte: 135 <-> mc: 75, class representative: 30 + 136, // quitte: 136 <-> mc: 68, class representative: 3 + 105, // quitte: 137 <-> mc: 69, class representative: 25 + 138, // quitte: 138 <-> mc: 70, class representative: 7 + 139, // quitte: 139 <-> mc: 71, class representative: 29 + 140, // quitte: 140 <-> mc: 76, class representative: 7 + 141, // quitte: 141 <-> mc: 77, class representative: 27 + 142, // quitte: 142 <-> mc: 78, class representative: 23 + 143, // quitte: 143 <-> mc: 79, class representative: 7 + 105, // quitte: 144 <-> mc: 80, class representative: 6 + 105, // quitte: 145 <-> mc: 81, class representative: 25 + 105, // quitte: 146 <-> mc: 82, class representative: 22 + 105, // quitte: 147 <-> mc: 83, class representative: 30 + 105, // quitte: 148 <-> mc: 88, class representative: 22 + 105, // quitte: 149 <-> mc: 89, class representative: 30 + 105, // quitte: 150 <-> mc: 90, class representative: 105 + 105, // quitte: 151 <-> mc: 91, class representative: 22 + 105, // quitte: 152 <-> mc: 84, class representative: 25 + 105, // quitte: 153 <-> mc: 85, class representative: 60 + 105, // quitte: 154 <-> mc: 86, class representative: 30 + 105, // quitte: 155 <-> mc: 87, class representative: 25 + 105, // quitte: 156 <-> mc: 92, class representative: 30 + 105, // quitte: 157 <-> mc: 93, class representative: 25 + 105, // quitte: 158 <-> mc: 94, class representative: 22 + 105, // quitte: 159 <-> mc: 95, class representative: 6 + 160, // quitte: 160 <-> mc: 96, class representative: 3 + 105, // quitte: 161 <-> mc: 97, class representative: 25 + 162, // quitte: 162 <-> mc: 98, class representative: 7 + 163, // quitte: 163 <-> mc: 99, class representative: 27 + 105, // quitte: 164 <-> mc: 104, class representative: 25 + 105, // quitte: 165 <-> mc: 105, class representative: 60 + 105, // quitte: 166 <-> mc: 106, class representative: 30 + 105, // quitte: 167 <-> mc: 107, class representative: 25 + 168, // quitte: 168 <-> mc: 100, class representative: 7 + 105, // quitte: 169 <-> mc: 101, class representative: 30 + 170, // quitte: 170 <-> mc: 102, class representative: 15 + 171, // quitte: 171 <-> mc: 103, class representative: 7 + 172, // quitte: 172 <-> mc: 108, class representative: 29 + 105, // quitte: 173 <-> mc: 109, class representative: 25 + 174, // quitte: 174 <-> mc: 110, class representative: 7 + 175, // quitte: 175 <-> mc: 111, class representative: 3 + 176, // quitte: 176 <-> mc: 112, class representative: 7 + 177, // quitte: 177 <-> mc: 113, class representative: 29 + 178, // quitte: 178 <-> mc: 114, class representative: 23 + 179, // quitte: 179 <-> mc: 115, class representative: 7 + 105, // quitte: 180 <-> mc: 120, class representative: 30 + 105, // quitte: 181 <-> mc: 121, class representative: 25 + 105, // quitte: 182 <-> mc: 122, class representative: 22 + 105, // quitte: 183 <-> mc: 123, class representative: 6 + 184, // quitte: 184 <-> mc: 116, class representative: 27 + 105, // quitte: 185 <-> mc: 117, class representative: 25 + 186, // quitte: 186 <-> mc: 118, class representative: 7 + 187, // quitte: 187 <-> mc: 119, class representative: 3 + 105, // quitte: 188 <-> mc: 124, class representative: 25 + 105, // quitte: 189 <-> mc: 125, class representative: 24 + 105, // quitte: 190 <-> mc: 126, class representative: 6 + 191, // quitte: 191 <-> mc: 127, class representative: 1 + 192, // quitte: 192 <-> mc: 192, class representative: 3 + 105, // quitte: 193 <-> mc: 193, class representative: 25 + 105, // quitte: 194 <-> mc: 194, class representative: 25 + 105, // quitte: 195 <-> mc: 195, class representative: 60 + 196, // quitte: 196 <-> mc: 200, class representative: 7 + 197, // quitte: 197 <-> mc: 201, class representative: 29 + 105, // quitte: 198 <-> mc: 202, class representative: 30 + 105, // quitte: 199 <-> mc: 203, class representative: 25 + 200, // quitte: 200 <-> mc: 196, class representative: 7 + 105, // quitte: 201 <-> mc: 197, class representative: 30 + 202, // quitte: 202 <-> mc: 198, class representative: 27 + 105, // quitte: 203 <-> mc: 199, class representative: 25 + 204, // quitte: 204 <-> mc: 204, class representative: 15 + 205, // quitte: 205 <-> mc: 205, class representative: 7 + 206, // quitte: 206 <-> mc: 206, class representative: 7 + 207, // quitte: 207 <-> mc: 207, class representative: 3 + 208, // quitte: 208 <-> mc: 208, class representative: 7 + 209, // quitte: 209 <-> mc: 209, class representative: 27 + 105, // quitte: 210 <-> mc: 210, class representative: 30 + 105, // quitte: 211 <-> mc: 211, class representative: 25 + 212, // quitte: 212 <-> mc: 216, class representative: 23 + 213, // quitte: 213 <-> mc: 217, class representative: 7 + 105, // quitte: 214 <-> mc: 218, class representative: 22 + 105, // quitte: 215 <-> mc: 219, class representative: 6 + 216, // quitte: 216 <-> mc: 212, class representative: 29 + 105, // quitte: 217 <-> mc: 213, class representative: 25 + 105, // quitte: 218 <-> mc: 214, class representative: 25 + 105, // quitte: 219 <-> mc: 215, class representative: 24 + 220, // quitte: 220 <-> mc: 220, class representative: 7 + 221, // quitte: 221 <-> mc: 221, class representative: 3 + 105, // quitte: 222 <-> mc: 222, class representative: 6 + 223, // quitte: 223 <-> mc: 223, class representative: 1 + 224, // quitte: 224 <-> mc: 224, class representative: 7 + 105, // quitte: 225 <-> mc: 225, class representative: 30 + 226, // quitte: 226 <-> mc: 226, class representative: 29 + 105, // quitte: 227 <-> mc: 227, class representative: 25 + 228, // quitte: 228 <-> mc: 232, class representative: 27 + 105, // quitte: 229 <-> mc: 233, class representative: 25 + 105, // quitte: 230 <-> mc: 234, class representative: 25 + 105, // quitte: 231 <-> mc: 235, class representative: 24 + 232, // quitte: 232 <-> mc: 228, class representative: 23 + 105, // quitte: 233 <-> mc: 229, class representative: 22 + 234, // quitte: 234 <-> mc: 230, class representative: 7 + 105, // quitte: 235 <-> mc: 231, class representative: 6 + 236, // quitte: 236 <-> mc: 236, class representative: 7 + 105, // quitte: 237 <-> mc: 237, class representative: 6 + 238, // quitte: 238 <-> mc: 238, class representative: 3 + 239, // quitte: 239 <-> mc: 239, class representative: 1 + 240, // quitte: 240 <-> mc: 240, class representative: 15 + 241, // quitte: 241 <-> mc: 241, class representative: 7 + 242, // quitte: 242 <-> mc: 242, class representative: 7 + 243, // quitte: 243 <-> mc: 243, class representative: 3 + 244, // quitte: 244 <-> mc: 248, class representative: 7 + 245, // quitte: 245 <-> mc: 249, class representative: 3 + 105, // quitte: 246 <-> mc: 250, class representative: 6 + 247, // quitte: 247 <-> mc: 251, class representative: 1 + 248, // quitte: 248 <-> mc: 244, class representative: 7 + 105, // quitte: 249 <-> mc: 245, class representative: 6 + 250, // quitte: 250 <-> mc: 246, class representative: 3 + 251, // quitte: 251 <-> mc: 247, class representative: 1 + 252, // quitte: 252 <-> mc: 252, class representative: 3 + 253, // quitte: 253 <-> mc: 253, class representative: 1 + 254, // quitte: 254 <-> mc: 254, class representative: 1 + 255 // quitte: 255 <-> mc: 255, class representative: 0 +}; + +} // namespace Cube_table + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_MARCHING_CUBES_3_INTERNAL_TABLES_H \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables_old.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables_old.h new file mode 100755 index 00000000000..09f2f5cb6f2 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tables_old.h @@ -0,0 +1,306 @@ +#pragma once +///////////////////////////////////////////////////// +// tables +///////////////////////////////////////////////////// + +// TODO: find better version with licence +// compare to mc paper +// maybe take from pierre / roberto + +// Polygonising a scalar field +// Also known as: "3D Contouring", "Marching Cubes", "Surface Reconstruction" +// Written by Paul Bourke +// May 1994 +// http://paulbourke.net/geometry/polygonise/ + +const static int edgeTable[256] = { + 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, + 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, + 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, + 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, + 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c, + 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, + 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac, + 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, + 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c, + 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc, + 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c, + 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc , + 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, + 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, + 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, + 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, + 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, + 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, + 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, + 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, + 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460, + 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, + 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0, + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, + 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230, + 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, + 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, + 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, + 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 }; + +const static int triTable[256][16] = +{ { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1 }, +{ 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1 }, +{ 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1 }, +{ 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1 }, +{ 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1 }, +{ 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1 }, +{ 10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1 }, +{ 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1 }, +{ 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1 }, +{ 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1 }, +{ 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1 }, +{ 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1 }, +{ 11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1 }, +{ 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1 }, +{ 11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1 }, +{ 11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1 }, +{ 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1 }, +{ 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1 }, +{ 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1 }, +{ 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1 }, +{ 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1 }, +{ 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1 }, +{ 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1 }, +{ 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1 }, +{ 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1 }, +{ 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1 }, +{ 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1 }, +{ 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1 }, +{ 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1 }, +{ 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1 }, +{ 10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1 }, +{ 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1 }, +{ 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1 }, +{ 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1 }, +{ 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1 }, +{ 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1 }, +{ 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1 }, +{ 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1 }, +{ 10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1 }, +{ 10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1 }, +{ 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1 }, +{ 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1 }, +{ 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1 }, +{ 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1 }, +{ 11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1 }, +{ 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1 }, +{ 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1 }, +{ 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1 }, +{ 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1 }, +{ 10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1 }, +{ 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1 }, +{ 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1 }, +{ 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1 }, +{ 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1 }, +{ 10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1 }, +{ 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1 }, +{ 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1 }, +{ 10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1 }, +{ 10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1 }, +{ 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1 }, +{ 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1 }, +{ 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1 }, +{ 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1 }, +{ 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1 }, +{ 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1 }, +{ 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1 }, +{ 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1 }, +{ 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1 }, +{ 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1 }, +{ 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1 }, +{ 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1 }, +{ 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1 }, +{ 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1 }, +{ 11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1 }, +{ 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1 }, +{ 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1 }, +{ 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1 }, +{ 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1 }, +{ 10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1 }, +{ 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1 }, +{ 11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1 }, +{ 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1 }, +{ 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1 }, +{ 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1 }, +{ 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1 }, +{ 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1 }, +{ 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1 }, +{ 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1 }, +{ 10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1 }, +{ 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1 }, +{ 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1 }, +{ 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1 }, +{ 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1 }, +{ 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1 }, +{ 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1 }, +{ 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1 }, +{ 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1 }, +{ 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1 }, +{ 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1 }, +{ 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1 }, +{ 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1 }, +{ 11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1 }, +{ 11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1 }, +{ 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1 }, +{ 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1 }, +{ 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1 }, +{ 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1 }, +{ 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1 }, +{ 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1 }, +{ 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1 }, +{ 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, +{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }; diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tmc_internal.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tmc_internal.h new file mode 100644 index 00000000000..35e04969c90 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tmc_internal.h @@ -0,0 +1,901 @@ +#ifndef CGAL_TMC_INTERNAL_TMC_H +#define CGAL_TMC_INTERNAL_TMC_H + +#include +#include + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +template +class TMC_functor { +private: + typedef Domain_ Domain; + typedef PointRange Point_range; + typedef PolygonRange Polygon_range; + + typedef typename Domain::FT FT; + typedef typename Domain::Point_3 Point; + typedef typename Domain::Vector_3 Vector; + typedef typename Domain::Vertex_handle Vertex_handle; + typedef typename Domain::Edge_handle Edge_handle; + typedef typename Domain::Cell_handle Cell_handle; + +public: + TMC_functor(const Domain& domain, const FT iso_value, Point_range& points, Polygon_range& polygons) + : domain(domain), iso_value(iso_value), points(points), polygons(polygons) {} + + void operator()(const Cell_handle& cell) { + + FT values[8]; + Point corners[8]; + const int i_case = get_cell_corners(domain, cell, iso_value, corners, values); + + const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1 + if (Cube_table::intersected_edges[i_case] == 0 || Cube_table::intersected_edges[i_case] == all_bits_set) { + return; + } + + // this is the only difference to mc + int tcm = (int)Cube_table::t_ambig[i_case]; + if (tcm == 105) { + p_slice(cell, iso_value, values, corners, i_case); + return; + } + + std::array vertices; + mc_construct_vertices(domain.cell_edges(cell), iso_value, i_case, corners, values, vertices); + + std::lock_guard lock(mutex); + mc_construct_triangles(i_case, vertices, points, polygons); + } + + void add_triangle(const std::size_t p0, const std::size_t p1, const std::size_t p2) { + polygons.push_back({}); + auto& triangle = polygons.back(); + + triangle.push_back(p0); + triangle.push_back(p1); + triangle.push_back(p2); + } + + void p_slice(const Cell_handle& cell, const double i0, FT* values, Point* corners, const int i_case) { + // there are 12 edges, assign to each vertex three edges, the global edge numbering + // consist of 3*global_vertex_id + edge_offset. + const unsigned long long gei_pattern_ = 670526590282893600ull; + + // code edge end vertices for each of the 12 edges + const unsigned char l_edges_[12] = {16, 49, 50, 32, 84, 117, 118, 100, 64, 81, 115, 98}; + auto get_edge_vertex = [](const int e, unsigned int& v0, unsigned int& v1, const unsigned char l_edges_[12]) { + v0 = (unsigned int)(l_edges_[e] & 0xF); + v1 = (unsigned int)(l_edges_[e] >> 4) & 0xF; + }; + + // A hexahedron has twelve edges, save the intersection of the isosurface with the edge + // save global edge and global vertex index of isosurface + std::vector vertices(12); + // save loca coordinate along the edge of intersection point + std::vector ecoord{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // collect vertices + ushort flag{1}; + for (int eg = 0; eg < 12; eg++) { + if (flag & Cube_table::intersected_edges[i_case]) { + // the edge global index is given by the vertex global index + the edge offset + // uint shift = 5 * eg; + // const int ix = i_index + (int)((gei_pattern_ >> shift) & 1); // global_edge_id[eg][0]; + // const int iy = j_index + (int)((gei_pattern_ >> (shift + 1)) & 1); // global_edge_id[eg][1]; + // const int iz = k_index + (int)((gei_pattern_ >> (shift + 2)) & 1); // global_edge_id[eg][2]; + // const int off_val = (int)((gei_pattern_ >> (shift + 3)) & 3); + + // int g_edg = int(m_cell_shift_factor * m_ugrid.global_index(ix, iy, iz) + off_val); + + // generate vertex here, do not care at this point if vertex already exist + uint v0, v1; + get_edge_vertex(eg, v0, v1, l_edges_); + + double l = (i0 - values[v0]) / (values[v1] - values[v0]); + ecoord[eg] = l; + // interpolate vertex + const FT px = (1 - l) * corners[v0][0] + l * corners[v1][0]; + const FT py = (1 - l) * corners[v0][1] + l * corners[v1][1]; + const FT pz = (1 - l) * corners[v0][2] + l * corners[v1][2]; + + // set vertex in map + // set vertex index + // auto s_index = m_vertices.find(vertices[eg].g_edg); + // if (s_index == m_vertices.end()) { + const int g_idx = (int)points.size(); + vertices[eg] = g_idx; + // m_vertices[vertices[eg].g_edg] = g_idx; + points.push_back(Point(px, py, pz)); + //} else { + // vertices[eg] = s_index->second; + //} + } + /*else { + e_set[eg] = false; + }*/ + // next edge + flag <<= 1; + } + + // compute oriented contours + // A countour consists of segment at the faces connecting the intersection of the + // iso-surface with the edges. For each edge we store the edge to which the segment + // is outgoing and the edge from which the segment in comming. Therefore a contour + // cab be reconstructed by connecting the edges in the direccion of the outgoing. + // The contour is oriented in such a way, that the positive vertices are outside. + // 1. build segments + // 2. connect segments + // build up segments + // set segments map + unsigned char segm_[12] = {0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF}; + auto set_segm = [](const int e, const int pos, const int val, unsigned char segm_[12]) { + if (pos == 0) { + segm_[e] &= 0xF0; + segm_[e] |= (unsigned char)val & 0xF; + } else if (pos == 1) { + segm_[e] &= 0xF; + segm_[e] |= val << 4; + } + }; + auto get_segm = [](const int e, const int pos, unsigned char segm_[12]) { + if (pos == 0) + return (int)(segm_[e] & 0xF); + else + return (int)((segm_[e] >> 4) & 0xF); + }; + auto is_segm_set = [](const int e, unsigned char segm_[12]) { return (segm_[e] != 0xFF); }; + auto unset_segm = [](const int e, unsigned char segm_[12]) { segm_[e] = 0xFF; }; + // In order to compute oriented segments, the hexahedron has to be flatten. + // The insides of the faces of the hexahedron have to be all at the same + // side of the flattend hexa. This requires changing the order of the + // edges when reading from the faces + // code edges at face + // unsigned short face_e_[6] = { 12816, 30292, 33936, 46754, 34739, 38305 }; + std::array e_face_{{291, 18277, 18696, 10859, 33719, 38305}}; + // code vertices at face + // unsigned short face_v_[6] = { 12816, 30292, 21520, 30258, 25632, 30001 }; + std::array v_face_{{12576, 25717, 5380, 29538, 8292, 30001}}; + + // reading edge from face + auto get_face_e = [e_face_](const int f, const int e) { return ((e_face_[f] >> (4 * e)) & 0xF); }; + auto get_face_v = [v_face_](const int f, const int e) { return ((v_face_[f] >> (4 * e)) & 0xF); }; + // compute oriented segments using the isoline scheme at the faces + const unsigned int BIT_1 = 1; + const unsigned int BIT_2 = 2; + const unsigned int BIT_3 = 4; + const unsigned int BIT_4 = 8; + auto asymptotic_decider = [](const double f0, const double f1, const double f2, const double f3) { + return (f0 * f3 - f1 * f2) / (f0 + f3 - f1 - f2); + }; + std::vector f_flag(6, false); + for (int f = 0; f < 6; f++) { + // classify face + unsigned int f_case{0}; + uint v0 = get_face_v(f, 0); + uint v1 = get_face_v(f, 1); + uint v2 = get_face_v(f, 2); + uint v3 = get_face_v(f, 3); + uint e0 = get_face_e(f, 0); + uint e1 = get_face_e(f, 1); + uint e2 = get_face_e(f, 2); + uint e3 = get_face_e(f, 3); + double f0 = values[v0]; + double f1 = values[v1]; + double f2 = values[v2]; + double f3 = values[v3]; + if (f0 >= i0) f_case |= BIT_1; + if (f1 >= i0) f_case |= BIT_2; + if (f2 >= i0) f_case |= BIT_3; + if (f3 >= i0) f_case |= BIT_4; + switch (f_case) { + case 1: + set_segm(e0, 0, e3, segm_); + set_segm(e3, 1, e0, segm_); + break; + case 2: + set_segm(e1, 0, e0, segm_); + set_segm(e0, 1, e1, segm_); + break; + case 3: + set_segm(e1, 0, e3, segm_); + set_segm(e3, 1, e1, segm_); + break; + case 4: + set_segm(e3, 0, e2, segm_); + set_segm(e2, 1, e3, segm_); + break; + case 5: + set_segm(e0, 0, e2, segm_); + set_segm(e2, 1, e0, segm_); + break; + case 6: { + const double val = asymptotic_decider(f0, f1, f2, f3); + if (val > i0) { + set_segm(e3, 0, e0, segm_); + set_segm(e0, 1, e3, segm_); + set_segm(e1, 0, e2, segm_); + set_segm(e2, 1, e1, segm_); + } else if (val < i0) { + set_segm(e1, 0, e0, segm_); + set_segm(e0, 1, e1, segm_); + set_segm(e3, 0, e2, segm_); + set_segm(e2, 1, e3, segm_); + } else { + f_flag[f] = true; + // singular case val == i0, there are no asymptotes + // check if there is a reasonable triangulation of the face + unsigned short e_flag = 0x218; + unsigned short bit_1 = 0x1; + unsigned short bit_2 = 0x2; + double ec0 = ecoord[e0]; + double ec1 = ecoord[e1]; + double ec2 = ecoord[e2]; + double ec3 = ecoord[e3]; + if ((e_flag >> (f * 2)) & bit_1) { + ec0 = 1 - ec0; + ec2 = 1 - ec2; + } + if ((e_flag >> (f * 2)) & bit_2) { + ec1 = 1 - ec1; + ec3 = 1 - ec3; + } + if (ec1 < ec3 && ec0 > ec2) { + set_segm(e1, 0, e0, segm_); + set_segm(e0, 1, e1, segm_); + set_segm(e3, 0, e2, segm_); + set_segm(e2, 1, e3, segm_); + } else if (ec1 > ec3 && ec0 < ec2) { + set_segm(e3, 0, e0, segm_); + set_segm(e0, 1, e3, segm_); + set_segm(e1, 0, e2, segm_); + set_segm(e2, 1, e1, segm_); + } else { + std::cerr << "ERROR: can't correctly triangulate cell's face\n"; + return; + } + } + } break; + case 7: + set_segm(e1, 0, e2, segm_); + set_segm(e2, 1, e1, segm_); + break; + case 8: + set_segm(e2, 0, e1, segm_); + set_segm(e1, 1, e2, segm_); + break; + case 9: { + const double val = asymptotic_decider(f0, f1, f2, f3); + if (val > i0) { + set_segm(e0, 0, e1, segm_); + set_segm(e1, 1, e0, segm_); + set_segm(e2, 0, e3, segm_); + set_segm(e3, 1, e2, segm_); + } else if (val < i0) { + set_segm(e0, 0, e3, segm_); + set_segm(e3, 1, e0, segm_); + set_segm(e2, 0, e1, segm_); + set_segm(e1, 1, e2, segm_); + } else { + f_flag[f] = true; + // singular case val == i0, there are no asymptotes + // check if there is a reasonable triangulation of the face + unsigned short e_flag = 0x218; + unsigned short bit_1 = 0x1; + unsigned short bit_2 = 0x2; + double ec0 = ecoord[e0]; + double ec1 = ecoord[e1]; + double ec2 = ecoord[e2]; + double ec3 = ecoord[e3]; + if ((e_flag >> (f * 2)) & bit_1) { + ec0 = 1 - ec0; + ec2 = 1 - ec2; + } + if ((e_flag >> (f * 2)) & bit_2) { + ec1 = 1 - ec1; + ec3 = 1 - ec3; + } + if (ec1 < ec3 && ec0 > ec2) { + set_segm(e0, 0, e1, segm_); + set_segm(e1, 1, e0, segm_); + set_segm(e2, 0, e3, segm_); + set_segm(e3, 1, e2, segm_); + } else if (ec1 > ec3 && ec0 < ec2) { + set_segm(e0, 0, e3, segm_); + set_segm(e3, 1, e0, segm_); + set_segm(e2, 0, e1, segm_); + set_segm(e1, 1, e2, segm_); + } else { + std::cerr << "ERROR: can't correctly triangulate cell's face\n"; + return; + } + } + } break; + case 10: + set_segm(e2, 0, e0, segm_); + set_segm(e0, 1, e2, segm_); + + break; + case 11: + set_segm(e2, 0, e3, segm_); + set_segm(e3, 1, e2, segm_); + + break; + case 12: + set_segm(e3, 0, e1, segm_); + set_segm(e1, 1, e3, segm_); + + break; + case 13: + set_segm(e0, 0, e1, segm_); + set_segm(e1, 1, e0, segm_); + + break; + case 14: + set_segm(e3, 0, e0, segm_); + set_segm(e0, 1, e3, segm_); + break; + default: + break; + } + } + + // connect oriented segments into oriented contours + // closed contours are coded in 64 bit unsigned long long + // 1) Each entry has 4 bits + // 2) The first 4 entries are reserved for the size of the contours + // 3) The next 12 entries are the indices of the edges constituting the contorus + // The indices are numbers from 0 to 12 + unsigned long long c_ = 0xFFFFFFFFFFFF0000; + // in the 4 first bits store size of contours + auto get_cnt_size = [](const int cnt, unsigned long long& c_) { + return (size_t)((c_ & (0xF << 4 * cnt)) >> 4 * cnt); + }; + auto set_cnt_size = [](const int cnt, const int size, unsigned long long& c_) { + // unset contour size + c_ &= ~(0xF << 4 * cnt); + c_ |= (size << 4 * cnt); + }; + // set corresponging edge + auto set_c = [](const int cnt, const int pos, const int val, unsigned long long& c_) { + const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; + const uint c_sz = c_ & mask[cnt]; + const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); + c_ &= ~(((unsigned long long)0xF) << e); + c_ |= (((unsigned long long)val) << e); + }; + // read edge from contour + auto get_c = [](const int cnt, const int pos, unsigned long long c_) { + const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; + const uint c_sz = (uint)(c_ & mask[cnt]); + const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); + return (int)((c_ >> e) & 0xF); + }; + + + // connect oriented contours + uint cnt_{0}; + for (uint e = 0; e < 12; e++) { + if (is_segm_set(e, segm_)) { + uint eTo = get_segm(e, 0, segm_); + uint eIn = get_segm(e, 1, segm_); + uint eStart = e; + uint pos = 0; + set_c(cnt_, pos, eStart, c_); + while (eTo != eStart) { + pos = pos + 1; + set_c(cnt_, pos, eTo, c_); + eIn = eTo; + eTo = get_segm(eIn, 0, segm_); + unset_segm(eIn, segm_); + } + // set contour length + set_cnt_size(cnt_, pos + 1, c_); + // update number of contours + cnt_ = cnt_ + 1; + } + } + + // compute intersection of opposite faces + // It is enough to compute a pair of solutions for one face + // The other solutions are obtained by evaluating the equations + // for the common variable + double ui[2]{}; + double vi[2]{}; + double wi[2]{}; + unsigned char q_sol{0}; + const double a = (values[0] - values[1]) * (-values[6] + values[7] + values[4] - values[5]) - + (values[4] - values[5]) * (-values[2] + values[3] + values[0] - values[1]); + const double b = (i0 - values[0]) * (-values[6] + values[7] + values[4] - values[5]) + + (values[0] - values[1]) * (values[6] - values[4]) - + (i0 - values[4]) * (-values[2] + values[3] + values[0] - values[1]) - + (values[4] - values[5]) * (values[2] - values[0]); + const double c = (i0 - values[0]) * (values[6] - values[4]) - (i0 - values[4]) * (values[2] - values[0]); + ; + double d = b * b - 4 * a * c; + if (d > 0) { + d = std::sqrt(d); + // compute u-coord of solutions + ui[0] = (-b - d) / (2 * a); + ui[1] = (-b + d) / (2 * a); + // compute v-coord of solutions + double g1 = values[0] * (1 - ui[0]) + values[1] * ui[0]; + double g2 = values[2] * (1 - ui[0]) + values[3] * ui[0]; + vi[0] = (i0 - g1) / (g2 - g1); + if (std::isnan(vi[0]) || std::isinf(vi[0])) vi[0] = -1.f; + g1 = values[0] * (1 - ui[1]) + values[1] * ui[1]; + g2 = values[2] * (1 - ui[1]) + values[3] * ui[1]; + vi[1] = (i0 - g1) / (g2 - g1); + if (std::isnan(vi[1]) || std::isinf(vi[1])) vi[1] = -1.f; + // compute w-coordinates of solutions + g1 = values[0] * (1 - ui[0]) + values[1] * ui[0]; + g2 = values[4] * (1 - ui[0]) + values[5] * ui[0]; + wi[0] = (i0 - g1) / (g2 - g1); + if (std::isnan(wi[0]) || std::isinf(wi[0])) wi[0] = -1.f; + g1 = values[0] * (1 - ui[1]) + values[1] * ui[1]; + g2 = values[4] * (1 - ui[1]) + values[5] * ui[1]; + wi[1] = (i0 - g1) / (g2 - g1); + if (std::isnan(wi[1]) || std::isinf(wi[1])) wi[1] = -1.f; + // correct values for roots of quadratic equations + // in case the asymptotic decider has failed + if (f_flag[0] == true) { // face 1, w = 0; + if (wi[0] < wi[1]) + wi[0] = 0; + else + wi[1] = 0; + } + if (f_flag[1] == true) { // face 2, w = 1 + if (wi[0] > wi[1]) + wi[1] = 1; + else + wi[1] = 1; + } + if (f_flag[2] == true) { // face 3, v = 0 + if (vi[0] < vi[1]) + vi[0] = 0; + else + vi[1] = 0; + } + if (f_flag[3] == true) { // face 4, v = 1 + if (vi[0] > vi[1]) + vi[0] = 1; + else + vi[1] = 1; + } + if (f_flag[4] == true) { // face 5, u = 0 + if (ui[0] < ui[1]) + ui[0] = 0; + else + ui[1] = 0; + } + if (f_flag[5] == true) { // face 6, u = 1 + if (ui[0] > ui[1]) + ui[0] = 1; + else + ui[1] = 1; + } + + // check solution intervals + if (0 < ui[0] && ui[0] < 1) { + q_sol |= 1; + } + if (0 < ui[1] && ui[1] < 1) { + q_sol |= 2; + } + if (0 < vi[0] && vi[0] < 1) { + q_sol |= 4; + } + if (0 < vi[1] && vi[1] < 1) { + q_sol |= 8; + } + if (0 < wi[0] && wi[0] < 1) { + q_sol |= 16; + } + if (0 < wi[1] && wi[1] < 1) { + q_sol |= 32; + } + } + + // + // count the number of set bits + auto numberOfSetBits = [](const unsigned char n) { + // C or C++: use uint32_t + uint b = (uint)n; + b = b - ((b >> 1) & 0x55555555); + b = (b & 0x33333333) + ((b >> 2) & 0x33333333); + return (((b + (b >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24; + }; + // compute the number of solutions to the quadratic equation for a given face + auto nrQSolFace = [](const uint f, const unsigned char n) { + uint nr{0}; + switch (f) { + case 0: + if ((n & 0x5) == 0x5) nr = nr + 1; + if ((n & 0xA) == 0xA) nr = nr + 1; + break; + case 1: + if ((n & 0x11) == 0x11) nr = nr + 1; + if ((n & 0x22) == 0x22) nr = nr + 1; + break; + case 2: + if ((n & 0x18) == 0x18) nr = nr + 1; + if ((n & 0x24) == 0x24) nr = nr + 1; + break; + } + return nr; + }; + + + // triangulate contours + // if all bits are set, then there are three pairs of nontrivial solutions + // to the quadratic equations. In this case, there is a tunnel or a contour + // with 12 vertices. If there are three contours, then there is a tunnel and + // one of the contorus with only three vertices is not part of it. + if (numberOfSetBits(q_sol) == 6) { + // there are at most three contours + // Possible cases: + // 1) a single contour with 12 vertices + // 2) two contours which build a tunnel + // 3) three contours, one has only 3 vertices and does not belong to the tunnel + + // construct the six vertices of the inner hexagon + double hvt[6][3]; + hvt[0][0] = ui[0]; + hvt[0][1] = vi[0]; + hvt[0][2] = wi[0]; + hvt[1][0] = ui[0]; + hvt[1][1] = vi[0]; + hvt[1][2] = wi[1]; + hvt[2][0] = ui[1]; + hvt[2][1] = vi[0]; + hvt[2][2] = wi[1]; + hvt[3][0] = ui[1]; + hvt[3][1] = vi[1]; + hvt[3][2] = wi[1]; + hvt[4][0] = ui[1]; + hvt[4][1] = vi[1]; + hvt[4][2] = wi[0]; + hvt[5][0] = ui[0]; + hvt[5][1] = vi[1]; + hvt[5][2] = wi[0]; + + // construct vertices at intersections with the edges + auto e_vert = [&ecoord](const int e, const int i) { + const unsigned int l_coord[3]{1324855, 5299420, 16733440}; + unsigned char flag = (l_coord[i] >> (2 * e)) & 3; + if (flag == 3) + return ecoord[e]; + else + return (FT)(flag); + }; + + // if there are three contours, then there is a tunnel and one + // of the contours is not part of it. + unsigned char _not_tunnel = 0xF; + if (cnt_ == 3) { + // loop over the contorus + // triangulate the contour which is not part of + // the tunnel + const double uc_min = (ui[0] < ui[1]) ? ui[0] : ui[1]; + const double uc_max = (ui[0] < ui[1]) ? ui[1] : ui[0]; + for (int t = 0; t < (int)cnt_; t++) { + if (get_cnt_size(t, c_) == 3) { + double umin = 2; + double umax = -2; + uint e0 = get_c(t, 0, c_); + uint e1 = get_c(t, 1, c_); + uint e2 = get_c(t, 2, c_); + const double u_e0 = e_vert(e0, 0); + const double u_e1 = e_vert(e1, 0); + const double u_e2 = e_vert(e2, 0); + umin = (u_e0 < umin) ? u_e0 : umin; + umin = (u_e1 < umin) ? u_e1 : umin; + umin = (u_e2 < umin) ? u_e2 : umin; + umax = (u_e0 > umax) ? u_e0 : umax; + umax = (u_e1 > umax) ? u_e1 : umax; + umax = (u_e2 > umax) ? u_e1 : umax; + if (uc_min > umax || uc_max < umin) { + // this contour is not part of the tunnel + _not_tunnel = t; + + add_triangle(vertices[e0], vertices[e1], vertices[e2]); + } + } + } + } + + // compute vertices of inner hexagon, save new vertices in list and compute and keep + // global vertice index to build triangle connectivity later on. + uint tg_idx[6]; + for (int i = 0; i < 6; i++) { + + const double u = hvt[i][0]; + const double v = hvt[i][1]; + const double w = hvt[i][2]; + const FT px = (1 - w) * ((1 - v) * (corners[0][0] + u * (corners[1][0] - corners[0][0])) + + v * (corners[2][0] + u * (corners[3][0] - corners[2][0]))) + + w * ((1 - v) * (corners[4][0] + u * (corners[5][0] - corners[4][0])) + + v * (corners[6][0] + u * (corners[7][0] - corners[6][0]))); + const FT py = (1 - w) * ((1 - v) * (corners[0][1] + u * (corners[1][1] - corners[0][1])) + + v * (corners[2][1] + u * (corners[3][1] - corners[2][1]))) + + w * ((1 - v) * (corners[4][1] + u * (corners[5][1] - corners[4][1])) + + v * (corners[6][1] + u * (corners[7][1] - corners[6][1]))); + const FT pz = (1 - w) * ((1 - v) * (corners[0][2] + u * (corners[1][2] - corners[0][2])) + + v * (corners[2][2] + u * (corners[3][2] - corners[2][2]))) + + w * ((1 - v) * (corners[4][2] + u * (corners[5][2] - corners[4][2])) + + v * (corners[6][2] + u * (corners[7][2] - corners[6][2]))); + + tg_idx[i] = (uint)points.size(); + points.push_back(Point(px, py, pz)); + } + + // triangulate contours with inner hexagon + unsigned char tcon_[12]; + for (int i = 0; i < (int)cnt_; i++) { + if (_not_tunnel != i) { // contour belongs to tunnel + const int cnt_sz = (int)get_cnt_size(i, c_); + for (int r = 0; r < cnt_sz; r++) { + uint index = -1; + double dist = 1000.; + uint ci = get_c(i, r, c_); + const double u_edge = e_vert(ci, 0); + const double v_edge = e_vert(ci, 1); + const double w_edge = e_vert(ci, 2); + for (int s = 0; s < 6; s++) { + const double uval = u_edge - hvt[s][0]; + const double vval = v_edge - hvt[s][1]; + const double wval = w_edge - hvt[s][2]; + double val = uval * uval + vval * vval + wval * wval; + if (dist > val) { + index = s; + dist = val; + } + } + tcon_[ci] = (unsigned char)index; + } + + // correspondence between vertices found + // create triangles + // needs some functions + auto distanceRingIntsModulo = [](const int d1, const int d2) { + const int r = (d1 - d2) < 0 ? d2 - d1 : d1 - d2; + return (r > 2 ? 6 - r : r); + }; + auto midpointRingIntModulo = [](const int d1, const int d2) { + const int dmax = (d1 > d2) ? d1 : d2; + const int dmin = (d1 < d2) ? d1 : d2; + return ((dmax + 2) % 6 == dmin) ? (dmax + 1) % 6 : (dmax + dmin) / 2; + }; + + for (int r = 0; r < cnt_sz; r++) { + const uint tid1 = get_c(i, r, c_); + const uint tid2 = get_c(i, ((r + 1) % cnt_sz), c_); + const uint cid1 = tcon_[tid1]; + const uint cid2 = tcon_[tid2]; + // compute index distance + const int dst = distanceRingIntsModulo(cid1, cid2); + switch (dst) { + case 0: { + add_triangle(vertices[tid1], vertices[tid2], tg_idx[cid1]); + } break; + case 1: { + // measure diagonals + // triangulate along shortest diagonal + double u_edge = e_vert(tid1, 0); + double v_edge = e_vert(tid1, 1); + double w_edge = e_vert(tid1, 2); + const double l1 = (u_edge - hvt[cid2][0]) * (u_edge - hvt[cid2][0]) + + (v_edge - hvt[cid2][1]) * (v_edge - hvt[cid2][1]) + + (w_edge - hvt[cid2][2]) * (w_edge - hvt[cid2][2]); + u_edge = e_vert(tid2, 0); + v_edge = e_vert(tid2, 1); + w_edge = e_vert(tid2, 2); + const double l2 = (u_edge - hvt[cid1][0]) * (u_edge - hvt[cid1][0]) + + (v_edge - hvt[cid1][1]) * (v_edge - hvt[cid1][1]) + + (w_edge - hvt[cid1][2]) * (w_edge - hvt[cid1][2]); + + if (l1 < l2) { + add_triangle(vertices[tid1], vertices[tid2], tg_idx[cid2]); + add_triangle(vertices[tid1], tg_idx[cid2], tg_idx[cid1]); + } else { + add_triangle(vertices[tid1], vertices[tid2], tg_idx[cid1]); + add_triangle(vertices[tid2], tg_idx[cid2], tg_idx[cid1]); + } + } break; + case 2: { + const int cidm = midpointRingIntModulo(cid1, cid2); + + add_triangle(vertices[tid1], vertices[tid2], tg_idx[cidm]); + add_triangle(vertices[tid1], tg_idx[cidm], tg_idx[cid1]); + add_triangle(vertices[tid2], tg_idx[cid2], tg_idx[cidm]); + } break; + } // switch + } // for loop over the vertices of the contour + } // if (_not_tunnel) + } // for loop over contours + if (cnt_ == 1) { + // there is a single contour + // triangulate and close inner hexagon + // triangle must have the correct orientation + // use asymptotic_decider() to see if positive vertices + // are separated, in thic case orientation must be changed + const bool s_ = (asymptotic_decider(values[0], values[1], values[2], values[3]) <= i0); + const bool of_ = (wi[1] < wi[0]) ? s_ : !s_; + + if (!of_) { + add_triangle(tg_idx[0], tg_idx[2], tg_idx[1]); + add_triangle(tg_idx[2], tg_idx[4], tg_idx[3]); + add_triangle(tg_idx[0], tg_idx[5], tg_idx[4]); + add_triangle(tg_idx[0], tg_idx[4], tg_idx[2]); + } else { + add_triangle(tg_idx[0], tg_idx[1], tg_idx[2]); + add_triangle(tg_idx[2], tg_idx[3], tg_idx[4]); + add_triangle(tg_idx[0], tg_idx[4], tg_idx[5]); + add_triangle(tg_idx[0], tg_idx[2], tg_idx[4]); + } + } + } else { + // there is no tunnel + // handle case with no saddle point as simple polygons with 3, 4, 5 or six vertices + const unsigned char nr_u{(unsigned char)nrQSolFace(0, q_sol)}; + const unsigned char nr_v{(unsigned char)nrQSolFace(1, q_sol)}; + const unsigned char nr_w{(unsigned char)nrQSolFace(2, q_sol)}; + const unsigned char nr_t{(unsigned char)(nr_u + nr_v + nr_w)}; + if (nr_t == nr_u || nr_t == nr_v || nr_t == nr_w) { + // loop over all contours + for (int i = 0; i < (int)cnt_; i++) { + switch (get_cnt_size(i, c_)) { + case 3: { + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 1, c_)], + vertices[get_c(i, 2, c_)]); + } break; + case 4: { + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 1, c_)], + vertices[get_c(i, 2, c_)]); + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 2, c_)], + vertices[get_c(i, 3, c_)]); + } break; + case 5: { + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 1, c_)], + vertices[get_c(i, 2, c_)]); + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 2, c_)], + vertices[get_c(i, 3, c_)]); + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 3, c_)], + vertices[get_c(i, 4, c_)]); + } break; + case 6: { + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 1, c_)], + vertices[get_c(i, 3, c_)]); + add_triangle(vertices[get_c(i, 1, c_)], vertices[get_c(i, 2, c_)], + vertices[get_c(i, 3, c_)]); + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 3, c_)], + vertices[get_c(i, 4, c_)]); + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 4, c_)], + vertices[get_c(i, 5, c_)]); + } break; + } // switch over size of contour + } // loop over contorus + } // thre are no saddle points + else { + // there are saddle points + // fc1 = fs(1, 1)*fs(2, 1) + fs(1, 2)*fs(2, 2); + // fc2 = fs(1, 1)*fs(3, 1) + fs(1, 2)*fs(3, 2); + // fc3 = fs(2, 1)*fs(3, 2) + fs(2, 2)*fs(3, 1); + typedef u_char uchar; // TODO + + unsigned char fs[3][2]{{(uchar)(q_sol & 1), (uchar)((q_sol >> 1) & 1)}, + {(uchar)((q_sol >> 2) & 1), (uchar)((q_sol >> 3) & 1)}, + {(uchar)((q_sol >> 4) & 1), (uchar)((q_sol >> 5) & 1)}}; + + const unsigned char fc1 = fs[0][0] * fs[1][0] + fs[0][1] * fs[1][1]; + const unsigned char fc2 = fs[0][0] * fs[2][0] + fs[0][1] * fs[2][1]; + const unsigned char fc3 = fs[1][0] * fs[2][1] + fs[1][1] * fs[2][0]; + const unsigned char c_faces = fc1 + fc2 + fc3; + double ucoord{}; + double vcoord{}; + double wcoord{}; + switch (c_faces) { + case 2: { + if (fc1 == 0) { + ucoord = fs[0][0] * ui[0] + fs[0][1] * ui[1]; + vcoord = fs[1][0] * vi[0] + fs[1][1] * vi[1]; + wcoord = fs[1][0] * wi[1] + fs[1][1] * wi[0]; + } else if (fc2 == 0) { + ucoord = fs[0][0] * ui[0] + fs[0][1] * ui[1]; + vcoord = fs[0][0] * vi[0] + fs[0][1] * vi[1]; + wcoord = fs[0][0] * wi[1] + fs[0][1] * wi[0]; + } else if (fc3 == 0) { + ucoord = fs[1][0] * ui[0] + fs[1][1] * ui[1]; + vcoord = fs[1][0] * vi[0] + fs[1][1] * vi[1]; + wcoord = fs[1][0] * wi[0] + fs[1][1] * wi[1]; + } + } break; + case 3: { + ucoord = (fs[0][0] * ui[0] + fs[0][1] * ui[1]) / (fs[0][0] + fs[0][1]); + vcoord = (fs[1][0] * vi[0] + fs[1][1] * vi[1]) / (fs[1][0] + fs[1][1]); + wcoord = (fs[2][0] * wi[0] + fs[2][1] * wi[1]) / (fs[2][0] + fs[2][1]); + } break; + case 4: { + const unsigned char nr_u = fs[0][0] + fs[0][1]; + const unsigned char nr_v = fs[1][0] + fs[1][1]; + const unsigned char nr_w = fs[2][0] + fs[2][1]; + if (nr_w == 1) { + ucoord = fs[2][0] * ui[0] + fs[2][1] * ui[1]; + vcoord = fs[2][1] * vi[0] + fs[2][0] * vi[1]; + wcoord = fs[2][0] * wi[0] + fs[2][1] * wi[1]; + } else if (nr_v == 1) { + ucoord = fs[1][0] * ui[0] + fs[1][1] * ui[1]; + vcoord = fs[1][0] * vi[0] + fs[1][1] * vi[1]; + wcoord = fs[1][1] * wi[0] + fs[1][0] * wi[1]; + } else if (nr_u == 1) { + ucoord = fs[0][0] * ui[0] + fs[0][1] * ui[1]; + vcoord = fs[0][0] * vi[0] + fs[0][1] * vi[1]; + wcoord = fs[0][0] * wi[0] + fs[0][1] * wi[1]; + } + } break; + } // switch(c_faces) + + // create inner vertex + const FT px = + (1 - wcoord) * ((1 - vcoord) * (corners[0][0] + ucoord * (corners[1][0] - corners[0][0])) + + vcoord * (corners[2][0] + ucoord * (corners[3][0] - corners[2][0]))) + + wcoord * ((1 - vcoord) * (corners[4][0] + ucoord * (corners[5][0] - corners[4][0])) + + vcoord * (corners[6][0] + ucoord * (corners[7][0] - corners[6][0]))); + const FT py = + (1 - wcoord) * ((1 - vcoord) * (corners[0][1] + ucoord * (corners[1][1] - corners[0][1])) + + vcoord * (corners[2][1] + ucoord * (corners[3][1] - corners[2][1]))) + + wcoord * ((1 - vcoord) * (corners[4][1] + ucoord * (corners[5][1] - corners[4][1])) + + vcoord * (corners[6][1] + ucoord * (corners[7][1] - corners[6][1]))); + const FT pz = + (1 - wcoord) * ((1 - vcoord) * (corners[0][2] + ucoord * (corners[1][2] - corners[0][2])) + + vcoord * (corners[2][2] + ucoord * (corners[3][2] - corners[2][2]))) + + wcoord * ((1 - vcoord) * (corners[4][2] + ucoord * (corners[5][2] - corners[4][2])) + + vcoord * (corners[6][2] + ucoord * (corners[7][2] - corners[6][2]))); + + const uint g_index = (uint)points.size(); + // loop over the contorus + bool pt_used = false; + for (int i = 0; i < (int)cnt_; i++) { + const unsigned char cnt_sz = (unsigned char)get_cnt_size(i, c_); + if (cnt_sz == 3) { + add_triangle(vertices[get_c(i, 0, c_)], vertices[get_c(i, 1, c_)], vertices[get_c(i, 2, c_)]); + } else { + pt_used = true; + for (int t = 0; t < cnt_sz; t++) { + add_triangle(vertices[get_c(i, t, c_)], vertices[get_c(i, (t + 1) % cnt_sz, c_)], g_index); + } + } + } + if (pt_used) { + points.push_back(Point(px, py, pz)); + } + } // else - there are saddle points + } + } + +private: + const Domain& domain; + FT iso_value; + + Point_range& points; + Polygon_range& polygons; + + // compute a unique global index for vertices + // use as key the unique edge number + std::map vertex_map; + + std::mutex mutex; +}; + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_TMC_INTERNAL_TMC_H \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Marching_cubes.h b/Isosurfacing_3/include/CGAL/Marching_cubes_3.h similarity index 70% rename from Isosurfacing_3/include/CGAL/Marching_cubes.h rename to Isosurfacing_3/include/CGAL/Marching_cubes_3.h index 59ca4f0f646..efbe5f0429e 100644 --- a/Isosurfacing_3/include/CGAL/Marching_cubes.h +++ b/Isosurfacing_3/include/CGAL/Marching_cubes_3.h @@ -1,7 +1,9 @@ -#ifndef CGAL_MARCHING_CUBES_H -#define CGAL_MARCHING_CUBES_H +#ifndef CGAL_MARCHING_CUBES_3_H +#define CGAL_MARCHING_CUBES_3_H -#include "Isosurfacing_3/internal/Marching_cubes_3_internal.h" +#include +#include +#include namespace CGAL { namespace Isosurfacing { @@ -25,15 +27,17 @@ namespace Isosurfacing { * \param points points making the polygons of the soup * \param polygons each element in the vector describes a polygon using the indices of the points in points */ -template +template void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typename Domain_::FT iso_value, PointRange& points, PolygonRange& polygons) { - internal::Marching_cubes_RG2 functor(domain, iso_value, points, polygons); - domain.iterate_cells(functor, ConcurrencyTag()); + // static_assert(Domain_::CELL_TYPE & CUBICAL_CELL); + + internal::Marching_cubes_functor functor(domain, iso_value, points, polygons); + domain.iterate_cells(functor, Concurrency_tag()); } } // namespace Isosurfacing } // namespace CGAL -#endif // CGAL_MARCHING_CUBES_H +#endif // CGAL_MARCHING_CUBES_3_H diff --git a/Isosurfacing_3/include/CGAL/Marching_cubes_3_old.h b/Isosurfacing_3/include/CGAL/Marching_cubes_3_old.h new file mode 100644 index 00000000000..5b5906ba481 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Marching_cubes_3_old.h @@ -0,0 +1,275 @@ +#ifndef CGAL_MARCHING_CUBES_3_OLD_H +#define CGAL_MARCHING_CUBES_3_OLD_H + +#include +#include + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template +void marching_cubes_cell_old(const std::size_t x, const std::size_t y, const std::size_t z, const Domain_& domain, + const typename Domain_::FT iso_value, PointRange& points, PolygonRange& polygons, + std::mutex& mutex) { + + typedef std::array Idx_3; + typedef typename Domain_::FT FT; + typedef typename Domain_::Point_3 Point_3; + + const Idx_3 idx0 = {x + 0, y + 1, z + 0}; + const Idx_3 idx1 = {x + 1, y + 1, z + 0}; + const Idx_3 idx2 = {x + 1, y + 0, z + 0}; + const Idx_3 idx3 = {x + 0, y + 0, z + 0}; + const Idx_3 idx4 = {x + 0, y + 1, z + 1}; + const Idx_3 idx5 = {x + 1, y + 1, z + 1}; + const Idx_3 idx6 = {x + 1, y + 0, z + 1}; + const Idx_3 idx7 = {x + 0, y + 0, z + 1}; + + const Point_3 pos0 = domain.position(idx0[0], idx0[1], idx0[2]); + const Point_3 pos1 = domain.position(idx1[0], idx1[1], idx1[2]); + const Point_3 pos2 = domain.position(idx2[0], idx2[1], idx2[2]); + const Point_3 pos3 = domain.position(idx3[0], idx3[1], idx3[2]); + const Point_3 pos4 = domain.position(idx4[0], idx4[1], idx4[2]); + const Point_3 pos5 = domain.position(idx5[0], idx5[1], idx5[2]); + const Point_3 pos6 = domain.position(idx6[0], idx6[1], idx6[2]); + const Point_3 pos7 = domain.position(idx7[0], idx7[1], idx7[2]); + + const FT dist0 = domain.value(idx0[0], idx0[1], idx0[2]); + const FT dist1 = domain.value(idx1[0], idx1[1], idx1[2]); + const FT dist2 = domain.value(idx2[0], idx2[1], idx2[2]); + const FT dist3 = domain.value(idx3[0], idx3[1], idx3[2]); + const FT dist4 = domain.value(idx4[0], idx4[1], idx4[2]); + const FT dist5 = domain.value(idx5[0], idx5[1], idx5[2]); + const FT dist6 = domain.value(idx6[0], idx6[1], idx6[2]); + const FT dist7 = domain.value(idx7[0], idx7[1], idx7[2]); + + unsigned int cubeindex = 0; + // set bit if corresponding corner is below iso + cubeindex |= (dist0 < iso_value) << 0; + cubeindex |= (dist1 < iso_value) << 1; + cubeindex |= (dist2 < iso_value) << 2; + cubeindex |= (dist3 < iso_value) << 3; + cubeindex |= (dist4 < iso_value) << 4; + cubeindex |= (dist5 < iso_value) << 5; + cubeindex |= (dist6 < iso_value) << 6; + cubeindex |= (dist7 < iso_value) << 7; + + Point_3 vertlist[12]; + // bitmap of edges that intersect the iso-plane(s) + const int edges = edgeTable[cubeindex]; + + if (edges == 0) { + return; + } + + // interpolation of vertices incident to the edge, according to diagram + if (edges & (1 << 0)) { + vertlist[0] = vertex_interpolation(pos0, pos1, dist0, dist1, iso_value); + } + if (edges & (1 << 1)) { + vertlist[1] = vertex_interpolation(pos1, pos2, dist1, dist2, iso_value); + } + if (edges & (1 << 2)) { + vertlist[2] = vertex_interpolation(pos2, pos3, dist2, dist3, iso_value); + } + if (edges & (1 << 3)) { + vertlist[3] = vertex_interpolation(pos3, pos0, dist3, dist0, iso_value); + } + if (edges & (1 << 4)) { + vertlist[4] = vertex_interpolation(pos4, pos5, dist4, dist5, iso_value); + } + if (edges & (1 << 5)) { + vertlist[5] = vertex_interpolation(pos5, pos6, dist5, dist6, iso_value); + } + if (edges & (1 << 6)) { + vertlist[6] = vertex_interpolation(pos6, pos7, dist6, dist7, iso_value); + } + if (edges & (1 << 7)) { + vertlist[7] = vertex_interpolation(pos7, pos4, dist7, dist4, iso_value); + } + if (edges & (1 << 8)) { + vertlist[8] = vertex_interpolation(pos0, pos4, dist0, dist4, iso_value); + } + if (edges & (1 << 9)) { + vertlist[9] = vertex_interpolation(pos1, pos5, dist1, dist5, iso_value); + } + if (edges & (1 << 10)) { + vertlist[10] = vertex_interpolation(pos2, pos6, dist2, dist6, iso_value); + } + if (edges & (1 << 11)) { + vertlist[11] = vertex_interpolation(pos3, pos7, dist3, dist7, iso_value); + } + + // std::lock_guard lock(mutex); + + for (int i = 0; triTable[cubeindex][i] != -1; i += 3) { + + const Point_3& p0 = vertlist[triTable[cubeindex][i + 0]]; + const Point_3& p1 = vertlist[triTable[cubeindex][i + 1]]; + const Point_3& p2 = vertlist[triTable[cubeindex][i + 2]]; + + const std::size_t p0_idx = points.size(); // TODO: not allowed + + points.push_back(p0); + points.push_back(p1); + points.push_back(p2); + + polygons.push_back({}); + auto& triangle = polygons.back(); + + triangle.push_back(p0_idx + 0); + triangle.push_back(p0_idx + 1); + triangle.push_back(p0_idx + 2); + } +} + +template +void marching_cubes_cell_RG(const std::size_t x, const std::size_t y, const std::size_t z, const Domain_& domain, + const typename Domain_::FT iso_value, PointRange& points, PolygonRange& polygons, + std::mutex& mutex, std::unordered_map& vertex_map) { + + typedef typename Domain_::FT FT; + typedef typename Domain_::Point_3 Point; + + /// The structure _Vertex_ represents a vertex by giving its unique global index and the unique index of the edge. + struct Vertex { + std::size_t g_idx; // vertices; + + // slice hex + // collect function values and build index + FT values[8]; + Point corners[8]; + + int vi = 0; + std::bitset<8> index = 0; + for (int kl = 0; kl <= 1; kl++) { + for (int jl = 0; jl <= 1; jl++) { + for (int il = 0; il <= 1; il++) { + // collect scalar values and computex index + corners[vi] = domain.position(x + il, y + jl, z + kl); + values[vi] = domain.value(x + il, y + jl, z + kl); + + if (values[vi] >= iso_value) { + // index.set(VertexMapping[vi]); + index.set(vi); + } + // next cell vertex + vi++; + } + } + } + + // collect edges from table and + // interpolate triangle vertex positon + const int i_case = int(index.to_ullong()); + // compute for this case the vertices + ushort flag = 1; + for (int edge = 0; edge < 12; edge++) { + if (flag & Cube_table::intersected_edges[i_case]) { + // the edge global index is given by the vertex global index + the edge offset + const std::size_t ix = x + Cube_table::global_edge_id[edge][0]; + const std::size_t iy = y + Cube_table::global_edge_id[edge][1]; + const std::size_t iz = z + Cube_table::global_edge_id[edge][2]; + const std::size_t global_index = (iz * domain.size_y() * domain.size_x() + iy * domain.size_x() + ix); + vertices[edge].g_edg = cell_shift_factor * global_index + Cube_table::global_edge_id[edge][3]; + // generate vertex here, do not care at this point if vertex already exist + // interpolation weight + const int v0 = Cube_table::edge_to_vertex[edge][0]; + const int v1 = Cube_table::edge_to_vertex[edge][1]; + const FT l = (iso_value - values[v0]) / (values[v1] - values[v0]); + // interpolate vertex + const FT px = (1 - l) * corners[v0][0] + l * corners[v1][0]; + const FT py = (1 - l) * corners[v0][1] + l * corners[v1][1]; + const FT pz = (1 - l) * corners[v0][2] + l * corners[v1][2]; + const Point position(px, py, pz); + + // std::lock_guard lock(mutex); + // set vertex index + // const auto s_index = vertex_map.find(vertices[edge].g_edg); + // if (s_index == vertex_map.end()) { + // index not found! Add index to hash map + const std::size_t g_idx = points.size(); + // vertex_map[vertices[edge].g_edg] = g_idx; + vertices[edge].g_idx = g_idx; + points.push_back(position); + //} else { + // vertices[edge].g_idx = s_index->second; // this is vertex global index g_idx + //} + } + flag <<= 1; + } + + // std::lock_guard lock(mutex); + // construct triangles + for (int t = 0; t < 16; t += 3) { + const int t_index = i_case * 16 + t; + // if (e_tris_list[t_index] == 0x7f) + if (Cube_table::triangle_cases[t_index] == -1) break; + + const int eg0 = Cube_table::triangle_cases[t_index]; + const int eg1 = Cube_table::triangle_cases[t_index + 1]; + const int eg2 = Cube_table::triangle_cases[t_index + 2]; + + // insert new triangle in list + polygons.push_back({}); + auto& triangle = polygons.back(); + + triangle.push_back(vertices[eg0].g_idx); + triangle.push_back(vertices[eg1].g_idx); + triangle.push_back(vertices[eg2].g_idx); + } +} + +template +void make_triangle_mesh_using_marching_cubes_old(const Domain_& domain, const typename Domain_::FT iso_value, + PointRange& points, PolygonRange& polygons) { + + std::mutex mutex; + + const std::size_t size_k = domain.size_z(); + const std::size_t size_j = domain.size_y(); + const std::size_t size_i = domain.size_x(); + + const std::size_t blocking_size = 100; + + // TODO: look at polygon mesh processing for tbb (also linking) + + // compute a unique global index for vertices + // use as key the unique edge number + std::unordered_map v_map; + + + for (std::size_t bj = 0; bj < size_j - 1; bj += blocking_size) { + //#pragma omp parallel for + for (std::size_t k = 0; k < size_k - 1; k++) { + + const std::size_t j_start = bj; + const std::size_t j_end = std::min(size_j - 1, bj + blocking_size); + + for (std::size_t j = j_start; j < j_end; j++) { + for (std::size_t i = 0; i < size_i - 1; i++) { + // internal::marching_cubes_cell_old(i, j, k, domain, iso_value, points, polygons, mutex); + internal::marching_cubes_cell_RG(i, j, k, domain, iso_value, points, polygons, mutex, v_map); + } + } + } + } +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_MARCHING_CUBES_3_OLD_H diff --git a/Isosurfacing_3/include/CGAL/Octree_domain.h b/Isosurfacing_3/include/CGAL/Octree_domain.h new file mode 100644 index 00000000000..b40d864a9fa --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Octree_domain.h @@ -0,0 +1,135 @@ +#ifndef CGAL_OCTREE_DOMAIN_H +#define CGAL_OCTREE_DOMAIN_H + +#include +#include +#include + +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Octree_domain { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; + + typedef Octree_wrapper Octree; + typedef typename Octree::Vertex_handle Vertex_handle; + typedef typename Octree::Edge_handle Edge_handle; + typedef typename Octree::Voxel_handle Cell_handle; + + static constexpr Cell_type CELL_TYPE = CUBICAL_CELL; + static constexpr std::size_t VERTICES_PER_CELL = 8; + static constexpr std::size_t EDGES_PER_CELL = 12; + + typedef std::array Edge_vertices; + typedef std::array Cells_incident_to_edge; // TODO: not alwayys 4 + typedef std::array Cell_vertices; + typedef std::array Cell_edges; + +public: + Octree_domain(const Octree& octree) : octree_(&octree) {} + + Point position(const Vertex_handle& v) const { + return octree_->point(v); + } + + Vector gradient(const Vertex_handle& v) const { + return octree_->gradient(v); + } + + FT value(const Vertex_handle& v) const { + return octree_->vertex_value(v); + } + + Edge_vertices edge_vertices(const Edge_handle& e_id) const { + return octree_->edge_vertices(e_id); + } + + Cells_incident_to_edge cells_incident_to_edge(const Edge_handle& e_id) const { + return octree_->edge_voxels(e_id); + } + + Cell_vertices cell_vertices(const Cell_handle& vox) const { + return octree_->voxel_vertices(vox); + } + + Cell_edges cell_edges(const Cell_handle& vox) const { + return octree_->voxel_edges(vox); + } + + template + void iterate_vertices(Functor& f, Sequential_tag = Sequential_tag()) const { + for (const Vertex_handle& v : octree_->leaf_vertices()) { + f(v); + } + } + + template + void iterate_edges(Functor& f, Sequential_tag = Sequential_tag()) const { + for (const Edge_handle& e : octree_->leaf_edges()) { + f(e); + } + } + + template + void iterate_cells(Functor& f, Sequential_tag = Sequential_tag()) const { + for (const Cell_handle& v : octree_->leaf_voxels()) { + f(v); + } + } + +#ifdef CGAL_LINKED_WITH_TBB + template + void iterate_vertices(Functor& f, Parallel_tag) const { + const auto& vertices = octree_->leaf_vertices(); + + auto iterator = [&](const tbb::blocked_range& r) { + for (std::size_t i = r.begin(); i != r.end(); i++) { + f(vertices[i]); + } + }; + + tbb::parallel_for(tbb::blocked_range(0, vertices.size()), iterator); + } + + template + void iterate_edges(Functor& f, Parallel_tag) const { + const auto& edges = octree_->leaf_edges(); + + auto iterator = [&](const tbb::blocked_range& r) { + for (std::size_t i = r.begin(); i != r.end(); i++) { + f(edges[i]); + } + }; + + tbb::parallel_for(tbb::blocked_range(0, edges.size()), iterator); + } + + template + void iterate_cells(Functor& f, Parallel_tag) const { + const auto& cells = octree_->leaf_voxels(); + + auto iterator = [&](const tbb::blocked_range& r) { + for (std::size_t i = r.begin(); i != r.end(); i++) { + f(cells[i]); + } + }; + + tbb::parallel_for(tbb::blocked_range(0, cells.size()), iterator); + } +#endif // CGAL_LINKED_WITH_TBB + +private: + const Octree* octree_; +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_OCTREE_DOMAIN_H diff --git a/Isosurfacing_3/include/CGAL/Octree_wrapper.h b/Isosurfacing_3/include/CGAL/Octree_wrapper.h new file mode 100644 index 00000000000..e3f5445995c --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Octree_wrapper.h @@ -0,0 +1,511 @@ +#pragma once + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Octree_wrapper { + /* + * Naming convention from "A parallel dual marching cubes approach to quad + * only surface reconstruction - Grosso & Zint" + * + * ^ y + * | + * v2------e2------v3 + * /| /| + * e11| e10| + * / e3 / e1 + * v6------e6------v7 | + * | | | | + * | v0------e0--|---v1 --> x + * e7 / e5 / + * | e8 | e9 + * |/ |/ + * v4------e4------v5 + * / + * < z + */ + +public: + typedef GeomTraits Kernel; + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Vector_3 Vector_3; + + typedef CGAL::Octree> Octree; + + typedef std::size_t Vertex_handle; + typedef std::tuple Edge_handle; + typedef std::size_t Voxel_handle; + + typedef typename Octree::Node Node; + typedef typename Node::Global_coordinates Uniform_coords; // coordinates on max depth level + +private: + std::size_t max_depth_ = 0; + + FT offset_x_; + FT offset_y_; + FT offset_z_; + + CGAL::Bbox_3 bbox_; + + std::size_t dim_ = 1; + + FT hx_ = 0; + + std::vector point_range_; + Octree octree_; + + // std::set leaf_node_uniform_coordinates_; + std::vector leaf_voxels_; + std::vector leaf_edges_; + std::vector leaf_vertices_; + std::map vertex_values_; + std::map vertex_gradients_; + +public: + Octree_wrapper(const CGAL::Bbox_3& bbox) + : offset_x_(bbox.xmin()), + offset_y_(bbox.ymin()), + offset_z_(bbox.zmin()), + bbox_(bbox), + point_range_({{bbox.xmin(), bbox.ymin(), bbox.zmin()}, {bbox.xmax(), bbox.ymax(), bbox.zmax()}}), + octree_(point_range_) {} + + template + void refine(const Split_predicate& split_predicate) { + namespace Tables = internal::Cube_table; + + octree_.refine(split_predicate); + + max_depth_ = octree_.depth(); + dim_ = std::size_t(1) << max_depth_; + hx_ = bbox_.x_span() / dim_; + + // store leaf elements in sets + initialize value maps + std::set leaf_voxels_set; + std::set leaf_edges_set; + std::set leaf_vertices_set; + for (Node node : octree_.traverse(CGAL::Orthtrees::Leaves_traversal())) { + const auto& coords_uniform = uniform_coordinates(node); + // write all leaf nodes in a set + leaf_voxels_set.insert(lex_index(coords_uniform[0], coords_uniform[1], coords_uniform[2], max_depth_)); + + // init vertex values + for (int i = 0; i < Tables::N_VERTICES; ++i) { + Uniform_coords vuc = vertex_uniform_coordinates(node, i); + const auto lex = lex_index(vuc[0], vuc[1], vuc[2], max_depth_); + leaf_vertices_set.insert(lex); + vertex_values_[lex] = 0; + } + + // write all leaf edges in a set + const auto& coords_global = node.global_coordinates(); + const auto& depth = node.depth(); + const auto& df = depth_factor(node.depth()); + for (const auto& edge_voxels : Tables::edge_to_voxel_neighbor) { + bool are_all_voxels_leafs = true; + for (const auto& node_ijk : edge_voxels) { + const std::size_t x = coords_uniform[0] + df * node_ijk[0]; + const std::size_t y = coords_uniform[1] + df * node_ijk[1]; + const std::size_t z = coords_uniform[2] + df * node_ijk[2]; + // check for overflow / ignore edges on boundary + if (x >= dim_ || y >= dim_ || z >= dim_) { + are_all_voxels_leafs = false; + break; + } + + const Node n = get_node(x, y, z); + if (n.depth() > depth) { + are_all_voxels_leafs = false; + break; + } + } + if (are_all_voxels_leafs) { + // add to leaf edge set + std::size_t e_gl = + e_glIndex(edge_voxels[0][3], coords_global[0], coords_global[1], coords_global[2], depth); + leaf_edges_set.insert({e_gl, depth}); + } + } + } + + leaf_voxels_ = std::vector(leaf_voxels_set.begin(), leaf_voxels_set.end()); + leaf_edges_ = std::vector(leaf_edges_set.begin(), leaf_edges_set.end()); + leaf_vertices_ = std::vector(leaf_vertices_set.begin(), leaf_vertices_set.end()); + } + + std::size_t dim() const { + return dim_; + } + FT hx() const { + return hx_; + } + FT offset_x() const { + return offset_x_; + } + FT offset_y() const { + return offset_y_; + } + FT offset_z() const { + return offset_z_; + } + std::size_t max_depth() const { + return max_depth_; + } + + const std::vector& leaf_edges() const { + return leaf_edges_; + } + const std::vector& leaf_vertices() const { + return leaf_vertices_; + } + const std::vector& leaf_voxels() const { + return leaf_voxels_; + } + + FT value(const Vertex_handle& v) const { + return vertex_values_.at(v); + } + FT& value(const Vertex_handle& v) { + return vertex_values_[v]; + } + + Vector_3 gradient(const Vertex_handle& v) const { + return vertex_gradients_.at(v); + } + Vector_3& gradient(const Vertex_handle& v) { + return vertex_gradients_[v]; + } + + std::size_t depth_factor(const std::size_t& depth) const { + return std::size_t(1) << (max_depth_ - depth); + } + + Uniform_coords uniform_coordinates(const Node& node) const { + auto coords = node.global_coordinates(); + const std::size_t df = depth_factor(node.depth()); + for (int i = 0; i < Node::Dimension::value; ++i) { + coords[i] *= df; + } + + return coords; + } + + std::array node_points(const Node& node) const { + auto coords = node.global_coordinates(); + const std::size_t df = depth_factor(node.depth()); + + const FT x0 = offset_x_ + coords[0] * df * hx_; + const FT y0 = offset_y_ + coords[1] * df * hx_; + const FT z0 = offset_z_ + coords[2] * df * hx_; + const FT x1 = offset_x_ + (coords[0] + 1) * df * hx_; + const FT y1 = offset_y_ + (coords[1] + 1) * df * hx_; + const FT z1 = offset_z_ + (coords[2] + 1) * df * hx_; + + std::array points; + points[0] = {x0, y0, z0}; + points[1] = {x1, y0, z0}; + points[2] = {x0, y1, z0}; + points[3] = {x1, y1, z0}; + + points[4] = {x0, y0, z1}; + points[5] = {x1, y0, z1}; + points[6] = {x0, y1, z1}; + points[7] = {x1, y1, z1}; + + return points; + } + + Point_3 point(const Uniform_coords& vertex_coordinates) const { + const FT x0 = offset_x_ + vertex_coordinates[0] * hx_; + const FT y0 = offset_y_ + vertex_coordinates[1] * hx_; + const FT z0 = offset_z_ + vertex_coordinates[2] * hx_; + return {x0, y0, z0}; + } + + Point_3 point(const Vertex_handle& v) const { + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(v, max_depth_); + + const FT x0 = offset_x_ + i * hx_; + const FT y0 = offset_y_ + j * hx_; + const FT z0 = offset_z_ + k * hx_; + return {x0, y0, z0}; + } + + Uniform_coords vertex_uniform_coordinates(const Node& node, + const typename Node::Local_coordinates local_coords) const { + const auto node_coords = node.global_coordinates(); + auto v_coords = node_coords; + for (int i = 0; i < Node::Dimension::value; ++i) { + v_coords[i] += std::size_t(local_coords[i]); + } + + const auto df = depth_factor(node.depth()); + for (int i = 0; i < Node::Dimension::value; ++i) { + v_coords[i] *= df; + } + + return v_coords; + } + + Node get_node(const std::size_t& i, const std::size_t& j, const std::size_t& k) const { + Node node = octree_.root(); + const std::size_t& x = i; + const std::size_t& y = j; + const std::size_t& z = k; + while (!node.is_leaf()) { + std::size_t dist_to_max = max_depth_ - node.depth() - 1; + typename Node::Local_coordinates loc; + if (x & (std::size_t(1) << dist_to_max)) { + loc[0] = true; + } + if (y & (std::size_t(1) << dist_to_max)) { + loc[1] = true; + } + if (z & (std::size_t(1) << dist_to_max)) { + loc[2] = true; + } + node = node[loc.to_ulong()]; + } + return node; + } + + Node get_node(const std::size_t lex_index) const { + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(lex_index, max_depth_); + return get_node(i, j, k); + } + + std::size_t lex_index(const std::size_t& i, const std::size_t& j, const std::size_t& k, + const std::size_t& depth) const { + std::size_t dim = (std::size_t(1) << depth) + 1; + return k * dim * dim + j * dim + i; + } + + std::size_t i_index(const std::size_t& lex_index, const std::size_t& depth) const { + std::size_t dim = (std::size_t(1) << depth) + 1; + return lex_index % dim; + } + std::size_t j_index(const std::size_t& lex_index, const std::size_t& depth) const { + std::size_t dim = (std::size_t(1) << depth) + 1; + return ((lex_index / dim) % dim); + } + std::size_t k_index(const std::size_t& lex_index, const std::size_t& depth) const { + std::size_t dim = (std::size_t(1) << depth) + 1; + return (lex_index / (dim * dim)); + } + + std::tuple ijk_index(const std::size_t& lex_index, + const std::size_t& depth) const { + return std::make_tuple(i_index(lex_index, depth), j_index(lex_index, depth), k_index(lex_index, depth)); + } + + /// + /// compute unique edge global index. + /// + /// local edge index + /// i-index of cell + /// j-index of cell + /// k-index of cell + /// + std::size_t e_glIndex(const std::size_t& e, const std::size_t& i_idx, const std::size_t& j_idx, + const std::size_t& k_idx, const std::size_t& depth) const { + const unsigned long long gei_pattern_ = 670526590282893600ull; + const size_t i = i_idx + (size_t)((gei_pattern_ >> 5 * e) & 1); // global_edge_id[eg][0]; + const size_t j = j_idx + (size_t)((gei_pattern_ >> (5 * e + 1)) & 1); // global_edge_id[eg][1]; + const size_t k = k_idx + (size_t)((gei_pattern_ >> (5 * e + 2)) & 1); // global_edge_id[eg][2]; + const size_t offs = (size_t)((gei_pattern_ >> (5 * e + 3)) & 3); + return (3 * lex_index(i, j, k, depth) + offs); + } + + std::array voxel_values(const Voxel_handle& vox) const { + namespace Tables = internal::Cube_table; + + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(vox, max_depth_); + Node node = get_node(i, j, k); + const auto& df = depth_factor(node.depth()); + + std::array v; + for (int v_id = 0; v_id < Tables::N_VERTICES; ++v_id) { + const auto& l = Tables::local_vertex_position[v_id]; + const auto lex = lex_index(i + df * l[0], j + df * l[1], k + df * l[2], max_depth_); + v[v_id] = lex; + } + + std::array s; + std::transform(v.begin(), v.end(), s.begin(), [this](const auto& e) { return this->vertex_values_.at(e); }); + + return s; + } + + FT vertex_value(const Vertex_handle& v) const { + return vertex_values_.at(v); + } + + std::array voxel_edges(const Voxel_handle& vox) const { + namespace Tables = internal::Cube_table; + + std::array v = voxel_vertices(vox); + + std::array edges; + for (int e_id = 0; e_id < Tables::N_EDGES; ++e_id) { + edges[e_id] = {v[Tables::edge_to_vertex[e_id][0]], v[Tables::edge_to_vertex[e_id][1]]}; + } + + return edges; + } + + std::array voxel_vertices(const Voxel_handle& vox) const { + namespace Tables = internal::Cube_table; + + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(vox, max_depth_); + Node node = get_node(i, j, k); + const auto& df = depth_factor(node.depth()); + + std::array v; + for (int v_id = 0; v_id < Tables::N_VERTICES; ++v_id) { + const auto& l = Tables::local_vertex_position[v_id]; + const auto lex = lex_index(i + df * l[0], j + df * l[1], k + df * l[2], max_depth_); + v[v_id] = lex; + } + + return v; + } + + std::array voxel_gradients(const Voxel_handle& vox) const { + namespace Tables = internal::Cube_table; + + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(vox, max_depth_); + Node node = get_node(i, j, k); + const auto& df = depth_factor(node.depth()); + + std::array v; + for (int v_id = 0; v_id < Tables::N_VERTICES; ++v_id) { + const auto& l = Tables::local_vertex_position[v_id]; + const auto lex = lex_index(i + df * l[0], j + df * l[1], k + df * l[2], max_depth_); + v[v_id] = lex; + } + + std::array s; + std::transform(v.begin(), v.end(), s.begin(), [this](const auto& e) { return this->vertex_gradients_.at(e); }); + + return s; + } + std::array voxel_vertex_positions(const Voxel_handle& vox) const { + Node node = get_node(vox); + return node_points(node); + } + + /// + /// Get the values at the incident two vertices. Vertices are sorted in + /// ascending order. + /// + /// + /// + std::array edge_values(const Edge_handle& e_id) const { + namespace Tables = internal::Cube_table; + + std::size_t e_global_id, depth; + std::tie(e_global_id, depth) = e_id; + const auto df = depth_factor(depth); + + const size_t v0_lex_index = e_global_id / 3; + std::size_t i0, j0, k0; + std::tie(i0, j0, k0) = ijk_index(v0_lex_index, depth); + + // v1 + const std::size_t e_local_index = Tables::edge_store_index[e_global_id % 3]; + const auto& v1_local = Tables::local_vertex_position[Tables::edge_to_vertex[e_local_index][1]]; + + const std::size_t i1 = i0 + v1_local[0]; + const std::size_t j1 = j0 + v1_local[1]; + const std::size_t k1 = k0 + v1_local[2]; + + const auto v0 = lex_index(df * i0, df * j0, df * k0, max_depth_); + const auto v1 = lex_index(df * i1, df * j1, df * k1, max_depth_); + + return {value(v0), value(v1)}; + } + + std::array edge_vertices(const Edge_handle& e_id) const { + namespace Tables = internal::Cube_table; + + std::size_t e_global_id, depth; + std::tie(e_global_id, depth) = e_id; + const auto df = depth_factor(depth); + + const size_t v0_lex_index = e_global_id / 3; + std::size_t i0, j0, k0; + std::tie(i0, j0, k0) = ijk_index(v0_lex_index, depth); + + // v1 + const std::size_t e_local_index = Tables::edge_store_index[e_global_id % 3]; + const auto& v1_local = Tables::local_vertex_position[Tables::edge_to_vertex[e_local_index][1]]; + + const std::size_t i1 = i0 + v1_local[0]; + const std::size_t j1 = j0 + v1_local[1]; + const std::size_t k1 = k0 + v1_local[2]; + + const auto v0 = lex_index(df * i0, df * j0, df * k0, max_depth_); + const auto v1 = lex_index(df * i1, df * j1, df * k1, max_depth_); + + return {v0, v1}; + } + + /// + /// Get the 4 voxels incident to an edge. If an edge has only three incident + /// voxels, one will appear twice. The voxels are given with the uniform + /// lexicographical index. + /// + /// + /// + std::array edge_voxels(const Edge_handle& e_id) const { + namespace Tables = internal::Cube_table; + + std::size_t e_global_id, depth; + std::tie(e_global_id, depth) = e_id; + const std::size_t e_local_index = Tables::edge_store_index[e_global_id % 3]; + + const auto df = depth_factor(depth); + + const size_t v0_lex_index = e_global_id / 3; + std::size_t i, j, k; + std::tie(i, j, k) = ijk_index(v0_lex_index, depth); + i *= df; + j *= df; + k *= df; + + const auto& voxel_neighbors = Tables::edge_to_voxel_neighbor[e_local_index]; + Node n0 = get_node(i + voxel_neighbors[0][0], j + voxel_neighbors[0][1], k + voxel_neighbors[0][2]); + Node n1 = get_node(i + voxel_neighbors[1][0], j + voxel_neighbors[1][1], k + voxel_neighbors[1][2]); + Node n2 = get_node(i + voxel_neighbors[2][0], j + voxel_neighbors[2][1], k + voxel_neighbors[2][2]); + Node n3 = get_node(i + voxel_neighbors[3][0], j + voxel_neighbors[3][1], k + voxel_neighbors[3][2]); + + const Uniform_coords n0_uniform_coords = uniform_coordinates(n0); + const Uniform_coords n1_uniform_coords = uniform_coordinates(n1); + const Uniform_coords n2_uniform_coords = uniform_coordinates(n2); + const Uniform_coords n3_uniform_coords = uniform_coordinates(n3); + + std::size_t n0_lex = lex_index(n0_uniform_coords[0], n0_uniform_coords[1], n0_uniform_coords[2], max_depth_); + std::size_t n1_lex = lex_index(n1_uniform_coords[0], n1_uniform_coords[1], n1_uniform_coords[2], max_depth_); + std::size_t n2_lex = lex_index(n2_uniform_coords[0], n2_uniform_coords[1], n2_uniform_coords[2], max_depth_); + std::size_t n3_lex = lex_index(n3_uniform_coords[0], n3_uniform_coords[1], n3_uniform_coords[2], max_depth_); + + return {n0_lex, n1_lex, n2_lex, n3_lex}; + + // return { value( i0, j0, k0 ), value( i1, j1, k1 ) }; + } +}; + +} // namespace Isosurfacing +} // namespace CGAL diff --git a/Isosurfacing_3/include/CGAL/TC_marching_cubes_3.h b/Isosurfacing_3/include/CGAL/TC_marching_cubes_3.h new file mode 100644 index 00000000000..3e3e80b14e6 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/TC_marching_cubes_3.h @@ -0,0 +1,24 @@ +#ifndef CGAL_TMC_3_H +#define CGAL_TMC_3_H + +#include +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template +void make_triangle_mesh_using_tmc(const Domain_& domain, const typename Domain_::FT iso_value, PointRange& points, + PolygonRange& polygons) { + + // static_assert(Domain_::CELL_TYPE & CUBICAL_CELL); + + internal::TMC_functor functor(domain, iso_value, points, polygons); + domain.iterate_cells(functor, Concurrency_tag()); +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_TMC_3_H