# HG changeset patch
# User Dmitriy Morozov <dmitriy@mrzv.org>
# Date 1511726869 28800
# Node ID 5d387c35465b8c0189723d2dee183efcf366bf03
# Parent 21b56c6c9512bbdd4dd3d99ff12c0715225536bb
Add Vanessa Robins' periodic alpha shapes 2D code
diff -r 21b56c6c9512 -r 5d387c35465b examples/alphashapes/CMakeLists.txt
--- a/examples/alphashapes/CMakeLists.txt Thu Sep 07 17:35:28 2017 -0700
+++ b/examples/alphashapes/CMakeLists.txt Sun Nov 26 12:07:49 2017 -0800
@@ -14,6 +14,7 @@
alphashapes2d
alphashapes3d-cohomology
alphashapes3d-periodic
+ alphashapes2d-periodic
#alpharadius
)
add_definitions (${CGAL_CXX_FLAGS_INIT})
diff -r 21b56c6c9512 -r 5d387c35465b examples/alphashapes/alphashapes2d-periodic.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d-periodic.cpp Sun Nov 26 12:07:49 2017 -0800
@@ -0,0 +1,120 @@
+#include <utilities/log.h>
+
+#include "alphashapes2d-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<AlphaSimplex2D> 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("error") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+#endif
+
+ //SetFrequency(GetCounter("filtration/pair"), 10000);
+ //SetTrigger(GetCounter("filtration/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;
+
+
+ // VR: The first two "points" from the input file must define the antipodal corners of the periodic domain.
+ // The periodic domain must be a square.
+ double a,b,c,d;
+ in >> a >> b >> c >> d;
+
+ Delaunay2D Dt(GT::Iso_rectangle_2(a,b,c,d));
+
+
+ while(in)
+ {
+ in >> x >> y;
+ // if (!in) break;
+ Point p(x,y);
+ 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: %i", af.size());
+
+ // Create the alpha-shape filtration
+ af.sort(AlphaSimplex2D::AlphaOrder());
+ rInfo("Filtration initialized");
+
+ Persistence p(af);
+ rInfo("Persistence initializaed");
+
+ p.pair_simplices();
+ rInfo("Simplices paired");
+
+ 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, AlphaSimplex2D::AlphaValueEvaluator()),
+ evaluate_through_map(m, AlphaSimplex2D::DimensionExtractor()));
+
+#if 0
+ std::cout << 0 << std::endl << dgms[0] << std::endl;
+ std::cout << 1 << std::endl << dgms[1] << std::endl;
+#endif
+
+ std::ofstream ofs(outfilename.c_str());
+ boost::archive::binary_oarchive oa(ofs);
+ oa << dgms;
+}
+
diff -r 21b56c6c9512 -r 5d387c35465b examples/alphashapes/alphashapes2d-periodic.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d-periodic.h Sun Nov 26 12:07:49 2017 -0800
@@ -0,0 +1,103 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ *
+ * Modified by Vanessa Robins, ANU, November 2017 to use
+ * Delaunay triangulations built with periodic boundary conditions
+ *
+ */
+
+#ifndef __ALPHASHAPES2D_H__
+#define __ALPHASHAPES2D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+// #include <CGAL/Delaunay_triangulation_2.h>
+// VR: replaced with
+#include <CGAL/Periodic_2_triangulation_traits_2.h>
+#include <CGAL/Periodic_2_Delaunay_triangulation_2.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_2<K> Delaunay2D;
+// VR: replaced with
+typedef CGAL::Periodic_2_triangulation_traits_2<K> GT;
+typedef CGAL::Periodic_2_Delaunay_triangulation_2<GT> Delaunay2D;
+
+typedef Delaunay2D::Point Point;
+typedef Delaunay2D::Vertex_handle Vertex_handle;
+typedef Delaunay2D::Face_handle Face_handle;
+
+typedef K::FT RealValue;
+
+// VR: had to change iterator typedefs from Finite_blah_iterator to Blah_iterator.
+//typedef Delaunay2D::Finite_vertices_iterator Vertex_iterator;
+//typedef Delaunay2D::Finite_edges_iterator Edge_iterator;
+//typedef Delaunay2D::Finite_faces_iterator Face_iterator;
+
+typedef Delaunay2D::Vertex_iterator Vertex_iterator;
+typedef Delaunay2D::Edge_iterator Edge_iterator;
+typedef Delaunay2D::Face_iterator Face_iterator;
+
+
+class AlphaSimplex2D: public Simplex<Vertex_handle>
+{
+ public:
+ typedef Simplex<Vertex_handle> Parent;
+ typedef std::set<AlphaSimplex2D, Parent::VertexComparison> SimplexSet;
+ typedef Parent::VertexContainer VertexSet;
+
+ public:
+ AlphaSimplex2D() {}
+ AlphaSimplex2D(const Parent& p):
+ Parent(p) {}
+ AlphaSimplex2D(const AlphaSimplex2D& s):
+ Parent(s) { attached_ = s.attached_; alpha_ = s.alpha_; }
+ AlphaSimplex2D(const Delaunay2D::Vertex& v);
+
+ AlphaSimplex2D(const Delaunay2D::Edge& e);
+ AlphaSimplex2D(const Delaunay2D::Edge& e, const SimplexSet& simplices, const Delaunay2D& Dt, const RealValue r);
+
+ AlphaSimplex2D(const Delaunay2D::Face& f);
+ AlphaSimplex2D(const Delaunay2D::Face& f, 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 AlphaSimplex2D& first, const AlphaSimplex2D& second) const; };
+
+ struct AlphaValueEvaluator
+ {
+ typedef AlphaSimplex2D first_argument_type;
+ typedef RealType result_type;
+
+ RealType operator()(const AlphaSimplex2D& s) const { return s.value(); }
+ };
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ RealValue alpha_;
+ bool attached_;
+};
+
+typedef std::vector<AlphaSimplex2D> AlphaSimplex2DVector;
+void fill_simplex_set(const Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices);
+template<class Filtration>
+void fill_complex(const Delaunay2D& Dt, Filtration& filtration);
+
+std::ostream& operator<<(std::ostream& out, const AlphaSimplex2D& s) { return s.operator<<(out); }
+
+#include "alphashapes2d-periodic.hpp"
+
+#endif // __ALPHASHAPES2D_H__
diff -r 21b56c6c9512 -r 5d387c35465b examples/alphashapes/alphashapes2d-periodic.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d-periodic.hpp Sun Nov 26 12:07:49 2017 -0800
@@ -0,0 +1,211 @@
+#include <utilities/log.h>
+#include <boost/foreach.hpp>
+
+AlphaSimplex2D::
+AlphaSimplex2D(const Delaunay2D::Vertex& v): alpha_(0), attached_(false)
+{
+ for (int i = 0; i < 3; ++i)
+ if (v.face()->vertex(i) != Vertex_handle() && v.face()->vertex(i)->point() == v.point())
+ Parent::add(v.face()->vertex(i));
+}
+
+AlphaSimplex2D::
+AlphaSimplex2D(const Delaunay2D::Edge& e): attached_(false)
+{
+ Face_handle f = e.first;
+ for (int i = 0; i < 3; ++i)
+ if (i != e.second)
+ Parent::add(f->vertex(i));
+}
+
+// AlphaSimplex2D::
+// AlphaSimplex2D(const Delaunay2D::Edge& e, const SimplexSet& simplices, const Delaunay2D& Dt): attached_(false)
+// {
+// Face_handle f = e.first;
+// for (int i = 0; i < 3; ++i)
+// if (i != e.second)
+// Parent::add(f->vertex(i));
+//
+// VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
+// const Point& p1 = (*v++)->point();
+// const Point& p2 = (*v)->point();
+//
+// Face_handle o = f->neighbor(e.second);
+// if (o == Face_handle())
+// {
+// alpha_ = CGAL::squared_radius(p1, p2);
+// return;
+// }
+// int oi = o->index(f);
+//
+// attached_ = false;
+// if (!Dt.is_infinite(f->vertex(e.second)) &&
+// CGAL::side_of_bounded_circle(p1, p2,
+// f->vertex(e.second)->point()) == CGAL::ON_BOUNDED_SIDE)
+// attached_ = true;
+// else if (!Dt.is_infinite(o->vertex(oi)) &&
+// CGAL::side_of_bounded_circle(p1, p2,
+// o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
+// attached_ = true;
+// else
+// alpha_ = CGAL::squared_radius(p1, p2);
+//
+// if (attached_)
+// {
+// if (Dt.is_infinite(f))
+// alpha_ = simplices.find(AlphaSimplex2D(*o))->alpha();
+// else if (Dt.is_infinite(o))
+// alpha_ = simplices.find(AlphaSimplex2D(*f))->alpha();
+// else
+// alpha_ = std::min(simplices.find(AlphaSimplex2D(*f))->alpha(),
+// simplices.find(AlphaSimplex2D(*o))->alpha());
+// }
+// }
+
+// VR modified for periodic DT
+AlphaSimplex2D::
+AlphaSimplex2D(const Delaunay2D::Edge& e, const SimplexSet& simplices, const Delaunay2D& Dt, const RealValue r): attached_(false)
+{
+ Face_handle f = e.first;
+ for (int i = 0; i < 3; ++i)
+ if (i != e.second)
+ Parent::add(f->vertex(i));
+
+ Face_handle o = f->neighbor(e.second);
+ int oi = o->index(f);
+ alpha_ = r;
+
+ attached_ = false;
+// First check if the circle centred on the edge contains its f-opposite vertex
+ Point fpt = Dt.point(Dt.periodic_point(e.first,e.second));
+ Delaunay2D::Segment fseg = Dt.segment(e.first,e.second);
+ if (CGAL::side_of_bounded_circle(fseg[0],fseg[1],fpt) == CGAL::ON_BOUNDED_SIDE)
+ attached_ = true;
+
+// Then check if the circle centred on the edge contains its o-opposite vertex.
+ Point opt = Dt.point(Dt.periodic_point(o,oi));
+ Delaunay2D::Segment oseg = Dt.segment(o,oi);
+ if (CGAL::side_of_bounded_circle(oseg[0],oseg[1],opt) == CGAL::ON_BOUNDED_SIDE)
+ attached_ = true;
+
+
+ if (attached_)
+ {
+ //std::cout << "found an attached edge" << std::endl;
+ alpha_ = std::min(simplices.find(AlphaSimplex2D(*f))->alpha(),
+ simplices.find(AlphaSimplex2D(*o))->alpha());
+ }
+}
+
+// AlphaSimplex2D::
+// AlphaSimplex2D(const Delaunay2D::Face& f): attached_(false)
+// {
+// for (int i = 0; i < 3; ++i)
+// Parent::add(f.vertex(i));
+// VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
+// Point p1 = (*v++)->point();
+// Point p2 = (*v++)->point();
+// Point p3 = (*v)->point();
+// alpha_ = CGAL::squared_radius(p1, p2, p3);
+// }
+
+// VR: new initialisation for periodic DT faces.
+
+AlphaSimplex2D::
+AlphaSimplex2D(const Delaunay2D::Face& f)
+{
+ for (int i = 0; i < 3; ++i)
+ Parent::add(f.vertex(i));
+}
+
+AlphaSimplex2D::
+AlphaSimplex2D(const Delaunay2D::Face& f, const RealValue r): attached_(false)
+{
+ for (int i = 0; i < 3; ++i)
+ Parent::add(f.vertex(i));
+
+ alpha_ = r;
+}
+
+
+bool
+AlphaSimplex2D::AlphaOrder::
+operator()(const AlphaSimplex2D& first, const AlphaSimplex2D& second) const
+{
+ if (first.alpha() == second.alpha())
+ return (first.dimension() < second.dimension());
+ else
+ return (first.alpha() < second.alpha());
+}
+
+std::ostream&
+AlphaSimplex2D::
+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 Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices)
+// {
+// for(Face_iterator cur = Dt.finite_faces_begin(); cur != Dt.finite_faces_end(); ++cur)
+// simplices.insert(AlphaSimplex2D(*cur));
+// rInfo("Faces inserted");
+// for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
+// simplices.insert(AlphaSimplex2D(*cur, simplices, Dt));
+// rInfo("Edges inserted");
+// for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
+// simplices.insert(AlphaSimplex2D(*cur));
+// rInfo("Vertices inserted");
+// }
+
+//VR: modified function for periodic DT
+
+void fill_simplex_set(const Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices)
+{
+ // Compute all simplices with their alpha values and attachment information
+ int counter = 0;
+ for(Face_iterator cur = Dt.faces_begin(); cur != Dt.faces_end(); ++cur){
+ //Delaunay2D::Periodic_triangle ptri = Dt.periodic_triangle(cur);
+ //Delaunay2D::Triangle tri = Dt.triangle(ptri);
+ Delaunay2D::Triangle tri = Dt.triangle(cur);
+ RealValue sqrad = CGAL::squared_radius(tri[0], tri[1], tri[2]);
+ simplices.insert(AlphaSimplex2D(*cur, sqrad));
+ counter++;
+ }
+ //std::cout << "Number of triangles inserted in filtration " << counter << std::endl;
+ rInfo("Facets inserted");
+ counter = 0;
+ for(Edge_iterator cur = Dt.edges_begin(); cur != Dt.edges_end(); ++cur){
+ //Delaunay2D::Periodic_segment pseg = Dt.periodic_segment(*cur);
+ //Delaunay2D::Segment seg = Dt.segment(pseg);
+ Delaunay2D::Segment seg = Dt.segment(cur);
+ RealValue sqrad = CGAL::squared_radius(seg[0],seg[1]);
+ simplices.insert(AlphaSimplex2D(*cur, simplices, Dt, sqrad));
+ counter++;
+ }
+ //std::cout << "Number of edges inserted in filtration " << counter << std::endl;
+ rInfo("Edges inserted");
+ for(Vertex_iterator cur = Dt.vertices_begin(); cur != Dt.vertices_end(); ++cur)
+ simplices.insert(AlphaSimplex2D(*cur));
+ rInfo("Vertices inserted");
+
+}
+
+
+
+template<class Filtration>
+void fill_complex(const Delaunay2D& Dt, Filtration& filtration)
+{
+ // Compute all simplices with their alpha values and attachment information
+ // TODO: this can be optimized; the new Filtration can act as a SimplexSet
+ AlphaSimplex2D::SimplexSet simplices;
+ fill_simplex_set(Dt, simplices);
+ BOOST_FOREACH(const AlphaSimplex2D& s, simplices)
+ filtration.push_back(s);
+}
+