--- 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);
+}
+