Merged geometry and include breakup with AR vineyard
authorDmitriy Morozov <morozov@cs.duke.edu>
Wed, 15 Aug 2007 11:39:34 -0400
changeset 22 c42734397d62
parent 21 565a3e60eb7e (diff)
parent 20 7bf6aa6b0ab6 (current diff)
child 23 cb700b407c0d
child 55 7e71f5984931
Merged geometry and include breakup with AR vineyard
examples/CMakeLists.txt
examples/alphashapes/alphashapes3d.cpp
examples/alphashapes/alphashapes3d.h
examples/ar-vineyard/ar-simplex3d.h
examples/ar-vineyard/ar-vineyard.cpp
examples/ar-vineyard/ar-vineyard.h
include/circular_list.h
include/conesimplex.h
include/conesimplex.hpp
include/consistencylist.h
include/consistencylist.hpp
include/counter.h
include/cycle.h
include/cycle.hpp
include/debug.h
include/filtration.h
include/filtration.hpp
include/filtrationcontainer.h
include/filtrationsimplex.h
include/lowerstarfiltration.h
include/lowerstarfiltration.hpp
include/orderlist.h
include/orderlist.hpp
include/simplex.h
include/simplex.hpp
include/sys.h
include/topology/conesimplex.h
include/topology/conesimplex.hpp
include/topology/filtrationsimplex.h
include/types.h
include/vineyard.h
include/vineyard.hpp
tests/test-consistencylist.cpp
tests/test-orderlist.cpp
--- a/examples/CMakeLists.txt	Tue Aug 14 17:15:08 2007 -0400
+++ b/examples/CMakeLists.txt	Wed Aug 15 11:39:34 2007 -0400
@@ -1,4 +1,5 @@
 add_subdirectory			(alphashapes)
+add_subdirectory			(ar-vineyard)
 add_subdirectory			(cech-complex)
 add_subdirectory			(grid)
 add_subdirectory			(triangle)
--- a/examples/alphashapes/alphashapes3d.cpp	Tue Aug 14 17:15:08 2007 -0400
+++ b/examples/alphashapes/alphashapes3d.cpp	Wed Aug 15 11:39:34 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	Tue Aug 14 17:15:08 2007 -0400
+++ b/examples/alphashapes/alphashapes3d.h	Wed Aug 15 11:39:34 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:39:34 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:39:34 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 <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>    		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:39:34 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:39:34 2007 -0400
@@ -0,0 +1,59 @@
+#include <utilities/sys.h>
+#include <utilities/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:39:34 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 "utilities/sys.h"
+#include "utilities/debug.h"
+
+#include "topology/conesimplex.h"
+#include "topology/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:39:34 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/topology/conesimplex.h	Wed Aug 15 11:39:34 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/topology/conesimplex.hpp	Wed Aug 15 11:39:34 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/topology/filtrationsimplex.h	Tue Aug 14 17:15:08 2007 -0400
+++ b/include/topology/filtrationsimplex.h	Wed Aug 15 11:39:34 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()								{}
 };
 
 /**