--- a/examples/CMakeLists.txt Wed Feb 07 14:00:39 2007 -0500
+++ b/examples/CMakeLists.txt Wed Aug 15 11:25:35 2007 -0400
@@ -1,3 +1,4 @@
add_subdirectory (alphashapes)
+add_subdirectory (ar-vineyard)
add_subdirectory (grid)
add_subdirectory (triangle)
--- a/examples/alphashapes/alphashapes3d.cpp Wed Feb 07 14:00:39 2007 -0500
+++ b/examples/alphashapes/alphashapes3d.cpp Wed Aug 15 11:25:35 2007 -0400
@@ -7,7 +7,6 @@
#include <fstream>
-typedef std::vector<AlphaSimplex3D> SimplexVector;
typedef Filtration<AlphaSimplex3D> AlphaFiltration;
int main(int argc, char** argv)
@@ -30,7 +29,7 @@
// Create the alpha-shape filtration
AlphaFiltration af;
- for (SimplexVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
+ for (AlphaSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
af.append(*cur);
af.fill_simplex_index_map();
af.pair_simplices(af.begin(), af.end());
--- a/examples/alphashapes/alphashapes3d.h Wed Feb 07 14:00:39 2007 -0500
+++ b/examples/alphashapes/alphashapes3d.h Wed Aug 15 11:25:35 2007 -0400
@@ -59,7 +59,7 @@
AlphaSimplex3D(const Cell& c);
- virtual RealType value() const { return CGAL::to_double(alpha_); }
+ RealType value() const { return CGAL::to_double(alpha_); }
RealValue alpha() const { return alpha_; }
bool attached() const { return attached_; }
Cycle boundary() const;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/CMakeLists.txt Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,7 @@
+set (targets
+ ar-vineyard)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp ${external_sources})
+ target_link_libraries (${t} ${libraries} ${cgal_libraries})
+endforeach (t)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-simplex3d.h Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,103 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ */
+
+#ifndef __ARSIMPLEX3D_H__
+#define __ARSIMPLEX3D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+//#include <CGAL/Kernel/global_functions_3.h> // FIXME: what do we include for circumcenter?
+#include <CGAL/squared_distance_3.h>
+
+#include <simplex.h>
+#include <types.h>
+
+#include <vector>
+#include <set>
+#include <iostream>
+
+struct K: CGAL::Exact_predicates_exact_constructions_kernel {};
+
+typedef CGAL::Delaunay_triangulation_3<K> Delaunay;
+typedef Delaunay::Point Point;
+typedef Delaunay::Vertex Vertex;
+typedef Delaunay::Vertex_handle Vertex_handle;
+typedef Delaunay::Edge Edge;
+typedef Delaunay::Facet Facet;
+typedef Delaunay::Cell Cell;
+typedef Delaunay::Cell_handle Cell_handle;
+typedef K::FT RealValue;
+
+typedef Delaunay::Finite_vertices_iterator Vertex_iterator;
+typedef Delaunay::Finite_edges_iterator Edge_iterator;
+typedef Delaunay::Finite_facets_iterator Facet_iterator;
+typedef Delaunay::Finite_cells_iterator Cell_iterator;
+typedef Delaunay::Facet_circulator Facet_circulator;
+
+
+class ARSimplex3D: public SimplexWithVertices<Vertex_handle>
+{
+ public:
+ typedef std::map<ARSimplex3D, RealValue> SimplexPhiMap;
+ typedef SimplexWithVertices<Vertex_handle> Parent;
+ typedef Parent::VertexContainer VertexSet;
+ typedef std::list<ARSimplex3D> Cycle;
+
+ public:
+ ARSimplex3D() {}
+ ARSimplex3D(const Parent& p):
+ Parent(p) {}
+ ARSimplex3D(const ARSimplex3D& s);
+
+ ARSimplex3D(const Parent::Vertex& v);
+ ARSimplex3D(const ::Vertex& v, const Point& z);
+
+ ARSimplex3D(const Edge& e);
+ ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices,
+ Facet_circulator facet_bg);
+
+ ARSimplex3D(const Facet& f);
+ ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices);
+
+ ARSimplex3D(const Cell& c, const Point& z);
+
+ RealType value() const { return CGAL::to_double(alpha_); }
+ RealValue alpha() const { return alpha_; }
+ RealValue phi_const() const { return phi_const_; }
+ RealValue rho() const { return rho_; }
+ RealValue s() const { return s_; }
+ RealValue v() const { return v_; }
+
+ bool attached() const { return attached_; }
+ Cycle boundary() const;
+
+ void set_phi_const(RealValue phi_const) { phi_const_ = phi_const; }
+
+ // Ordering
+ struct AlphaOrder
+ { bool operator()(const ARSimplex3D& first, const ARSimplex3D& second) const; };
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ RealValue alpha_; // alpha_ is the squared radius of the smallest _empty_ circumsphere
+ RealValue rho_; // rho_ is the squared radius of the smallest circumsphere
+ RealValue s_; // s_ is the squared distance from z to the affine hull of the simplex
+ RealValue v_; // v_ is the squared distance from z to the affine hull of the dual Voronoi cell
+ RealValue phi_const_; // see LHI paper, Appendix 1
+ bool attached_;
+
+};
+
+typedef std::vector<ARSimplex3D> ARSimplex3DVector;
+
+void update_simplex_phi_map(const ARSimplex3D& s, ARSimplex3D::SimplexPhiMap& simplices);
+void fill_alpha_order(const Delaunay& Dt, const Point& z, ARSimplex3DVector& alpha_order);
+
+std::ostream& operator<<(std::ostream& out, const ARSimplex3D& s) { return s.operator<<(out); }
+
+#include "ar-simplex3d.hpp"
+
+#endif // __ARSIMPLEX3D_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,253 @@
+ARSimplex3D::
+ARSimplex3D(const ARSimplex3D& s): Parent(s)
+{
+ attached_ = s.attached_;
+ alpha_ = s.alpha_;
+ rho_ = s.rho_;
+ s_ = s.s_;
+ v_ = s.v_;
+ phi_const_ = s.phi_const_;
+}
+
+ARSimplex3D::
+ARSimplex3D(const Parent::Vertex& v)
+{
+ Parent::add(v);
+}
+
+ARSimplex3D::
+ARSimplex3D(const ::Vertex& v, const Point& z): alpha_(0), rho_(0), v_(0), attached_(false)
+{
+ for (int i = 0; i < 4; ++i)
+ if (v.cell()->vertex(i)->point() == v.point())
+ Parent::add(v.cell()->vertex(i));
+ s_ = CGAL::squared_distance((*Parent::vertices().begin())->point(), z);
+
+ // phi_const_ will be set by an edge
+}
+
+ARSimplex3D::
+ARSimplex3D(const Edge& e)
+{
+ Cell_handle c = e.first;
+ Parent::add(c->vertex(e.second));
+ Parent::add(c->vertex(e.third));
+}
+
+ARSimplex3D::
+ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, Facet_circulator facet_bg)
+{
+ Cell_handle c = e.first;
+ Parent::add(c->vertex(e.second));
+ Parent::add(c->vertex(e.third));
+
+ Facet_circulator cur = facet_bg;
+ SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
+ while (cur_iter == simplices.end())
+ {
+ ++cur;
+ cur_iter = simplices.find(ARSimplex3D(*cur));
+ }
+ RealValue min = cur_iter->first.alpha();
+ RealValue phi_const_min = cur_iter->first.phi_const();
+
+ VertexSet::const_iterator v2 = Parent::vertices().begin();
+ VertexSet::const_iterator v1 = v2++;
+ const Point& p1 = (*v1)->point();
+ const Point& p2 = (*v2)->point();
+ attached_ = false;
+ rho_ = CGAL::squared_radius(p1, p2);
+
+ if (facet_bg != 0) do
+ {
+ int i0 = (*cur).first->index(*v1);
+ int i1 = (*cur).first->index(*v2);
+ int i = 6 - i0 - i1 - (*cur).second;
+ Point p3 = (*cur).first->vertex(i)->point();
+
+ cur_iter = simplices.find(ARSimplex3D(*cur));
+ if (cur_iter == simplices.end()) // cur is infinite
+ {
+ ++cur; continue;
+ // FIXME: what do we do with infinite cofaces (i.e., what
+ // phi_const does a simplex get if its dual Voronoi cell is
+ // infinite?
+ }
+
+ if (CGAL::side_of_bounded_sphere(p1, p2, p3) == CGAL::ON_BOUNDED_SIDE)
+ attached_ = true;
+
+ RealValue val = cur_iter->first.alpha();
+ if (val < min) min = val;
+ RealValue phi_const_val = cur_iter->first.phi_const();
+ if (phi_const_val < phi_const_min) phi_const_min = phi_const_val;
+
+ ++cur;
+ } while (cur != facet_bg);
+
+ if (attached_)
+ alpha_ = min;
+ else
+ alpha_ = rho_;
+ phi_const_ = phi_const_min;
+
+ // update phi_const_ for v1 and v2 if necessary
+ SimplexPhiMap::iterator v1_iter = simplices.find(ARSimplex3D(*v1));
+ SimplexPhiMap::iterator v2_iter = simplices.find(ARSimplex3D(*v2));
+ if (phi_const_ < v1_iter->second) v1_iter->second = phi_const_;
+ if (phi_const_ < v2_iter->second) v2_iter->second = phi_const_;
+
+
+ s_ = CGAL::squared_distance(z, CGAL::Segment_3<K>(p1,p2).supporting_line());
+ Point origin;
+ Point cc = origin + ((p1 - origin) + (p2 - origin))/2; // CGAL is funny
+ v_ = CGAL::squared_distance(z, cc) - s_;
+}
+
+ARSimplex3D::
+ARSimplex3D(const Facet& f)
+{
+ Cell_handle c = f.first;
+ for (int i = 0; i < 4; ++i)
+ if (i != f.second)
+ Parent::add(c->vertex(i));
+}
+
+ARSimplex3D::
+ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices)
+{
+ 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 v3 = Parent::vertices().begin();
+ VertexSet::const_iterator v1 = v3++;
+ VertexSet::const_iterator v2 = v3++;
+ const Point& p1 = (*v1)->point();
+ const Point& p2 = (*v2)->point();
+ const Point& p3 = (*v3)->point();
+ rho_ = squared_radius(p1, p2, p3);
+
+ attached_ = false;
+ if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ c->vertex(f.second)->point()) == CGAL::ON_BOUNDED_SIDE)
+ attached_ = true;
+ else if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
+ attached_ = true;
+ else
+ alpha_ = rho_;
+
+ SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
+ SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
+ if (c_iter == simplices.end()) // c is infinite
+ {
+ if (attached_) alpha_ = o_iter->first.alpha();
+ phi_const_ = o_iter->first.phi_const(); // FIXME: it's probably the other way around
+ }
+ else if (o_iter == simplices.end()) // o is infinite
+ {
+ if (attached_) alpha_ = c_iter->first.alpha();
+ phi_const_ = c_iter->first.phi_const(); // FIXME: it's probably the other way around
+ }
+ else
+ {
+ if (attached_) alpha_ = std::min(c_iter->first.alpha(), o_iter->first.alpha());
+ phi_const_ = std::min(c_iter->first.phi_const(), o_iter->first.phi_const());
+ }
+
+ Point cc = CGAL::circumcenter(p1, p2, p3);
+ v_ = CGAL::squared_distance<K>(z, CGAL::Plane_3<K>(p1,p2,p3));
+ s_ = CGAL::squared_distance(z, cc) - v_;
+}
+
+ARSimplex3D::
+ARSimplex3D(const Cell& c, const Point& z): attached_(false)
+{
+ for (int i = 0; i < 4; ++i)
+ Parent::add(c.vertex(i));
+ VertexSet::const_iterator v = Parent::vertices().begin();
+ Point p1 = (*v++)->point();
+ Point p2 = (*v++)->point();
+ Point p3 = (*v++)->point();
+ Point p4 = (*v)->point();
+ rho_ = alpha_ = CGAL::squared_radius(p1, p2, p3, p4);
+
+ s_ = 0;
+ v_ = CGAL::squared_distance(z, CGAL::circumcenter(p1, p2, p3, p4));
+
+ phi_const_ = rho_ - v_;
+}
+
+ARSimplex3D::Cycle
+ARSimplex3D::boundary() const
+{
+ Cycle bdry;
+ Parent::Cycle pbdry = Parent::boundary();
+ for (Parent::Cycle::const_iterator cur = pbdry.begin(); cur != pbdry.end(); ++cur)
+ bdry.push_back(*cur);
+ return bdry;
+}
+
+
+bool
+ARSimplex3D::AlphaOrder::
+operator()(const ARSimplex3D& first, const ARSimplex3D& second) const
+{
+ if (first.alpha() == second.alpha())
+ return (first.dimension() < second.dimension());
+ else
+ return (first.alpha() < second.alpha());
+}
+
+std::ostream&
+ARSimplex3D::
+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 update_simplex_phi_map(const ARSimplex3D& s, ARSimplex3D::SimplexPhiMap& simplices)
+{
+ simplices[s] = s.phi_const();
+}
+
+void fill_alpha_order(const Delaunay& Dt, const Point& z, ARSimplex3DVector& alpha_order)
+{
+ ARSimplex3D::SimplexPhiMap simplices;
+ for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
+ update_simplex_phi_map(ARSimplex3D(*cur, z), simplices);
+ std::cout << "Cells inserted" << std::endl;
+ for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
+ simplices[ARSimplex3D(*cur, z)] = 0; // only one tetrahedron can have non-negative phi_const value
+ // (namely the one containing z); all other simplices will have a
+ // negative phi_const value, so 0 is safe
+ std::cout << "Vertices inserted" << std::endl;
+
+ for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices), simplices);
+ std::cout << "Facets inserted" << std::endl;
+ for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt.incident_facets(*cur)), simplices);
+ std::cout << "Edges inserted" << std::endl;
+
+ // Sort simplices by their alpha values
+ alpha_order.resize(simplices.size()); ARSimplex3DVector::iterator out = alpha_order.begin();
+ for (ARSimplex3D::SimplexPhiMap::const_iterator in = simplices.begin(); in != simplices.end(); ++in, ++out)
+ *out = in->first;
+ std::sort(alpha_order.begin(), alpha_order.end(), ARSimplex3D::AlphaOrder());
+
+ // Update phi_const for vertices
+ for (ARSimplex3DVector::iterator cur = alpha_order.begin(); cur != alpha_order.end(); ++cur)
+ if (cur->dimension() == 0) cur->set_phi_const(simplices[*cur]);
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-vineyard.cpp Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,59 @@
+#include <sys.h>
+#include <debug.h>
+
+#include "ar-vineyard.h"
+
+#include <iostream>
+#include <fstream>
+
+
+int main(int argc, char** argv)
+{
+#ifdef CWDEBUG
+ Debug(dc::filtration.off());
+ Debug(dc::cycle.off());
+ Debug(dc::vineyard.off());
+ Debug(dc::transpositions.off());
+ Debug(dc::lsfiltration.off());
+
+ dionysus::debug::init();
+#endif
+
+ // Read command-line arguments
+ if (argc < 6)
+ {
+ std::cout << "Usage: ar-vineyard POINTS X Y Z OUTFILENAME" << std::endl;
+ std::cout << " POINTS - filename containing points" << std::endl;
+ std::cout << " X,Y,Z - center-point z at which to compute the vineyard" << std::endl;
+ std::cout << " OUTFILENAME - filename for the resulting vineyard" << std::endl;
+ std::cout << std::endl;
+ std::cout << "Computes an (alpha,r)-vineyard of the given pointset around the given point." << std::endl;
+ exit(0);
+ }
+ std::string infilename = argv[1];
+ double zx,zy,zz; std::istringstream(argv[2]) >> zx;
+ std::istringstream(argv[3]) >> zy; std::istringstream(argv[4]) >> zz;
+ std::string outfilename = argv[5];
+
+
+ // Read in the point set and compute its Delaunay triangulation
+ std::ifstream in(infilename.c_str());
+ double x,y,z;
+ ARVineyard::PointList points;
+ while(in)
+ {
+ in >> x >> y >> z;
+ points.push_back(Point(x,y,z));
+ }
+
+
+ // Setup vineyard and compute initial pairing
+ ARVineyard arv(points, Point(zx,zy,zz));
+ arv.compute_pairing();
+
+ // Compute vineyard
+ arv.compute_vineyard(true);
+ std::cout << "Vineyard computed" << std::endl;
+ arv.vineyard()->save_edges(outfilename);
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-vineyard.h Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,191 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ */
+
+#ifndef __AR_VINEYARD_H__
+#define __AR_VINEYARD_H__
+
+#include "sys.h"
+#include "debug.h"
+
+#include "conesimplex.h"
+#include "filtration.h"
+
+#include <CGAL/Kinetic/Inexact_simulation_traits_1.h>
+#include <CGAL/Kinetic/Sort.h>
+#include <CGAL/Kinetic/Sort_visitor_base.h>
+
+#include <list>
+#include "ar-simplex3d.h"
+
+#include <vector>
+
+
+class ARVineyardBase
+{
+ public:
+ /// \name CGAL Kinetic Sort types
+ /// @{
+ class SortVisitor;
+ typedef CGAL::Kinetic::Inexact_simulation_traits_1 Traits;
+ typedef CGAL::Kinetic::Sort<Traits, SortVisitor> Sort;
+ typedef Traits::Simulator Simulator;
+ typedef Traits::Active_points_1_table ActivePointsTable;
+ typedef ActivePointsTable::Key Key;
+
+ typedef Traits::Kinetic_kernel::
+ Function_kernel::Construct_function CF;
+ typedef Traits::Kinetic_kernel::Motion_function F;
+ /// @}
+
+ class ARConeSimplex;
+ class MembershipFunctionChangeEvent;
+};
+
+class ARVineyardBase::ARConeSimplex: public ConeSimplex<ARSimplex3D>
+{
+ public:
+ typedef ConeSimplex<ARSimplex3D> Parent;
+ typedef ARSimplex3D ARSimplex3D;
+
+ ARConeSimplex(const ARSimplex3D& s, bool coned = false):
+ Parent(s, coned) {}
+
+ Key kinetic_key() const { return key_; }
+ void set_kinetic_key(Key k) { key_ = k; }
+
+ private:
+ Key key_;
+};
+
+
+class ARVineyard: public ARVineyardBase
+{
+ public:
+ typedef ARVineyard Self;
+
+ typedef Filtration<ARConeSimplex> ARFiltration;
+ typedef ARFiltration::Simplex Simplex;
+ typedef ARFiltration::Index Index;
+ typedef ARFiltration::Vineyard Vineyard;
+ typedef Vineyard::Evaluator Evaluator;
+ typedef std::map<Key, Index> KeyIndexMap;
+
+ typedef std::list<Point> PointList;
+
+ class StaticEvaluator;
+ class KineticEvaluator;
+
+ public:
+ ARVineyard(const PointList& points, const Point& z);
+ ~ARVineyard();
+
+ void compute_pairing();
+ void compute_vineyard(bool explicit_events = false);
+
+ const ARFiltration* filtration() const { return filtration_; }
+ const Vineyard* vineyard() const { return vineyard_; }
+
+ public:
+ // For Kinetic Sort
+ void swap(Key a, Key b);
+
+ private:
+ void add_simplices();
+ void change_evaluator(Evaluator* eval);
+
+ private:
+ ARFiltration* filtration_;
+ Vineyard* vineyard_;
+ Evaluator* evaluator_;
+
+ KeyIndexMap kinetic_map_;
+
+ Point z_;
+ Delaunay dt_;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ ARVineyard() {}
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type )
+ {
+ // FIXME
+ };
+#endif
+};
+
+//BOOST_CLASS_EXPORT(ARVineyard)
+
+
+class ARVineyardBase::MembershipFunctionChangeEvent
+{
+ public:
+ MembershipFunctionChangeEvent(Key k, F function,
+ ActivePointsTable::Handle apt):
+ key_(k), function_(function), apt_(apt) {}
+
+ void process(Simulator::Time t) const;
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ Key key_;
+ F function_;
+ ActivePointsTable::Handle apt_;
+};
+
+std::ostream& operator<<(std::ostream& out, const ARVineyardBase::MembershipFunctionChangeEvent& e)
+{ return e.operator<<(out); }
+
+class ARVineyard::StaticEvaluator: public Evaluator
+{
+ public:
+ StaticEvaluator(RealType t): time_(t) {}
+
+ virtual RealType time() const { return time_; }
+ virtual RealType value(const Simplex& s) const { return s.value(); }
+
+ private:
+ RealType time_;
+};
+
+class ARVineyard::KineticEvaluator: public Evaluator
+{
+ public:
+ KineticEvaluator(Simulator::Handle sp,
+ ActivePointsTable::Handle apt,
+ RealType time_offset):
+ sp_(sp), apt_(apt) {}
+
+ virtual RealType time() const { return CGAL::to_double(get_time()); }
+ virtual RealType value(const Simplex& s) const { return CGAL::to_double(apt_->at(s.kinetic_key()).x()(get_time())); }
+
+ private:
+ Simulator::Time get_time() const { return sp_->current_time(); }
+
+ Simulator::Handle sp_;
+ ActivePointsTable::Handle apt_;
+};
+
+
+class ARVineyardBase::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
+{
+ public:
+ SortVisitor(ARVineyard* arv): arv_(arv) {}
+
+ template<class Vertex_handle>
+ void before_swap(Vertex_handle a, Vertex_handle b) const;
+
+ private:
+ ARVineyard* arv_;
+};
+
+
+#include "ar-vineyard.hpp"
+
+#endif // __AR_VINEYARD_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-vineyard.hpp Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,170 @@
+/* Implementation */
+
+ARVineyard::
+ARVineyard(const PointList& points, const Point& z): z_(z)
+{
+ for (PointList::const_iterator cur = points.begin(); cur != points.end(); ++cur)
+ dt_.insert(*cur);
+ std::cout << "Delaunay triangulation computed" << std::endl;
+
+ ARSimplex3DVector alpha_ordering;
+ fill_alpha_order(dt_, z_, alpha_ordering);
+ std::cout << "Delaunay simplices: " << alpha_ordering.size() << std::endl;
+
+ evaluator_ = new StaticEvaluator(0);
+ vineyard_ = new Vineyard(evaluator_);
+
+ filtration_ = new ARFiltration(vineyard_);
+ for (ARSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
+ {
+ filtration_->append(*cur); // Delaunay simplex
+ filtration_->append(ARConeSimplex(*cur, true)); // Coned off delaunay simplex
+ }
+}
+
+ARVineyard::
+~ARVineyard()
+{
+ delete filtration_;
+ delete vineyard_;
+ delete evaluator_;
+}
+
+void
+ARVineyard::
+compute_pairing()
+{
+ filtration_->fill_simplex_index_map();
+ filtration_->pair_simplices(filtration_->begin(), filtration_->end());
+ vineyard_->start_vines(filtration_->begin(), filtration_->end());
+ std::cout << "Simplices paired" << std::endl;
+}
+
+void
+ARVineyard::
+compute_vineyard(bool explicit_events)
+{
+ AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
+
+ typedef Traits::Kinetic_kernel::Point_1 Point_1;
+ typedef Simulator::Time Time;
+
+ Traits tr;
+ Simulator::Handle sp = tr.simulator_handle();
+ ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
+ Sort sort(tr, SortVisitor(this));
+
+ // Setup the kinetic sort and membership changes
+ std::cout << "Setting up the kinetic sort and membership events" << std::endl;
+ CF cf;
+ kinetic_map_.clear();
+ for (Index cur = filtration_->begin(); cur != filtration_->end(); ++cur)
+ {
+ F x = cf(F::NT(CGAL::to_double(cur->alpha())));
+ Point_1 p(x);
+ cur->set_kinetic_key(apt->insert(p));
+ kinetic_map_[cur->kinetic_key()] = cur;
+
+ if (!cur->coned()) continue; // non-coned simplices stay put, so we are done
+
+ Time lambda_alpha = CGAL::to_double((cur->alpha() - cur->rho())); // when lambda becomes greater than alpha
+ lambda_alpha += 2*CGAL::sqrt(CGAL::to_double(cur->s()*lambda_alpha));
+ lambda_alpha += CGAL::to_double(cur->s() + cur->v());
+
+ Time phi_alpha = CGAL::to_double(cur->alpha() - cur->phi_const());
+
+ Time phi_lambda = CGAL::to_double(cur->rho() + cur->s() - cur->v() - cur->phi_const());
+ phi_lambda *= phi_lambda;
+ phi_lambda /= CGAL::to_double(4*cur->s());
+ phi_lambda += CGAL::to_double(cur->v());
+
+ Time sv = CGAL::to_double(cur->s() + cur->v());
+
+ if (true || phi_lambda < sv || phi_lambda < phi_alpha) // FIXME: remove true
+ {
+ sp->new_event(Time(phi_alpha),
+ MembershipFunctionChangeEvent(cur->kinetic_key(),
+ cf(F::NT(CGAL::to_double(cur->phi_const())), 1),
+ apt)); // \phi^2 = r^2 + \phi_c^2
+ std::cout << "Scheduled" << std::endl;
+ } else
+ std::cout << "Not scheduled" << std::endl;
+
+
+ //sp->new_event(Time(...), MembershipFunctionChangeEvent(cur->kinetic_key()));
+
+ std::cout << *cur << std::endl;
+ std::cout << "lambda_alpha: " << lambda_alpha << std::endl;
+ std::cout << "phi_alpha: " << phi_alpha << std::endl;
+ std::cout << "phi_lambda: " << phi_lambda << std::endl;
+ std::cout << "s^2 + v^2: " << sv << std::endl;
+ std::cout << std::endl;
+ }
+
+ // Process all the events (compute the vineyard in the process)
+ // FIXME: the time should not be 1, but something like twice the radius of
+ // the pointset as seen from z
+ change_evaluator(new KineticEvaluator(sp, apt, 0));
+ if (explicit_events)
+ {
+ while (sp->next_event_time() < 1)
+ {
+ std::cout << "Next event time: " << sp->next_event_time() << std::endl;
+ sp->set_current_event_number(sp->current_event_number() + 1);
+ std::cout << "Processed event" << std::endl;
+ }
+ } else
+ sp->set_current_time(1.0);
+ std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
+
+ //change_evaluator(new StaticEvaluator(1));
+ vineyard_->record_diagram(filtration_->begin(), filtration_->end());
+}
+
+void
+ARVineyard::
+swap(Key a, Key b)
+{
+ Index ao = kinetic_map_[a], bo = kinetic_map_[b];
+ AssertMsg(filtration_->get_trails_cmp()(ao, bo), "In swap(a,b), a must precede b");
+ filtration_->transpose(ao);
+ AssertMsg(filtration_->get_trails_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
+}
+
+void
+ARVineyard::
+change_evaluator(Evaluator* eval)
+{
+ AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
+
+ delete evaluator_;
+ evaluator_ = eval;
+ vineyard_->set_evaluator(evaluator_);
+}
+
+void
+ARVineyardBase::MembershipFunctionChangeEvent::
+process(Simulator::Time t) const
+{
+ apt_->set(key_, function_);
+ std::cout << "Updated for phi's dominance" << std::endl;
+}
+
+
+template<class Vertex_handle>
+void
+ARVineyardBase::SortVisitor::
+before_swap(Vertex_handle a, Vertex_handle b) const
+{
+ std::cout << "Swapping elements" << *a << " and " << *b << std::endl;
+ arv_->swap(*a,*b);
+}
+
+
+std::ostream&
+ARVineyardBase::MembershipFunctionChangeEvent::
+operator<<(std::ostream& out) const
+{
+ return out << "Membership change" << std::endl;
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/conesimplex.h Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,39 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ */
+
+#ifndef __CONESIMPLEX_H__
+#define __CONESIMPLEX_H__
+
+#include <list>
+#include <iostream>
+
+template<class S>
+class ConeSimplex: public S
+{
+ public:
+ typedef S Parent;
+ typedef ConeSimplex<S> Self;
+ typedef std::list<Self> Cycle;
+
+ public:
+ ConeSimplex(const Parent& parent,
+ bool coned = false):
+ Parent(parent), coned_(coned) {}
+
+ Cycle boundary() const;
+ bool coned() const { return coned_; }
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ bool coned_;
+};
+
+template<class S>
+std::ostream& operator<<(std::ostream& out, const ConeSimplex<S>& s) { return s.operator<<(out); }
+
+#include "conesimplex.hpp"
+
+#endif // __CONESIMPLEX_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/conesimplex.hpp Wed Aug 15 11:25:35 2007 -0400
@@ -0,0 +1,24 @@
+template<class S>
+typename ConeSimplex<S>::Cycle
+ConeSimplex<S>::boundary() const
+{
+ Cycle bdry;
+ typename Parent::Cycle pbdry = Parent::boundary();
+
+ for (typename Parent::Cycle::const_iterator cur = pbdry.begin(); cur != pbdry.end(); ++cur)
+ bdry.push_back(Self(*cur, coned_));
+
+ if (coned_)
+ bdry.push_back(Self(*this, false));
+
+ return bdry;
+}
+
+template<class S>
+std::ostream&
+ConeSimplex<S>::operator<<(std::ostream& out) const
+{
+ Parent::operator<<(out) << ' ';
+ if (coned_) out << "[coned]";
+ return out;
+}
--- a/include/filtrationsimplex.h Wed Feb 07 14:00:39 2007 -0500
+++ b/include/filtrationsimplex.h Wed Aug 15 11:25:35 2007 -0400
@@ -31,8 +31,10 @@
public:
typedef Smplx Simplex;
- virtual RealType time() { return 0; }
- virtual RealType value(const Simplex& s) { return 0; }
+ virtual RealType time() const { return 0; }
+ virtual RealType value(const Simplex& s) const { return 0; }
+
+ virtual ~Evaluator() {}
};
/**