summaryrefslogtreecommitdiff
path: root/addons
diff options
context:
space:
mode:
authorulrich.bauer@gmail.com <ulrich.bauer@gmail.com@8e3bb3c2-eed4-f18f-5264-0b6c94e6926d>2013-04-29 11:41:20 +0000
committerulrich.bauer@gmail.com <ulrich.bauer@gmail.com@8e3bb3c2-eed4-f18f-5264-0b6c94e6926d>2013-04-29 11:41:20 +0000
commit89afd70e622d03bca4e4189ada72d602b7070992 (patch)
tree9e702745bf6dc2103968a5296ac5dee8b8b8f586 /addons
parent5b3c289e9a16218061a0ff609c1bbf6f32ee7d4d (diff)
alpha shapes
git-svn-id: https://phat.googlecode.com/svn/trunk@46 8e3bb3c2-eed4-f18f-5264-0b6c94e6926d
Diffstat (limited to 'addons')
-rw-r--r--addons/alpha_3.cpp317
1 files changed, 317 insertions, 0 deletions
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 <phat/boundary_matrix.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Triangulation_hierarchy_3.h>
+
+#include <CGAL/Triangulation_vertex_base_with_info_3.h>
+#include <CGAL/Triangulation_cell_base_with_info_3.h>
+
+#include <fstream>
+#include <list>
+#include <cassert>
+#include <sstream>
+
+template<typename Index_>
+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> index_;
+};
+
+
+template<typename Index_>
+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<Index> edge_index_[6];
+ boost::optional<Index> facet_index_[4];
+
+};
+
+typedef phat::dimension Dim;
+typedef std::vector<Dim> Dim_container;
+typedef phat::index Index;
+typedef phat::column Column;
+typedef std::vector<std::vector<Index> > Matrix;
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Gt;
+typedef Gt::FT FT;
+
+typedef CGAL::Triangulation_cell_base_with_info_3<Cell_info_3<Index>,Gt> Cb;
+typedef CGAL::Triangulation_vertex_base_with_info_3<Vertex_info_3<Index>,Gt> Vb;
+typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb> Vbh;
+typedef CGAL::Triangulation_data_structure_3<Vbh,Cb> Tds;
+typedef CGAL::Delaunay_triangulation_3<Gt,Tds> Triangulation_3;
+typedef CGAL::Triangulation_hierarchy_3<Triangulation_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<typename Triple>
+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<Point> 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<Point>::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<FT,int,CGAL::Object> Triple;
+
+ std::vector<Triple > 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<Triple>());
+
+ 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<Triple>::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;
+}