diff options
author | salinasd <salinasd@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-02-10 17:03:43 +0000 |
---|---|---|
committer | salinasd <salinasd@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-02-10 17:03:43 +0000 |
commit | def05e2da637a43c02e439af8faaf789214395cd (patch) | |
tree | b2839ae762693a2a09c4876c976772332fc1342a | |
parent | e50f918673baf48d35c2e2ffa7bfc0e23206612f (diff) |
skbl correction bug constructor + contraction add garland heckbert example
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@466 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 1789e903625df3c4e7c689fa4888bebd86e616eb
18 files changed, 2371 insertions, 1509 deletions
diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt index ba411bd4..6d65de73 100644 --- a/src/Contraction/example/CMakeLists.txt +++ b/src/Contraction/example/CMakeLists.txt @@ -1,7 +1,10 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIskbl) + add_executable(RipsContraction Rips_contraction.cpp) +add_executable(GarlandHeckbert Garland_heckbert.cpp) + target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) add_test(RipsContraction.sphere.0.2 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/sphere3D_2646.off 0.2) diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp new file mode 100644 index 00000000..912199a4 --- /dev/null +++ b/src/Contraction/example/Garland_heckbert.cpp @@ -0,0 +1,115 @@ +/* + * Garland_heckbert.h + * Created on: Feb 10, 2015 + * This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + */ + + +#ifndef GARLAND_HECKBERT_H_ +#define GARLAND_HECKBERT_H_ + + +#include <iostream> +#include "Garland_heckbert/Point.h" +#include "Garland_heckbert/Error_quadric.h" +#include "Garland_heckbert/Garland_heckbert_policies.h" + +#include "gudhi/Edge_contraction.h" +#include "gudhi/Skeleton_blocker.h" +#include "gudhi/Off_reader.h" + +using namespace std; +using namespace Gudhi; +using namespace skbl; +using namespace contraction; + + +struct Geometry_trait{ + typedef Point_d Point; +}; + +struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits<Geometry_trait> { + +public: + // the vertex stored in the complex contains a quadric + struct Garland_heckbert_vertex : public Simple_geometric_vertex{ + Error_quadric<Geometry_trait::Point> quadric; + }; + typedef Garland_heckbert_vertex Graph_vertex; + +}; + +typedef Skeleton_blocker_geometric_complex< Garland_heckbert_traits > Complex; +typedef Edge_profile<Complex> Profile; +typedef Skeleton_blocker_contractor<Complex> Complex_contractor; + + + + +int main(int argc, char *argv[]){ + if (argc!=3){ + std::cerr << "Usage "<<argv[0]<<" ../../../data/meshes/test.off N to load the file ../../data/test.off and contract N edges.\n"; + return -1; + } + + Complex complex; + + // load the points + Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex); + if(!off_reader.is_valid()){ + std::cerr << "Unable to read file:"<<argv[1]<<std::endl; + return EXIT_FAILURE; + } + + int num_contractions = atoi(argv[2]); + + std::cout << "Initial complex has "<< + complex.num_vertices()<<" vertices, "<< + complex.num_blockers()<<" blockers, "<< + complex.num_edges()<<" edges and" << + complex.num_triangles()<<" triangles."; + + Complex_contractor contractor(complex, + new GH_cost<Profile>(complex), + new GH_placement<Profile>(complex), + contraction::make_link_valid_contraction<Profile>(), + new GH_visitor<Profile>(complex) + ); + + std::cout<<"Contract "<<num_contractions<<" edges\n"; + contractor.contract_edges(num_contractions); + + std::cout << "Final complex has "<< + complex.num_vertices()<<" vertices, "<< + complex.num_edges()<<" edges and" << + complex.num_triangles()<<" triangles."; + + Skeleton_blocker_off_writer<Complex> off_writer("reduced.off",complex); + + + return EXIT_SUCCESS; +} + + + +#endif /* GARLAND_HECKBERT_H_ */ diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h new file mode 100755 index 00000000..7a30c35a --- /dev/null +++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h @@ -0,0 +1,164 @@ +/*
+ * Error_quadric.h
+ *
+ * Created on: 24 janv. 2014
+ * Author: dsalinas
+ */
+
+#ifndef ERROR_QUADRIC_H_
+#define ERROR_QUADRIC_H_
+
+#include <vector>
+#include <utility>
+#include <boost/optional/optional.hpp>
+
+
+template <typename Point> class Error_quadric{
+private :
+ double coeff[10];
+
+public :
+ Error_quadric(){
+ clear();
+ }
+
+ /**
+ * Quadric corresponding to the L2 distance to the plane.
+ *
+ * According to the notation of Garland Heckbert, they
+ * denote a quadric symetric matrix as :
+ * Q = [ q11 q12 q13 q14]
+ * [ q12 q22 q23 q24]
+ * [ q13 q23 q33 q34]
+ * [ q14 q24 q34 q44]
+ *
+ * which is represented by a vector with 10 elts that
+ * are denoted ci for clarity with :
+ * Q = [ c0 c1 c2 c3 ]
+ * [ c1 c4 c5 c6 ]
+ * [ c2 c5 c7 c8 ]
+ * [ c3 c6 c8 c9 ]
+ *
+ * The constructor return the quadrics that represents
+ * the squared distance to the plane defined by triangle p0,p1,p2
+ * times the area of triangle p0,p1,p2.
+ */
+ Error_quadric(const Point & p0,const Point & p1,const Point & p2){
+
+ Point normal(unit_normal(p0,p1,p2));
+ double a=normal[0];
+ double b=normal[1];
+ double c=normal[2];
+ double d= -a*p0[0]-b*p0[1]-c*p0[2];
+ coeff[0] = a*a ;
+ coeff[1] = a*b ;
+ coeff[2] = a*c ;
+ coeff[3] = a*d ;
+ coeff[4] = b*b ;
+ coeff[5] = b*c ;
+ coeff[6] = b*d ;
+ coeff[7] = c*c ;
+ coeff[8] = c*d ;
+ coeff[9] = d*d ;
+
+ double area_p0p1p2 = std::sqrt(squared_area(p0,p1,p2));
+ for(auto& x : coeff)
+ x*= area_p0p1p2;
+ }
+
+
+ inline double squared_area(const Point& p0,const Point& p1,const Point& p2) {
+ //if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0
+ //then the squared area is = (u^2+v^2+w^2)/4
+ //with: u = x2 * y3 - x3 * y2;
+ // v = x3 * y1 - x1 * y3;
+ // w = x1 * y2 - x2 * y1;
+ Point p0p1(p1-p0);
+ Point p0p2(p2-p0);
+ double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1];
+ double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2];
+ double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0];
+ return 1./4. * (A*A+B*B+C*C);
+ }
+
+
+ void clear(){
+ for(auto& x:coeff)
+ x=0;
+ }
+
+ Error_quadric& operator+=(const Error_quadric& other){
+ if(this!=&other)
+ for(int i = 0 ; i < 10; ++i)
+ coeff[i] += other.coeff[i];
+ return *this;
+ }
+
+ /**
+ * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert
+ */
+ inline double cost(const Point& point) const{
+ double cost =
+ coeff[0]*point.x()*point.x()+coeff[4]*point.y()*point.y()+coeff[7]*point.z()*point.z()
+ +2*(coeff[1]*point.x()*point.y()+coeff[5]*point.y()*point.z()+coeff[2]*point.z()*point.x())
+ +2*(coeff[3]*point.x()+coeff[6]*point.y()+coeff[8]*point.z())
+ +coeff[9];
+ if(cost<0) return 0;
+ else {
+ return cost;
+ }
+ }
+
+ inline double grad_determinant() const{
+ return
+ coeff[0] * coeff[4] * coeff[7]
+ - coeff[0] * coeff[5] * coeff[5]
+ - coeff[1] * coeff[1] * coeff[7]
+ +2*coeff[1] * coeff[5] * coeff[2]
+ - coeff[4] * coeff[2] * coeff[2];
+ }
+
+ /**
+ * Return the point such that it minimizes the gradient of the quadric.
+ * Det must be passed with the determinant value of the gradient (should be non zero).
+ */
+ inline Point solve_linear_gradient(double det = grad_determinant()) const{
+ return Point(
+ (-coeff[1]*coeff[5]*coeff[8]+coeff[1]*coeff[7]*coeff[6]+coeff[2]*coeff[8]*coeff[4]-coeff[2]*coeff[5]*coeff[6]-coeff[3]*coeff[4]*coeff[7]+coeff[3]*coeff[5]*coeff[5])/ det,
+ (coeff[0]*coeff[5]*coeff[8]-coeff[0]*coeff[7]*coeff[6]-coeff[5]*coeff[2]*coeff[3]-coeff[1]*coeff[2]*coeff[8]+coeff[6]*coeff[2]*coeff[2]+coeff[1]*coeff[3]*coeff[7])/det,
+ (-coeff[8]*coeff[0]*coeff[4]+coeff[8]*coeff[1]*coeff[1]+coeff[2]*coeff[3]*coeff[4]+coeff[5]*coeff[0]*coeff[6]-coeff[5]*coeff[1]*coeff[3]-coeff[1]*coeff[2]*coeff[6])/det
+ );
+ }
+
+
+ /**
+ * returns the point that minimizes the quadric.
+ * It inverses the quadric if its determinant is higher that a given threshold .
+ * If the determinant is lower than this value the returned value is uninitialized.
+ */
+ boost::optional<Point> min_cost(double scale=1) const{
+ // const double min_determinant = 1e-4 * scale*scale;
+ const double min_determinant = 1e-5;
+ boost::optional<Point> pt_res;
+ double det = grad_determinant();
+ if (std::abs(det)>min_determinant)
+ pt_res = solve_linear_gradient(det);
+ return pt_res;
+ }
+
+ friend std::ostream& operator<< (std::ostream& stream, const Error_quadric& quadric) {
+ stream << "\n[ "<<quadric.coeff[0]<<","<<quadric.coeff[1]<<","<<quadric.coeff[2]<<","<<quadric.coeff[3]<<";\n";
+ stream << " "<<quadric.coeff[1]<<","<<quadric.coeff[4]<<","<<quadric.coeff[5]<<","<<quadric.coeff[6]<<";\n";
+ stream << " "<<quadric.coeff[2]<<","<<quadric.coeff[5]<<","<<quadric.coeff[7]<<","<<quadric.coeff[8]<<";\n";
+ stream << " "<<quadric.coeff[3]<<","<<quadric.coeff[6]<<","<<quadric.coeff[8]<<","<<quadric.coeff[9]<<"]";
+ return stream;
+ }
+
+
+};
+
+
+
+
+#endif /* ERROR_QUADRIC_H_ */
+
diff --git a/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h b/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h new file mode 100644 index 00000000..c792b5b3 --- /dev/null +++ b/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h @@ -0,0 +1,115 @@ +/* + * Garland_heckbert_policies.h + * Created on: Feb 10, 2015 + * This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + */ + + +#ifndef GARLAND_HECKBERT_POLICIES_H_ +#define GARLAND_HECKBERT_POLICIES_H_ + +#include "gudhi/Edge_contraction.h" +#include "Error_quadric.h" + +template<typename EdgeProfile> +class GH_visitor: public Gudhi::contraction::Contraction_visitor<EdgeProfile> { + typedef typename EdgeProfile::Complex Complex; + typedef typename Complex::Point Point; + Complex& complex_; +public: + GH_visitor(Complex& complex):complex_(complex){} + + + void on_started(Complex & complex) override{ + //Compute quadrics for every vertex v + //The quadric of v consists in the sum of quadric + //of every triangles passing through v weighted by its area + for(auto v : complex.vertex_range()){ + auto & quadric_v(complex[v].quadric); + for(auto t : complex.triangle_range(v)){ + auto t_it = t.begin(); + const auto& p0(complex.point(*t_it++)); + const auto& p1(complex.point(*t_it++)); + const auto& p2(complex.point(*t_it++)); + quadric_v+=Error_quadric<Point>(p0,p1,p2); + } + } + } + + /** + * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. + */ + void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement) + override{ + profile.v0().quadric += profile.v1().quadric; + } +}; + + +template<typename EdgeProfile> +class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile>{ + typedef typename EdgeProfile::Complex Complex; + Complex& complex_; +public: + typedef typename EdgeProfile::Point Point; + typedef typename Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; + + GH_placement(Complex& complex):complex_(complex){} + + Placement_type operator()(const EdgeProfile& profile) const override{ + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + + boost::optional<Point> min_quadric_pt(sum_quad.min_cost()); + if (min_quadric_pt) + return Placement_type(*min_quadric_pt); + else + return profile.p0(); + } +}; + +template<typename EdgeProfile> +class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile>{ + typedef typename EdgeProfile::Complex Complex; + Complex& complex_; +public: + + typedef typename Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; + typedef typename Complex::Point Point; + + GH_cost(Complex& complex):complex_(complex){} + + Cost_type operator()( EdgeProfile const& profile, boost::optional<Point> const& new_point ) const override { + Cost_type res; + if (new_point){ + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + res = sum_quad.cost(*new_point); + } + return res; + } +}; + + + +#endif /* GARLAND_HECKBERT_POLICIES_H_ */ diff --git a/src/Contraction/example/Garland_heckbert/Point.h b/src/Contraction/example/Garland_heckbert/Point.h new file mode 100644 index 00000000..eb250347 --- /dev/null +++ b/src/Contraction/example/Garland_heckbert/Point.h @@ -0,0 +1,176 @@ +/* + * Basic_geometry.h + * Created on: Feb 10, 2015 + * This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + */ + + +#ifndef BASIC_GEOMETRY_H_ +#define BASIC_GEOMETRY_H_ + +#include <cmath> +#include <vector> +#include <cassert> +#include <cstddef> + +class Point_d{ +public: + Point_d():coords_(3,0){} + Point_d(size_t dim):coords_(dim,0){} + Point_d(const Point_d& other):coords_(other.coords_){} + Point_d(double x,double y,double z):coords_({x,y,z}){} + + template<typename CoordsIt> + Point_d(CoordsIt begin,CoordsIt end):coords_(begin,end){} + + size_t dimension() const{ + return coords_.size(); + } + + double x() const{ + return coords_[0]; + } + + double y() const{ + return coords_[1]; + } + + double z() const{ + return coords_[2]; + } + + double& x(){ + return coords_[0]; + } + + double& y(){ + return coords_[1]; + } + + double& z(){ + return coords_[2]; + } + + std::vector<double>::const_iterator begin() const{ + return coords_.begin(); + } + + std::vector<double>::const_iterator end() const{ + return coords_.end(); + } + + double& operator[](unsigned i){ + return coords_[i]; + } + const double& operator[](unsigned i) const{ + return coords_[i]; + } + + double squared_norm() const{ + double res = 0; + for(auto x : coords_) + res+= x*x; + return res; + } + + double squared_dist(const Point_d& other) const{ + assert(dimension()==other.dimension()); + double res = 0; + for(unsigned i = 0; i < coords_.size(); ++i) + res+= coords_[i]*coords_[i] + other[i]*other[i]; + return res; + } + + /** + * dot product + */ + double operator*(const Point_d& other) const{ + assert(dimension()==other.dimension()); + double res = 0; + for(unsigned i = 0; i < coords_.size(); ++i) + res+= coords_[i]*other[i]; + return res; + } + + /** + * only if points have dimension 3 + */ + Point_d cross_product(const Point_d& other){ + assert(dimension()==3 && other.dimension()==3); + Point_d res(3); + res[0] = (*this)[1] * other[2] - (*this)[2] * other[1]; + res[1] = (*this)[2] * other[0] - (*this)[0] * other[2]; + res[2] = (*this)[0] * other[1] - (*this)[1] * other[0]; + return res; + } + + Point_d operator+(const Point_d& other) const{ + assert(dimension()==other.dimension()); + Point_d res(dimension()); + for(unsigned i = 0; i < coords_.size(); ++i) + res[i] = (*this)[i] + other[i]; + return res; + } + + Point_d operator*(double lambda) const{ + Point_d res(dimension()); + for(unsigned i = 0; i < coords_.size(); ++i) + res[i] = (*this)[i] * lambda; + return res; + } + + Point_d operator/(double lambda) const{ + Point_d res(dimension()); + for(unsigned i = 0; i < coords_.size(); ++i) + res[i] = (*this)[i] / lambda; + return res; + } + + Point_d operator-(const Point_d& other) const{ + assert(dimension()==other.dimension()); + Point_d res(dimension()); + for(unsigned i = 0; i < coords_.size(); ++i) + res[i] = (*this)[i] - other[i]; + return res; + } + + friend Point_d unit_normal(const Point_d& p1,const Point_d& p2,const Point_d& p3){ + assert(p1.dimension()==3); + assert(p2.dimension()==3); + assert(p3.dimension()==3); + Point_d p1p2 = p2 - p1; + Point_d p1p3 = p3 - p1; + Point_d res(p1p2.cross_product(p1p3)); + return res / std::sqrt(res.squared_norm()); + } + + +private: + std::vector<double> coords_; +}; + + + + + +#endif /* BASIC_GEOMETRY_H_ */ diff --git a/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h index c7941bfd..ba2138a1 100644 --- a/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h +++ b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h @@ -391,9 +391,7 @@ public: } } else - { DBG("uncomputable cost"); - } } if(contraction_visitor_) contraction_visitor_->on_stop_condition_reached(); } diff --git a/src/Skeleton_blocker/example/CMakeLists.txt b/src/Skeleton_blocker/example/CMakeLists.txt index d9eab4e1..d1228526 100644 --- a/src/Skeleton_blocker/example/CMakeLists.txt +++ b/src/Skeleton_blocker/example/CMakeLists.txt @@ -3,6 +3,8 @@ project(GUDHIskbl) add_executable(SkeletonBlockerFromSimplices Skeleton_blocker_from_simplices.cpp) - add_executable(SkeletonBlockerIteration Skeleton_blocker_iteration.cpp) +add_executable(SkeletonBlockerLink Skeleton_blocker_link.cpp) + + target_link_libraries(SkeletonBlockerIteration ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) diff --git a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp index 221c7881..92fa17f3 100644 --- a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp +++ b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp @@ -41,7 +41,7 @@ typedef Complex::Simplex_handle Simplex; Complex build_complete_complex(int n){ - // build a full complex with 10 vertices and 2^n-1 simplices + // build a full complex with n vertices and 2^n-1 simplices Complex complex; for(int i=0;i<n;i++) complex.add_vertex(); @@ -57,7 +57,7 @@ int main (int argc, char *argv[]){ const int n = 15; - // build a full complex with 10 vertices and 2^n-1 simplices + // build a full complex with n vertices and 2^n-1 simplices Complex complex(build_complete_complex(n)); // this is just to illustrate iterators, to count number of vertices diff --git a/src/Skeleton_blocker/example/Skeleton_blocker_link.cpp b/src/Skeleton_blocker/example/Skeleton_blocker_link.cpp new file mode 100644 index 00000000..0987cc89 --- /dev/null +++ b/src/Skeleton_blocker/example/Skeleton_blocker_link.cpp @@ -0,0 +1,71 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + + +#include <stdio.h> +#include <stdlib.h> +#include <string> +#include <fstream> +#include <sstream> + + +#include "gudhi/Skeleton_blocker.h" + +using namespace std; +using namespace Gudhi; +using namespace skbl; + +typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Complex::Vertex_handle Vertex_handle; +typedef Complex::Root_vertex_handle Root_vertex_handle; +typedef Complex::Simplex_handle Simplex; + + +int main (int argc, char *argv[]){ + // build a full complex with 4 vertices and 2^4-1 simplices + // Initial vertices are (0,1,2,3,4) + Simplex tetrahedron(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + Complex complex; + complex.add_simplex(tetrahedron); + + cout<<"complex:"<<complex.to_string()<<endl; + + //build the link of vertex 1, eg a triangle {0,2,3} + auto link = complex.link(Vertex_handle(1)); + cout<<"link:"<<link.to_string()<<endl; + + //Internally link is a subcomplex of 'complex' and its vertices are stored in a vector. + //They can be accessed via Vertex_handle(x) where x is an index of the vector. + //In that example, link has three vertices and thus it contains only + // Vertex_handle(0),Vertex_handle(1) and Vertex_handle(2) are). + for(int i = 0; i<5; ++i) + cout << "link.contains_vertex(Vertex_handle("<<i<<")):"<<link.contains_vertex(Vertex_handle(i))<<endl; + cout<<endl; + + //To access to the initial vertices eg (0,1,2,3,4), Root_vertex_handle must be used. + //For instance, to test if the link contains the vertex that was labeled i: + for(int i = 0; i<5; ++i) + cout << "link.contains_vertex(Root_vertex_handle("<<i<<")):"<<link.contains_vertex(Root_vertex_handle(i))<<endl; + + return EXIT_SUCCESS; +} + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h index bc832d5d..aaa682c3 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h @@ -1,24 +1,24 @@ - /* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): David Salinas - * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) - * - * This program is free software: 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. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ #ifndef SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ #define SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ @@ -35,72 +35,173 @@ namespace skbl { *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex. */ template<typename Complex> -class Skeleton_blocker_off_visitor_reader { - Complex& complex_; - typedef typename Complex::Vertex_handle Vertex_handle; - typedef typename Complex::Point Point; +class Skeleton_blocker_off_flag_visitor_reader { + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Point Point; - const bool load_only_points_; + const bool load_only_points_; - public: - explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false): - complex_(complex), - load_only_points_(load_only_points) {} +public: + explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false): + complex_(complex), + load_only_points_(load_only_points) {} - void init(int dim, int num_vertices, int num_faces, int num_edges) { - // todo do an assert to check that this number are correctly read - // todo reserve size for vector points - } + void init(int dim, int num_vertices, int num_faces, int num_edges) { + // todo do an assert to check that this number are correctly read + // todo reserve size for vector points + } - void point(const std::vector<double>& point) { - complex_.add_vertex(point); - } + void point(const std::vector<double>& point) { + complex_.add_vertex(Point(point.begin(),point.end())); + } - void maximal_face(const std::vector<int>& face) { - if (!load_only_points_) { - for (size_t i = 0; i < face.size(); ++i) - for (size_t j = i+1; j < face.size(); ++j) { - complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j])); - } - } - } + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + for (size_t i = 0; i < face.size(); ++i) + for (size_t j = i+1; j < face.size(); ++j) { + complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j])); + } + } + } - void done() {} + void done() {} }; /** -*@brief Class that allows to load a Skeleton_blocker_complex from an off file. -*/ + *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex. + */ +template<typename Complex> +class Skeleton_blocker_off_visitor_reader { + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Simplex_handle Simplex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + std::vector<Point> points_; + std::vector<Simplex_handle> maximal_faces_; + +public: + explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false): + complex_(complex), + load_only_points_(load_only_points) {} + + + void init(int dim, int num_vertices, int num_faces, int num_edges){ + maximal_faces_.reserve(num_faces); + points_.reserve(num_vertices); + } + + + void point(const std::vector<double>& point) { + points_.emplace_back(point.begin(),point.end()); + } + + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + Simplex_handle s; + for(auto x : face) + s.add_vertex(Vertex_handle(x)); + maximal_faces_.emplace_back(s); + } + } + + void done() { + complex_ = make_complex_from_top_faces( + maximal_faces_.begin(), + maximal_faces_.end(), + points_.begin(), + points_.end() + ); + } +}; + + +/** + *@brief Class that allows to load a Skeleton_blocker_complex from an off file. + */ template<typename Complex> class Skeleton_blocker_off_reader { - public: - /** - * name_file : file to read - * read_complex : complex that will receive the file content - * read_only_points : specify true if only the points must be read - */ - Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false):valid_(false) { - std::ifstream stream(name_file); - if (stream.is_open()) { - Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points); - Off_reader off_reader(stream); - valid_ = off_reader.read(off_visitor); - } - } - - /** - * return true iff reading did not meet problems. - */ - bool is_valid() const { - return valid_; - } - - private: - bool valid_; +public: + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,bool is_flag = false):valid_(false) { + std::ifstream stream(name_file); + if (stream.is_open()) { + if(is_flag){ + Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + else{ + Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + } + } + + /** + * return true iff reading did not meet problems. + */ + bool is_valid() const { + return valid_; + } + +private: + bool valid_; +}; + + +template<typename Complex> +class Skeleton_blocker_off_writer{ +public: + /** + * name_file : file where the off will be written + * save_complex : complex that be outputted in the file + * for now only save triangles. + */ + Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex){ + typedef typename Complex::Vertex_handle Vertex_handle; + + std::ofstream stream(name_file); + if (stream.is_open()) { + stream<<"OFF\n"; + size_t num_triangles = std::distance(save_complex.triangle_range().begin(),save_complex.triangle_range().end()); + stream<<save_complex.num_vertices()<<" "<<num_triangles<< " 0 \n"; + + //in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance + //we compute a map from 0 2 5 7 to 0 1 2 3 + std::map<Vertex_handle,size_t> vertex_num; + size_t current_vertex=0; + + for(auto v : save_complex.vertex_range()){ + vertex_num[v]=current_vertex++; + const auto& pt(save_complex.point(v)); + for(auto it = pt.begin();it!=pt.end();++it) + stream<<*it<<" "; + stream<<std::endl; + } + + for(const auto & t : save_complex.triangle_range()){ + stream<<"3 "; + for(auto x : t) + stream<<vertex_num[x]<<" "; + stream<<std::endl; + } + stream.close(); + } else std::cerr<<"could not open file "<<name_file<<std::endl; + } + }; + } // namespace skbl diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h index 717b6052..434538fe 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h @@ -64,9 +64,6 @@ class Skeleton_blocker_simplex { //@{ // Skeleton_blocker_simplex():simplex_set() {} - /** - * Clear the simplex - */ void clear() { simplex_set.clear(); } @@ -121,7 +118,6 @@ class Skeleton_blocker_simplex { * Add the vertex v to the simplex: * Note that adding two times the same vertex is * the same that adding it once. - * \f$ (*this) \leftarrow (*this) \cup \{ v \} \f$ */ void add_vertex(T v) { simplex_set.insert(v); @@ -129,15 +125,13 @@ class Skeleton_blocker_simplex { /** * Remove the vertex v from the simplex: - * \f$ (*this) \leftarrow (*this) \setminus \{ v \} \f$ */ void remove_vertex(T v) { simplex_set.erase(v); } /** - * Intersects the simplex with the simplex a: - * \f$ (*this) \leftarrow (*this) \cap a \f$ + * Intersects the simplex with the simplex a. */ void intersection(const Skeleton_blocker_simplex & a) { std::vector<T> v; @@ -152,8 +146,7 @@ class Skeleton_blocker_simplex { } /** - * Substracts a from the simplex: - * \f$ (*this) \leftarrow (*this) \setminus a \f$ + * Substracts a from the simplex. */ void difference(const Skeleton_blocker_simplex & a) { std::vector<T> v; @@ -169,8 +162,7 @@ class Skeleton_blocker_simplex { } /** - * Add vertices of a to the simplex: - * \f$ (*this) \leftarrow (*this) \cup a \f$ + * Add vertices of a to the simplex. */ void union_vertices(const Skeleton_blocker_simplex & a) { std::vector<T> v; @@ -234,7 +226,7 @@ class Skeleton_blocker_simplex { return *(simplex_set.rbegin()); } /** - * @return true iff the simplex contains the simplex a, i.e. iff \f$ a \subset (*this) \f$. + * @return true iff the simplex contains the simplex a. */ bool contains(const Skeleton_blocker_simplex & a) const { return includes(simplex_set.cbegin(), simplex_set.cend(), @@ -278,7 +270,7 @@ class Skeleton_blocker_simplex { auto last2 = a.end(); while (first2 != last2) { - // we ignore vertices of b + // we ignore vertices x if (x == *first2) { ++first2; } else { @@ -293,7 +285,7 @@ class Skeleton_blocker_simplex { } /** - * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$. + * @return true iff the simplex contains the difference \f$ a \setminus \{ x,y \} \f$. */ bool contains_difference(const Skeleton_blocker_simplex& a, T x, T y) const { auto first1 = begin(); @@ -303,7 +295,7 @@ class Skeleton_blocker_simplex { auto last2 = a.end(); while (first2 != last2) { - // we ignore vertices of b + // we ignore vertices of x,y if (x == *first2 || y == *first2) { ++first2; } else { @@ -317,16 +309,10 @@ class Skeleton_blocker_simplex { return true; } - /** - * @return true iff the simplex contains the vertex v, i.e. iff \f$ v \in (*this) \f$. - */ bool contains(T v) const { return (simplex_set.find(v) != simplex_set.end()); } - /** - * @return \f$ (*this) \cap a = \emptyset \f$. - */ bool disjoint(const Skeleton_blocker_simplex& a) const { std::vector<T> v; v.reserve(std::min(simplex_set.size(), a.simplex_set.size())); @@ -354,9 +340,6 @@ class Skeleton_blocker_simplex { //@} - /** - * Display a simplex - */ friend std::ostream& operator <<(std::ostream& o, const Skeleton_blocker_simplex & sigma) { bool first = true; diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h index 7b7ac7a9..81fc28cb 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h @@ -63,977 +63,988 @@ namespace skbl { */ template<class SkeletonBlockerDS> class Skeleton_blocker_complex { - template<class ComplexType> friend class Complex_vertex_iterator; - template<class ComplexType> friend class Complex_neighbors_vertices_iterator; - template<class ComplexType> friend class Complex_edge_iterator; - template<class ComplexType> friend class Complex_edge_around_vertex_iterator; - - template<class ComplexType> friend class Skeleton_blocker_link_complex; - template<class ComplexType> friend class Skeleton_blocker_link_superior; - template<class ComplexType> friend class Skeleton_blocker_sub_complex; - - public: - /** - * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS - */ - typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; - - /** - * @brief The type of stored edge node, specified by the template SkeletonBlockerDS - */ - typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; - - typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; - - /** - * @brief The type of an handle to a vertex of the complex. - */ - typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; - typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; - - /** - * @brief A ordered set of integers that represents a simplex. - */ - typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; - typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; - - /** - * @brief Handle to a blocker of the complex. - */ - typedef Simplex_handle* Blocker_handle; - - typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; - typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; - - protected: - typedef typename boost::adjacency_list<boost::setS, // edges - boost::vecS, // vertices - boost::undirectedS, Graph_vertex, Graph_edge> Graph; - // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink - // (eg Link with greater vertices) - // that burdens simplex iteration / complex initialization via list of simplices. - // to avoid that, one should modify the graph by storing two lists of adjacency for every - // vertex, the one with superior and the one with lower vertices, that way there is - // no more extra cost for computation of SuperiorLink - typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; - typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; - - protected: - typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; - - public: - /** - * @brief Handle to an edge of the complex. - */ - typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; - - protected: - typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator; - - protected: - int num_vertices_; - int num_blockers_; - - typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; - // typedef Visitor* Visitor_ptr; - Visitor* visitor; - - /** - * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. - * - * This quantity is updated when adding/removing edge. - * - * This is useful because the operation - * list.size() is done in linear time. - */ - std::vector<boost_vertex_handle> degree_; - Graph skeleton; /** 1-skeleton of the simplicial complex. */ - - /** Each vertex can access to the blockers passing through it. */ - BlockerMap blocker_map_; - - public: - ///////////////////////////////////////////////////////////////////////////// - /** @name Constructors, Destructors - */ - //@{ - - /** - *@brief constructs a simplicial complex with a given number of vertices and a visitor. - */ - explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) - : visitor(visitor_) { - clear(); - for (int i = 0; i < num_vertices_; ++i) { - add_vertex(); - } - } - - private: - typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie; - - template<typename SimpleHandleOutputIterator> - /** - * return a vector of simplex trie for every vertex - */ - std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) { - std::vector<STrie*> cofaces(num_vertices(), 0); - for (auto i = 0u; i < num_vertices(); ++i) - cofaces[i] = new STrie(Vertex_handle(i)); - for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { - if (s_it->dimension() >= 1) - cofaces[s_it->first_vertex()]->add_simplex(*s_it); - } - return cofaces; - } - - - public: - /** - * @brief Constructor with a list of simplices. - * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). - */ - template<typename SimpleHandleOutputIterator> - Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL) - : num_vertices_(0), - num_blockers_(0), - visitor(visitor_) { - std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; - - // first pass to add vertices and edges - int num_vertex = -1; - for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { - if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex); - if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex()); - } - while (num_vertex-- >= 0) add_vertex(); - - for (const auto& e : edges) - add_edge(e.first, e.second); - - if (!is_flag_complex) { - // need to compute blockers - std::vector<STrie*> cofaces(compute_cofaces(simplex_begin, simplex_end)); - std::vector<Simplex_handle> max_faces; - - for (auto i = 0u; i < num_vertices(); ++i) { - auto max_faces_i = cofaces[i]->maximal_faces(); - max_faces.insert(max_faces.end(), max_faces_i.begin(), max_faces_i.end()); - } - - // all time here - - // for every maximal face s, it picks its first vertex v and check for all nv \in N(v) - // if s union nv is in the complex, if not it is a blocker. - for (auto max_face : max_faces) { - if (max_face.dimension() > 0) { - auto first_v = max_face.first_vertex(); - for (auto nv : vertex_range(first_v)) { - if (!max_face.contains(nv) && max_face.first_vertex() < nv) { - // check that all edges in vertices(max_face)\cup nv are here - // since max_face is a simplex, we only need to check that edges - // (x nv) where x \in max_face are present - bool all_edges_here = true; - for (auto x : max_face) - if (!contains_edge(x, nv)) { - all_edges_here = false; - break; - } - if (all_edges_here) { // eg this->contains(max_face) - max_face.add_vertex(nv); - if (!cofaces[first_v]->contains(max_face)) { - // if there exists a blocker included in max_face, we remove it - // as it is not a minimum missing face - // the other alternative would be to check to all properfaces - // are in the complex before adding a blocker but that - // would be more expensive if there are few blockers - std::vector<Blocker_handle> blockers_to_remove; - for (auto b : blocker_range(first_v)) - if (b->contains(max_face)) - blockers_to_remove.push_back(b); - for (auto b : blockers_to_remove) - this->delete_blocker(b); - - add_blocker(max_face); - } - max_face.remove_vertex(nv); - } - } - } - } - } - for (auto i = 0u; i < num_vertices(); ++i) - delete(cofaces[i]); - } - } - - // We cannot use the default copy constructor since we need - // to make a copy of each of the blockers - - Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) { - visitor = NULL; - degree_ = copy.degree_; - skeleton = Graph(copy.skeleton); - num_vertices_ = copy.num_vertices_; - - num_blockers_ = 0; - // we copy the blockers - for (auto blocker : copy.const_blocker_range()) { - add_blocker(*blocker); - } - } - - /** - */ - Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) { - clear(); - visitor = NULL; - degree_ = copy.degree_; - skeleton = Graph(copy.skeleton); - num_vertices_ = copy.num_vertices_; - - num_blockers_ = 0; - // we copy the blockers - for (auto blocker : copy.const_blocker_range()) - add_blocker(*blocker); - return *this; - } - - - /** - * return true if both complexes have the same simplices. - */ - bool operator==(const Skeleton_blocker_complex& other) const{ - if(other.num_vertices() != num_vertices()) return false; - if(other.num_edges() != num_edges()) return false; - if(other.num_blockers() != num_blockers()) return false; - - for(auto v : vertex_range()) - if(!other.contains_vertex(v)) return false; - - for(auto e : edge_range()) - if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false; - - for(const auto b : const_blocker_range()) - if(!other.contains_blocker(*b)) return false; - - return true; - } - - bool operator!=(const Skeleton_blocker_complex& other) const{ - return !(*this == other); - } - - /** - * The destructor delete all blockers allocated. - */ - virtual ~Skeleton_blocker_complex() { - clear(); - } - - /** - * @details Clears the simplicial complex. After a call to this function, - * blockers are destroyed. The 1-skeleton and the set of blockers - * are both empty. - */ - virtual void clear() { - // xxx for now the responsabilty of freeing the visitor is for - // the user - visitor = NULL; - - degree_.clear(); - num_vertices_ = 0; - - remove_blockers(); - - skeleton.clear(); - } - - /** - *@brief allows to change the visitor. - */ - void set_visitor(Visitor* other_visitor) { - visitor = other_visitor; - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Vertices operations - */ - //@{ - public: - /** - * @brief Return a local Vertex_handle of a vertex given a global one. - * @remark Assume that the vertex is present in the complex. - */ - Vertex_handle operator[](Root_vertex_handle global) const { - auto local(get_address(global)); - assert(local); - return *local; - } - - /** - * @brief Return the vertex node associated to local Vertex_handle. - * @remark Assume that the vertex is present in the complex. - */ - Graph_vertex& operator[](Vertex_handle address) { - assert( - 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); - return skeleton[address.vertex]; - } - - /** - * @brief Return the vertex node associated to local Vertex_handle. - * @remark Assume that the vertex is present in the complex. - */ - const Graph_vertex& operator[](Vertex_handle address) const { - assert( - 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); - return skeleton[address.vertex]; - } - - /** - * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. - */ - Vertex_handle add_vertex() { - Vertex_handle address(boost::add_vertex(skeleton)); - num_vertices_++; - (*this)[address].activate(); - // safe since we now that we are in the root complex and the field 'address' and 'id' - // are identical for every vertices - (*this)[address].set_id(Root_vertex_handle(address.vertex)); - degree_.push_back(0); - if (visitor) - visitor->on_add_vertex(address); - return address; - } - - /** - * @brief Remove a vertex from the simplicial complex - * @remark It just deactivates the vertex with a boolean flag but does not - * remove it from vertices from complexity issues. - */ - void remove_vertex(Vertex_handle address) { - assert(contains_vertex(address)); - // We remove b - boost::clear_vertex(address.vertex, skeleton); - (*this)[address].deactivate(); - num_vertices_--; - degree_[address.vertex] = -1; - if (visitor) - visitor->on_remove_vertex(address); - } - - /** - */ - bool contains_vertex(Vertex_handle u) const { - if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton)) - return false; - return (*this)[u].is_active(); - } - - /** - */ - bool contains_vertex(Root_vertex_handle u) const { - boost::optional<Vertex_handle> address = get_address(u); - return address && (*this)[*address].is_active(); - } - - /** - * @return true iff the simplicial complex contains all vertices - * of simplex sigma - */ - bool contains_vertices(const Simplex_handle & sigma) const { - for (auto vertex : sigma) - if (!contains_vertex(vertex)) - return false; - return true; - } - - /** - * @brief Given an Id return the address of the vertex having this Id in the complex. - * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex. - */ - virtual boost::optional<Vertex_handle> get_address( - Root_vertex_handle id) const { - boost::optional<Vertex_handle> res; - if (id.vertex < boost::num_vertices(skeleton)) - res = Vertex_handle(id.vertex); // xxx - return res; - } - - /** - * return the id of a vertex of adress local present in the graph - */ - Root_vertex_handle get_id(Vertex_handle local) const { - assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); - return (*this)[local].get_id(); - } - - /** - * @brief Convert an address of a vertex of a complex to the address in - * the current complex. - * @details - * If the current complex is a sub (or sup) complex of 'other', it converts - * the address of a vertex v expressed in 'other' to the address of the vertex - * v in the current one. - * @remark this methods uses Root_vertex_handle to identify the vertex and - * assumes the vertex is present in the current complex. - */ - Vertex_handle convert_handle_from_another_complex( - const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const { - auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); - assert(vh_in_current_complex); - return *vh_in_current_complex; - } - - /** - * @brief return the graph degree of a vertex. - */ - int degree(Vertex_handle local) const { - assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); - return degree_[local.vertex]; - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Edges operations - */ - //@{ - public: - /** - * @brief return an edge handle if the two vertices forms - * an edge in the complex - */ - boost::optional<Edge_handle> operator[]( - const std::pair<Vertex_handle, Vertex_handle>& ab) const { - boost::optional<Edge_handle> res; - std::pair<Edge_handle, bool> edge_pair( - boost::edge(ab.first.vertex, ab.second.vertex, skeleton)); - if (edge_pair.second) - res = edge_pair.first; - return res; - } - - /** - * @brief returns the stored node associated to an edge - */ - Graph_edge& operator[](Edge_handle edge_handle) { - return skeleton[edge_handle]; - } - - /** - * @brief returns the stored node associated to an edge - */ - const Graph_edge& operator[](Edge_handle edge_handle) const { - return skeleton[edge_handle]; - } - - /** - * @brief returns the first vertex of an edge - * @details it assumes that the edge is present in the complex - */ - Vertex_handle first_vertex(Edge_handle edge_handle) const { - return static_cast<Vertex_handle> (source(edge_handle, skeleton)); - } - - /** - * @brief returns the first vertex of an edge - * @details it assumes that the edge is present in the complex - */ - Vertex_handle second_vertex(Edge_handle edge_handle) const { - return static_cast<Vertex_handle> (target(edge_handle, skeleton)); - } - - /** - * @brief returns the simplex made with the two vertices of the edge - * @details it assumes that the edge is present in the complex - - */ - Simplex_handle get_vertices(Edge_handle edge_handle) const { - auto edge((*this)[edge_handle]); - return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]); - } - - /** - * @brief Adds an edge between vertices a and b and all its cofaces. - */ - Edge_handle add_edge(Vertex_handle a, Vertex_handle b) { - assert(contains_vertex(a) && contains_vertex(b)); - assert(a != b); - - auto edge_handle((*this)[std::make_pair(a, b)]); - // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)]; - // Edge_handle edge_descr; - // bool edge_present = pair_descr_bool.second; - if (!edge_handle) { - edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first; - (*this)[*edge_handle].setId(get_id(a), get_id(b)); - degree_[a.vertex]++; - degree_[b.vertex]++; - if (visitor) - visitor->on_add_edge(a, b); - } - return *edge_handle; - } - - /** - * @brief Adds all edges and their cofaces of a simplex to the simplicial complex. - */ - void add_edges(const Simplex_handle & sigma) { - Simplex_handle_iterator i, j; - for (i = sigma.begin(); i != sigma.end(); ++i) - for (j = i, j++; j != sigma.end(); ++j) - add_edge(*i, *j); - } - - /** - * @brief Removes an edge from the simplicial complex and all its cofaces. - * @details returns the former Edge_handle representing the edge - */ - virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) { - bool found; - Edge_handle edge; - tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton); - if (found) { - if (visitor) - visitor->on_remove_edge(a, b); - // if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge); - boost::remove_edge(a.vertex, b.vertex, skeleton); - degree_[a.vertex]--; - degree_[b.vertex]--; - } - return edge; - } - - /** - * @brief Removes edge and its cofaces from the simplicial complex. - */ - void remove_edge(Edge_handle edge) { - assert(contains_vertex(first_vertex(edge))); - assert(contains_vertex(second_vertex(edge))); - remove_edge(first_vertex(edge), second_vertex(edge)); - } - - /** - * @brief The complex is reduced to its set of vertices. - * All the edges and blockers are removed. - */ - void keep_only_vertices() { - remove_blockers(); - - for (auto u : vertex_range()) { - while (this->degree(u) > 0) { - Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); - this->remove_edge(u, v); - } - } - } - - /** - * @return true iff the simplicial complex contains an edge between - * vertices a and b - */ - bool contains_edge(Vertex_handle a, Vertex_handle b) const { - // if (a.vertex<0 || b.vertex <0) return false; - return boost::edge(a.vertex, b.vertex, skeleton).second; - } - - /** - * @return true iff the simplicial complex contains all vertices - * and all edges of simplex sigma - */ - bool contains_edges(const Simplex_handle & sigma) const { - for (auto i = sigma.begin(); i != sigma.end(); ++i) { - if (!contains_vertex(*i)) - return false; - for (auto j = i; ++j != sigma.end();) { - if (!contains_edge(*i, *j)) - return false; - } - } - return true; - } - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Blockers operations - */ - //@{ - - /** - * @brief Adds the simplex to the set of blockers and - * returns a Blocker_handle toward it if was not present before and 0 otherwise. - */ - Blocker_handle add_blocker(const Simplex_handle& blocker) { - assert(blocker.dimension() > 1); - if (contains_blocker(blocker)) { - // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; - return 0; - } else { - if (visitor) - visitor->on_add_blocker(blocker); - Blocker_handle blocker_pt = new Simplex_handle(blocker); - num_blockers_++; - auto vertex = blocker_pt->begin(); - while (vertex != blocker_pt->end()) { - blocker_map_.insert(BlockerPair(*vertex, blocker_pt)); - ++vertex; - } - return blocker_pt; - } - } - - protected: - /** - * @brief Adds the simplex to the set of blockers - */ - void add_blocker(Blocker_handle blocker) { - if (contains_blocker(*blocker)) { - // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; - return; - } else { - if (visitor) - visitor->on_add_blocker(*blocker); - num_blockers_++; - auto vertex = blocker->begin(); - while (vertex != blocker->end()) { - blocker_map_.insert(BlockerPair(*vertex, blocker)); - ++vertex; - } - } - } - - protected: - /** - * Removes sigma from the blocker map of vertex v - */ - void remove_blocker(const Blocker_handle sigma, Vertex_handle v) { - Complex_blocker_around_vertex_iterator blocker; - for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end(); - ++blocker) { - if (*blocker == sigma) - break; - } - if (*blocker != sigma) { - std::cerr - << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; - assert(false); - } else { - blocker_map_.erase(blocker.current_position()); - } - } - - public: - /** - * @brief Removes the simplex from the set of blockers. - * @remark sigma has to belongs to the set of blockers - */ - void remove_blocker(const Blocker_handle sigma) { - for (auto vertex : *sigma) - remove_blocker(sigma, vertex); - num_blockers_--; - } - - /** - * @brief Remove all blockers, in other words, it expand the simplicial - * complex to the smallest flag complex that contains it. - */ - void remove_blockers() { - // Desallocate the blockers - while (!blocker_map_.empty()) { - delete_blocker(blocker_map_.begin()->second); - } - num_blockers_ = 0; - blocker_map_.clear(); - } - - protected: - /** - * Removes the simplex sigma from the set of blockers. - * sigma has to belongs to the set of blockers - * - * @remark contrarily to delete_blockers does not call the destructor - */ - void remove_blocker(const Simplex_handle& sigma) { - assert(contains_blocker(sigma)); - for (auto vertex : sigma) - remove_blocker(sigma, vertex); - num_blockers_--; - } - - public: - /** - * Removes the simplex s from the set of blockers - * and desallocate s. - */ - void delete_blocker(Blocker_handle sigma) { - if (visitor) - visitor->on_delete_blocker(sigma); - remove_blocker(sigma); - delete sigma; - } - - /** - * @return true iff s is a blocker of the simplicial complex - */ - bool contains_blocker(const Blocker_handle s) const { - if (s->dimension() < 2) - return false; - - Vertex_handle a = s->first_vertex(); - - for (const auto blocker : const_blocker_range(a)) { - if (s == *blocker) - return true; - } - return false; - } - - /** - * @return true iff s is a blocker of the simplicial complex - */ - bool contains_blocker(const Simplex_handle & s) const { - if (s.dimension() < 2) - return false; - - Vertex_handle a = s.first_vertex(); - - for (auto blocker : const_blocker_range(a)) { - if (s == *blocker) - return true; - } - return false; - } - - private: - /** - * @return true iff a blocker of the simplicial complex - * is a face of sigma. - */ - bool blocks(const Simplex_handle & sigma) const { - for (auto blocker : const_blocker_range()) { - if (sigma.contains(*blocker)) - return true; - } - return false; - } - - //@} - - protected: - /** - * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$. - * If keep_only_superior is true then only vertices that are greater than v are added. - */ - virtual void add_neighbours(Vertex_handle v, Simplex_handle & n, - bool keep_only_superior = false) const { - boost_adjacency_iterator ai, ai_end; - for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; - ++ai) { - if (keep_only_superior) { - if (*ai > v.vertex) { - n.add_vertex(Vertex_handle(*ai)); - } - } else { - n.add_vertex(Vertex_handle(*ai)); - } - } - } - - /** - * @details Add to simplex res all vertices which are - * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. - * - * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. - * todo revoir - * - */ - virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res, - bool keep_only_superior = false) const { - res.clear(); - auto alpha_vertex = alpha.begin(); - add_neighbours(*alpha_vertex, res, keep_only_superior); - for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end(); - ++alpha_vertex) - keep_neighbours(*alpha_vertex, res, keep_only_superior); - } - - /** - * @details Remove from simplex n all vertices which are - * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$. - * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. - */ - virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res, - bool keep_only_superior = false) const { - Simplex_handle nv; - add_neighbours(v, nv, keep_only_superior); - res.intersection(nv); - } - - /** - * @details Remove from simplex all vertices which are - * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$. - * If 'keep_only_superior' is true then only vertices that are greater than v are added. - */ - virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res, - bool keep_only_superior = false) const { - Simplex_handle nv; - add_neighbours(v, nv, keep_only_superior); - res.difference(nv); - } - - public: - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex; - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(Vertex_handle v) const { - return Link_complex(*this, Simplex_handle(v)); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(Edge_handle edge) const { - return Link_complex(*this, edge); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(const Simplex_handle& simplex) const { - return Link_complex(*this, simplex); - } - - /** - * @brief Compute the local vertices of 's' in the current complex - * If one of them is not present in the complex then the return value is uninitialized. - * - * - */ - // xxx rename get_address et place un using dans sub_complex - - boost::optional<Simplex_handle> get_simplex_address( - const Root_simplex_handle& s) const { - boost::optional<Simplex_handle> res; - - Simplex_handle s_address; - // Root_simplex_const_iterator i; - for (auto i = s.begin(); i != s.end(); ++i) { - boost::optional<Vertex_handle> address = get_address(*i); - if (!address) - return res; - else - s_address.add_vertex(*address); - } - res = s_address; - return res; - } - - /** - * @brief returns a simplex with vertices which are the id of vertices of the - * argument. - */ - Root_simplex_handle get_id(const Simplex_handle& local_simplex) const { - Root_simplex_handle global_simplex; - for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) { - global_simplex.add_vertex(get_id(*x)); - } - return global_simplex; - } - - /** - * @brief returns true iff the simplex s belongs to the simplicial - * complex. - */ - virtual bool contains(const Simplex_handle & s) const { - if (s.dimension() == -1) { - return false; - } else if (s.dimension() == 0) { - return contains_vertex(s.first_vertex()); - } else { - return (contains_edges(s) && !blocks(s)); - } - } - - /* - * @brief returnrs true iff the complex is empty. - */ - bool empty() const { - return num_vertices() == 0; - } - - /* - * @brief returns the number of vertices in the complex. - */ - int num_vertices() const { - // remark boost::num_vertices(skeleton) counts deactivated vertices - return num_vertices_; - } - - /* - * @brief returns the number of edges in the complex. - * @details currently in O(n) - */ - // todo cache the value - - int num_edges() const { - return boost::num_edges(skeleton); - } - - /* - * @brief returns the number of blockers in the complex. - */ - int num_blockers() const { - return num_blockers_; - } - - /* - * @brief returns true iff the graph of the 1-skeleton of the complex is complete. - */ - bool complete() const { - return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges(); - } - - /** - * @brief returns the number of connected components in the graph of the 1-skeleton. - */ - int num_connected_components() const { - int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); - std::vector<int> component(skeleton.vertex_set().size()); - return boost::connected_components(this->skeleton, &component[0]) - - num_vert_collapsed; - } - - /** - * @brief %Test if the complex is a cone. - * @details Runs in O(n) where n is the number of vertices. - */ - bool is_cone() const { - if (num_vertices() == 0) - return false; - if (num_vertices() == 1) - return true; - for (auto vi : vertex_range()) { - // xxx todo faire une methode bool is_in_blocker(Vertex_handle) - if (blocker_map_.find(vi) == blocker_map_.end()) { - // no blocker passes through the vertex, we just need to - // check if the current vertex is linked to all others vertices of the complex - if (degree_[vi.vertex] == num_vertices() - 1) - return true; - } - } - return false; - } - - //@} + template<class ComplexType> friend class Complex_vertex_iterator; + template<class ComplexType> friend class Complex_neighbors_vertices_iterator; + template<class ComplexType> friend class Complex_edge_iterator; + template<class ComplexType> friend class Complex_edge_around_vertex_iterator; + + template<class ComplexType> friend class Skeleton_blocker_link_complex; + template<class ComplexType> friend class Skeleton_blocker_link_superior; + template<class ComplexType> friend class Skeleton_blocker_sub_complex; + +public: + /** + * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; + + /** + * @brief The type of stored edge node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; + + typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; + + /** + * @brief The type of an handle to a vertex of the complex. + */ + typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; + typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; + + /** + * @brief A ordered set of integers that represents a simplex. + */ + typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; + typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; + + /** + * @brief Handle to a blocker of the complex. + */ + typedef Simplex_handle* Blocker_handle; + + typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; + typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; + +protected: + typedef typename boost::adjacency_list<boost::setS, // edges + boost::vecS, // vertices + boost::undirectedS, Graph_vertex, Graph_edge> Graph; + // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink + // (eg Link with greater vertices) + // that burdens simplex iteration / complex initialization via list of simplices. + // to avoid that, one should modify the graph by storing two lists of adjacency for every + // vertex, the one with superior and the one with lower vertices, that way there is + // no more extra cost for computation of SuperiorLink + typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; + typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; + +protected: + typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; + +public: + /** + * @brief Handle to an edge of the complex. + */ + typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; + +protected: + typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator; + +protected: + int num_vertices_; + int num_blockers_; + + typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; + // typedef Visitor* Visitor_ptr; + Visitor* visitor; + + /** + * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. + * + * This quantity is updated when adding/removing edge. + * + * This is useful because the operation + * list.size() is done in linear time. + */ + std::vector<boost_vertex_handle> degree_; + Graph skeleton; /** 1-skeleton of the simplicial complex. */ + + /** Each vertex can access to the blockers passing through it. */ + BlockerMap blocker_map_; + +public: + ///////////////////////////////////////////////////////////////////////////// + /** @name Constructors, Destructors + */ + //@{ + + /** + *@brief constructs a simplicial complex with a given number of vertices and a visitor. + */ + explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) + : visitor(visitor_) { + clear(); + for (int i = 0; i < num_vertices_; ++i) { + add_vertex(); + } + } + +private: + typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie; + + template<typename SimpleHandleOutputIterator> + /** + * return a vector of simplex trie for every vertex + */ + std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) { + std::vector<STrie*> cofaces(num_vertices(), 0); + for (auto i = 0u; i < num_vertices(); ++i) + cofaces[i] = new STrie(Vertex_handle(i)); + for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { + if (s_it->dimension() >= 1) + cofaces[s_it->first_vertex()]->add_simplex(*s_it); + } + return cofaces; + } + + +public: + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL) + : num_vertices_(0), + num_blockers_(0), + visitor(visitor_) { + std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; + // first pass to add vertices and edges + int num_vertex = -1; + for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { + if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex); + if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex()); + } + while (num_vertex-- >= 0) add_vertex(); + + for (const auto& e : edges) + add_edge(e.first, e.second); + + if (!is_flag_complex) { + // need to compute blockers + std::vector<STrie*> cofaces(compute_cofaces(simplex_begin, simplex_end)); + std::vector<Simplex_handle> max_faces; + + for (auto i = 0u; i < num_vertices(); ++i) { + auto max_faces_i = cofaces[i]->maximal_faces(); + max_faces.insert(max_faces.end(), max_faces_i.begin(), max_faces_i.end()); + } + + // all time here + + // for every maximal face s, it picks its first vertex v and check for all nv \in N(v) + // if s union nv is in the complex, if not it is a blocker. + for (auto max_face : max_faces) { + if (max_face.dimension() > 0) { + auto first_v = max_face.first_vertex(); + for (auto nv : vertex_range(first_v)) { + if (!max_face.contains(nv) && max_face.first_vertex() < nv) { + // check that all edges in vertices(max_face)\cup nv are here + // since max_face is a simplex, we only need to check that edges + // (x nv) where x \in max_face are present + bool all_edges_here = true; + for (auto x : max_face) + if (!contains_edge(x, nv)) { + all_edges_here = false; + break; + } + if (all_edges_here) { // eg this->contains(max_face) + max_face.add_vertex(nv); + if (!cofaces[first_v]->contains(max_face)) { + // if there exists a blocker included in max_face, we remove it + // as it is not a minimum missing face + // the other alternative would be to check to all properfaces + // are in the complex before adding a blocker but that + // would be more expensive if there are few blockers + std::vector<Blocker_handle> blockers_to_remove; + bool is_blocker = true; + for (auto b : blocker_range(first_v)) + if (b->contains(max_face)){ + blockers_to_remove.push_back(b); + } + else + if (max_face.contains(*b)) { + is_blocker = false; + break; + } + for (auto b : blockers_to_remove) + this->delete_blocker(b); + + if(is_blocker) + add_blocker(max_face); + } + max_face.remove_vertex(nv); + } + } + } + } + } + for (auto i = 0u; i < num_vertices(); ++i) + delete(cofaces[i]); + } + } + + // We cannot use the default copy constructor since we need + // to make a copy of each of the blockers + + Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) { + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) { + add_blocker(*blocker); + } + } + + /** + */ + Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) { + clear(); + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) + add_blocker(*blocker); + return *this; + } + + + /** + * return true if both complexes have the same simplices. + */ + bool operator==(const Skeleton_blocker_complex& other) const{ + if(other.num_vertices() != num_vertices()) return false; + if(other.num_edges() != num_edges()) return false; + if(other.num_blockers() != num_blockers()) return false; + + for(auto v : vertex_range()) + if(!other.contains_vertex(v)) return false; + + for(auto e : edge_range()) + if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false; + + for(const auto b : const_blocker_range()) + if(!other.contains_blocker(*b)) return false; + + return true; + } + + bool operator!=(const Skeleton_blocker_complex& other) const{ + return !(*this == other); + } + + /** + * The destructor delete all blockers allocated. + */ + virtual ~Skeleton_blocker_complex() { + clear(); + } + + /** + * @details Clears the simplicial complex. After a call to this function, + * blockers are destroyed. The 1-skeleton and the set of blockers + * are both empty. + */ + virtual void clear() { + // xxx for now the responsabilty of freeing the visitor is for + // the user + visitor = NULL; + + degree_.clear(); + num_vertices_ = 0; + + remove_blockers(); + + skeleton.clear(); + } + + /** + *@brief allows to change the visitor. + */ + void set_visitor(Visitor* other_visitor) { + visitor = other_visitor; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertices operations + */ + //@{ +public: + /** + * @brief Return a local Vertex_handle of a vertex given a global one. + * @remark Assume that the vertex is present in the complex. + */ + Vertex_handle operator[](Root_vertex_handle global) const { + auto local(get_address(global)); + assert(local); + return *local; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + Graph_vertex& operator[](Vertex_handle address) { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + const Graph_vertex& operator[](Vertex_handle address) const { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. + */ + Vertex_handle add_vertex() { + Vertex_handle address(boost::add_vertex(skeleton)); + num_vertices_++; + (*this)[address].activate(); + // safe since we now that we are in the root complex and the field 'address' and 'id' + // are identical for every vertices + (*this)[address].set_id(Root_vertex_handle(address.vertex)); + degree_.push_back(0); + if (visitor) + visitor->on_add_vertex(address); + return address; + } + + /** + * @brief Remove a vertex from the simplicial complex + * @remark It just deactivates the vertex with a boolean flag but does not + * remove it from vertices from complexity issues. + */ + void remove_vertex(Vertex_handle address) { + assert(contains_vertex(address)); + // We remove b + boost::clear_vertex(address.vertex, skeleton); + (*this)[address].deactivate(); + num_vertices_--; + degree_[address.vertex] = -1; + if (visitor) + visitor->on_remove_vertex(address); + } + + /** + */ + bool contains_vertex(Vertex_handle u) const { + if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton)) + return false; + return (*this)[u].is_active(); + } + + /** + */ + bool contains_vertex(Root_vertex_handle u) const { + boost::optional<Vertex_handle> address = get_address(u); + return address && (*this)[*address].is_active(); + } + + /** + * @return true iff the simplicial complex contains all vertices + * of simplex sigma + */ + bool contains_vertices(const Simplex_handle & sigma) const { + for (auto vertex : sigma) + if (!contains_vertex(vertex)) + return false; + return true; + } + + /** + * @brief Given an Id return the address of the vertex having this Id in the complex. + * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex. + */ + virtual boost::optional<Vertex_handle> get_address( + Root_vertex_handle id) const { + boost::optional<Vertex_handle> res; + if (id.vertex < boost::num_vertices(skeleton)) + res = Vertex_handle(id.vertex); // xxx + return res; + } + + /** + * return the id of a vertex of adress local present in the graph + */ + Root_vertex_handle get_id(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return (*this)[local].get_id(); + } + + /** + * @brief Convert an address of a vertex of a complex to the address in + * the current complex. + * @details + * If the current complex is a sub (or sup) complex of 'other', it converts + * the address of a vertex v expressed in 'other' to the address of the vertex + * v in the current one. + * @remark this methods uses Root_vertex_handle to identify the vertex and + * assumes the vertex is present in the current complex. + */ + Vertex_handle convert_handle_from_another_complex( + const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const { + auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); + assert(vh_in_current_complex); + return *vh_in_current_complex; + } + + /** + * @brief return the graph degree of a vertex. + */ + int degree(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return degree_[local.vertex]; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Edges operations + */ + //@{ +public: + /** + * @brief return an edge handle if the two vertices forms + * an edge in the complex + */ + boost::optional<Edge_handle> operator[]( + const std::pair<Vertex_handle, Vertex_handle>& ab) const { + boost::optional<Edge_handle> res; + std::pair<Edge_handle, bool> edge_pair( + boost::edge(ab.first.vertex, ab.second.vertex, skeleton)); + if (edge_pair.second) + res = edge_pair.first; + return res; + } + + /** + * @brief returns the stored node associated to an edge + */ + Graph_edge& operator[](Edge_handle edge_handle) { + return skeleton[edge_handle]; + } + + /** + * @brief returns the stored node associated to an edge + */ + const Graph_edge& operator[](Edge_handle edge_handle) const { + return skeleton[edge_handle]; + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle first_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (source(edge_handle, skeleton)); + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle second_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (target(edge_handle, skeleton)); + } + + /** + * @brief returns the simplex made with the two vertices of the edge + * @details it assumes that the edge is present in the complex + + */ + Simplex_handle get_vertices(Edge_handle edge_handle) const { + auto edge((*this)[edge_handle]); + return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]); + } + + /** + * @brief Adds an edge between vertices a and b and all its cofaces. + */ + Edge_handle add_edge(Vertex_handle a, Vertex_handle b) { + assert(contains_vertex(a) && contains_vertex(b)); + assert(a != b); + + auto edge_handle((*this)[std::make_pair(a, b)]); + // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)]; + // Edge_handle edge_descr; + // bool edge_present = pair_descr_bool.second; + if (!edge_handle) { + edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first; + (*this)[*edge_handle].setId(get_id(a), get_id(b)); + degree_[a.vertex]++; + degree_[b.vertex]++; + if (visitor) + visitor->on_add_edge(a, b); + } + return *edge_handle; + } + + /** + * @brief Adds all edges and their cofaces of a simplex to the simplicial complex. + */ + void add_edges(const Simplex_handle & sigma) { + Simplex_handle_iterator i, j; + for (i = sigma.begin(); i != sigma.end(); ++i) + for (j = i, j++; j != sigma.end(); ++j) + add_edge(*i, *j); + } + + /** + * @brief Removes an edge from the simplicial complex and all its cofaces. + * @details returns the former Edge_handle representing the edge + */ + virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) { + bool found; + Edge_handle edge; + tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton); + if (found) { + if (visitor) + visitor->on_remove_edge(a, b); + // if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge); + boost::remove_edge(a.vertex, b.vertex, skeleton); + degree_[a.vertex]--; + degree_[b.vertex]--; + } + return edge; + } + + /** + * @brief Removes edge and its cofaces from the simplicial complex. + */ + void remove_edge(Edge_handle edge) { + assert(contains_vertex(first_vertex(edge))); + assert(contains_vertex(second_vertex(edge))); + remove_edge(first_vertex(edge), second_vertex(edge)); + } + + /** + * @brief The complex is reduced to its set of vertices. + * All the edges and blockers are removed. + */ + void keep_only_vertices() { + remove_blockers(); + + for (auto u : vertex_range()) { + while (this->degree(u) > 0) { + Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); + this->remove_edge(u, v); + } + } + } + + /** + * @return true iff the simplicial complex contains an edge between + * vertices a and b + */ + bool contains_edge(Vertex_handle a, Vertex_handle b) const { + // if (a.vertex<0 || b.vertex <0) return false; + return boost::edge(a.vertex, b.vertex, skeleton).second; + } + + /** + * @return true iff the simplicial complex contains all vertices + * and all edges of simplex sigma + */ + bool contains_edges(const Simplex_handle & sigma) const { + for (auto i = sigma.begin(); i != sigma.end(); ++i) { + if (!contains_vertex(*i)) + return false; + for (auto j = i; ++j != sigma.end();) { + if (!contains_edge(*i, *j)) + return false; + } + } + return true; + } + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Blockers operations + */ + //@{ + + /** + * @brief Adds the simplex to the set of blockers and + * returns a Blocker_handle toward it if was not present before and 0 otherwise. + */ + Blocker_handle add_blocker(const Simplex_handle& blocker) { + assert(blocker.dimension() > 1); + if (contains_blocker(blocker)) { + // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return 0; + } else { + if (visitor) + visitor->on_add_blocker(blocker); + Blocker_handle blocker_pt = new Simplex_handle(blocker); + num_blockers_++; + auto vertex = blocker_pt->begin(); + while (vertex != blocker_pt->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker_pt)); + ++vertex; + } + return blocker_pt; + } + } + +protected: + /** + * @brief Adds the simplex to the set of blockers + */ + void add_blocker(Blocker_handle blocker) { + if (contains_blocker(*blocker)) { + // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return; + } else { + if (visitor) + visitor->on_add_blocker(*blocker); + num_blockers_++; + auto vertex = blocker->begin(); + while (vertex != blocker->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker)); + ++vertex; + } + } + } + +protected: + /** + * Removes sigma from the blocker map of vertex v + */ + void remove_blocker(const Blocker_handle sigma, Vertex_handle v) { + Complex_blocker_around_vertex_iterator blocker; + for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end(); + ++blocker) { + if (*blocker == sigma) + break; + } + if (*blocker != sigma) { + std::cerr + << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; + assert(false); + } else { + blocker_map_.erase(blocker.current_position()); + } + } + +public: + /** + * @brief Removes the simplex from the set of blockers. + * @remark sigma has to belongs to the set of blockers + */ + void remove_blocker(const Blocker_handle sigma) { + for (auto vertex : *sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + + /** + * @brief Remove all blockers, in other words, it expand the simplicial + * complex to the smallest flag complex that contains it. + */ + void remove_blockers() { + // Desallocate the blockers + while (!blocker_map_.empty()) { + delete_blocker(blocker_map_.begin()->second); + } + num_blockers_ = 0; + blocker_map_.clear(); + } + +protected: + /** + * Removes the simplex sigma from the set of blockers. + * sigma has to belongs to the set of blockers + * + * @remark contrarily to delete_blockers does not call the destructor + */ + void remove_blocker(const Simplex_handle& sigma) { + assert(contains_blocker(sigma)); + for (auto vertex : sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + +public: + /** + * Removes the simplex s from the set of blockers + * and desallocate s. + */ + void delete_blocker(Blocker_handle sigma) { + if (visitor) + visitor->on_delete_blocker(sigma); + remove_blocker(sigma); + delete sigma; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Blocker_handle s) const { + if (s->dimension() < 2) + return false; + + Vertex_handle a = s->first_vertex(); + + for (const auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Simplex_handle & s) const { + if (s.dimension() < 2) + return false; + + Vertex_handle a = s.first_vertex(); + + for (auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + +private: + /** + * @return true iff a blocker of the simplicial complex + * is a face of sigma. + */ + bool blocks(const Simplex_handle & sigma) const { + for (auto blocker : const_blocker_range()) { + if (sigma.contains(*blocker)) + return true; + } + return false; + } + + //@} + +protected: + /** + * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$. + * If keep_only_superior is true then only vertices that are greater than v are added. + */ + virtual void add_neighbours(Vertex_handle v, Simplex_handle & n, + bool keep_only_superior = false) const { + boost_adjacency_iterator ai, ai_end; + for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; + ++ai) { + if (keep_only_superior) { + if (*ai > v.vertex) { + n.add_vertex(Vertex_handle(*ai)); + } + } else { + n.add_vertex(Vertex_handle(*ai)); + } + } + } + + /** + * @details Add to simplex res all vertices which are + * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. + * todo revoir + * + */ + virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res, + bool keep_only_superior = false) const { + res.clear(); + auto alpha_vertex = alpha.begin(); + add_neighbours(*alpha_vertex, res, keep_only_superior); + for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end(); + ++alpha_vertex) + keep_neighbours(*alpha_vertex, res, keep_only_superior); + } + + /** + * @details Remove from simplex n all vertices which are + * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. + */ + virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res, + bool keep_only_superior = false) const { + Simplex_handle nv; + add_neighbours(v, nv, keep_only_superior); + res.intersection(nv); + } + + /** + * @details Remove from simplex all vertices which are + * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are added. + */ + virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res, + bool keep_only_superior = false) const { + Simplex_handle nv; + add_neighbours(v, nv, keep_only_superior); + res.difference(nv); + } + +public: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Vertex_handle v) const { + return Link_complex(*this, Simplex_handle(v)); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Edge_handle edge) const { + return Link_complex(*this, edge); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(const Simplex_handle& simplex) const { + return Link_complex(*this, simplex); + } + + /** + * @brief Compute the local vertices of 's' in the current complex + * If one of them is not present in the complex then the return value is uninitialized. + * + * + */ + // xxx rename get_address et place un using dans sub_complex + + boost::optional<Simplex_handle> get_simplex_address( + const Root_simplex_handle& s) const { + boost::optional<Simplex_handle> res; + + Simplex_handle s_address; + // Root_simplex_const_iterator i; + for (auto i = s.begin(); i != s.end(); ++i) { + boost::optional<Vertex_handle> address = get_address(*i); + if (!address) + return res; + else + s_address.add_vertex(*address); + } + res = s_address; + return res; + } + + /** + * @brief returns a simplex with vertices which are the id of vertices of the + * argument. + */ + Root_simplex_handle get_id(const Simplex_handle& local_simplex) const { + Root_simplex_handle global_simplex; + for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) { + global_simplex.add_vertex(get_id(*x)); + } + return global_simplex; + } + + /** + * @brief returns true iff the simplex s belongs to the simplicial + * complex. + */ + virtual bool contains(const Simplex_handle & s) const { + if (s.dimension() == -1) { + return false; + } else if (s.dimension() == 0) { + return contains_vertex(s.first_vertex()); + } else { + return (contains_edges(s) && !blocks(s)); + } + } + + /* + * @brief returnrs true iff the complex is empty. + */ + bool empty() const { + return num_vertices() == 0; + } + + /* + * @brief returns the number of vertices in the complex. + */ + int num_vertices() const { + // remark boost::num_vertices(skeleton) counts deactivated vertices + return num_vertices_; + } + + /* + * @brief returns the number of edges in the complex. + * @details currently in O(n) + */ + // todo cache the value + + int num_edges() const { + return boost::num_edges(skeleton); + } + + int num_triangles() const{ + auto triangles = triangle_range(); + return std::distance(triangles.begin(),triangles.end()); + } + /* + * @brief returns the number of blockers in the complex. + */ + int num_blockers() const { + return num_blockers_; + } + + /* + * @brief returns true iff the graph of the 1-skeleton of the complex is complete. + */ + bool complete() const { + return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges(); + } + + /** + * @brief returns the number of connected components in the graph of the 1-skeleton. + */ + int num_connected_components() const { + int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); + std::vector<int> component(skeleton.vertex_set().size()); + return boost::connected_components(this->skeleton, &component[0]) + - num_vert_collapsed; + } + + /** + * @brief %Test if the complex is a cone. + * @details Runs in O(n) where n is the number of vertices. + */ + bool is_cone() const { + if (num_vertices() == 0) + return false; + if (num_vertices() == 1) + return true; + for (auto vi : vertex_range()) { + // xxx todo faire une methode bool is_in_blocker(Vertex_handle) + if (blocker_map_.find(vi) == blocker_map_.end()) { + // no blocker passes through the vertex, we just need to + // check if the current vertex is linked to all others vertices of the complex + if (degree_[vi.vertex] == num_vertices() - 1) + return true; + } + } + return false; + } + + //@} /** @Simplification operations */ //@{ @@ -1244,279 +1255,279 @@ private: //@} public: - ///////////////////////////////////////////////////////////////////////////// - /** @name Vertex iterators - */ - //@{ - typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI; // todo rename - - // - // @brief Range over the vertices of the simplicial complex. - // Methods .begin() and .end() return a Complex_vertex_iterator. - // - typedef boost::iterator_range< - Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range; - - /** - * @brief Returns a Complex_vertex_range over all vertices of the complex - */ - Complex_vertex_range vertex_range() const { - auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this); - auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this, 0); - return Complex_vertex_range(begin, end); - } - - typedef boost::iterator_range< - Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v - */ - Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const { - auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>( - this, v); - auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>( - this, v, 0); - return Complex_neighbors_vertices_range(begin, end); - } - - //@} - - /** @name Edge iterators - */ - //@{ - - typedef boost::iterator_range< - Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> Complex_edge_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex - */ - Complex_edge_range edge_range() const { - auto begin = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this); - auto end = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, 0); - return Complex_edge_range(begin, end); - } - - typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> - Complex_edge_around_vertex_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes - * through 'v' - */ - Complex_edge_around_vertex_range edge_range(Vertex_handle v) const { - auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v); - auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v, 0); - return Complex_edge_around_vertex_range(begin, end); - } - - //@} - - /** @name Triangles iterators - */ - //@{ - private: - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; - typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; - - public: - typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator; - - typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range; - - /** - * @brief Range over triangles around a vertex of the simplicial complex. - * Methods .begin() and .end() return a Triangle_around_vertex_iterator. - * - */ - Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const { - auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v); - auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0); - return Complex_triangle_around_vertex_range(begin, end); - } - - typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; - - /** - * @brief Range over triangles of the simplicial complex. - * Methods .begin() and .end() return a Triangle_around_vertex_iterator. - * - */ - Complex_triangle_range triangle_range() const { - auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0); - if (empty()) { - return Complex_triangle_range(end, end); - } else { - auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); - return Complex_triangle_range(begin, end); - } - } - - //@} - - /** @name Simplices iterators - */ - //@{ - typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator; - - /** - * @brief Range over the simplices of the simplicial complex around a vertex. - * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. - */ - typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; - - /** - * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex - */ - Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const { - assert(contains_vertex(v)); - return Complex_simplex_around_vertex_range( - Complex_simplex_around_vertex_iterator(this, v), - Complex_simplex_around_vertex_iterator(this, v, true)); - } - - // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; - typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; - - typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; - - /** - * @brief Returns a Complex_simplex_range over all the simplices of the complex - */ - Complex_simplex_range simplex_range() const { - Complex_simplex_iterator end(this, true); - if (empty()) { - return Complex_simplex_range(end, end); - } else { - Complex_simplex_iterator begin(this); - return Complex_simplex_range(begin, end); - } - } - - //@} - - /** @name Blockers iterators - */ - //@{ - private: - /** - * @brief Iterator over the blockers adjacent to a vertex - */ - typedef Blocker_iterator_around_vertex_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, - Blocker_handle> - Complex_blocker_around_vertex_iterator; - - /** - * @brief Iterator over (constant) blockers adjacent to a vertex - */ - typedef Blocker_iterator_around_vertex_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, - const Blocker_handle> - Const_complex_blocker_around_vertex_iterator; - - typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; - typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range; - - public: - /** - * @brief Returns a range of the blockers of the complex passing through a vertex - */ - Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) { - auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); - auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); - return Complex_blocker_around_vertex_range(begin, end); - } - - /** - * @brief Returns a range of the blockers of the complex passing through a vertex - */ - Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const { - auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); - auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); - return Const_complex_blocker_around_vertex_range(begin, end); - } - - private: - /** - * @brief Iterator over the blockers. - */ - typedef Blocker_iterator_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, - Blocker_handle> - Complex_blocker_iterator; - - /** - * @brief Iterator over the (constant) blockers. - */ - typedef Blocker_iterator_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, - const Blocker_handle> - Const_complex_blocker_iterator; - - typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; - typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; - - public: - /** - * @brief Returns a range of the blockers of the complex - */ - Complex_blocker_range blocker_range() { - auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); - auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); - return Complex_blocker_range(begin, end); - } - - /** - * @brief Returns a range of the blockers of the complex - */ - Const_complex_blocker_range const_blocker_range() const { - auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); - auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); - return Const_complex_blocker_range(begin, end); - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Print and IO methods - */ - //@{ - public: - std::string to_string() const { - std::ostringstream stream; - stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl; - stream << num_edges() << " edges:\n" << edges_to_string() << std::endl; - stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl; - return stream.str(); - } - - std::string vertices_to_string() const { - std::ostringstream stream; - for (auto vertex : vertex_range()) { - stream << "{" << (*this)[vertex].get_id() << "} "; - } - stream << std::endl; - return stream.str(); - } - - std::string edges_to_string() const { - std::ostringstream stream; - for (auto edge : edge_range()) - stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} "; - stream << std::endl; - return stream.str(); - } - - std::string blockers_to_string() const { - std::ostringstream stream; - - for (auto b : const_blocker_range()) - stream << *b << std::endl; - return stream.str(); - } - //@} + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertex iterators + */ + //@{ + typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI; // todo rename + + // + // @brief Range over the vertices of the simplicial complex. + // Methods .begin() and .end() return a Complex_vertex_iterator. + // + typedef boost::iterator_range< + Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range; + + /** + * @brief Returns a Complex_vertex_range over all vertices of the complex + */ + Complex_vertex_range vertex_range() const { + auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this); + auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this, 0); + return Complex_vertex_range(begin, end); + } + + typedef boost::iterator_range< + Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v + */ + Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const { + auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>( + this, v); + auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>( + this, v, 0); + return Complex_neighbors_vertices_range(begin, end); + } + + //@} + + /** @name Edge iterators + */ + //@{ + + typedef boost::iterator_range< + Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> Complex_edge_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex + */ + Complex_edge_range edge_range() const { + auto begin = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this); + auto end = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, 0); + return Complex_edge_range(begin, end); + } + + typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> + Complex_edge_around_vertex_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes + * through 'v' + */ + Complex_edge_around_vertex_range edge_range(Vertex_handle v) const { + auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v); + auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v, 0); + return Complex_edge_around_vertex_range(begin, end); + } + + //@} + + /** @name Triangles iterators + */ + //@{ +private: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; + typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; + +public: + typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator; + + typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range; + + /** + * @brief Range over triangles around a vertex of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const { + auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v); + auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0); + return Complex_triangle_around_vertex_range(begin, end); + } + + typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; + + /** + * @brief Range over triangles of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_range triangle_range() const { + auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0); + if (empty()) { + return Complex_triangle_range(end, end); + } else { + auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); + return Complex_triangle_range(begin, end); + } + } + + //@} + + /** @name Simplices iterators + */ + //@{ + typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator; + + /** + * @brief Range over the simplices of the simplicial complex around a vertex. + * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. + */ + typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; + + /** + * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex + */ + Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const { + assert(contains_vertex(v)); + return Complex_simplex_around_vertex_range( + Complex_simplex_around_vertex_iterator(this, v), + Complex_simplex_around_vertex_iterator(this, v, true)); + } + + // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; + typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; + + typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; + + /** + * @brief Returns a Complex_simplex_range over all the simplices of the complex + */ + Complex_simplex_range simplex_range() const { + Complex_simplex_iterator end(this, true); + if (empty()) { + return Complex_simplex_range(end, end); + } else { + Complex_simplex_iterator begin(this); + return Complex_simplex_range(begin, end); + } + } + + //@} + + /** @name Blockers iterators + */ + //@{ +private: + /** + * @brief Iterator over the blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_around_vertex_iterator; + + /** + * @brief Iterator over (constant) blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_around_vertex_iterator; + + typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; + typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range; + +public: + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) { + auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Complex_blocker_around_vertex_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const { + auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Const_complex_blocker_around_vertex_range(begin, end); + } + +private: + /** + * @brief Iterator over the blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_iterator; + + /** + * @brief Iterator over the (constant) blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_iterator; + + typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; + typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; + +public: + /** + * @brief Returns a range of the blockers of the complex + */ + Complex_blocker_range blocker_range() { + auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Complex_blocker_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex + */ + Const_complex_blocker_range const_blocker_range() const { + auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Const_complex_blocker_range(begin, end); + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Print and IO methods + */ + //@{ +public: + std::string to_string() const { + std::ostringstream stream; + stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl; + stream << num_edges() << " edges:\n" << edges_to_string() << std::endl; + stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl; + return stream.str(); + } + + std::string vertices_to_string() const { + std::ostringstream stream; + for (auto vertex : vertex_range()) { + stream << "{" << (*this)[vertex].get_id() << "} "; + } + stream << std::endl; + return stream.str(); + } + + std::string edges_to_string() const { + std::ostringstream stream; + for (auto edge : edge_range()) + stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} "; + stream << std::endl; + return stream.str(); + } + + std::string blockers_to_string() const { + std::ostringstream stream; + + for (auto b : const_blocker_range()) + stream << *b << std::endl; + return stream.str(); + } + //@} @@ -1533,78 +1544,13 @@ public: */ template<typename Complex, typename SimplexHandleIterator> Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, bool is_flag_complex = false) { - typedef typename Complex::Simplex_handle Simplex_handle; - typedef typename Complex::Blocker_handle Blocker_handle; - typedef typename Complex::Vertex_handle Vertex_handle; - Complex complex; - - complex.clear(); - - std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; - // first pass to add vertices and edges - for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { - // if meet simplex 9 12 15, need to have at least 16 vertices - int max_vertex = 0; - for (auto vh : *s_it) - max_vertex = (std::max)(max_vertex, static_cast<int>(vh)); - while (complex.num_vertices() <= max_vertex) - complex.add_vertex(); - - // for all pairs in s, add an edge - for (auto u_it = s_it->begin(); u_it != s_it->end(); ++u_it) - for (auto v_it = u_it; ++v_it != s_it->end(); /**/) - complex.add_edge(*u_it, *v_it); - } - - if (!is_flag_complex) { - // need to compute blockers - - // store a structure to decide faster if a simplex is in the complex defined by the max faces - - std::vector<std::set < Simplex_handle >> vertex_to_maxfaces(complex.num_vertices()); - for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) - vertex_to_maxfaces[s_it->first_vertex()].insert(*s_it); - - // for every maximal face s, it picks its first vertex v and check for all nv \in N(v) - // if s union nv is in the complex, if not it is a blocker. - for (auto max_face = simplex_begin; max_face != simplex_end; ++max_face) { - if (max_face->dimension() > 0) { - auto first_v = max_face->first_vertex(); - for (auto nv : complex.vertex_range(first_v)) { - if (!max_face->contains(nv) && max_face->first_vertex() < nv) { - // check that all edges in vertices(max_face)\cup nv are here - // since max_face is a simplex, we only need to check that edges - // (x nv) where x \in max_face are present - bool all_edges_here = true; - for (auto x : *max_face) - if (!complex.contains_edge(x, nv)) { - all_edges_here = false; - break; - } - if (all_edges_here) { // eg this->contains(max_face) - max_face->add_vertex(nv); - if (vertex_to_maxfaces[first_v].find(*max_face) == vertex_to_maxfaces[first_v].end()) { // xxxx - // if there exists a blocker included in max_face, we remove it - // as it is not a minimum missing face - // the other alternative would be to check to all properfaces - // are in the complex before adding a blocker but that - // would be more expensive if there are few blockers - std::vector<Blocker_handle> blockers_to_remove; - for (auto b : complex.blocker_range(first_v)) - if (b->contains(*max_face)) - blockers_to_remove.push_back(b); - for (auto b : blockers_to_remove) - complex.delete_blocker(b); - complex.add_blocker(*max_face); - } - max_face->remove_vertex(nv); - } - } - } - } - } - } - return complex; + typedef typename Complex::Simplex_handle Simplex_handle; + std::vector<Simplex_handle> simplices; + for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) { + auto subfaces_topface = subfaces(*top_face); + simplices.insert(simplices.end(),subfaces_topface.begin(),subfaces_topface.end()); + } + return Complex(simplices.begin(),simplices.end(),is_flag_complex); } } // namespace skbl diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h index f2c1e447..bb58d0dc 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h @@ -37,98 +37,162 @@ namespace skbl { */ template<typename SkeletonBlockerGeometricDS> class Skeleton_blocker_geometric_complex : -public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> { - public: - typedef typename SkeletonBlockerGeometricDS::GT GT; - - typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; - - typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; - typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; - typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; - typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle; - - typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; - - typedef typename SkeletonBlockerGeometricDS::Point Point; - - /** - * @brief Add a vertex to the complex with its associated point. - */ - Vertex_handle add_vertex(const Point& point) { - Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); - (*this)[ad].point() = point; - return ad; - } - - /** - * @brief Returns the Point associated to the vertex v. - */ - const Point& point(Vertex_handle v) const { - assert(this->contains_vertex(v)); - return (*this)[v].point(); - } - - /** - * @brief Returns the Point associated to the vertex v. - */ - Point& point(Vertex_handle v) { - assert(this->contains_vertex(v)); - return (*this)[v].point(); - } - - const Point& point(Root_vertex_handle global_v) const { - Vertex_handle local_v((*this)[global_v]); - assert(this->contains_vertex(local_v)); - return (*this)[local_v].point(); - } - - Point& point(Root_vertex_handle global_v) { - Vertex_handle local_v((*this)[global_v]); - assert(this->contains_vertex(local_v)); - return (*this)[local_v].point(); - } - - typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(Vertex_handle v) const { - Geometric_link link(*this, Simplex_handle(v)); - // we now add the point info - add_points_to_link(link); - return link; - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(const Simplex_handle& simplex) const { - Geometric_link link(*this, simplex); - // we now add the point info - add_points_to_link(link); - return link; - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(Edge_handle edge) const { - Geometric_link link(*this, edge); - // we now add the point info - add_points_to_link(link); - return link; - } - - private: - - void add_points_to_link(Geometric_link& link) const { - for (Vertex_handle v : link.vertex_range()) { - Root_vertex_handle v_root(link.get_id(v)); - link.point(v) = (*this).point(v_root); - } - } + public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> { + public: + typedef typename SkeletonBlockerGeometricDS::GT GT; + + typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; + + typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; + typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; + typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; + typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle; + + typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; + + typedef typename SkeletonBlockerGeometricDS::Point Point; + + + + Skeleton_blocker_geometric_complex(){ + } + + + + /** + * constructor given a list of points + */ + template<typename PointIterator> + explicit Skeleton_blocker_geometric_complex(int num_vertices,PointIterator begin,PointIterator end){ + for(auto point = begin; point != end; ++point) + add_vertex(*point); + } + + + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator,typename PointIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + PointIterator points_begin,PointIterator points_end, + bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){ + unsigned current = 0; + for(auto point = points_begin; point != points_end; ++point) + (*this)[Vertex_handle(current++)].point() = Point(point->begin(),point->end()); + } + + template<typename SimpleHandleOutputIterator,typename PointIterator> + friend Skeleton_blocker_geometric_complex make_complex_from_top_faces( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + PointIterator points_begin,PointIterator points_end, + bool is_flag_complex = false){ + Skeleton_blocker_geometric_complex complex; + unsigned current = 0; + complex=make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin,simplex_end,is_flag_complex); + for(auto point = points_begin; point != points_end; ++point) + complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end()); + return complex; + } + + + /** + * @brief Constructor with a list of simplices. + * Points of every vertex are the point constructed with default constructor. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){ + } + + + /** + * @brief Add a vertex to the complex with a default constructed associated point. + */ + Vertex_handle add_vertex() { + return SimplifiableSkeletonblocker::add_vertex(); + } + + /** + * @brief Add a vertex to the complex with its associated point. + */ + Vertex_handle add_vertex(const Point& point) { + Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); + (*this)[ad].point() = point; + return ad; + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + const Point& point(Vertex_handle v) const { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + Point& point(Vertex_handle v) { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + const Point& point(Root_vertex_handle global_v) const { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + Point& point(Root_vertex_handle global_v) { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Vertex_handle v) const { + Geometric_link link(*this, Simplex_handle(v)); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(const Simplex_handle& simplex) const { + Geometric_link link(*this, simplex); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Edge_handle edge) const { + Geometric_link link(*this, edge); + // we now add the point info + add_points_to_link(link); + return link; + } + + private: + + void add_points_to_link(Geometric_link& link) const { + for (Vertex_handle v : link.vertex_range()) { + Root_vertex_handle v_root(link.get_id(v)); + link.point(v) = (*this).point(v_root); + } + } }; } // namespace skbl diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h index 163fb7e3..9eab7f1e 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -468,6 +468,17 @@ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_handle& sigma) { assert(!this->contains(sigma)); assert(sigma.dimension() > 1); + + int num_vertex_to_add = 0; + for(auto v : sigma) + if(!contains_vertex(v)) ++num_vertex_to_add; + while(num_vertex_to_add--) add_vertex(); + + for(auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it) + for(auto v_it = u_it; ++v_it != sigma.end(); /**/){ + std::cout <<"add edge"<<*u_it<<" "<<*v_it<<std::endl; + add_edge(*u_it,*v_it); + } remove_blocker_include_in_simplex(sigma); } diff --git a/src/Skeleton_blocker/test/CMakeLists.txt b/src/Skeleton_blocker/test/CMakeLists.txt index d70b5d7f..c69bfec7 100644 --- a/src/Skeleton_blocker/test/CMakeLists.txt +++ b/src/Skeleton_blocker/test/CMakeLists.txt @@ -11,9 +11,12 @@ endif() add_executable(TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp) add_executable(TestSimplifiable TestSimplifiable.cpp) +add_executable(TestGeometricComplex TestGeometricComplex.cpp) add_test(TestSkeletonBlockerComplex ${CMAKE_CURRENT_BINARY_DIR}/TestSkeletonBlockerComplex) add_test(TestSimplifiable ${CMAKE_CURRENT_BINARY_DIR}/TestSimplifiable) +add_test(TestGeometricComplex ${CMAKE_CURRENT_BINARY_DIR}/TestGeometricComplex) + if (LCOV_PATH) # Lcov code coverage of unitary test diff --git a/src/Skeleton_blocker/test/TestGeometricComplex.cpp b/src/Skeleton_blocker/test/TestGeometricComplex.cpp new file mode 100644 index 00000000..1a81b7f7 --- /dev/null +++ b/src/Skeleton_blocker/test/TestGeometricComplex.cpp @@ -0,0 +1,85 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * + * This program is free software: 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. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string> +#include <fstream> +#include <sstream> +#include "gudhi/Test.h" +#include "gudhi/Skeleton_blocker.h" + + +using namespace std; +using namespace Gudhi; +using namespace skbl; + +struct Geometry_trait{ + typedef std::vector<double> Point; +}; + +typedef Geometry_trait::Point Point; +typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits; +typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex; + + + +bool test_constructor1(){ + Complex complex; + Skeleton_blocker_off_reader<Complex> off_reader("test2.off",complex); + if(!off_reader.is_valid()){ + std::cerr << "Unable to read file"<<std::endl; + return false; + } + + + std::cout << "complex has "<< + complex.num_vertices()<<" vertices, "<< + complex.num_blockers()<<" blockers, "<< + complex.num_edges()<<" edges and" << + complex.num_triangles()<<" triangles."; + + if(complex.num_vertices()!=7 || complex.num_edges()!=12 || complex.num_triangles() !=6) + return false; + + Skeleton_blocker_off_writer<Complex> off_writer("tmp.off",complex); + Complex same; + Skeleton_blocker_off_reader<Complex> off_reader2("tmp.off",same); + + std::cout<<"\ncomplex:"<<complex.to_string()<<endl; + std::cout<<"\nsame:"<<same.to_string()<<endl; + + return (complex==same); +} + + + +int main (int argc, char *argv[]) +{ + Tests tests_geometric_complex; + tests_geometric_complex.add("Test constructor 1",test_constructor1); + + if(tests_geometric_complex.run()) + return EXIT_SUCCESS; + else + return EXIT_FAILURE; +} diff --git a/src/Skeleton_blocker/test/test.off b/src/Skeleton_blocker/test/test.off new file mode 100644 index 00000000..61ef2d2b --- /dev/null +++ b/src/Skeleton_blocker/test/test.off @@ -0,0 +1,10 @@ +OFF +4 4 0 +0 0 0 +1 0 0 +0 1 0 +.5 .5 1 +3 0 1 2 +3 0 1 3 +3 0 2 3 +3 1 2 3 diff --git a/src/Skeleton_blocker/test/test2.off b/src/Skeleton_blocker/test/test2.off new file mode 100644 index 00000000..b12ce572 --- /dev/null +++ b/src/Skeleton_blocker/test/test2.off @@ -0,0 +1,15 @@ +OFF +7 6 0 +-1.5 2 0 +0 2 0 +1.5 2 0 +0 1.5 0 +-1 1 0 +1 1 0 +0 0 0 +3 0 1 4 +3 1 3 4 +3 1 3 5 +3 1 2 5 +3 3 4 5 +3 4 5 6 |