songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Add alpha TC + point set protection + another test to solve inconsistencies 

Test: see check_and_solve_inconsistencies_by_filtering_simplices_out()
This commit is contained in:
Clement Jamin 2016-01-19 15:54:31 +01:00
parent 4016d0d247
commit 018ee66557
2 changed files with 599 additions and 221 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
Tangential_complex/benchmark/Tangential_complex/benchmark_tc.cpp
View File

@ -6,6 +6,9 @@
# define TBB_USE_THREADING_TOOL
#endif
#include <CGAL/Tangential_complex/protected_sets.h> // CJTODO TEST
#include <CGAL/assertions_behaviour.h>
#include <CGAL/Epick_d.h>
#include <CGAL/Tangential_complex.h>
@ -36,14 +39,17 @@ const char * const BENCHMARK_SCRIPT_FILENAME = "benchmark_script.txt";
typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_d Point;
typedef Kernel::Vector_d Vector;
typedef CGAL::Tangential_complex<
Kernel, CGAL::Dynamic_dimension_tag,
CGAL::Parallel_tag> TC;
//#define CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
//#define TC_PROTECT_POINT_SET_DELTA 0.003
//#define JUST_BENCHMARK_SPATIAL_SEARCH // CJTODO: test
//#define CHECK_IF_ALL_SIMPLICES_ARE_IN_THE_AMBIENT_DELAUNAY
#define TC_INPUT_STRIDES 3 // only take one point every TC_INPUT_STRIDES points
#define TC_NO_EXPORT
//#define TC_INPUT_STRIDES 10 // only take one point every TC_INPUT_STRIDES points
//#define TC_NO_EXPORT
#ifdef JUST_BENCHMARK_SPATIAL_SEARCH
std::ofstream spatial_search_csv_file("benchmark_spatial_search.csv");
@ -217,8 +223,10 @@ bool export_to_off(
}
void make_tc(std::vector<Point> &points,
TC::TS_container const& tangent_spaces, // can be empty
int intrinsic_dim,
double sparsity = 0.,
bool sparsify = true,
double sparsity = 0.01,
bool perturb = true,
bool add_high_dim_simpl = false,
bool collapse = false,
@ -226,6 +234,12 @@ void make_tc(std::vector<Point> &points,
const char *input_name = "tc")
{
Kernel k;
if (sparsify && !tangent_spaces.empty())
{
std::cerr << "Error: cannot sparsify point set with pre-computed normals.\n";
return;
}
// CJTODO TEMP TEST
//TC::Simplicial_complex compl;
@ -275,14 +289,14 @@ void make_tc(std::vector<Point> &points,
CGAL_TC_SET_PERFORMANCE_DATA("Num_points_in_input", points.size());
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
std::vector<Point> points_not_sparse = points;
#endif
//===========================================================================
// Sparsify point set if requested
//===========================================================================
if (sparsity != 0.)
if (sparsify)
{
std::size_t num_points_before = points.size();
points = sparsify_point_set(k, points, sparsity*sparsity);
@ -290,6 +304,27 @@ void make_tc(std::vector<Point> &points,
<< num_points_before << " / " << points.size() << std::endl;
}
#ifdef TC_PROTECT_POINT_SET_DELTA
// CJTODO TEST
# ifdef CGAL_TC_PROFILING
Wall_clock_timer t_protection;
# endif
std::vector<Point> points2;
std::vector<int> dummy;
std::vector<std::vector<int> > dummy2;
landmark_choice_protected_delaunay(
points, points.size(), points2, dummy, TC_PROTECT_POINT_SET_DELTA, dummy2, false, true);
points = points2;
# ifdef CGAL_TC_PROFILING
std::cerr << "Point set protected in " << t_protection.elapsed()
<< " seconds." << std::endl;
# endif
std::cerr << "Number of points after PROTECTION: " << points.size() << "\n";
#endif
CGAL_TC_SET_PERFORMANCE_DATA("Sparsity", sparsity);
CGAL_TC_SET_PERFORMANCE_DATA("Num_points", points.size());
@ -297,13 +332,18 @@ void make_tc(std::vector<Point> &points,
// Compute Tangential Complex
//===========================================================================
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
TC tc(points.begin(), points.end(), sparsity, intrinsic_dim,
points_not_sparse.begin(), points_not_sparse.end(), k);
#else
TC tc(points.begin(), points.end(), sparsity, intrinsic_dim, k);
#endif
if (!tangent_spaces.empty())
{
tc.set_tangent_planes(tangent_spaces);
}
double init_time = t.elapsed(); t.reset();
tc.compute_tangential_complex();
@ -356,7 +396,7 @@ void make_tc(std::vector<Point> &points,
//===========================================================================
t.reset();
double export_before_time =
(export_to_off(tc, input_name_stripped, "_BEFORE_FIX") ? t.elapsed() : -1);
(export_to_off(tc, input_name_stripped, "_INITIAL_TC") ? t.elapsed() : -1);
t.reset();
unsigned int num_perturb_steps = 0;
@ -406,6 +446,9 @@ void make_tc(std::vector<Point> &points,
CGAL_TC_SET_PERFORMANCE_DATA("Final_num_inconsistent_local_tr", "N/A");
}
// CJTODO TEST
//tc.check_and_solve_inconsistencies_by_filtering_simplices_out();
int max_dim = -1;
double fix2_time = -1;
double export_after_fix2_time = -1.;
@ -447,14 +490,14 @@ void make_tc(std::vector<Point> &points,
complex.display_stats();
// Export to OFF with higher-dim simplices colored
std::set<std::set<std::size_t> > higher_dim_simplices;
// CJTODO TEMP: Export to OFF with higher-dim simplices colored
/*std::set<std::set<std::size_t> > higher_dim_simplices;
complex.get_simplices_matching_test(
Test_dim(intrinsic_dim + 1),
std::inserter(higher_dim_simplices, higher_dim_simplices.begin()));
export_to_off(
tc, input_name_stripped, "_BEFORE_COLLAPSE", false, &complex,
&higher_dim_simplices);
&higher_dim_simplices);*/
//===========================================================================
// Collapse
@ -484,18 +527,23 @@ void make_tc(std::vector<Point> &points,
//===========================================================================
// Export to OFF
//===========================================================================
t.reset();
bool exported = export_to_off(
tc, input_name_stripped, "_AFTER_COLLAPSE", false, &complex,
&wrong_dim_simplices, &wrong_number_of_cofaces_simplices,
&unconnected_stars_simplices);
std::cerr
<< " OFF colors:" << std::endl
<< " * Red: wrong dim simplices" << std::endl
<< " * Green: wrong number of cofaces simplices" << std::endl
<< " * Blue: not-connected stars" << std::endl;
double export_after_collapse_time = (exported ? t.elapsed() : -1);
t.reset();
double export_after_collapse_time = -1.;
if (collapse)
{
t.reset();
bool exported = export_to_off(
tc, input_name_stripped, "_AFTER_COLLAPSE", false, &complex,
&wrong_dim_simplices, &wrong_number_of_cofaces_simplices,
&unconnected_stars_simplices);
std::cerr
<< " OFF colors:" << std::endl
<< " * Red: wrong dim simplices" << std::endl
<< " * Green: wrong number of cofaces simplices" << std::endl
<< " * Blue: not-connected stars" << std::endl;
double export_after_collapse_time = (exported ? t.elapsed() : -1.);
t.reset();
}
//===========================================================================
// Display info
@ -550,7 +598,7 @@ int main()
# ifdef _DEBUG
int num_threads = 1;
# else
int num_threads = 10;
int num_threads = 8;
# endif
#endif
@ -608,6 +656,7 @@ int main()
std::size_t num_points;
int ambient_dim;
int intrinsic_dim;
char sparsify;
double sparsity;
char perturb, add_high_dim_simpl, collapse;
double time_limit_for_perturb;
@ -619,6 +668,7 @@ int main()
sstr >> num_points;
sstr >> ambient_dim;
sstr >> intrinsic_dim;
sstr >> sparsify;
sstr >> sparsity;
sstr >> perturb;
sstr >> add_high_dim_simpl;
@ -657,6 +707,7 @@ int main()
#endif
std::vector<Point> points;
TC::TS_container tangent_spaces;
if (input == "generate_moment_curve")
{
@ -728,8 +779,9 @@ int main()
}
else
{
load_points_from_file<Point>(
input, std::back_inserter(points)/*, 600*/);
load_points_from_file<Kernel, typename TC::Tangent_space_basis>(
input, std::back_inserter(points),
std::back_inserter(tangent_spaces)/*, 600*/);
}
#ifdef CGAL_TC_PROFILING
@ -747,9 +799,9 @@ int main()
<< " points.\n"
<< "****************************************\n";
#endif
make_tc(points, intrinsic_dim, sparsity, // strided: CJTODO TEMP
perturb=='Y', add_high_dim_simpl=='Y', collapse=='Y',
time_limit_for_perturb, input.c_str());
make_tc(points, tangent_spaces, intrinsic_dim, sparsify=='Y',
sparsity, perturb=='Y', add_high_dim_simpl=='Y',
collapse=='Y', time_limit_for_perturb, input.c_str());
std::cerr << "TC #" << i++ << " done." << std::endl;
std::cerr << std::endl << "---------------------------------"

712
Tangential_complex/include/CGAL/Tangential_complex.h
View File

@ -37,6 +37,8 @@
#include <CGAL/Delaunay_triangulation.h>
#include <CGAL/Combination_enumerator.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
# include <CGAL/Mesh_3/Profiling_tools.h>
@ -47,6 +49,8 @@
#include <boost/optional.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/range/adaptor/transformed.hpp>
#include <boost/range/counting_range.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <vector>
@ -65,6 +69,12 @@
# include <tbb/mutex.h>
#endif
// choose exact integral type for QP solver
// (Gmpzf is not thread-safe)
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
//#define CGAL_QP_NO_ASSERTIONS // CJTODO: NECESSARY? http://doc.cgal.org/latest/QP_solver/group__PkgQPSolverFunctions.html#ga1fefbd0436aca0e281f88e8e6cd8eb74
//#define CGAL_TC_EXPORT_NORMALS // Only for 3D surfaces (k=2, d=3)
//static std::ofstream csv_stream("output/stats.csv"); // CJTODO TEMP
@ -98,7 +108,7 @@ private:
/// The class Tangential_complex represents a tangential complex
template <
typename Kernel_, // ambiant dimension
typename Kernel_, // ambiant kernel
typename DimensionTag, // intrinsic dimension
typename Concurrency_tag = CGAL::Parallel_tag,
#ifdef CGAL_ALPHA_TC
@ -149,9 +159,6 @@ class Tangential_complex
typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle;
typedef typename Tr_traits::Vector_d Tr_vector;
typedef Basis<K> Tangent_space_basis;
typedef Basis<K> Orthogonal_space_basis;
typedef std::vector<Point> Points;
#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_TC_GLOBAL_REFRESH)
typedef tbb::mutex Mutex_for_perturb;
@ -212,11 +219,14 @@ class Tangential_complex
};
public:
typedef typename std::vector<Tangent_space_basis> TS_container;
typedef typename std::vector<Orthogonal_space_basis> OS_container;
typedef Basis<K> Tangent_space_basis;
typedef Basis<K> Orthogonal_space_basis;
typedef std::vector<Tangent_space_basis> TS_container;
typedef std::vector<Orthogonal_space_basis> OS_container;
private:
typedef typename std::vector<Tr_and_VH> Tr_container;
typedef typename std::vector<Vector> Vectors;
typedef std::vector<Tr_and_VH> Tr_container;
typedef std::vector<Vector> Vectors;
// An Incident_simplex is the list of the vertex indices
// except the center vertex
@ -262,6 +272,22 @@ private:
{
return vh->point();
}
struct First_of_pair
{
template<typename> struct result;
template <typename F, typename Pair>
struct result<F(Pair)>
{
typedef typename boost::remove_reference<Pair>::type::first_type const& type;
};
template <typename Pair>
typename Pair::first_type const& operator()(Pair const& pair) const
{
return pair.first;
}
};
public:
typedef Tangential_complex_::Simplicial_complex Simplicial_complex;
@ -270,7 +296,7 @@ public:
template <typename InputIterator>
Tangential_complex(InputIterator first, InputIterator last,
double sparsity, int intrinsic_dimension,
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
InputIterator first_for_tse, InputIterator last_for_tse,
#endif
const K &k = K()
@ -295,7 +321,7 @@ public:
#ifdef CGAL_TC_EXPORT_NORMALS
, m_orth_spaces(m_points.size(), Orthogonal_space_basis())
#endif
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
, m_points_for_tse(first_for_tse, last_for_tse)
, m_points_ds_for_tse(m_points_for_tse)
#endif
@ -581,9 +607,10 @@ public:
<< stats_before.second
<< " (" << 100. * stats_before.second / stats_before.first << "%)"
<< std::endl
<< " * Num inconsistent stars: "
/*<< " * Num inconsistent stars: "
<< num_inconsistent_local_tr
<< " (" << 100. * num_inconsistent_local_tr / m_points.size() << "%)"
<< std::endl*/
<< std::endl
<< " * AFTER fix_inconsistencies_using_perturbation:" << std::endl
<< " - Total number of simplices in stars (incl. duplicates): "
@ -752,7 +779,46 @@ public:
<< 100. * stats_after.second / stats_after.first << "%"
<< std::endl;
}
// Returns true if some inconsistencies were found
bool check_and_solve_inconsistencies_by_filtering_simplices_out()
{
// CJTODO TEMP
std::pair<std::size_t, std::size_t> stats_before =
number_of_inconsistent_simplices(false);
std::cerr << "BEFORE check_and_solve_inconsistencies_by_filtering_simplices_out():\n"
<< " - Total number of simplices in stars (incl. duplicates): "
<< stats_before.first << std::endl
<< " - Num inconsistent simplices in stars (incl. duplicates): "
<< stats_before.second << std::endl
<< " - Percentage of inconsistencies: "
<< 100. * stats_before.second / stats_before.first << "%"
<< std::endl;
bool inconsistencies_found = false;
// CJTODO: parallel_for???
for (std::size_t idx = 0 ; idx < m_triangulations.size() ; ++idx)
{
if (filter_inconsistent_simplices_in_a_local_triangulation(idx))
inconsistencies_found = true;
}
// CJTODO TEMP
std::pair<std::size_t, std::size_t> stats_after =
number_of_inconsistent_simplices(false);
std::cerr << "AFTER check_and_solve_inconsistencies_by_filtering_simplices_out():\n"
<< " - Total number of simplices in stars (incl. duplicates): "
<< stats_after.first << std::endl
<< " - Num inconsistent simplices in stars (incl. duplicates): "
<< stats_after.second << std::endl
<< " - Percentage of inconsistencies: "
<< 100. * stats_after.second / stats_after.first << "%"
<< std::endl;
return inconsistencies_found;
}
#ifdef CGAL_ALPHA_TC
private:
@ -804,21 +870,71 @@ private:
// star(p) does not contain "s"
std::size_t p = *it_p;
// Compute the intersection between aff(Voronoi_cell(s)) and Tq
boost::optional<Tr_bare_point> intersection =
compute_aff_of_voronoi_face_and_tangent_subspace_intersection(
boost::optional<Tr_bare_point> intersection;
// If we know the intersection between Tq and voronoi(s) is a point
// I.e. dim(Tq) + dim(V(full_simplex)) = D
// <=> dim(Tq) = dim(full_simplex)
if (q_tr.current_dimension() == s.size())
{
// Compute the intersection between aff(Voronoi_cell(s)) and Tq
// Note that the computation in done in the sub-space defined by Tq
intersection =
compute_aff_of_voronoi_face_and_tangent_subspace_intersection(
q_tr.current_dimension(),
project_points_and_compute_weights(
full_simplex, m_tangent_spaces[i], q_tr.geom_traits()),
m_tangent_spaces[i],
q_tr.geom_traits());
// CJTODO: replace with an assertion?
if (!intersection)
{
std::cerr << "ERROR fill_pqueues_for_alpha_tc: "
"aff(Voronoi_cell(s)) and Tq do not intersect.\n";
continue;
}
}
else
{
std::cerr << "***************** Intersection is not a point ********************\n"; // CJTODO TEMP
// Compute the intersection between Voronoi_cell(s) and Tq
// We need Q, i.e. the vertices which are common neighbors of all
// vertices of full_simplex in the Delaunay triangulation, but we
// don't have, so we use the 5^d nearest neighbors instead
int num_neighbors_for_Q_approx = std::pow(
BASE_VALUE_FOR_ALPHA_TC_NEIGHBORHOOD, m_intrinsic_dim);
const Point center_pt = compute_perturbed_point(p);
KNS_range ins_range = m_points_ds.query_ANN(
center_pt, num_neighbors_for_Q_approx);
// CJTODO: optimize that, use a vector!
std::set<std::size_t> Q(
boost::make_counting_iterator(std::size_t(0)),
boost::make_counting_iterator(m_points.size()));
for (auto i : full_simplex)
Q.erase(i);
intersection =
compute_voronoi_face_and_tangent_subspace_intersection(
q_tr.current_dimension(),
project_points_and_compute_weights(
full_simplex, m_tangent_spaces[i], q_tr.geom_traits()),
full_simplex, m_tangent_spaces[i], q_tr.geom_traits()),
project_points_and_compute_weights(
Q,
//boost::counting_range(std::size_t(0), m_points.size() - 1),
//boost::adaptors::transform(ins_range, First_of_pair()),
m_tangent_spaces[i], q_tr.geom_traits()),
m_tangent_spaces[i],
q_tr.geom_traits());
// CJTODO: replace with an assertion?
if (!intersection)
{
std::cerr << "ERROR fill_pqueues_for_alpha_tc: "
"aff(Voronoi_cell(s)) and Tq do not intersect.\n";
continue;
// CJTODO: replace with an assertion?
if (!intersection)
{
std::cerr << "ERROR fill_pqueues_for_alpha_tc: "
"Voronoi_cell(s) and Tq do not intersect.\n";
continue;
}
}
// The following computations are done in the Euclidian space
@ -1436,6 +1552,205 @@ private:
return *s.rbegin() == std::numeric_limits<std::size_t>::max();
}
// Output: "triangulation" is a Regular Triangulation containing at least the
// star of "center_pt"
// Returns the handle of the center vertex
Tr_vertex_handle compute_star(
std::size_t i, const Point &center_pt, const Tangent_space_basis &tsb,
Triangulation &triangulation, bool verbose = false)
{
int tangent_space_dim = tsb.dimension();
const Tr_traits &local_tr_traits = triangulation.geom_traits();
Tr_vertex_handle center_vertex;
// Kernel functor & objects
typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
// Triangulation's traits functor & objects
typename Tr_traits::Point_weight_d point_weight =
local_tr_traits.point_weight_d_object();
typename Tr_traits::Power_center_d power_center =
local_tr_traits.power_center_d_object();
//***************************************************
// Build a minimal triangulation in the tangent space
// (we only need the star of p)
//***************************************************
// Insert p
Tr_point proj_wp;
if (i == tsb.origin())
{
// Insert {(0, 0, 0...), m_weights[i]}
proj_wp = local_tr_traits.construct_weighted_point_d_object()(
local_tr_traits.construct_point_d_object()(tangent_space_dim, ORIGIN),
m_weights[i]);
}
else
{
const Weighted_point& wp = compute_perturbed_weighted_point(i);
proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
}
center_vertex = triangulation.insert(proj_wp);
center_vertex->data() = i;
if (verbose)
std::cerr << "* Inserted point #" << i << std::endl;
#ifdef CGAL_TC_VERY_VERBOSE
std::size_t num_attempts_to_insert_points = 1;
std::size_t num_inserted_points = 1;
#endif
//const int NUM_NEIGHBORS = 150;
//KNS_range ins_range = m_points_ds.query_ANN(center_pt, NUM_NEIGHBORS);
INS_range ins_range = m_points_ds.query_incremental_ANN(center_pt);
// While building the local triangulation, we keep the radius
// of the sphere "star sphere" centered at "center_vertex"
// and which contains all the
// circumspheres of the star of "center_vertex"
boost::optional<FT> squared_star_sphere_radius_plus_margin;
#ifdef CGAL_ALPHA_TC
/*FT max_absolute_alpha = tsb.max_absolute_alpha();
// "2*m_half_sparsity" because both points can be perturbed
FT max_sqdist_to_tangent_space = (max_absolute_alpha == FT(0) ?
FT(0) : CGAL::square(2*max_absolute_alpha + 2*m_half_sparsity));*/
std::size_t number_of_attempts_to_insert_points = 0;
const std::size_t MAX_NUM_INSERTED_POINTS =
tsb.num_thickening_vectors() > 0 ?
static_cast<std::size_t>(std::pow(BASE_VALUE_FOR_ALPHA_TC_NEIGHBORHOOD, /*tsb.dimension()*/m_intrinsic_dim))
: std::numeric_limits<std::size_t>::max();
#endif
// Insert points until we find a point which is outside "star shere"
for (INS_iterator nn_it = ins_range.begin() ;
nn_it != ins_range.end() ;
++nn_it)
{
#ifdef CGAL_ALPHA_TC
++number_of_attempts_to_insert_points;
/*if (number_of_attempts_to_insert_points > MAX_NUM_INSERTED_POINTS)
break;*/
#endif
std::size_t neighbor_point_idx = nn_it->first;
// ith point = p, which is already inserted
if (neighbor_point_idx != i)
{
// No need to lock the Mutex_for_perturb here since this will not be
// called while other threads are perturbing the positions
Point neighbor_pt;
FT neighbor_weight;
compute_perturbed_weighted_point(
neighbor_point_idx, neighbor_pt, neighbor_weight);
if (squared_star_sphere_radius_plus_margin
&& k_sqdist(center_pt, neighbor_pt)
> *squared_star_sphere_radius_plus_margin)
break;
Tr_point proj_pt = project_point_and_compute_weight(
neighbor_pt, neighbor_weight, tsb,
local_tr_traits);
#ifdef CGAL_ALPHA_TC
/*for (int i = 0 ; i < tsb.num_thickening_vectors() ; ++i)
{
FT c = coord(proj_pt, m_intrinsic_dim + i);
if (c > 2.5*tsb.alpha_plus(i) + 2*m_half_sparsity //CJTODO: il faut 2*tsb.alpha_plus(i)
|| c < 2.5*tsb.alpha_minus(i) - 2*m_half_sparsity) //CJTODO: il faut 2*tsb.alpha_minus(i)
{
goto end_of_insert_loop; // "continue" in n-2 for-loop
}
}*/
/*if (max_sqdist_to_tangent_space != FT(0)
&& center_to_nbor_sqdist > max_sqdist_to_tangent_space)
break;*/
#endif
#ifdef CGAL_TC_VERY_VERBOSE
++num_attempts_to_insert_points;
#endif
Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex);
//Tr_vertex_handle vh = triangulation.insert(proj_pt);
if (vh != Tr_vertex_handle())
{
#ifdef CGAL_TC_VERY_VERBOSE
++num_inserted_points;
#endif
if (verbose)
std::cerr << "* Inserted point #" << neighbor_point_idx << std::endl;
vh->data() = neighbor_point_idx;
// Let's recompute squared_star_sphere_radius_plus_margin
if (triangulation.current_dimension() >= tangent_space_dim)
{
squared_star_sphere_radius_plus_margin = boost::none;
// Get the incident cells and look for the biggest circumsphere
std::vector<Tr_full_cell_handle> incident_cells;
triangulation.incident_full_cells(
center_vertex,
std::back_inserter(incident_cells));
for (typename std::vector<Tr_full_cell_handle>::iterator cit =
incident_cells.begin(); cit != incident_cells.end(); ++cit)
{
Tr_full_cell_handle cell = *cit;
if (triangulation.is_infinite(cell))
{
squared_star_sphere_radius_plus_margin = boost::none;
break;
}
else
{
Tr_point c = power_center(
boost::make_transform_iterator(
cell->vertices_begin(),
vertex_handle_to_point<Tr_point, Tr_vertex_handle>),
boost::make_transform_iterator(
cell->vertices_end(),
vertex_handle_to_point<Tr_point, Tr_vertex_handle>));
FT sq_power_sphere_diam = 4 * point_weight(c);
if (!squared_star_sphere_radius_plus_margin
|| sq_power_sphere_diam >
*squared_star_sphere_radius_plus_margin)
{
squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
}
}
}
// Let's add the margin, now
// The value depends on whether we perturb weight or position
if (squared_star_sphere_radius_plus_margin)
{
#ifdef CGAL_TC_PERTURB_WEIGHT
// "4*m_sq_half_sparsity" because both points can be perturbed
squared_star_sphere_radius_plus_margin =
*squared_star_sphere_radius_plus_margin + 4 * m_sq_half_sparsity;
#else
// "2*m_half_sparsity" because both points can be perturbed
squared_star_sphere_radius_plus_margin = CGAL::square(
CGAL::sqrt(*squared_star_sphere_radius_plus_margin)
+ 2 * m_half_sparsity);
#endif
}
}
}
}
}
return center_vertex;
}
void compute_tangent_triangulation(std::size_t i, bool verbose = false)
{
if (verbose)
@ -1484,176 +1799,9 @@ private:
int tangent_space_dim = tangent_basis_dim(i);
Triangulation &local_tr =
m_triangulations[i].construct_triangulation(tangent_space_dim);
const Tr_traits &local_tr_traits = local_tr.geom_traits();
Tr_vertex_handle &center_vertex = m_triangulations[i].center_vertex();
// Kernel functor & objects
typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
// Triangulation's traits functor & objects
typename Tr_traits::Point_weight_d point_weight =
local_tr_traits.point_weight_d_object();
typename Tr_traits::Power_center_d power_center =
local_tr_traits.power_center_d_object();
//***************************************************
// Build a minimal triangulation in the tangent space
// (we only need the star of p)
//***************************************************
// Insert p
Tr_point proj_wp;
if(i == tsb.origin())
{
proj_wp = local_tr_traits.construct_weighted_point_d_object()(
local_tr_traits.construct_point_d_object()(tangent_space_dim, ORIGIN),
m_weights[i]);
}
else
{
const Weighted_point& wp = compute_perturbed_weighted_point(i);
proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
}
center_vertex = local_tr.insert(proj_wp);
center_vertex->data() = i;
if (verbose)
std::cerr << "* Inserted point #" << i << std::endl;
#ifdef CGAL_TC_VERY_VERBOSE
std::size_t num_attempts_to_insert_points = 1;
std::size_t num_inserted_points = 1;
#endif
//const int NUM_NEIGHBORS = 150;
//KNS_range ins_range = m_points_ds.query_ANN(center_pt, NUM_NEIGHBORS);
INS_range ins_range = m_points_ds.query_incremental_ANN(center_pt);
// While building the local triangulation, we keep the radius
// of the sphere "star sphere" centered at "center_vertex"
// and which contains all the
// circumspheres of the star of "center_vertex"
boost::optional<FT> squared_star_sphere_radius_plus_margin;
#ifdef CGAL_ALPHA_TC
/*FT max_absolute_alpha = tsb.max_absolute_alpha();
// "2*m_half_sparsity" because both points can be perturbed
FT max_sqdist_to_tangent_space = (max_absolute_alpha == FT(0) ?
FT(0) : CGAL::square(2*max_absolute_alpha + 2*m_half_sparsity));*/
std::size_t number_of_attempts_to_insert_points = 0;
const std::size_t MAX_NUM_INSERTED_POINTS =
tsb.num_thickening_vectors() > 0 ?
static_cast<std::size_t>(std::pow(4, /*tsb.dimension()*/m_intrinsic_dim))
: std::numeric_limits<std::size_t>::max();
#endif
// Insert points until we find a point which is outside "star shere"
for (INS_iterator nn_it = ins_range.begin() ;
nn_it != ins_range.end() ;
++nn_it)
{
#ifdef CGAL_ALPHA_TC
++number_of_attempts_to_insert_points;
if (number_of_attempts_to_insert_points > MAX_NUM_INSERTED_POINTS)
break;
#endif
std::size_t neighbor_point_idx = nn_it->first;
// ith point = p, which is already inserted
if (neighbor_point_idx != i)
{
// No need to lock the Mutex_for_perturb here since this will not be
// called while other threads are perturbing the positions
Point neighbor_pt;
FT neighbor_weight;
compute_perturbed_weighted_point(
neighbor_point_idx, neighbor_pt, neighbor_weight);
if (squared_star_sphere_radius_plus_margin
&& k_sqdist(center_pt, neighbor_pt)
> *squared_star_sphere_radius_plus_margin)
break;
Tr_point proj_pt = project_point_and_compute_weight(
neighbor_pt, neighbor_weight, tsb,
local_tr_traits);
#ifdef CGAL_TC_VERY_VERBOSE
++num_attempts_to_insert_points;
#endif
Tr_vertex_handle vh = local_tr.insert_if_in_star(proj_pt, center_vertex);
//Tr_vertex_handle vh = local_tr.insert(proj_pt);
if (vh != Tr_vertex_handle())
{
#ifdef CGAL_TC_VERY_VERBOSE
++num_inserted_points;
#endif
if (verbose)
std::cerr << "* Inserted point #" << neighbor_point_idx << std::endl;
vh->data() = neighbor_point_idx;
// Let's recompute squared_star_sphere_radius_plus_margin
if (local_tr.current_dimension() >= tangent_space_dim)
{
squared_star_sphere_radius_plus_margin = boost::none;
// Get the incident cells and look for the biggest circumsphere
std::vector<Tr_full_cell_handle> incident_cells;
local_tr.incident_full_cells(
center_vertex,
std::back_inserter(incident_cells));
for (typename std::vector<Tr_full_cell_handle>::iterator cit =
incident_cells.begin(); cit != incident_cells.end(); ++cit)
{
Tr_full_cell_handle cell = *cit;
if (local_tr.is_infinite(cell))
{
squared_star_sphere_radius_plus_margin = boost::none;
break;
}
else
{
Tr_point c = power_center(
boost::make_transform_iterator(
cell->vertices_begin(),
vertex_handle_to_point<Tr_point, Tr_vertex_handle>),
boost::make_transform_iterator(
cell->vertices_end(),
vertex_handle_to_point<Tr_point, Tr_vertex_handle>));
FT sq_power_sphere_diam = 4*point_weight(c);
if (!squared_star_sphere_radius_plus_margin
|| sq_power_sphere_diam >
*squared_star_sphere_radius_plus_margin)
{
squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
}
}
}
// Let's add the margin, now
// The value depends on whether we perturb weight or position
if (squared_star_sphere_radius_plus_margin)
{
#ifdef CGAL_TC_PERTURB_WEIGHT
// "4*m_sq_half_sparsity" because both points can be perturbed
squared_star_sphere_radius_plus_margin =
*squared_star_sphere_radius_plus_margin + 4*m_sq_half_sparsity;
#else
// "2*m_half_sparsity" because both points can be perturbed
squared_star_sphere_radius_plus_margin = CGAL::square(
CGAL::sqrt(*squared_star_sphere_radius_plus_margin)
+ 2*m_half_sparsity);
#endif
}
}
}
}
}
m_triangulations[i].center_vertex() =
compute_star(i, center_pt, tsb, local_tr, verbose);
#if defined(CGAL_TC_PROFILING) && defined(CGAL_TC_VERY_VERBOSE)
std::cerr << " - triangulation construction: " << t.elapsed() << " s.\n";
@ -2003,7 +2151,7 @@ next_face:
//typename K::Translated_point_d transl =
// m_k.translated_point_d_object();
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
KNS_range kns_range = m_points_ds_for_tse.query_ANN(
p, num_points_for_pca, false);
const Points &points_for_pca = m_points_for_tse;
@ -2127,6 +2275,95 @@ next_face:
*/
}
// Compute the space tangent to a simplex (p1, p2, ... pn)
// CJTODO: Improve this?
// Basically, it takes all the neighbor points to p1, p2... pn and runs PCA
// on it. Note that most points are duplicated.
Tangent_space_basis compute_tangent_space(
const Indexed_simplex &s, bool normalize_basis = true)
{
unsigned int num_points_for_pca = static_cast<unsigned int>(
std::pow(BASE_VALUE_FOR_PCA, m_intrinsic_dim));
// Kernel functors
typename K::Construct_vector_d constr_vec =
m_k.construct_vector_d_object();
typename K::Compute_coordinate_d coord =
m_k.compute_coordinate_d_object();
typename K::Squared_length_d sqlen =
m_k.squared_length_d_object();
typename K::Scaled_vector_d scaled_vec =
m_k.scaled_vector_d_object();
typename K::Scalar_product_d scalar_pdct =
m_k.scalar_product_d_object();
typename K::Difference_of_vectors_d diff_vec =
m_k.difference_of_vectors_d_object();
//typename K::Translated_point_d transl =
// m_k.translated_point_d_object();
// One row = one point
Eigen::MatrixXd mat_points(s.size()*num_points_for_pca, m_ambient_dim);
unsigned int current_row = 0;
for (Indexed_simplex::const_iterator it_index = s.begin();
it_index != s.end() ; ++it_index)
{
const Point &p = m_points[*it_index];
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
KNS_range kns_range = m_points_ds_for_tse.query_ANN(
p, num_points_for_pca, false);
const Points &points_for_pca = m_points_for_tse;
#else
KNS_range kns_range = m_points_ds.query_ANN(p, num_points_for_pca, false);
const Points &points_for_pca = m_points;
#endif
KNS_iterator nn_it = kns_range.begin();
for ( ;
current_row < num_points_for_pca && nn_it != kns_range.end() ;
++current_row, ++nn_it)
{
for (int i = 0 ; i < m_ambient_dim ; ++i)
{
//const Point p = transl(
// points_for_pca[nn_it->first], m_translations[nn_it->first]);
mat_points(current_row, i) =
CGAL::to_double(coord(points_for_pca[nn_it->first], i));
}
}
}
Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
Eigen::MatrixXd cov = centered.adjoint() * centered;
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
Tangent_space_basis tsb;
// The eigenvectors are sorted in increasing order of their corresponding
// eigenvalues
for (int j = m_ambient_dim - 1 ;
j >= m_ambient_dim - m_intrinsic_dim ;
--j)
{
if (normalize_basis)
{
Vector v = constr_vec(m_ambient_dim,
eig.eigenvectors().col(j).data(),
eig.eigenvectors().col(j).data() + m_ambient_dim);
tsb.push_back(normalize_vector(v, m_k));
}
else
{
tsb.push_back(constr_vec(
m_ambient_dim,
eig.eigenvectors().col(j).data(),
eig.eigenvectors().col(j).data() + m_ambient_dim));
}
}
return tsb;
}
// Returns the dimension of the ith local triangulation
// This is particularly useful for the alpha-TC
int tangent_basis_dim(std::size_t i) const
@ -2830,6 +3067,95 @@ next_face:
return is_inconsistent;
}
bool filter_inconsistent_simplices_in_a_local_triangulation(
std::size_t tr_index)
{
bool is_inconsistent = false;
Star &star = m_stars[tr_index];
Triangulation const& tr = m_triangulations[tr_index].tr();
Tr_vertex_handle center_vh = m_triangulations[tr_index].center_vertex();
const Tr_traits &local_tr_traits = tr.geom_traits();
int cur_dim = tr.current_dimension();
// For each incident simplex
Star::iterator it_inc_simplex = star.begin();
while (it_inc_simplex != star.end())
{
const Incident_simplex &incident_simplex = *it_inc_simplex;
// Don't check infinite cells
if (is_infinite(incident_simplex))
{
++it_inc_simplex;
continue;
}
Indexed_simplex c = incident_simplex;
c.insert(tr_index); // Add the missing index
// Inconsistent?
if (!is_simplex_consistent(c))
{
is_inconsistent = true;
// Compute T, the space tangent to "c"
Tangent_space_basis tsb = compute_tangent_space(c);
std::size_t center_index = *(c.begin());
// Compute the star of "center_index" in the RT restricted to T
tsb.set_origin(center_index);
Triangulation tr(m_intrinsic_dim);
Tr_vertex_handle center_vh = compute_star(
center_index, compute_perturbed_point(center_index), tsb, tr);
// Get incident cells
std::vector<Tr_full_cell_handle> incident_cells;
tr.incident_full_cells(
center_vh, std::back_inserter(incident_cells));
// Check if "c" is among them
typename std::vector<Tr_full_cell_handle>::const_iterator it_c = incident_cells.begin();
typename std::vector<Tr_full_cell_handle>::const_iterator it_c_end = incident_cells.end();
bool found = false;
// For each cell
for (; !found && it_c != it_c_end ; ++it_c)
{
Incident_simplex incident_simplex;
for (int j = 0 ; j < tr.current_dimension() + 1 ; ++j)
{
std::size_t index = (*it_c)->vertex(j)->data();
incident_simplex.insert(index);
}
if (incident_simplex == c)
found = true;
}
// If the simplex is not in the new star, we remove it
if (!found)
{
if (it_inc_simplex != star.end() - 1)
{
*it_inc_simplex = star.back();
star.pop_back();
}
else
{
star.pop_back();
it_inc_simplex = star.end(); // end the loop
}
}
else
++it_inc_simplex;
}
// Consistent
else
++it_inc_simplex;
}
return is_inconsistent;
}
std::ostream &export_vertices_to_off(
std::ostream & os, std::size_t &num_vertices,
bool use_perturbed_points = false) const
@ -4105,7 +4431,7 @@ private:
//std::vector<Tr_mutex> m_tr_mutexes;
#endif
#ifdef USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
#ifdef CGAL_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
Points m_points_for_tse;
Points_ds m_points_ds_for_tse;
#endif