--- a/examples/SConscript	Fri Dec 22 10:25:26 2006 -0500
+++ b/examples/SConscript	Fri Dec 22 12:52:35 2006 -0500
@@ -1,6 +1,6 @@
 Import('*')
 
 SConscript				(['triangle/SConscript',
-#						  'alphashapes/SConscript',
+						  'alphashapes/SConscript',
 						  'grid/SConscript'],
 						 exports=['env', 'external_sources'])

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/SConscript	Fri Dec 22 12:52:35 2006 -0500
@@ -0,0 +1,7 @@
+Import('*')
+
+# Sources
+sources = 				['alphashapes3d.cpp', 'alpharadius.cpp']
+
+for s in sources:
+	Default				(env.Program([s] + external_sources))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alpharadius.cpp	Fri Dec 22 12:52:35 2006 -0500
@@ -0,0 +1,113 @@
+#include <sys.h>
+#include <debug.h>
+
+#include "alphashapes3d.h"
+#include <filtration.h>
+#include <iostream>
+#include <fstream>
+
+
+struct SimplexWithDistance: public AlphaSimplex3D
+{
+	SimplexWithDistance(const SimplexWithDistance& s): AlphaSimplex3D(s)				{ distance = s.distance; }
+	SimplexWithDistance(const AlphaSimplex3D& s): AlphaSimplex3D(s)						{ }
+	SimplexWithDistance(const AlphaSimplex3D::Parent& s): AlphaSimplex3D(s)				{ }
+	SimplexWithDistance(const AlphaSimplex3D& s, const Point& p): AlphaSimplex3D(s)		{ set_distance(p); }
+
+	void set_distance(const Point& p)
+	{
+		AlphaSimplex3D::VertexContainer::const_iterator cur = vertices().begin();
+		CGAL::Vector_3<K> v = ((*cur)->point() - p);
+		RealValue min_distance = v*v;
+	
+		while (cur != vertices().end())
+		{
+			v = ((*cur)->point() - p);
+			min_distance = std::min(v*v, min_distance);
+			++cur;
+		}
+	
+		distance = min_distance;
+	}
+
+	RealValue distance;
+};
+typedef std::vector<SimplexWithDistance> 		SimplexVector;
+typedef Filtration<SimplexWithDistance> 		RadiusFiltration;
+
+struct RadiusOrder
+{
+	bool operator()(const SimplexWithDistance& first, const SimplexWithDistance& second) const
+	{
+		if (first.distance == second.distance)
+			return (first.dimension() < second.dimension());
+		else
+			return (first.distance > second.distance); 
+	}
+};
+
+int main(int argc, char** argv) 
+{
+#ifdef CWDEBUG
+	Debug(dc::filtration.on());
+	//Debug(dc::cycle.on());
+
+	dionysus::debug::init();
+#endif
+
+	std::istream& in = std::cin;
+	std::ofstream point_vineyard("point_vineyard.vrt");
+
+	double x,y,z;
+	Delaunay Dt;
+	while(in)
+	{
+		in >> x >> y >> z;
+		Point p(x,y,z);
+		Dt.insert(p);
+	}
+	std::cout << "Delaunay triangulation computed" << std::endl;
+ 
+	AlphaSimplex3DVector alpha_ordering;
+	fill_alpha_order(Dt, alpha_ordering);
+    
+	// Compute r-filtration for each distinct alpha
+	Point p(0,0,0);
+	RadiusOrder ro;
+	SimplexVector radius_ordering;
+	for (AlphaSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
+	{
+		radius_ordering.push_back(*cur);
+		radius_ordering.back().set_distance(p);
+		if (boost::next(cur)->alpha() == cur->alpha())
+			continue;
+		
+		double current_alpha = CGAL::to_double(cur->value());
+		std::cout << "Current alpha: " << current_alpha << std::endl;
+		std::sort(radius_ordering.begin(), radius_ordering.end(), ro);
+		std::cout << "Radius ordering size: " << radius_ordering.size() << std::endl;
+
+		RadiusFiltration rf;
+		for (SimplexVector::const_iterator cur = radius_ordering.begin(); cur != radius_ordering.end(); ++cur)
+			rf.append(*cur);
+		rf.fill_simplex_index_map();
+		std::cout << "Simplex index map filled" << std::endl;
+		rf.pair_simplices(rf.begin(), rf.end());
+		std::cout << "Pairing computed" << std::endl;
+	
+		for (RadiusFiltration::const_Index cur = rf.begin(); cur != rf.end(); ++cur)
+		{
+			if (!cur->sign()) continue;
+	
+			RealValue d1 = cur->distance;
+			//if (cur == cur->pair())
+			//	std::cout << "Unpaired " << cur->dimension() << ' ' << CGAL::to_double(d1) << std::endl;
+			
+			RealValue d2 = cur->pair()->distance;
+			if (d1 == d2)	continue;
+			point_vineyard << CGAL::to_double(d1) << ' ' << CGAL::to_double(d2) << ' ' << current_alpha << std::endl;
+		}
+	}
+
+	point_vineyard.close();
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d.cpp	Fri Dec 22 12:52:35 2006 -0500
@@ -0,0 +1,39 @@
+#include <sys.h>
+#include <debug.h>
+
+#include "alphashapes3d.h"
+#include <filtration.h>
+#include <iostream>
+#include <fstream>
+
+
+typedef std::vector<AlphaSimplex3D> 			SimplexVector;
+typedef Filtration<AlphaSimplex3D>				AlphaFiltration;
+
+int main(int argc, char** argv) 
+{
+	// Read in the point set and compute its Delaunay triangulation
+	std::istream& in = std::cin;
+	double x,y,z;
+	Delaunay Dt;
+	while(in)
+	{
+		in >> x >> y >> z;
+		Point p(x,y,z);
+		Dt.insert(p);
+	}
+	std::cout << "Delaunay triangulation computed" << std::endl;
+   
+	AlphaSimplex3DVector alpha_ordering;
+	fill_alpha_order(Dt, alpha_ordering);
+	std::cout << "Simplices: " << alpha_ordering.size() << std::endl;
+
+	// Create the alpha-shape filtration
+	AlphaFiltration af;
+	for (SimplexVector::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());
+	std::cout << "Simplices paired" << std::endl;
+}
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d.h	Fri Dec 22 12:52:35 2006 -0500
@@ -0,0 +1,84 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2006
+ */
+
+#ifndef __ALPHASHAPES3D_H__
+#define __ALPHASHAPES3D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_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 AlphaSimplex3D: public SimplexWithVertices<Vertex_handle>
+{
+	public:
+		typedef 	std::set<AlphaSimplex3D>							SimplexSet;
+		typedef		SimplexWithVertices<Vertex_handle>					Parent;
+		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 ::Vertex& v);
+		
+								    AlphaSimplex3D(const Edge& e);
+								    AlphaSimplex3D(const Edge& e, const SimplexSet& simplices, Facet_circulator facet_bg);
+		
+								    AlphaSimplex3D(const Facet& f);
+								    AlphaSimplex3D(const Facet& f, const SimplexSet& simplices);
+	    
+									AlphaSimplex3D(const Cell& c);
+	    
+		virtual 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; };
+		
+		std::ostream& 				operator<<(std::ostream& out) const;
+		
+	private:
+		RealValue 					alpha_;
+		bool 						attached_;
+};
+
+typedef 			std::vector<AlphaSimplex3D>								AlphaSimplex3DVector;
+void 				fill_alpha_order(const Delaunay& Dt, 
+									 AlphaSimplex3DVector& alpha_order);
+
+std::ostream& 		operator<<(std::ostream& out, const AlphaSimplex3D& s)	{ return s.operator<<(out); }
+
+#include "alphashapes3d.hpp"
+
+#endif // __ALPHASHAPES3D_H__

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d.hpp	Fri Dec 22 12:52:35 2006 -0500
@@ -0,0 +1,171 @@
+AlphaSimplex3D::	    
+AlphaSimplex3D(const ::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));
+}
+
+AlphaSimplex3D::	    
+AlphaSimplex3D(const Edge& e)
+{
+    Cell_handle c = e.first;
+	Parent::add(c->vertex(e.second));
+	Parent::add(c->vertex(e.third));
+}
+
+AlphaSimplex3D::	    
+AlphaSimplex3D(const Edge& e, const SimplexSet& 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;
+	SimplexSet::const_iterator cur_iter = simplices.find(AlphaSimplex3D(*cur));
+	while (cur_iter == simplices.end())
+	{
+		++cur; 
+		cur_iter = simplices.find(AlphaSimplex3D(*cur));
+	}
+	RealValue min = cur_iter->alpha();
+	
+	VertexSet::const_iterator v = Parent::vertices().begin();
+	const Point& p1 = (*v++)->point();
+	const Point& p2 = (*v)->point();
+	attached_ = false;
+
+	if (facet_bg != 0) do
+	{
+		VertexSet::const_iterator v = Parent::vertices().begin();
+		int i0 = (*cur).first->index(*v++);
+		int i1 = (*cur).first->index(*v);
+		int i = 6 - i0 - i1 - (*cur).second;
+		Point p3 = (*cur).first->vertex(i)->point();
+
+		cur_iter = simplices.find(AlphaSimplex3D(*cur));
+		if (cur_iter == simplices.end())			// cur is infinite
+		{
+			++cur; continue;
+		}
+		
+		if (CGAL::side_of_bounded_sphere(p1, p2, p3) == CGAL::ON_BOUNDED_SIDE)
+			attached_ = true;
+		RealValue val = cur_iter->alpha();
+		if (val < min)
+			min = val;
+		++cur;
+	} while (cur != facet_bg);
+
+	if (attached_)
+		alpha_ = min;
+	else
+		alpha_ = CGAL::squared_radius(p1, p2);
+}
+
+AlphaSimplex3D::	    
+AlphaSimplex3D(const 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 Facet& f, const SimplexSet& 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 v = Parent::vertices().begin();
+	const Point& p1 = (*v++)->point();
+	const Point& p2 = (*v++)->point();
+	const Point& p3 = (*v)->point();
+	
+	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_ = squared_radius(p1, p2, p3);
+	
+	if (attached_)
+	{
+		SimplexSet::const_iterator c_iter = simplices.find(AlphaSimplex3D(*c));
+		SimplexSet::const_iterator o_iter = simplices.find(AlphaSimplex3D(*o));
+		if (c_iter == simplices.end())			// c is infinite
+			alpha_ = o_iter->alpha();
+		else if (o_iter == simplices.end())		// o is infinite
+			alpha_ = c_iter->alpha();
+		else
+			alpha_ = std::min(c_iter->alpha(), o_iter->alpha());
+	}
+}
+
+AlphaSimplex3D::	    
+AlphaSimplex3D(const Cell& c): 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();
+	alpha_ = CGAL::squared_radius(p1, p2, p3, p4);
+}
+
+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_alpha_order(const Delaunay& Dt, AlphaSimplex3DVector& alpha_order)
+{
+	// Compute all simplices with their alpha values and attachment information
+	AlphaSimplex3D::SimplexSet simplices;
+	for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
+		simplices.insert(AlphaSimplex3D(*cur));
+	std::cout << "Cells inserted" << std::endl;
+	for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
+		simplices.insert(AlphaSimplex3D(*cur, simplices));
+	std::cout << "Facets inserted" << std::endl;
+	for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
+		simplices.insert(AlphaSimplex3D(*cur, simplices, Dt.incident_facets(*cur)));
+	std::cout << "Edges inserted" << std::endl;
+	for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
+		simplices.insert(AlphaSimplex3D(*cur));
+	std::cout << "Vertices inserted" << std::endl;
+    
+	// Sort simplices by their alpha values
+	alpha_order.resize(simplices.size());
+	std::copy(simplices.begin(), simplices.end(), alpha_order.begin());
+	std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex3D::AlphaOrder());
+}
+