summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--addons/alpha_3.cpp210
-rw-r--r--src/benchmark.cpp6
-rw-r--r--src/phat.cpp8
3 files changed, 178 insertions, 46 deletions
diff --git a/addons/alpha_3.cpp b/addons/alpha_3.cpp
index 6df251c..5ce855a 100644
--- a/addons/alpha_3.cpp
+++ b/addons/alpha_3.cpp
@@ -39,10 +39,11 @@ private:
};
-template<typename Index_>
+template<typename Index_, typename FT_>
struct Cell_info_3 {
- typedef Index_ Index;
+ typedef Index_ Index;
+ typedef FT_ FT;
Cell_info_3() {
for(std::size_t i=0; i<6; i++) {
@@ -92,8 +93,13 @@ struct Cell_info_3 {
facet_index_[i]=I;
}
+ FT filtration_value;
+ FT filtration_value_of_facet[4];
+
private:
+
+
boost::optional<Index> edge_index_[6];
boost::optional<Index> facet_index_[4];
@@ -108,29 +114,29 @@ 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_cell_base_with_info_3<Cell_info_3<Index,FT>,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 CGAL::Triangulation_hierarchy_vertex_base_3<Vb> Vbh;
+typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
+typedef CGAL::Delaunay_triangulation_3<Gt,Tds> DT_3;
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 DT_3::Vertex_handle Vertex_handle;
+typedef DT_3::Edge Edge;
+typedef DT_3::Facet Facet;
+typedef DT_3::Cell_handle Cell_handle;
-typedef Triangulation_3::Cell_circulator Cell_circulator;
+typedef DT_3::Cell_circulator Cell_circulator;
+typedef DT_3::Facet_circulator Facet_circulator;
-void set_index_of_edge(const Triangulation_3& T, const Edge& e, Index I) {
+
+void set_index_of_edge(const DT_3& dt, 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=dt.incident_cells(e);
Cell_circulator ch_start=ch;
int count=0;
do {
@@ -141,10 +147,10 @@ void set_index_of_edge(const Triangulation_3& T, const Edge& e, Index I) {
//std::cout << "Did " << count << " updates" << std::endl;
}
-void set_index_of_facet(const Triangulation_3& T, const Facet& f, Index I) {
+void set_index_of_facet(const DT_3& dt, const Facet& f, Index I) {
f.first->info().set_facet_index(f.second,I);
- Facet mf = T.mirror_facet(f);
+ Facet mf = dt.mirror_facet(f);
mf.first->info().set_facet_index(mf.second,I);
}
@@ -160,6 +166,58 @@ struct Sort_triples {
};
+// Returns whether the facet facet is attached to the cell
+// that is used to represent facet
+bool edge_attached_to(const DT_3& dt, const Edge& edge, const Facet& facet) {
+
+ if(dt.is_infinite(facet)) {
+ return false;
+ }
+
+ Vertex_handle v1 = edge.first->vertex(edge.second);
+ Vertex_handle v2 = edge.first->vertex(edge.third);
+
+ Cell_handle cell = facet.first;
+ int i1 = facet.second;
+ int i2 = cell->index(v1);
+ int i3 = cell->index(v2);
+ CGAL_assertion(i1!=i2);
+ CGAL_assertion(i1!=i3);
+ CGAL_assertion(i2!=i3);
+ int j = 0;
+
+ while(j==i1 || j==i2 || j==i3) {
+ j++;
+ }
+ // j is the index of the third point of the facet
+ Vertex_handle w = cell->vertex(j);
+
+ return CGAL::side_of_bounded_sphere(v1->point(),v2->point(),w->point())==CGAL::ON_BOUNDED_SIDE;
+
+}
+
+
+// Returns whether the facet facet is attached to the cell
+// that is used to represent facet
+bool triangle_attached_to(const DT_3& dt, const Facet& facet) {
+ Cell_handle cell = facet.first;
+ int index = facet.second;
+
+ if(dt.is_infinite(cell)) {
+ return false;
+ }
+
+ Vertex_handle v1 = cell->vertex((index+1)%4);
+ Vertex_handle v2 = cell->vertex((index+2)%4);
+ Vertex_handle v3 = cell->vertex((index+3)%4);
+
+ Vertex_handle w = cell->vertex(facet.second);
+
+ return CGAL::side_of_bounded_sphere(v1->point(),v2->point(),v3->point(),w->point())==CGAL::ON_BOUNDED_SIDE;
+}
+
+
+
int main(int argc, char** argv)
{
@@ -183,8 +241,10 @@ int main(int argc, char** argv)
}
}
+ std::cerr << "Read " << lp.size() << " points" << std::endl;
+
std::cerr << "Compute Delaunay triangulation..." << std::endl;
- Triangulation_hierarchy dt(lp.begin(),lp.end());
+ DT_3 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
@@ -195,33 +255,105 @@ int main(int argc, char** argv)
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();
+ for (DT_3::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)));
+ FT circumradius = CGAL::squared_radius(v1->point(),v2->point(),v3->point(),v4->point());
+ circumradii.push_back(CGAL::make_triple(circumradius,3,CGAL::make_object(ch)));
+ ch->info().filtration_value = circumradius;
+ }
+
+ for (DT_3::Finite_facets_iterator f = dt.finite_facets_begin();
+ f != dt.finite_facets_end(); f++) {
+
+ // Check Gabrielness (attachment) of the triangle to its incident tetrahedra
+ bool is_attached = false;
+ FT filtration_value;
+
+ {
+ Cell_handle cell = f->first;
+ int index = f->second;
+ if( triangle_attached_to(dt,*f) ) {
+ filtration_value = cell->info().filtration_value;
+ is_attached=true;
+ }
+ }
+
+ {
+ if(! is_attached) {
+ Facet g = dt.mirror_facet(*f);
+ Cell_handle cell = g.first;
+ int index = g.second;
+ if( triangle_attached_to(dt,g) ) {
+ filtration_value = cell->info().filtration_value;
+ is_attached=true;
+ }
+ }
+ }
+ if(! is_attached) {
+ Cell_handle cell = f->first;
+ int index = f->second;
+ filtration_value = CGAL::squared_radius(cell->vertex((index+1)%4)->point(),
+ cell->vertex((index+2)%4)->point(),
+ cell->vertex((index+3)%4)->point());
+ }
+ circumradii.push_back(CGAL::make_triple(filtration_value,2,CGAL::make_object(*f)));
+ // Store filtration values in info object for further use
+ f->first->info().filtration_value_of_facet[f->second] = filtration_value;
+ Facet g = dt.mirror_facet(*f);
+ g.first->info().filtration_value_of_facet[g.second] = filtration_value;
}
+
+
+ for (DT_3::Finite_edges_iterator edge = dt.finite_edges_begin();
+ edge != dt.finite_edges_end(); edge++) {
+
+ // Check Gabrielness of the edge
+ bool is_attached = false;
+ FT filtration_value;
+
+ Facet_circulator f_start = dt.incident_facets(*edge);
+ Facet_circulator f_it = f_start;
+ f_it++;
+ while( f_it!=f_start ) {
+ if(edge_attached_to(dt,*edge,*f_it)) {
+ Cell_handle cell = f_it->first;
+ int index = f_it->second;
+ FT new_filtration_value = cell->info().filtration_value_of_facet[index];
+ // Now, maybe the edge is already attached to another triangle (can this really happen?)
+ if(is_attached) {
+ filtration_value = CGAL::min(filtration_value,new_filtration_value);
+ } else {
+ filtration_value = new_filtration_value;
+ }
+ is_attached=true;
+ }
+ f_it++;
+ }
+
+
+ if(! is_attached) {
+ Vertex_handle v1 = edge->first->vertex(edge->second);
+ Vertex_handle v2 = edge->first->vertex(edge->third);
+ filtration_value = CGAL::squared_radius(v1->point(),v2->point());
+ }
+
+ circumradii.push_back(CGAL::make_triple(filtration_value,1,CGAL::make_object(*edge)));
+ }
+
+
+
+
+ for (DT_3::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)));
+ }
+
std::cerr << "Sort circumradii..." << std::endl;
@@ -308,7 +440,7 @@ int main(int argc, char** argv)
curr_index++;
}
- boundary_matrix.save_binary("alpha_filtration.bin");
+ boundary_matrix.save_ascii("alpha_filtration.bin");
//phat::write(std::cout,M,dim_container);
diff --git a/src/benchmark.cpp b/src/benchmark.cpp
index 0f2a38e..b552d33 100644
--- a/src/benchmark.cpp
+++ b/src/benchmark.cpp
@@ -148,14 +148,14 @@ void benchmark( std::string input_filename, bool use_binary, Ansatz_type ansatz
dualize( matrix );
double dualization_time = omp_get_wtime() - dualization_timer;
double dualization_time_rounded = floor( dualization_time * 10.0 + 0.5 ) / 10.0;
- std::cout << " Dualization time: " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << dualization_time_rounded <<"s,";
+ std::cout << " Dualization time: " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << dualization_time_rounded <<"s,";
reduction_timer = omp_get_wtime();
reduction_algorithm( matrix );
}
double running_time = omp_get_wtime() - reduction_timer;
double running_time_rounded = floor( running_time * 10.0 + 0.5 ) / 10.0;
- std::cout << " Reduction time: " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << running_time_rounded <<"s" << std::endl;
+ std::cout << " Reduction time: " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << running_time_rounded <<"s" << std::endl;
}
template<typename Representation, typename Algorithm>
@@ -186,7 +186,7 @@ void benchmark_latex( std::string input_filename, bool use_binary, Ansatz_type a
//double running_time = omp_get_wtime() - reduction_timer + dualization_time;
double running_time = omp_get_wtime( ) - reduction_timer;
double running_time_rounded = floor( running_time * 10.0 + 0.5 ) / 10.0;
- std::cout << " && " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << std::setw( 12 ) << running_time_rounded << std::setw( 1 );
+ std::cout << " && " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << std::setw( 12 ) << running_time_rounded << std::setw( 1 );
}
#define COMPUTE(Representation) \
diff --git a/src/phat.cpp b/src/phat.cpp
index b099b59..3770cf3 100644
--- a/src/phat.cpp
+++ b/src/phat.cpp
@@ -109,7 +109,7 @@ void compute_pairing( std::string input_filename, std::string output_filename, b
}
double read_time = omp_get_wtime() - read_timer;
double read_time_rounded = floor( read_time * 10.0 + 0.5 ) / 10.0;
- LOG( "Reading input file took " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << read_time_rounded <<"s" )
+ LOG( "Reading input file took " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << read_time_rounded <<"s" )
if( !read_successful ) {
std::cerr << "Error opening file " << input_filename << std::endl;
@@ -124,7 +124,7 @@ void compute_pairing( std::string input_filename, std::string output_filename, b
phat::dualize ( matrix );
double dualize_time = omp_get_wtime() - dualize_timer;
double dualize_time_rounded = floor( dualize_time * 10.0 + 0.5 ) / 10.0;
- LOG( "Dualizing took " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << dualize_time_rounded <<"s" )
+ LOG( "Dualizing took " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << dualize_time_rounded <<"s" )
}
double pairs_timer = omp_get_wtime();
@@ -133,7 +133,7 @@ void compute_pairing( std::string input_filename, std::string output_filename, b
phat::compute_persistence_pairs < Algorithm > ( pairs, matrix );
double pairs_time = omp_get_wtime() - pairs_timer;
double pairs_time_rounded = floor( pairs_time * 10.0 + 0.5 ) / 10.0;
- LOG( "Computing persistence pairs took " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << pairs_time_rounded <<"s" )
+ LOG( "Computing persistence pairs took " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << pairs_time_rounded <<"s" )
if( dualize ) dualize_persistence_pairs( pairs, num_cols );
@@ -148,7 +148,7 @@ void compute_pairing( std::string input_filename, std::string output_filename, b
}
double write_time = omp_get_wtime() - write_timer;
double write_time_rounded = floor( write_time * 10.0 + 0.5 ) / 10.0;
- LOG( "Writing output file took " << setiosflags( std::ios::fixed ) << setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << write_time_rounded <<"s" )
+ LOG( "Writing output file took " << std::setiosflags( std::ios::fixed ) << std::setiosflags( std::ios::showpoint ) << std::setprecision( 1 ) << write_time_rounded <<"s" )
}
#define COMPUTE_PAIRING(Representation) \
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-18 06:40:42 +0000