Add Vanessa Robins' periodic alpha shapes code
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu, 07 Sep 2017 17:35:28 -0700
changeset 290 21b56c6c9512
parent 289 f7dbbfdac250
child 291 5d387c35465b
Add Vanessa Robins' periodic alpha shapes code
examples/alphashapes/CMakeLists.txt
examples/alphashapes/alphashapes3d-periodic.cpp
examples/alphashapes/alphashapes3d-periodic.h
examples/alphashapes/alphashapes3d-periodic.hpp
--- a/examples/alphashapes/CMakeLists.txt	Fri Jul 14 14:07:05 2017 -0700
+++ b/examples/alphashapes/CMakeLists.txt	Thu Sep 07 17:35:28 2017 -0700
@@ -13,6 +13,7 @@
     set                         (targets                        alphashapes3d
                                                                 alphashapes2d
                                                                 alphashapes3d-cohomology
+                                                                alphashapes3d-periodic
                                                                 #alpharadius
                                 )
     add_definitions             (${CGAL_CXX_FLAGS_INIT})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d-periodic.cpp	Thu Sep 07 17:35:28 2017 -0700
@@ -0,0 +1,120 @@
+#include <utilities/log.h>
+#include <utilities/timer.h>
+
+#include "alphashapes3d-periodic.h"
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/persistence-diagram.h>
+#include <iostream>
+
+#include <fstream>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/serialization/map.hpp>
+
+#include <boost/program_options.hpp>
+
+
+typedef Filtration<AlphaSimplex3D>              AlphaFiltration;
+typedef StaticPersistence<>                     Persistence;
+typedef PersistenceDiagram<>                    PDgm;
+
+namespace po = boost::program_options;
+
+int main(int argc, char** argv) 
+{
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+    stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+    //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/persistence") );
+    //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/chain") );
+#endif
+
+    // SetFrequency(GetCounter("persistence/pair"), 10000);
+    // SetTrigger(GetCounter("persistence/pair"), GetCounter(""));
+
+    std::string     infilename, outfilename;
+
+    po::options_description hidden("Hidden options");
+    hidden.add_options()
+        ("input-file",   po::value<std::string>(&infilename),     "Point set whose alpha shape filtration and persistence we want to compute")
+        ("output-file",  po::value<std::string>(&outfilename),    "Where to write the collection of persistence diagrams");
+
+    po::positional_options_description pos;
+    pos.add("input-file", 1);
+    pos.add("output-file", 2);
+    
+    po::options_description all; all.add(hidden);
+
+    po::variables_map vm;
+    po::store(po::command_line_parser(argc, argv).
+                  options(all).positional(pos).run(), vm);
+    po::notify(vm);
+
+    if (!vm.count("input-file") || !vm.count("output-file"))
+    { 
+        std::cout << "Usage: " << argv[0] << " input-file output-file" << std::endl;
+        std::cout << hidden << std::endl; 
+        return 1; 
+    }
+
+
+    // Read in the point set and compute its Delaunay triangulation
+    std::ifstream in(infilename.c_str());
+    double x,y,z;
+    
+    // VR: The first two "points" from the input file must define the antipodal corners of the periodic domain.
+    // The periodic domain must be a cube.
+    double a,b,c,d,e,f;
+    in >> a >> b >> c >> d >> e >> f;
+
+    Delaunay3D Dt(GT::Iso_cuboid_3(a,b,c,d,e,f));
+    
+    while(in)
+    {
+        in >> x >> y >> z;
+        Point p(x,y,z);
+        Dt.insert(p);
+    }
+    // VR: Switch to 1-sheeted cover if possible
+    if (Dt.is_triangulation_in_1_sheet())
+        Dt.convert_to_1_sheeted_covering();
+    else
+        std::cout << "Periodic Triangulation not in one sheet. There may be problems ahead" << std::endl;
+    
+    rInfo("Delaunay triangulation computed");
+    
+   
+    AlphaFiltration  af;
+    fill_complex(Dt, af);
+    rInfo("Simplices: %d", af.size());
+
+    // Create the alpha-shape filtration
+    af.sort(AlphaSimplex3D::AlphaOrder());
+    rInfo("Filtration initialized");
+    
+    Persistence p(af);
+    rInfo("Persistence initialized");
+
+    Timer persistence_timer; persistence_timer.start();
+    p.pair_simplices();
+    persistence_timer.stop();
+    rInfo("Simplices paired");
+    persistence_timer.check("Persistence timer");
+
+    Persistence::SimplexMap<AlphaFiltration>    m       = p.make_simplex_map(af);
+    std::map<Dimension, PDgm> dgms;
+    init_diagrams(dgms, p.begin(), p.end(), 
+                  evaluate_through_map(m, AlphaSimplex3D::AlphaValueEvaluator()), 
+                  evaluate_through_map(m, AlphaSimplex3D::DimensionExtractor()));
+#if 0
+    std::cout << 0 << std::endl << dgms[0] << std::endl;
+    std::cout << 1 << std::endl << dgms[1] << std::endl;
+    std::cout << 2 << std::endl << dgms[2] << std::endl;
+#endif
+
+    std::ofstream ofs(outfilename.c_str());
+    boost::archive::binary_oarchive oa(ofs);
+    oa << dgms;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d-periodic.h	Thu Sep 07 17:35:28 2017 -0700
@@ -0,0 +1,103 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2006
+ * 
+ * Modified by Vanessa Robins August 2017 to use
+ * Delaunay triangulations built with periodic boundary conditions
+ *
+ */
+
+#ifndef __ALPHASHAPES3D_H__
+#define __ALPHASHAPES3D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+// #include <CGAL/Delaunay_triangulation_3.h>
+// VR: replaced with
+#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
+
+
+#include <topology/simplex.h>
+#include <utilities/types.h>
+
+#include <vector>
+#include <set>
+#include <iostream>
+
+struct K: CGAL::Exact_predicates_exact_constructions_kernel {};
+
+//typedef CGAL::Delaunay_triangulation_3<K>           Delaunay3D;
+// VR: replaced with
+typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> GT;
+typedef CGAL::Periodic_3_Delaunay_triangulation_3<GT>   Delaunay3D;
+
+typedef Delaunay3D::Point                           Point;
+typedef Delaunay3D::Offset                          Offset;
+//typedef Delaunay3D::Periodic_point                  Periodic_point;
+typedef Delaunay3D::Vertex_handle                   Vertex_handle;
+typedef Delaunay3D::Cell_handle                     Cell_handle;
+typedef K::FT                                       RealValue;
+
+typedef Delaunay3D::Vertex_iterator         Vertex_iterator;
+typedef Delaunay3D::Edge_iterator           Edge_iterator;
+typedef Delaunay3D::Facet_iterator          Facet_iterator;
+typedef Delaunay3D::Cell_iterator           Cell_iterator;
+typedef Delaunay3D::Facet_circulator        Facet_circulator;
+
+
+class AlphaSimplex3D: public Simplex<Vertex_handle>
+{
+	public:
+		typedef		Simplex<Vertex_handle>					            Parent;
+		typedef 	std::set<AlphaSimplex3D, Parent::VertexComparison>  SimplexSet;
+		typedef		Parent::VertexContainer								VertexSet;
+
+    public:
+									AlphaSimplex3D()					{}
+									AlphaSimplex3D(const Parent& p): 
+											Parent(p) 					{}
+									AlphaSimplex3D(const AlphaSimplex3D& s): 
+											Parent(s) 					{ attached_ = s.attached_; alpha_ = s.alpha_; }
+	    							AlphaSimplex3D(const Delaunay3D::Vertex& v);
+		
+								    AlphaSimplex3D(const Delaunay3D::Edge& e);
+								    AlphaSimplex3D(const Delaunay3D::Edge& e, const SimplexSet& simplices, const Delaunay3D& Dt, Facet_circulator facet_bg, const RealValue r);
+		
+								    AlphaSimplex3D(const Delaunay3D::Facet& f);
+								    AlphaSimplex3D(const Delaunay3D::Facet& f, const SimplexSet& simplices, const Delaunay3D& Dt, const RealValue r);
+	    
+									AlphaSimplex3D(const Delaunay3D::Cell& c);
+                                    AlphaSimplex3D(const Delaunay3D::Cell& c, const RealValue r);
+	    
+		RealType					value() const						{ return CGAL::to_double(alpha_); }
+		RealValue					alpha() const						{ return alpha_; }
+		bool						attached() const					{ return attached_; }
+
+		// Ordering
+		struct AlphaOrder
+		{ bool operator()(const AlphaSimplex3D& first, const AlphaSimplex3D& second) const; };
+	
+        struct AlphaValueEvaluator
+        { 
+            typedef                 AlphaSimplex3D                                  first_argument_type;
+            typedef                 RealType                                        result_type;
+
+            RealType                operator()(const AlphaSimplex3D& s) const       { return s.value(); }
+        };
+
+		std::ostream& 				operator<<(std::ostream& out) const;
+		
+	private:
+		RealValue 					alpha_;
+		bool 						attached_;
+};
+
+void 				fill_simplex_set(const Delaunay3D& Dt, AlphaSimplex3D::SimplexSet& simplices);
+template<class Filtration>
+void 				fill_complex(const Delaunay3D& Dt,     Filtration& filtration);
+
+std::ostream& 		operator<<(std::ostream& out, const AlphaSimplex3D& s)	{ return s.operator<<(out); }
+
+#include "alphashapes3d-periodic.hpp"
+
+#endif // __ALPHASHAPES3D_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d-periodic.hpp	Thu Sep 07 17:35:28 2017 -0700
@@ -0,0 +1,258 @@
+#include <utilities/log.h>
+#include <boost/foreach.hpp>
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Vertex& v): alpha_(0), attached_(false)
+{
+    for (int i = 0; i < 4; ++i)
+        if (v.cell()->vertex(i)->point() == v.point())
+            Parent::add(v.cell()->vertex(i));
+}
+// VR:  I'm not sure why there is the extra step of checking the cell that the vertex belongs to
+// It's possible that in the periodic version that the point locations might be different, i.e. differ by an offset.
+// Or perhaps more problematically, that the same point might appear twice in a single cell (but with different offsets)
+
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Edge& e)
+{
+    Cell_handle c = e.first;
+    Parent::add(c->vertex(e.second));
+    Parent::add(c->vertex(e.third));
+}
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Edge& e, const SimplexSet& simplices, const Delaunay3D& Dt, Facet_circulator facet_bg,  const RealValue r)
+{
+    Cell_handle c = e.first;
+    Parent::add(c->vertex(e.second));
+    Parent::add(c->vertex(e.third));
+
+    Facet_circulator cur = facet_bg;
+    // VR: Don't need to worry about infinite vertex in a periodic triangulation
+    //    while (Dt.is_infinite(*cur))    ++cur;
+    SimplexSet::const_iterator cur_iter = simplices.find(AlphaSimplex3D(*cur));
+    RealValue min = cur_iter->alpha();
+
+    //VR: I've rewritten this function substantially to make use of the CGAL is_Gabriel function
+    if (Dt.is_Gabriel(e)){
+        attached_ = false;
+        alpha_ = r;
+    }
+    else {
+        attached_ = true;
+        // VR: now we just have to figure out the smallest circumradius of the facets adjacent to the given edge
+        if (facet_bg != 0) do
+        {
+            cur_iter = simplices.find(AlphaSimplex3D(*cur));
+            RealValue val = cur_iter->alpha();
+            if (val < min)
+                min = val;
+            ++cur;
+        } while (cur != facet_bg);
+        alpha_ = min;
+    }
+}
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Facet& f)
+{
+    Cell_handle c = f.first;
+    for (int i = 0; i < 4; ++i)
+        if (i != f.second)
+            Parent::add(c->vertex(i));
+}
+
+/*
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Facet& f, const SimplexSet& simplices, const Delaunay3D& Dt)
+{
+    Cell_handle c = f.first;
+    for (int i = 0; i < 4; ++i)
+        if (i != f.second)
+            Parent::add(c->vertex(i));
+
+    Cell_handle o = c->neighbor(f.second);
+    int oi = o->index(c);
+
+    VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
+    const Point& p1 = (*v)->point();
+    const Offset& o1 = (*v++)->offset();
+    const Point& p2 = (*v)->point();
+    const Offset& o2 = (*v++)->offset();
+    const Point& p3 = (*v)->point();
+    const Offset& o3 = (*v)->offset();
+
+    //VR: I've rewritten this function substantially to make use of the CGAL is_Gabriel function
+    if (Dt.is_Gabriel(f)){
+        attached_ = false;
+        //alpha_ = CGAL::squared_radius(p1+o1, p2+o2, p3+o3);
+        alpha_ = CGAL::squared_radius(p1, p2, p3);
+    }
+    else {
+        attached_ = true;
+        alpha_ = std::min(simplices.find(AlphaSimplex3D(*c))->alpha(),
+                          simplices.find(AlphaSimplex3D(*o))->alpha());
+    }
+}
+*/
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Facet& f, const SimplexSet& simplices, const Delaunay3D& Dt, const RealValue r)
+{
+    Cell_handle c = f.first;
+    for (int i = 0; i < 4; ++i)
+        if (i != f.second)
+            Parent::add(c->vertex(i));
+    
+    Cell_handle o = c->neighbor(f.second);
+    int oi = o->index(c);
+    
+    //VR: I've rewritten this function substantially to make use of the CGAL is_Gabriel function
+    if (Dt.is_Gabriel(f)){
+        attached_ = false;
+        //alpha_ = CGAL::squared_radius(p1+o1, p2+o2, p3+o3);
+        alpha_ = r;
+    }
+    else {
+        attached_ = true;
+        alpha_ = std::min(simplices.find(AlphaSimplex3D(*c))->alpha(),
+                          simplices.find(AlphaSimplex3D(*o))->alpha());
+    }
+}
+
+
+/*
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Cell& c): attached_(false)
+{
+    for (int i = 0; i < 4; ++i)
+        Parent::add(c.vertex(i));
+        
+    VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
+    Point p1 = (*v)->point();
+    Offset o1 = (*v++)->offset();
+    Point p2 = (*v)->point();
+    Offset o2 = (*v++)->offset();
+    Point p3 = (*v)->point();
+    Offset o3 = (*v++)->offset();
+    Point p4 = (*v)->point();
+    Offset o4 = (*v)->offset();
+    alpha_ = CGAL::squared_radius(p1+o1, p2+o2, p3+o3, p4+o4);
+}
+*/
+
+/*
+ AlphaSimplex3D::
+ AlphaSimplex3D(const Delaunay3D::Cell& c): attached_(false)
+ {
+ for (int i = 0; i < 4; ++i)
+ Parent::add(c.vertex(i));
+ 
+ Point p1 = c.vertex(0)->point();
+ int o1 = c.offset(0);
+ Point p2 = c.vertex(1)->point();
+ int o2 = c.offset(1);
+ Point p3 = c.vertex(2)->point();
+ int o3 = c.offset(2);
+ Point p4 = c.vertex(3)->point();
+ int o4 = c.offset(3);
+ // can't do the following because the offsets aren't the right type
+ //alpha_ = CGAL::squared_radius(p1+o1, p2+o2, p3+o3, p4+o4);
+ }
+*/
+
+
+/* Okay so here's the problem with the above code:
+ The vertex base class has offsets defined as triplets of integers (using the Offset type) but
+ the cell base class returns an offset defined as a single integer that is manipulated in a fancy way (it seems). 
+ 
+ If I access the vertex and its offset (as commented out above) I get a null offset, because we're in the 1-sheeted cover, and the vertex is always inside the original domain.
+ To get the proper geometry of the cell, I should access the vertex offsets stored with the cell.  
+ I can do this but the offsets are returned as a single integer.   So I guess I would need to insert the CGAL code to convert these integers into triples.
+ 
+ I also tried to use the higher-level geometric access functions but wasn't accessing them in the right way and got a seg fault.
+ The problem here is that I pass a Cell& to the AlphaSimplex3D function, but then need to use a Cell_handle to get the info I want. I don't know how to convert between the two. 
+ 
+ Hooray! Problem fixed by putting the geometric calculations into the fill_simplex_set function instead of the AlphaSimplex3D creation function.
+*/
+
+
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Cell& c)
+{
+    for (int i = 0; i < 4; ++i)
+        Parent::add(c.vertex(i));
+}
+
+AlphaSimplex3D::
+AlphaSimplex3D(const Delaunay3D::Cell& c, const RealValue r): attached_(false)
+{
+    for (int i = 0; i < 4; ++i)
+        Parent::add(c.vertex(i));
+    
+    alpha_ = r;
+}
+
+
+
+bool
+AlphaSimplex3D::AlphaOrder::
+operator()(const AlphaSimplex3D& first, const AlphaSimplex3D& second) const
+{
+    if (first.alpha() == second.alpha())
+        return (first.dimension() < second.dimension());
+    else
+        return (first.alpha() < second.alpha());
+}
+
+std::ostream&
+AlphaSimplex3D::
+operator<<(std::ostream& out) const
+{
+    for (VertexSet::const_iterator cur = Parent::vertices().begin(); cur != Parent::vertices().end(); ++cur)
+        out << **cur << ", ";
+    out << "value = " << value();
+
+    return out;
+}
+
+void fill_simplex_set(const Delaunay3D& Dt, AlphaSimplex3D::SimplexSet& simplices)
+{
+    // Compute all simplices with their alpha values and attachment information
+    for(Cell_iterator cur = Dt.cells_begin(); cur != Dt.cells_end(); ++cur){
+        Delaunay3D::Periodic_tetrahedron ptet = Dt.periodic_tetrahedron(cur);
+        Delaunay3D::Tetrahedron tet = Dt.tetrahedron(ptet);
+        RealValue sqrad = CGAL::squared_radius(tet[0], tet[1], tet[2], tet[3]);
+        simplices.insert(AlphaSimplex3D(*cur, sqrad));
+    }
+    rInfo("Cells inserted");
+    for(Facet_iterator cur = Dt.facets_begin(); cur != Dt.facets_end(); ++cur){
+        Delaunay3D::Periodic_triangle ptri = Dt.periodic_triangle(*cur);
+        Delaunay3D::Triangle tri = Dt.triangle(ptri);
+        RealValue sqrad = CGAL::squared_radius(tri[0], tri[1], tri[2]);
+        simplices.insert(AlphaSimplex3D(*cur, simplices, Dt, sqrad));
+    }
+    rInfo("Facets inserted");
+    for(Edge_iterator cur = Dt.edges_begin(); cur != Dt.edges_end(); ++cur){
+        Delaunay3D::Periodic_segment pseg = Dt.periodic_segment(*cur);
+        Delaunay3D::Segment seg = Dt.segment(pseg);
+        RealValue sqrad = CGAL::squared_radius(seg[0],seg[1]);
+        simplices.insert(AlphaSimplex3D(*cur, simplices, Dt, Dt.incident_facets(*cur), sqrad));
+    }
+    rInfo("Edges inserted");
+    for(Vertex_iterator cur = Dt.vertices_begin(); cur != Dt.vertices_end(); ++cur)
+        simplices.insert(AlphaSimplex3D(*cur));
+    rInfo("Vertices inserted");
+}
+
+template<class Filtration>
+void fill_complex(const Delaunay3D& Dt, Filtration& filtration)
+{
+    AlphaSimplex3D::SimplexSet simplices;
+    fill_simplex_set(Dt, simplices);
+    BOOST_FOREACH(const AlphaSimplex3D& s, simplices)
+        filtration.push_back(s);
+}
+