add a marching cubes mesher in the demo (for comparison purpose)

This commit is contained in:
Gaël Guennebaud 2009-06-22 15:06:41 +00:00
parent 53c4fe7531
commit 20c6046a43
7 changed files with 893 additions and 60 deletions

2
.gitattributes vendored
View File

@ -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

View File

@ -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)

View File

@ -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)
{

View File

@ -1,126 +1,166 @@
<ui version="4.0" >
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>ApssDialog</class>
<widget class="QDialog" name="ApssDialog" >
<property name="geometry" >
<widget class="QDialog" name="ApssDialog">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>324</width>
<height>190</height>
<width>300</width>
<height>210</height>
</rect>
</property>
<property name="windowTitle" >
<property name="windowTitle">
<string>Dialog</string>
</property>
<layout class="QGridLayout" >
<item row="0" column="0" >
<widget class="QLabel" name="label" >
<property name="text" >
<layout class="QGridLayout" name="gridLayout">
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Min triangle angle:</string>
</property>
</widget>
</item>
<item row="0" column="1" >
<widget class="QDoubleSpinBox" name="m_inputAngle" >
<property name="suffix" >
<item row="0" column="1" colspan="2">
<widget class="QDoubleSpinBox" name="m_inputAngle">
<property name="suffix">
<string> °</string>
</property>
<property name="minimum" >
<property name="minimum">
<double>1.000000000000000</double>
</property>
<property name="maximum" >
<property name="maximum">
<double>30.000000000000000</double>
</property>
<property name="value" >
<property name="value">
<double>20.000000000000000</double>
</property>
</widget>
</item>
<item row="1" column="0" >
<widget class="QLabel" name="label_2" >
<property name="text" >
<item row="1" column="0">
<widget class="QLabel" name="label_2">
<property name="text">
<string>Max triangle size:</string>
</property>
</widget>
</item>
<item row="1" column="1" >
<widget class="QDoubleSpinBox" name="m_inputRadius" >
<property name="suffix" >
<item row="1" column="1" colspan="2">
<widget class="QDoubleSpinBox" name="m_inputRadius">
<property name="suffix">
<string> %</string>
</property>
<property name="decimals" >
<property name="decimals">
<number>2</number>
</property>
<property name="minimum" >
<property name="minimum">
<double>1.000000000000000</double>
</property>
<property name="maximum" >
<property name="maximum">
<double>100.000000000000000</double>
</property>
<property name="singleStep" >
<property name="singleStep">
<double>1.000000000000000</double>
</property>
<property name="value" >
<property name="value">
<double>10.000000000000000</double>
</property>
</widget>
</item>
<item row="2" column="0" >
<widget class="QLabel" name="label_3" >
<property name="text" >
<item row="2" column="0">
<widget class="QLabel" name="label_3">
<property name="text">
<string>Approximation error:</string>
</property>
</widget>
</item>
<item row="2" column="1" >
<widget class="QDoubleSpinBox" name="m_inputDistance" >
<property name="decimals" >
<item row="2" column="1" colspan="2">
<widget class="QDoubleSpinBox" name="m_inputDistance">
<property name="decimals">
<number>6</number>
</property>
<property name="minimum" >
<property name="minimum">
<double>0.001000000000000</double>
</property>
<property name="maximum" >
<property name="maximum">
<double>0.010000000000000</double>
</property>
<property name="singleStep" >
<property name="singleStep">
<double>0.001000000000000</double>
</property>
<property name="value" >
<property name="value">
<double>0.003000000000000</double>
</property>
</widget>
</item>
<item row="3" column="0" >
<widget class="QLabel" name="label_4" >
<property name="text" >
<item row="3" column="0">
<widget class="QLabel" name="label_4">
<property name="text">
<string>Smoothness:</string>
</property>
</widget>
</item>
<item row="3" column="1" >
<widget class="QDoubleSpinBox" name="m_inputSmoothness" >
<property name="minimum" >
<item row="3" column="1" colspan="2">
<widget class="QDoubleSpinBox" name="m_inputSmoothness">
<property name="minimum">
<double>1.000000000000000</double>
</property>
<property name="maximum" >
<property name="maximum">
<double>100.000000000000000</double>
</property>
<property name="singleStep" >
<property name="singleStep">
<double>0.200000000000000</double>
</property>
<property name="value" >
<property name="value">
<double>2.000000000000000</double>
</property>
</widget>
</item>
<item row="4" column="1" >
<widget class="QDialogButtonBox" name="buttonBox" >
<property name="orientation" >
<item row="4" column="0" colspan="3">
<widget class="QCheckBox" name="m_inputUseMC">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>Use marching cubes</string>
</property>
</widget>
</item>
<item row="5" column="0" colspan="2">
<widget class="QLabel" name="label_6">
<property name="enabled">
<bool>false</bool>
</property>
<property name="text">
<string>Marching cubes resolution:</string>
</property>
</widget>
</item>
<item row="5" column="2">
<widget class="QSpinBox" name="m_inputMcGridSize">
<property name="enabled">
<bool>false</bool>
</property>
<property name="minimum">
<number>1</number>
</property>
<property name="maximum">
<number>2000</number>
</property>
<property name="singleStep">
<number>50</number>
</property>
<property name="value">
<number>300</number>
</property>
</widget>
</item>
<item row="6" column="1" colspan="2">
<widget class="QDialogButtonBox" name="buttonBox">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="standardButtons" >
<set>QDialogButtonBox::Cancel|QDialogButtonBox::NoButton|QDialogButtonBox::Ok</set>
<property name="standardButtons">
<set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
</property>
</widget>
</item>
@ -134,11 +174,11 @@
<receiver>ApssDialog</receiver>
<slot>accept()</slot>
<hints>
<hint type="sourcelabel" >
<hint type="sourcelabel">
<x>248</x>
<y>254</y>
</hint>
<hint type="destinationlabel" >
<hint type="destinationlabel">
<x>157</x>
<y>274</y>
</hint>
@ -150,15 +190,95 @@
<receiver>ApssDialog</receiver>
<slot>reject()</slot>
<hints>
<hint type="sourcelabel" >
<hint type="sourcelabel">
<x>316</x>
<y>260</y>
</hint>
<hint type="destinationlabel" >
<hint type="destinationlabel">
<x>286</x>
<y>274</y>
</hint>
</hints>
</connection>
<connection>
<sender>m_inputUseMC</sender>
<signal>toggled(bool)</signal>
<receiver>m_inputMcGridSize</receiver>
<slot>setEnabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>131</x>
<y>261</y>
</hint>
<hint type="destinationlabel">
<x>206</x>
<y>312</y>
</hint>
</hints>
</connection>
<connection>
<sender>m_inputUseMC</sender>
<signal>toggled(bool)</signal>
<receiver>label_6</receiver>
<slot>setEnabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>149</x>
<y>135</y>
</hint>
<hint type="destinationlabel">
<x>106</x>
<y>163</y>
</hint>
</hints>
</connection>
<connection>
<sender>m_inputUseMC</sender>
<signal>toggled(bool)</signal>
<receiver>m_inputAngle</receiver>
<slot>setDisabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>149</x>
<y>135</y>
</hint>
<hint type="destinationlabel">
<x>210</x>
<y>16</y>
</hint>
</hints>
</connection>
<connection>
<sender>m_inputUseMC</sender>
<signal>toggled(bool)</signal>
<receiver>m_inputDistance</receiver>
<slot>setDisabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>149</x>
<y>135</y>
</hint>
<hint type="destinationlabel">
<x>210</x>
<y>76</y>
</hint>
</hints>
</connection>
<connection>
<sender>m_inputUseMC</sender>
<signal>toggled(bool)</signal>
<receiver>m_inputRadius</receiver>
<slot>setDisabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>149</x>
<y>135</y>
</hint>
<hint type="destinationlabel">
<x>210</x>
<y>46</y>
</hint>
</hints>
</connection>
</connections>
</ui>

View File

@ -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 <CGAL/Timer.h>
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
// This package
#include <CGAL/APSS_reconstruction_function.h>
#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
// APSS implicit function
typedef CGAL::APSS_reconstruction_function<Kernel> 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;
}

View File

@ -0,0 +1,574 @@
#ifndef MARCHINGCUBES_H
#define MARCHINGCUBES_H
#include <CGAL/Polyhedron_incremental_builder_3.h>
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 Surface, class HDS>
class Marching_cubes : public CGAL::Modifier_base<HDS>
{
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<GridElement> grid(maxBlockSize*maxBlockSize*maxBlockSize);
// start a new 3D mesh
int vertex_count = 0;
Polyhedron_incremental_builder_3<HDS> B(hds, true); // true means verbose ??
B.begin_surface(10,10,10);
typedef std::pair<int,int> 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]<nbBlocks[2] ; ++bi[2])
for(bi[1]=0 ; bi[1]<nbBlocks[1] ; ++bi[1])
for(bi[0]=0 ; bi[0]<nbBlocks[0] ; ++bi[0])
{
// record all vertices
std::map<EdgeKey, int> unique_vertex_map;
// compute the size of the local grid
uint gridSize[3];
for (uint k=0 ; k<3 ; ++k)
{
gridSize[k] = std::min<int>(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]<gridSize[0] ; ++ci[0])
for(ci[1]=0 ; ci[1]<gridSize[1] ; ++ci[1])
for(ci[2]=0 ; ci[2]<gridSize[2] ; ++ci[2])
{
GridElement& el = grid[(ci[2]*maxBlockSize + ci[1])*maxBlockSize + ci[0]];
el.position = origin + step*Vector_3(ci[0],ci[1],ci[2]);
bool ok;
el.value = surface(el.position, &ok);
if (!ok)
el.markInvalid();
}
// polygonize the grid (marching cube)
// for each cell...
for(ci[0]=0 ; ci[0]<gridSize[0]-1 ; ++ci[0])
for(ci[1]=0 ; ci[1]<gridSize[1]-1 ; ++ci[1])
for(ci[2]=0 ; ci[2]<gridSize[2]-1 ; ++ci[2])
{
uint cellId = ci[0]+maxBlockSize*(ci[1]+maxBlockSize*ci[2]);
// FIXME check if one corner is outside the surface definition domain
bool out = grid[cellId+offsets[0]].isInvalid() || grid[cellId+offsets[1]].isInvalid()
|| grid[cellId+offsets[2]].isInvalid() || grid[cellId+offsets[3]].isInvalid()
|| grid[cellId+offsets[4]].isInvalid() || grid[cellId+offsets[5]].isInvalid()
|| grid[cellId+offsets[6]].isInvalid() || grid[cellId+offsets[7]].isInvalid();
for (int k=0; k<8 && (!out); ++k)
out = out || (!is_finite(grid[cellId+offsets[k]].value));
if (!out)
{
// compute the mask
int mask = 0;
if (grid[cellId+offsets[0]].value <= iso_value) mask |= 1;
if (grid[cellId+offsets[1]].value <= iso_value) mask |= 2;
if (grid[cellId+offsets[2]].value <= iso_value) mask |= 4;
if (grid[cellId+offsets[3]].value <= iso_value) mask |= 8;
if (grid[cellId+offsets[4]].value <= iso_value) mask |= 16;
if (grid[cellId+offsets[5]].value <= iso_value) mask |= 32;
if (grid[cellId+offsets[6]].value <= iso_value) mask |= 64;
if (grid[cellId+offsets[7]].value <= iso_value) mask |= 128;
if (::CGALi::msEdgeTable[mask] != 0)
{
Point_3 edges[12];
if (::CGALi::msEdgeTable[mask] & 1)
edges[0] = interpolEdge(grid[cellId+offsets[0]],grid[cellId+offsets[1]]);
if (::CGALi::msEdgeTable[mask] & 2)
edges[1] = interpolEdge(grid[cellId+offsets[1]],grid[cellId+offsets[2]]);
if (::CGALi::msEdgeTable[mask] & 4)
edges[2] = interpolEdge(grid[cellId+offsets[2]],grid[cellId+offsets[3]]);
if (::CGALi::msEdgeTable[mask] & 8)
edges[3] = interpolEdge(grid[cellId+offsets[3]],grid[cellId+offsets[0]]);
if (::CGALi::msEdgeTable[mask] & 16)
edges[4] = interpolEdge(grid[cellId+offsets[4]],grid[cellId+offsets[5]]);
if (::CGALi::msEdgeTable[mask] & 32)
edges[5] = interpolEdge(grid[cellId+offsets[5]],grid[cellId+offsets[6]]);
if (::CGALi::msEdgeTable[mask] & 64)
edges[6] = interpolEdge(grid[cellId+offsets[6]],grid[cellId+offsets[7]]);
if (::CGALi::msEdgeTable[mask] & 128)
edges[7] = interpolEdge(grid[cellId+offsets[7]],grid[cellId+offsets[4]]);
if (::CGALi::msEdgeTable[mask] & 256)
edges[8] = interpolEdge(grid[cellId+offsets[0]],grid[cellId+offsets[4]]);
if (::CGALi::msEdgeTable[mask] & 512)
edges[9] = interpolEdge(grid[cellId+offsets[1]],grid[cellId+offsets[5]]);
if (::CGALi::msEdgeTable[mask] & 1024)
edges[10] = interpolEdge(grid[cellId+offsets[2]],grid[cellId+offsets[6]]);
if (::CGALi::msEdgeTable[mask] & 2048)
edges[11] = interpolEdge(grid[cellId+offsets[3]],grid[cellId+offsets[7]]);
for (int i=0 ; ::CGALi::msTriTable[mask][i]!=-1 ; i+=3)
{
int auxId[3];
int countAddedVertex = 0;
for (int j=0;j<3;++j)
{
int local_edge_id = ::CGALi::msTriTable[mask][i+j];
int v0 = cellId + offsets[edge_to_vertex[local_edge_id][0]];
int v1 = cellId + offsets[edge_to_vertex[local_edge_id][1]];
if (v0>v1)
std::swap(v0,v1);
EdgeKey key(v0,v1);
std::map<EdgeKey, int>::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<FT>::max()) && (x<=std::numeric_limits<FT>::max());
}
protected:
const Surface surface;
int grid_size;
};
template<typename Surface, typename Polyhedron>
void marching_cubes(const Surface& surface, int grid_size, Polyhedron& target)
{
Marching_cubes<Surface, typename Polyhedron::HalfedgeDS> mc(surface, grid_size);
target.delegate(mc);
}
} // end of namespace CGAL
#endif // MARCHINGCUBES_H

View File

@ -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<Gt> > 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<KdTreeElement> treeElements;
std::vector<FT> 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