From 20c6046a43d23854a2ba2fa4f79f5e14be679aa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ga=C3=ABl=20Guennebaud?= Date: Mon, 22 Jun 2009 15:06:41 +0000 Subject: [PATCH] add a marching cubes mesher in the demo (for comparison purpose) --- .gitattributes | 2 + .../Point_set_demo/CMakeLists.txt | 6 +- ...nt_set_demo_APSS_reconstruction_plugin.cpp | 18 +- ...int_set_demo_APSS_reconstruction_plugin.ui | 230 +++++-- ...emo_APSS_reconstruction_plugin_cgal_mc.cpp | 96 +++ .../Point_set_demo/marching_cubes.h | 574 ++++++++++++++++++ .../CGAL/APSS_reconstruction_function.h | 27 +- 7 files changed, 893 insertions(+), 60 deletions(-) create mode 100644 Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp create mode 100644 Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/marching_cubes.h diff --git a/.gitattributes b/.gitattributes index dfc4bb3726b..5fb290e4d3e 100644 --- a/.gitattributes +++ b/.gitattributes @@ -3968,6 +3968,7 @@ Surface_mesher/test/Surface_mesher/combined_spheres.cin -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/MainWindow.ui -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo.qrc -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.ui -text +Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_cleaning_plugin.cpp -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_local_spacing_plugin.cpp -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_normal_estimation_plugin.cpp -text @@ -3977,6 +3978,7 @@ Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_d Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/cgal_test_with_cmake -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/cgal_test_with_cmake.bat eol=crlf Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/data/ChineseDragon-points.off -text svneol=unset#application/octet-stream +Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/marching_cubes.h -text Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/resources/about.html svneol=native#text/html Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/resources/check-off.png -text svneol=unset#image/png Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/resources/check-off.svg -text diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt index c375df831cc..25b4367a0f6 100644 --- a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt @@ -236,7 +236,11 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) endmacro(polyhedron_demo_plugin) qt4_wrap_ui( APSS_UI_FILES Point_set_demo_APSS_reconstruction_plugin.ui) - polyhedron_demo_plugin(APSS_reconstruction_plugin Point_set_demo_APSS_reconstruction_plugin Point_set_demo_APSS_reconstruction_plugin_cgal_code.cpp ${APSS_UI_FILES}) + polyhedron_demo_plugin(APSS_reconstruction_plugin + Point_set_demo_APSS_reconstruction_plugin + Point_set_demo_APSS_reconstruction_plugin_cgal_code.cpp + Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp + ${APSS_UI_FILES}) target_link_libraries(APSS_reconstruction_plugin scene_polyhedron_item point_set) polyhedron_demo_plugin(inside_out_plugin Polyhedron_demo_inside_out_plugin) diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.cpp b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.cpp index a414cc90aea..ac964a18d4b 100644 --- a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.cpp +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.cpp @@ -21,6 +21,14 @@ Polyhedron* APSS_reconstruct(const Point_set& points, FT sm_distance, // Approximation error w.r.t. p.s.r.. For APSS: 0.015=fast, 0.003=smooth. FT smoothness = 2); // Smoothness factor. In the range 2 (clean datasets) and 8 (noisy datasets). +// same but using a marching cube +Polyhedron* APSS_reconstruct_mc(const Point_set& points, + FT sm_angle, // Min triangle angle (degrees). 20=fast, 30 guaranties convergence. + FT sm_radius, // Max triangle radius w.r.t. point set radius. 0.1 is fine. + FT sm_distance, // Approximation error w.r.t. p.s.r.. For APSS: 0.015=fast, 0.003=smooth. + FT smoothness, // Smoothness factor. In the range 2 (clean datasets) and 8 (noisy datasets). + int grid_size); // size of the grid + class Point_set_demo_APSS_reconstruction_plugin : public QObject, protected Polyhedron_demo_plugin_helper @@ -65,11 +73,12 @@ class Point_set_demo_APSS_reconstruction_plugin_dialog : public QDialog, private double triangleRadius() const { return m_inputRadius->value() * 0.01; } double triangleError() const { return m_inputDistance->value(); } double mlsSmoothness() const { return m_inputSmoothness->value(); } + bool useMarchingCube() const { return m_inputUseMC->isChecked(); } + int mcGridSize() const { return m_inputMcGridSize->value(); } private slots: void on_buttonBox_accepted() { -// std::cerr << "m_inputAngle.value() = " << m_inputAngle->value() << "\n"; } }; @@ -100,7 +109,12 @@ void Point_set_demo_APSS_reconstruction_plugin::reconstruct() QApplication::setOverrideCursor(Qt::WaitCursor); // Reconstruct point set as a polyhedron - Polyhedron* pRemesh = APSS_reconstruct(*points, sm_angle, sm_radius, sm_distance, smoothness); + Polyhedron* pRemesh = 0; + + if (dialog.useMarchingCube()) + pRemesh = APSS_reconstruct_mc(*points, sm_angle, sm_radius, sm_distance, smoothness, dialog.mcGridSize()); + else + pRemesh = APSS_reconstruct(*points, sm_angle, sm_radius, sm_distance, smoothness); if(pRemesh) { diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.ui b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.ui index 2e102e11560..cbc53af41eb 100644 --- a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.ui +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin.ui @@ -1,126 +1,166 @@ - + + ApssDialog - - + + 0 0 - 324 - 190 + 300 + 210 - + Dialog - - - - + + + + Min triangle angle: - - - + + + ° - + 1.000000000000000 - + 30.000000000000000 - + 20.000000000000000 - - - + + + Max triangle size: - - - + + + % - + 2 - + 1.000000000000000 - + 100.000000000000000 - + 1.000000000000000 - + 10.000000000000000 - - - + + + Approximation error: - - - + + + 6 - + 0.001000000000000 - + 0.010000000000000 - + 0.001000000000000 - + 0.003000000000000 - - - + + + Smoothness: - - - + + + 1.000000000000000 - + 100.000000000000000 - + 0.200000000000000 - + 2.000000000000000 - - - + + + + true + + + Use marching cubes + + + + + + + false + + + Marching cubes resolution: + + + + + + + false + + + 1 + + + 2000 + + + 50 + + + 300 + + + + + + Qt::Horizontal - - QDialogButtonBox::Cancel|QDialogButtonBox::NoButton|QDialogButtonBox::Ok + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok @@ -134,11 +174,11 @@ ApssDialog accept() - + 248 254 - + 157 274 @@ -150,15 +190,95 @@ ApssDialog reject() - + 316 260 - + 286 274 + + m_inputUseMC + toggled(bool) + m_inputMcGridSize + setEnabled(bool) + + + 131 + 261 + + + 206 + 312 + + + + + m_inputUseMC + toggled(bool) + label_6 + setEnabled(bool) + + + 149 + 135 + + + 106 + 163 + + + + + m_inputUseMC + toggled(bool) + m_inputAngle + setDisabled(bool) + + + 149 + 135 + + + 210 + 16 + + + + + m_inputUseMC + toggled(bool) + m_inputDistance + setDisabled(bool) + + + 149 + 135 + + + 210 + 76 + + + + + m_inputUseMC + toggled(bool) + m_inputRadius + setDisabled(bool) + + + 149 + 135 + + + 210 + 46 + + + diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp new file mode 100644 index 00000000000..93129390469 --- /dev/null +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/Point_set_demo_APSS_reconstruction_plugin_cgal_mc.cpp @@ -0,0 +1,96 @@ +//---------------------------------------------------------- +// APSS reconstruction method: +// Reconstruct a surface mesh from a point set and return it as a polyhedron. +//---------------------------------------------------------- + +#include "Kernel_type.h" +#include "Polyhedron_type.h" +#include "Point_set_scene_item.h" +#include "marching_cubes.h" + +// CGAL +#include +#include +#include +#include + +// This package +#include +#include + + +// APSS implicit function +typedef CGAL::APSS_reconstruction_function APSS_reconstruction_function; + +// APSS reconstruction method: +// Reconstruct a surface mesh from a point set and return it as a polyhedron. +Polyhedron* APSS_reconstruct_mc(const Point_set& points, + FT sm_angle, // Min triangle angle (degrees). 20 = fast, 30 guaranties convergence. + FT sm_radius, // Max triangle radius w.r.t. point set radius. 0.1 is fine. + FT sm_distance, // Approximation error w.r.t. p.s.r.. For APSS: 0.015 = fast, 0.003 = smooth. + FT smoothness, // smoothness factor + int grid_size) +{ + CGAL::Timer task_timer; task_timer.start(); + + //*************************************** + // Check requirements + //*************************************** + + int nb_points = points.size(); + if (nb_points == 0) + { + std::cerr << "Error: empty file" << std::endl; + return NULL; + } + + bool points_have_normals = (points.begin()->normal() != CGAL::NULL_VECTOR); + if ( ! points_have_normals ) + { + std::cerr << "Input point set not supported: this reconstruction method requires oriented normals" << std::endl; + return NULL; + } + + //*************************************** + // Creates implicit function + //*************************************** + + std::cerr << "Creates APSS implicit function (smoothness=" << smoothness << ")...\n"; + + // Creates implicit function + // Creates implicit function + APSS_reconstruction_function implicit_function(points.begin(), points.end(), + CGAL::make_normal_of_point_with_normal_pmap(points.begin()), + smoothness); + + // Prints status + std::cerr << "Creates implicit function: " << task_timer.time() << " seconds\n"; + task_timer.reset(); + + //*************************************** + // Surface mesh generation + //*************************************** + + std::cerr << "Marching cubes...\n"; + + Polyhedron* output_mesh = new Polyhedron; + marching_cubes(implicit_function, grid_size, *output_mesh); + + //*************************************** + // Erases small connected components + //*************************************** + +// std::cerr << "Erases small connected components...\n"; + +// unsigned int nb_erased_components = +// CGAL::keep_largest_connected_components(*output_mesh, 1/* keep largest component only*/); + + // Prints status +// std::cerr << "Erases small connected components: " << task_timer.time() << " seconds, " +// << nb_erased_components << " components erased" +// << std::endl; + task_timer.reset(); + + return output_mesh; +} + diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/marching_cubes.h b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/marching_cubes.h new file mode 100644 index 00000000000..8f24a0acd96 --- /dev/null +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/marching_cubes.h @@ -0,0 +1,574 @@ + +#ifndef MARCHINGCUBES_H +#define MARCHINGCUBES_H + +#include + +namespace CGALi +{ + +//************************************************************************ +// Static precalculated table +//************************************************************************ +int msEdgeTable[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}; + +int msTriTable[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}}; + +} // end of namespace CGALi + +namespace CGAL { + +template +class Marching_cubes : public CGAL::Modifier_base +{ + typedef typename Surface::Geom_traits Traits; + typedef typename Traits::FT FT; + typedef typename Traits::Point_3 Point_3; + typedef typename Traits::Vector_3 Vector_3; + typedef typename Surface::BBox BBox; + + // represent a vertices of the lattice + struct GridElement + { + Point_3 position; + FT value; + inline bool isInvalid() const { return value == Invalid; } + inline void markInvalid() { value = Invalid; } + }; + static const FT Invalid = FT(987654321); + + public: + + Marching_cubes(const Surface& s, int gs) : surface(s), grid_size(gs) {} + + void operator() (HDS& hds) + { + const FT iso_value = FT(0); + + // size of the blocks + static const int maxBlockSize = 128; + + // precomputed offsets to access the cell corners + static const int offsets[8] = { + 0, + 1, + 1+maxBlockSize*maxBlockSize, + maxBlockSize*maxBlockSize, + maxBlockSize, + 1+maxBlockSize, + 1+maxBlockSize+maxBlockSize*maxBlockSize, + maxBlockSize+maxBlockSize*maxBlockSize}; + + // returns the local extremity indices of an edge from its local index + static const int edge_to_vertex[12][2] = { + {0,1}, // 0 + {1,2}, // 1 + {2,3}, // 2 + {0,3}, // 3 + {4,5}, // 4 + {5,6}, // 5 + {6,7}, // 6 + {4,7}, // 7 + {0,4}, // 8 + {1,5}, // 9 + {2,6}, // 10 + {3,7}}; // 11 + + // get the bouding box: + BBox bbox = surface.bounding_box(); + + Vector_3 diag = Vector_3(bbox.xmax() - bbox.xmin(), bbox.ymax() - bbox.ymin(), bbox.zmax() - bbox.zmin()); + + if ( (diag.x()<=0.) + || (diag.y()<=0.) + || (diag.z()<=0.) + || (grid_size<=0)) + { + return; + } + + std::vector grid(maxBlockSize*maxBlockSize*maxBlockSize); + + // start a new 3D mesh + int vertex_count = 0; + Polyhedron_incremental_builder_3 B(hds, true); // true means verbose ?? + B.begin_surface(10,10,10); + + typedef std::pair EdgeKey; + + FT step = std::max( std::max(diag.x(), diag.y()), diag.z())/FT(grid_size); + + unsigned int nbCells[3]; + unsigned int nbBlocks[3]; + + for (int k=0 ; k<3 ; ++k) + { + nbCells[k] = int(diag[k]/step)+2; + nbBlocks[k] = nbCells[k]/maxBlockSize + ( (nbCells[k]%maxBlockSize)==0 ? 0 : 1); + } + + // for each macro block + uint bi[3]; // block id + for(bi[2]=0 ; bi[2] unique_vertex_map; + + // compute the size of the local grid + uint gridSize[3]; + for (uint k=0 ; k<3 ; ++k) + { + gridSize[k] = std::min(maxBlockSize, nbCells[k]-(maxBlockSize-1)*bi[k]); + } + Point_3 origin = Point_3(bbox.xmin(),bbox.ymin(),bbox.zmin()) + step * (maxBlockSize-1) * Vector_3(bi[0],bi[1],bi[2]); + + // fill the grid + uint ci[3]; // local cell id + + // for each vertex... + for(ci[0]=0 ; ci[0]v1) + std::swap(v0,v1); + + EdgeKey key(v0,v1); + std::map::iterator it = unique_vertex_map.find(key); + + if (it!=unique_vertex_map.end()) + { + //int count + // the vertex already exist + auxId[j] = it->second; + } + else + { + const Point_3& p = edges[::CGALi::msTriTable[mask][i+j]]; + // add a new vertex + auxId[j] = vertex_count; + countAddedVertex++; + B.add_vertex(p); + unique_vertex_map[key] = vertex_count; + vertex_count++; + } + } + if (auxId[0]!=auxId[1] && auxId[1]!=auxId[2] && auxId[2]!=auxId[0]) + { + B.begin_facet(); + for (uint j=0;j<3;++j) + B.add_vertex_to_facet(auxId[j]); + B.end_facet(); + } + } + } + } + } + } + B.end_surface(); + } + + protected: + inline Point_3 interpolEdge(const GridElement& v1, const GridElement& v2) + { + FT epsilon = 1e-6; // FIXME not very nice ;) + + if (std::abs(v1.value) < epsilon) + return v1.position; + if (std::abs(v2.value) < epsilon) + return v2.position; + if (std::abs(v1.value-v2.value) < epsilon) + return v1.position + FT(0.5) * (v2.position - v1.position); + + FT a = (-v1.value) / (v2.value - v1.value); + return v1.position + a * (v2.position - v1.position); + } + + inline bool is_finite(FT x) + { + return (x>=-std::numeric_limits::max()) && (x<=std::numeric_limits::max()); + } + + protected: + const Surface surface; + int grid_size; +}; + + + + +template +void marching_cubes(const Surface& surface, int grid_size, Polyhedron& target) +{ + Marching_cubes mc(surface, grid_size); + target.delegate(mc); +} + +} // end of namespace CGAL + +#endif // MARCHINGCUBES_H diff --git a/Surface_reconstruction_points_3/include/CGAL/APSS_reconstruction_function.h b/Surface_reconstruction_points_3/include/CGAL/APSS_reconstruction_function.h index 7573514e843..5a86f33b00b 100644 --- a/Surface_reconstruction_points_3/include/CGAL/APSS_reconstruction_function.h +++ b/Surface_reconstruction_points_3/include/CGAL/APSS_reconstruction_function.h @@ -33,7 +33,7 @@ CGAL_BEGIN_NAMESPACE -/// APSS_reconstruction_function computes an implicit function +/// APSS_reconstruction_function computes a signed distance field /// that defines a Point Set Surface (PSS) based on /// moving least squares (MLS) fitting of algebraic spheres. /// @@ -80,6 +80,7 @@ public: typedef typename Geom_traits::Point_3 Point; ///< typedef to Geom_traits::Point_3 typedef typename Geom_traits::Vector_3 Vector; ///< typedef to Geom_traits::Vector_3 typedef typename Geom_traits::Sphere_3 Sphere; ///< typedef to Geom_traits::Sphere_3 + typedef CGAL::Bbox_3 BBox; ///< typedef to CGAL::Bbox_3 // Private types private: @@ -186,6 +187,20 @@ private: Min_sphere_d< CGAL::Optimisation_d_traits_3 > ms3(first, beyond); m->bounding_sphere = Sphere(ms3.center(), ms3.squared_radius()); + // Compute axis aligned bounding box + // FIXME: I cannot believe there is nothing better to compute the bbox: + { + Point p = get(point_pmap,first); + m->bounding_box = BBox(p.x(), p.y(), p.z(), p.x(), p.y(), p.z()); + } + for (InputIterator it=first ; it != beyond ; ++it) + { + Point p = get(point_pmap,it); + m->bounding_box = m->bounding_box + BBox(p.x(), p.y(), p.z(), p.x(), p.y(), p.z()); + } + std::cerr << m->bounding_box.xmin() << " " << m->bounding_box.ymin() << " " << m->bounding_box.zmin() << " " + << m->bounding_box.xmax() << " " << m->bounding_box.ymax() << " " << m->bounding_box.zmax() << "\n"; + // Find a point inside the surface. find_inner_point(); @@ -298,6 +313,11 @@ public: return m->bounding_sphere; } + BBox bounding_box() const + { + return m->bounding_box; + } + private: /** Fit an algebraic sphere on a set of neigbors in a Moving Least Square sense. @@ -352,8 +372,10 @@ public: // Implementation note: this function is called a large number of times, // thus us heavily optimized. The bottleneck is Neighbor_search's constructor, // which we try to avoid calling. - FT operator()(const Point& p) const + FT operator()(const Point& p, bool* ok = 0) const { + if (ok) + *ok = true; // Is 'p' close to the surface? // Optimization: test first if 'p' is close to one of the neighbors // computed during the previous call. @@ -688,6 +710,7 @@ private: std::vector treeElements; std::vector radii; Sphere bounding_sphere; // Points' bounding sphere + BBox bounding_box; // Points' bounding box FT sqError; // Dichotomy error when projecting point (squared) unsigned int nofNeighbors; // Number of nearest neighbors Point inner_point; // Point inside the surface