diff --git a/Tangential_complex/dont_submit b/Tangential_complex/dont_submit new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Tangential_complex/include/CGAL/Tangential_complex.h b/Tangential_complex/include/CGAL/Tangential_complex.h new file mode 100644 index 00000000000..146f6990f17 --- /dev/null +++ b/Tangential_complex/include/CGAL/Tangential_complex.h @@ -0,0 +1,404 @@ +// Copyright (c) 2014 INRIA Sophia-Antipolis (France) +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL: $ +// $Id: $ +// +// +// Author(s) : Clement Jamin + + +#ifndef TANGENTIAL_COMPLEX_H +#define TANGENTIAL_COMPLEX_H + +#include +#include + +#include +#include +#include + +#include // CJTODO TEMP + +#include +#include +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +# include +#endif + +namespace CGAL { + +/// The class Tangential_complex represents a tangential complex +template < + typename Kernel, + int Intrinsic_dimension, + typename Concurrency_tag = CGAL::Parallel_tag, + typename Tr = Regular_triangulation + < + Regular_triangulation_euclidean_traits< + Epick_d > >, + + Triangulation_data_structure + < + typename Regular_triangulation_euclidean_traits< + Epick_d > >::Dimension, + Triangulation_vertex > >, std::size_t >, + Triangulation_full_cell > > > + > + > +> +class Tangential_complex +{ + typedef typename Kernel::FT FT; + typedef typename Kernel::Point_d Point; + typedef typename Kernel::Vector_d Vector; + + typedef Tr Triangulation; + typedef typename Triangulation::Geom_traits Tr_traits; + typedef typename Triangulation::Point Tr_point; + typedef typename Triangulation::Bare_point Tr_bare_point; + typedef typename Triangulation::Vertex_handle Tr_vertex_handle; + typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle; + + typedef typename std::vector Tangent_space_basis; + + typedef std::pair Tr_and_VH; + typedef typename std::vector Point_container; + typedef typename std::vector Tr_container; + typedef typename std::vector TS_container; + + // Stores the index of the original Point in the ambient space + /*struct Tr_point_with_index + : public Tr_point + { + Tr_point_with_index(const Tr_point &p, std::size_t i) + : Tr_point(p), index(i) {} + + std::size_t index; + };*/ + +public: + /// Constructor + Tangential_complex(const Kernel &k = Kernel()) + : m_k(k){} + + /// Constructor for a range of points + template + Tangential_complex(InputIterator first, InputIterator last, + const Kernel &k = Kernel()) + : m_k(k), m_points(first, last) {} + + /// Destructor + ~Tangential_complex() {} + + void compute_tangential_complex() + { + // We need to do that because we don't want the container to copy the + // already-computed triangulations (while resizing) since it would + // invalidate the vertex handles stored beside the triangulations + m_triangulations.resize( + m_points.size(), + std::make_pair((Triangulation*)NULL, Tr_vertex_handle())); + m_tangent_spaces.resize(m_points.size()); + +#ifdef CGAL_LINKED_WITH_TBB + // Parallel + if (boost::is_convertible::value) + { + // Apply moves in triangulation + tbb::parallel_for(tbb::blocked_range(0, m_points.size()), + Compute_tangent_triangulation(*this) + ); + } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB + { + for (std::size_t i = 0 ; i < m_points.size() ; ++i) + compute_tangent_triangulation(i); + } + } + + std::ostream &export_to_off(std::ostream & os) + { + const int ambient_dim = Ambient_dimension::value; + if (ambient_dim < 2 || ambient_dim > 3) + { + std::cerr << "Error: export_to_off => ambient dimension should be 2 or 3."; + os << "Error: export_to_off => ambient dimension should be 2 or 3."; + return os; + } + + if (Intrinsic_dimension < 1 || Intrinsic_dimension > 3) + { + std::cerr << "Error: export_to_off => intrinsic dimension should be between 1 and 3."; + os << "Error: export_to_off => intrinsic dimension should be between 1 and 3."; + return os; + } + + std::stringstream output; + + //******** VERTICES ************ + + Point_container::const_iterator it_p = m_points.begin(); + Point_container::const_iterator it_p_end = m_points.end(); + // For each point p + for ( ; it_p != it_p_end ; ++it_p) + { + int i = 0; + for ( ; i < ambient_dim ; ++i) + output << (*it_p)[i] << " "; + if (i == 2) + output << "0"; + output << std::endl; + } + + //******** CELLS ************ + + std::size_t num_cells = 0; + Tr_container::const_iterator it_tr = m_triangulations.begin(); + Tr_container::const_iterator it_tr_end = m_triangulations.end(); + // For each triangulation + for ( ; it_tr != it_tr_end ; ++it_tr) + { + const Triangulation &tr = *it_tr->first; + Tr_vertex_handle center_vh = it_tr->second; + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + std::vector::const_iterator it_c = incident_cells.begin(); + std::vector::const_iterator it_c_end= incident_cells.end(); + // For each triangulation + for ( ; it_c != it_c_end ; ++it_c) + { + output << Intrinsic_dimension + 1 << " "; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + output << (*it_c)->vertex(i)->data() << " "; + output << std::endl; + ++num_cells; + } + } + + os << "OFF \n" + << m_points.size() << " " + << num_cells << " " + << "0 \n" + << output.str(); + + return os; + } + +private: + +#ifdef CGAL_LINKED_WITH_TBB + // Functor for compute_tangential_complex function + class Compute_tangent_triangulation + { + Tangential_complex & m_tc; + + public: + // Constructor + Compute_tangent_triangulation(Tangential_complex &tc) + : m_tc(tc) + {} + + // Constructor + Compute_tangent_triangulation(const Compute_tangent_triangulation &ctt) + : m_tc(ctt.m_tc) + {} + + // operator() + void operator()( const tbb::blocked_range& r ) const + { + for( size_t i = r.begin() ; i != r.end() ; ++i) + m_tc.compute_tangent_triangulation(i); + } + }; +#endif // CGAL_LINKED_WITH_TBB + + void compute_tangent_triangulation(std::size_t i) + { + Triangulation *p_local_tr = + m_triangulations[i].first = + new Triangulation(Intrinsic_dimension); + const Tr_traits &local_tr_traits = p_local_tr->geom_traits(); + Tr_vertex_handle ¢er_vertex = m_triangulations[i].second; + + // Estimate the tangent space + const Point ¢er_pt = m_points[i]; + m_tangent_spaces[i] = compute_tangent_space(center_pt); + + //*************************************************** + // Build a minimal triangulation in the tangent space + // (we only need the star of p) + //*************************************************** + + // First, compute the projected points + std::vector projected_points; + FT max_squared_weight = 0; + projected_points.reserve(m_points.size() - 1); + Point_container::const_iterator it_p = m_points.begin(); + Point_container::const_iterator it_p_end = m_points.end(); + for (std::size_t j = 0 ; it_p != it_p_end ; ++it_p, ++j) + { + // ith point = p, which is already inserted + if (j != i) + { + Tr_point wp = project_point(*it_p, center_pt, m_tangent_spaces[i]); + projected_points.push_back(wp); + FT w = local_tr_traits.point_weight_d_object()(wp); + if (w > max_squared_weight) + max_squared_weight = w; + } + } + + // Now we can insert the points + + // Insert p + Tr_point wp = local_tr_traits.construct_weighted_point_d_object()( + local_tr_traits.construct_point_d_object()(0, 0), + CGAL::sqrt(max_squared_weight)); + center_vertex = p_local_tr->insert(wp); + center_vertex->data() = i; + //std::cerr << "Inserted CENTER POINT of weight " << CGAL::sqrt(max_squared_weight) << std::endl; + + /*std::cerr << 0 << " " + << 0 << " " + << CGAL::sqrt(max_squared_weight) << std::endl;*/ + + // Insert the other points + std::vector::const_iterator it_wp = projected_points.begin(); + it_p = m_points.begin(); + for (std::size_t j = 0 ; it_p != it_p_end ; ++it_p, ++j) + { + // ith point = p, which is already inserted + if (j != i) + { + FT squared_dist_to_tangent_plane = + local_tr_traits.point_weight_d_object()(*it_wp); + FT w = CGAL::sqrt(max_squared_weight - squared_dist_to_tangent_plane); + Tr_point wp = local_tr_traits.construct_weighted_point_d_object()( + local_tr_traits.point_drop_weight_d_object()(*it_wp), + w); + /*Tr_bare_point bp = traits.point_drop_weight_d_object()(*it_wp); + Tr_point wp(traits.point_drop_weight_d_object()(*it_wp), w);*/ + + Tr_vertex_handle vh = p_local_tr->insert_if_in_star(wp, center_vertex); + //Tr_vertex_handle vh = p_local_tr->insert(wp); + if (vh != Tr_vertex_handle()) + { + /*std::cerr << traits.point_drop_weight_d_object()(*it_wp)[0] << " " + << traits.point_drop_weight_d_object()(*it_wp)[1] << " " + << w << std::endl;*/ + vh->data() = j; + } + ++it_wp; + } + } + + // CJTODO DEBUG + //std::cerr << "\nChecking topology and geometry..." + // << (p_local_tr->is_valid(true) ? "OK.\n" : "Error.\n"); + // DEBUG: output the local mesh into an OFF file + //std::stringstream sstr; + //sstr << "data/local_tri_" << i << ".off"; + //std::ofstream off_stream_tr(sstr.str()); + //CGAL::export_triangulation_to_off(off_stream_tr, *p_local_tr); + } + + Tangent_space_basis compute_tangent_space(const Point &p) const + { + Tangent_space_basis ts; + ts.reserve(Intrinsic_dimension); + // CJTODO: this is only for a sphere in R^3 + Vector t1(-p[1] - p[2], p[0], p[0]); + Vector t2(p[1] * t1[2] - p[2] * t1[1], + p[2] * t1[0] - p[0] * t1[2], + p[0] * t1[1] - p[1] * t1[0]); + // CJTODO: this is for a plane (test) + //Vector t1(1, 0, 0); + //Vector t2(0, 1, 0); + + // Normalize t1 and t2 + Get_functor::type sqlen(m_k); + + // CJTODO: use this when Scaled_vector is fixed + //Get_functor::type scale(m_k); + //ts.push_back(scale(t1, 1./CGAL::sqrt(sqlen(t1)))); + //ts.push_back(scale(t2, 1./CGAL::sqrt(sqlen(t2)))); + + // ****** Temporary code ******* + FT t1_len = CGAL::sqrt(sqlen(t1)); + FT t2_len = CGAL::sqrt(sqlen(t2)); + for (int i = 0 ; i < Ambient_dimension::value ; ++i) + { + t1[i] /= t1_len; + t2[i] /= t2_len; + } + ts.push_back(t1); + ts.push_back(t2); + // ****** /Temporary code ******* + + return ts; + } + + // Project the point in the tangent space + // The weight will be the squared distance between p and the projection of p + Tr_point project_point(const Point &p, const Point &origin, + const Tangent_space_basis &ts) const + { + Get_functor::type inner_pdct(m_k); + Get_functor::type diff_points(m_k); + + std::vector coords; + // Ambiant-space coords of the projected point + std::vector p_proj(origin.cartesian_begin(), origin.cartesian_end()); + coords.reserve(Intrinsic_dimension); + for (std::size_t i = 0 ; i < Intrinsic_dimension ; ++i) + { + // Compute the inner product p * ts[i] + Vector v = diff_points(p, origin); + FT coord = inner_pdct(v, ts[i]); + coords.push_back(coord); + + // p_proj += coord * v; + for (int j = 0 ; j < Ambient_dimension::value ; ++j) + p_proj[i] += coord * ts[i][j]; + } + + Point projected_pt(Ambient_dimension::value, + p_proj.begin(), p_proj.end()); + return Tr_point( + Tr_bare_point(Intrinsic_dimension, coords.begin(), coords.end()), + m_k.squared_distance_d_object()(p, projected_pt)); + } + +private: + const Kernel m_k; + Point_container m_points; + TS_container m_tangent_spaces; + Tr_container m_triangulations; // Contains the triangulations + // and their center vertex + +}; // /class Tangential_complex + +} // end namespace CGAL + +#endif // TANGENTIAL_COMPLEX_H diff --git a/Tangential_complex/package_info/Tangential_complex/copyright b/Tangential_complex/package_info/Tangential_complex/copyright new file mode 100644 index 00000000000..8932b3233d2 --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/copyright @@ -0,0 +1,2 @@ +INRIA Sophia-Antipolis (France) + diff --git a/Tangential_complex/package_info/Tangential_complex/copyright.txt b/Tangential_complex/package_info/Tangential_complex/copyright.txt new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Tangential_complex/package_info/Tangential_complex/description.txt b/Tangential_complex/package_info/Tangential_complex/description.txt new file mode 100644 index 00000000000..b8969f7fc3e --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/description.txt @@ -0,0 +1,3 @@ +Tangential complex + +This CGAL component provides an implementation of the tangential complex. diff --git a/Tangential_complex/package_info/Tangential_complex/license.txt b/Tangential_complex/package_info/Tangential_complex/license.txt new file mode 100644 index 00000000000..8bb8efcb72b --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Tangential_complex/package_info/Tangential_complex/maintainer b/Tangential_complex/package_info/Tangential_complex/maintainer new file mode 100644 index 00000000000..769c1668e20 --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/maintainer @@ -0,0 +1 @@ +ClĂ©ment Jamin \ No newline at end of file diff --git a/Tangential_complex/test/Tangential_complex/CMakeLists.txt b/Tangential_complex/test/Tangential_complex/CMakeLists.txt new file mode 100644 index 00000000000..8d58f3631fe --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/CMakeLists.txt @@ -0,0 +1,46 @@ +# Created by the script cgal_create_cmake_script +# This is the CMake script for compiling a CGAL application. + + +project( Tangential_complex_test ) + +cmake_minimum_required(VERSION 2.6.2) +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6) + if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3) + cmake_policy(VERSION 2.8.4) + else() + cmake_policy(VERSION 2.6) + endif() +endif() + +find_package(CGAL QUIET COMPONENTS Core ) + +if ( CGAL_FOUND ) + include( ${CGAL_USE_FILE} ) + + find_package( TBB QUIET ) + + if( TBB_FOUND ) + include(${TBB_USE_FILE}) + list(APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES}) + endif() + + + include( CGAL_CreateSingleSourceCGALProgram ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + include_directories (BEFORE "include") + + create_single_source_cgal_program( "test_tangential_complex.cpp" ) + + else() + message(STATUS "NOTICE: Some of the executables in this directory need Eigen 3.1 (or greater) and will not be compiled.") + endif() + +else() + message(STATUS "This program requires the CGAL library, and will not be compiled.") +endif() + diff --git a/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp b/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp new file mode 100644 index 00000000000..2fd7fe8a1c1 --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp @@ -0,0 +1,55 @@ +// Without TBB_USE_THREADING_TOOL Intel Inspector XE will report false positives in Intel TBB +// (http://software.intel.com/en-us/articles/compiler-settings-for-threading-error-analysis-in-intel-inspector-xe/) +#ifdef _DEBUG +# define TBB_USE_THREADING_TOOL +#endif + +#include +#include +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +# include +#endif + +int main() +{ + typedef CGAL::Epick_d > Kernel; + typedef Kernel::Point_d Point; + const int INTRINSIC_DIMENSION = 2; + +#ifdef CGAL_LINKED_WITH_TBB + tbb::task_scheduler_init init(10); +#endif + +#ifdef _DEBUG + const int NUM_POINTS = 50; +#else + const int NUM_POINTS = 5000; +#endif + CGAL::default_random = CGAL::Random(0); // NO RANDOM + CGAL::Random_points_on_sphere_3 generator(3.0); + std::vector points; + points.reserve(NUM_POINTS); + for (int i = 0 ; i != NUM_POINTS ; ++i) + { + points.push_back(*generator++); + //points.push_back(Point((double)(rand()%10000)/5000, (double)(rand()%10000)/5000, 0.)); // CJTODO : plane + } + + CGAL::Tangential_complex< + Kernel, + INTRINSIC_DIMENSION, + CGAL::Parallel_tag> tc(points.begin(), points.end()); + + tc.compute_tangential_complex(); + + std::stringstream output_filename; + output_filename << "data/test_tc_" << INTRINSIC_DIMENSION << ".off"; + std::ofstream off_stream(output_filename.str()); + tc.export_to_off(off_stream); + + return 0; +} \ No newline at end of file diff --git a/Triangulation/include/CGAL/Regular_triangulation.h b/Triangulation/include/CGAL/Regular_triangulation.h index d0552b6d7eb..c5e127532c8 100644 --- a/Triangulation/include/CGAL/Regular_triangulation.h +++ b/Triangulation/include/CGAL/Regular_triangulation.h @@ -1029,7 +1029,7 @@ Regular_triangulation } else { - Orientation_d ori = geom_traits().orientation_d_object(); + Orientation_d ori = geom_traits().orientation_d_object(); // CJTODO: create member variables for this? Power_test_d side = geom_traits().power_test_d_object(); Conflict_pred_in_fullspace c(*this, p, ori, side); return c(s);