From 89afd70e622d03bca4e4189ada72d602b7070992 Mon Sep 17 00:00:00 2001 From: "ulrich.bauer@gmail.com" Date: Mon, 29 Apr 2013 11:41:20 +0000 Subject: alpha shapes git-svn-id: https://phat.googlecode.com/svn/trunk@46 8e3bb3c2-eed4-f18f-5264-0b6c94e6926d --- addons/alpha_3.cpp | 317 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 addons/alpha_3.cpp (limited to 'addons') diff --git a/addons/alpha_3.cpp b/addons/alpha_3.cpp new file mode 100644 index 0000000..548500a --- /dev/null +++ b/addons/alpha_3.cpp @@ -0,0 +1,317 @@ +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +template +struct Vertex_info_3 { + + typedef Index_ Index; + + Vertex_info_3() { + index_=boost::none; + } + + bool has_index() { + return index_; + } + + Index index() { + CGAL_assertion(has_index()); + return index_.get(); + } + + void set_index(Index I) { + index_=I; + } + +private: + boost::optional index_; +}; + + +template +struct Cell_info_3 { + + typedef Index_ Index; + + Cell_info_3() { + for(std::size_t i=0; i<6; i++) { + edge_index_[i] = boost::none; + } + for(std::size_t i=0; i<4; i++) { + facet_index_[i] = boost::none; + } + } + + int edge_conv(int i, int j) { + if(i>j) std::swap(i,j); + if(i==0 && j==1) return 0; + if(i==0 && j==2) return 1; + if(i==0 && j==3) return 2; + if(i==1 && j==2) return 3; + if(i==1 && j==3) return 4; + if(i==2 && j==3) return 5; + } + + + bool has_edge_index(int i, int j) { + return edge_index_[edge_conv(i,j)]; + } + + Index edge_index(int i, int j) { + CGAL_assertion(has_edge_index(i,j)); + int k = edge_conv(i,j); + return edge_index_[k].get(); + } + + bool has_facet_index(int i) { + CGAL_assertion(i>=0 && i<4); + return facet_index_[i]; + } + + Index facet_index(int i) { + CGAL_assertion(has_facet_index(i)); + return facet_index_[i].get(); + } + + void set_edge_index(int i, int j, Index I) { + edge_index_[edge_conv(i,j)]=I; + } + + void set_facet_index(int i, Index I) { + facet_index_[i]=I; + } + +private: + + boost::optional edge_index_[6]; + boost::optional facet_index_[4]; + +}; + +typedef phat::dimension Dim; +typedef std::vector Dim_container; +typedef phat::index Index; +typedef phat::column Column; +typedef std::vector > Matrix; + +typedef CGAL::Exact_predicates_exact_constructions_kernel Gt; +typedef Gt::FT FT; + +typedef CGAL::Triangulation_cell_base_with_info_3,Gt> Cb; +typedef CGAL::Triangulation_vertex_base_with_info_3,Gt> Vb; +typedef CGAL::Triangulation_hierarchy_vertex_base_3 Vbh; +typedef CGAL::Triangulation_data_structure_3 Tds; +typedef CGAL::Delaunay_triangulation_3 Triangulation_3; +typedef CGAL::Triangulation_hierarchy_3 Triangulation_hierarchy; + +typedef Gt::Point_3 Point; + + +typedef Triangulation_3::Vertex_handle Vertex_handle; +typedef Triangulation_3::Edge Edge; +typedef Triangulation_3::Facet Facet; +typedef Triangulation_3::Cell_handle Cell_handle; + +typedef Triangulation_3::Cell_circulator Cell_circulator; + +void set_index_of_edge(const Triangulation_3& T, const Edge& e, Index I) { + + Vertex_handle v1 = e.first->vertex(e.second); + Vertex_handle v2 = e.first->vertex(e.third); + + Cell_circulator ch=T.incident_cells(e); + Cell_circulator ch_start=ch; + int count=0; + do { + ch->info().set_edge_index(ch->index(v1),ch->index(v2),I); + ch++; + count++; + } while(ch!=ch_start); + //std::cout << "Did " << count << " updates" << std::endl; +} + +void set_index_of_facet(const Triangulation_3& T, const Facet& f, Index I) { + + f.first->info().set_facet_index(f.second,I); + Facet mf = T.mirror_facet(f); + mf.first->info().set_facet_index(mf.second,I); + +} + +template +struct Sort_triples { + + bool operator() (const Triple& a, const Triple& b) { + if(a.first < b.first) return true; + if(a.first > b.first) return false; + return a.second < b.second; + } + +}; + +int main(int argc, char** argv) +{ + + if(argc<2) { + std::cerr << "Need input file" << std::endl; + std::exit(1); + } + + std::list lp; + Point p; + + //read input + std::ifstream is(argv[1]); + std::string next_line; + while( getline( is, next_line ) ) { + if( next_line != "" && next_line[ 0 ] != '#' ) { + std::cerr << next_line; + std::stringstream sstr(next_line); + sstr >> p; + lp.push_back(p); + } + } + + std::cerr << "Compute Delaunay triangulation..." << std::endl; + Triangulation_hierarchy dt(lp.begin(),lp.end()); +// for (std::list::iterator it = lp.begin(); it != lp.end(); ++it) { +// dt.insert(*it); +// }/Users/uli/Downloads/neptune-raw.off/782_neptune-raw.pts.txt + + std::cerr << "Compute circumradii..." << std::endl; + + typedef CGAL::Triple Triple; + + std::vector circumradii; + + for (Triangulation_hierarchy::Finite_vertices_iterator vertex = dt.finite_vertices_begin(); + vertex != dt.finite_vertices_end(); vertex++) { + Vertex_handle vh(vertex); + circumradii.push_back(CGAL::make_triple(FT(0),0,CGAL::make_object(vh))); + } + for (Triangulation_hierarchy::Finite_edges_iterator edge = dt.finite_edges_begin(); + edge != dt.finite_edges_end(); edge++) { + Vertex_handle v1 = edge->first->vertex(edge->second); + Vertex_handle v2 = edge->first->vertex(edge->third); + circumradii.push_back(CGAL::make_triple(CGAL::squared_radius(v1->point(),v2->point()),1,CGAL::make_object(*edge))); + } + for (Triangulation_hierarchy::Finite_facets_iterator f = dt.finite_facets_begin(); + f != dt.finite_facets_end(); f++) { + Vertex_handle v1 = f->first->vertex((f->second+1)%4); + Vertex_handle v2 = f->first->vertex((f->second+2)%4); + Vertex_handle v3 = f->first->vertex((f->second+3)%4); + circumradii.push_back(CGAL::make_triple(CGAL::squared_radius(v1->point(),v2->point(),v3->point()),2,CGAL::make_object(*f))); + } + for (Triangulation_hierarchy::Finite_cells_iterator cit = dt.finite_cells_begin(); + cit != dt.finite_cells_end(); cit++) { + Cell_handle ch(cit); + Vertex_handle v1 = cit->vertex(0); + Vertex_handle v2 = cit->vertex(1); + Vertex_handle v3 = cit->vertex(2); + Vertex_handle v4 = cit->vertex(3); + circumradii.push_back(CGAL::make_triple(CGAL::squared_radius(v1->point(),v2->point(),v3->point(),v4->point()),3,CGAL::make_object(ch))); + } + + std::cerr << "Sort circumradii..." << std::endl; + + std::sort(circumradii.begin(),circumradii.end(),Sort_triples()); + + std::cerr << "Filtration of size " << circumradii.size() << std::endl; + + Matrix M; + Dim_container dim_container; + + Vertex_handle v; + Edge e; + Facet f; + Cell_handle c; + + std::size_t filtration_index = 0; + std::size_t filtration_size = circumradii.size(); + + Index curr_index = 0; + + for(std::vector::const_iterator it = circumradii.begin(); + it != circumradii.end(); it++) { + + if(filtration_index % 1000 == 0) { + std::cerr << filtration_index << " of " << filtration_size + << std::endl; + } + filtration_index++; + + const CGAL::Object& obj = it->third; + + Column col; + + if(CGAL::assign(v,obj)) { + //std::cout << "Vertex " << it->second << std::endl; + dim_container.push_back(0); + v->info().set_index(curr_index); + M.push_back(col); + } + if(CGAL::assign(e,obj)) { + //std::cout << "Edge " << it->second << std::endl; + dim_container.push_back(1); + Vertex_handle v1 = e.first->vertex(e.second); + CGAL_assertion(v1->info().has_index()); + Vertex_handle v2 = e.first->vertex(e.third); + CGAL_assertion(v2->info().has_index()); + Index i1 = v1->info().index(); + Index i2 = v2->info().index(); + if(i1>i2) { + std::swap(v1,v2); + std::swap(i1,i2); + } + col.push_back(i1); + col.push_back(i2); + M.push_back(col); + set_index_of_edge(dt,e,curr_index); + } + if(CGAL::assign(f,obj)) { + //std::cout << "Facet " << it->second << std::endl; + dim_container.push_back(2); + + Index i1= f.first->info().edge_index( (f.second+1)%4, (f.second+2)%4 ); + col.push_back(i1); + Index i2= f.first->info().edge_index( (f.second+1)%4, (f.second+3)%4 ); + col.push_back(i2); + Index i3= f.first->info().edge_index( (f.second+2)%4, (f.second+3)%4 ); + col.push_back(i3); + std::sort(col.begin(),col.end()); + M.push_back(col); + set_index_of_facet(dt,f,curr_index); + + } + if(CGAL::assign(c,obj)) { + //std::cout << "Cell " << it->second << std::endl; + dim_container.push_back(3); + col.push_back(c->info().facet_index(0)); + col.push_back(c->info().facet_index(1)); + col.push_back(c->info().facet_index(2)); + col.push_back(c->info().facet_index(3)); + std::sort(col.begin(),col.end()); + M.push_back(col); + } + curr_index++; + } + + phat::boundary_matrix< phat::vector_vector > boundary_matrix; + boundary_matrix.init(M, dim_container); + boundary_matrix.save_binary("alpha_filtration.bin"); + + //phat::write(std::cout,M,dim_container); + + return 0; +} -- cgit v1.2.3