Added 2D alpha shape filtration example
authorDmitriy Morozov <morozov@cs.duke.edu>
Thu, 11 Oct 2007 11:46:00 -0400
changeset 39 b85ded6b1530
parent 38 aee8f549d324
child 40 122e4a1fa117
Added 2D alpha shape filtration example
examples/alphashapes/CMakeLists.txt
examples/alphashapes/alphashapes2d.cpp
examples/alphashapes/alphashapes2d.h
examples/alphashapes/alphashapes2d.hpp
include/topology/filtration.hpp
--- a/examples/alphashapes/CMakeLists.txt	Thu Sep 20 08:02:49 2007 -0400
+++ b/examples/alphashapes/CMakeLists.txt	Thu Oct 11 11:46:00 2007 -0400
@@ -1,5 +1,6 @@
 set							(targets						
 							 alphashapes3d
+							 alphashapes2d
 							 alpharadius)
 							 
 add_definitions				(${cgal_cxxflags})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.cpp	Thu Oct 11 11:46:00 2007 -0400
@@ -0,0 +1,58 @@
+#include <utilities/log.h>
+
+#include "alphashapes2d.h"
+#include <topology/filtration.h>
+#include <iostream>
+#include <fstream>
+
+
+typedef Filtration<AlphaSimplex2D>				AlphaFiltration;
+
+int main(int argc, char** argv) 
+{
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
+
+	stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+#endif
+
+	SetFrequency(GetCounter("filtration/pair"), 10000);
+	SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+
+	// Read in the point set and compute its Delaunay triangulation
+	std::istream& in = std::cin;
+	double x,y;
+	Delaunay Dt;
+	while(in)
+	{
+		in >> x >> y;
+		Point p(x,y);
+		Dt.insert(p);
+	}
+	std::cout << "Delaunay triangulation computed" << std::endl;
+   
+	AlphaSimplex2DVector alpha_ordering;
+	fill_alpha_order(Dt, alpha_ordering);
+	std::cout << "Simplices: " << alpha_ordering.size() << std::endl;
+
+	// Create the alpha-shape filtration
+	AlphaFiltration af;
+	for (AlphaSimplex2DVector::const_iterator cur = alpha_ordering.begin(); 
+											  cur != alpha_ordering.end(); ++cur)
+		af.append(*cur);
+	af.fill_simplex_index_map();
+	std::cout << "Filled simplex-index map" << std::endl;
+	af.pair_simplices(af.begin(), af.end(), false);
+	std::cout << "Simplices paired" << std::endl;
+
+	for (AlphaFiltration::Index i = af.begin(); i != af.end(); ++i)
+		if (i->is_paired())
+		{
+			if (i->sign())
+				std::cout << i->dimension() << " " << i->value() << " " << i->pair()->value() << std::endl;
+		} //else std::cout << i->value() << std::endl;
+
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.h	Thu Oct 11 11:46:00 2007 -0400
@@ -0,0 +1,80 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ */
+
+#ifndef __ALPHASHAPES2D_H__
+#define __ALPHASHAPES2D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/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>    		Delaunay;
+typedef Delaunay::Point                				Point;
+typedef Delaunay::Vertex            				Vertex;
+typedef Delaunay::Vertex_handle            			Vertex_handle;
+typedef Delaunay::Edge								Edge;
+typedef Delaunay::Face								Face;
+typedef Delaunay::Face_handle						Face_handle;
+typedef K::FT										RealValue;
+
+typedef Delaunay::Finite_vertices_iterator    		Vertex_iterator;
+typedef Delaunay::Finite_edges_iterator        		Edge_iterator;
+typedef Delaunay::Finite_faces_iterator        		Face_iterator;
+
+
+class AlphaSimplex2D: public SimplexWithVertices<Vertex_handle>
+{
+	public:
+		typedef 	std::set<AlphaSimplex2D>							SimplexSet;
+		typedef		SimplexWithVertices<Vertex_handle>					Parent;
+		typedef		Parent::VertexContainer								VertexSet;
+		typedef		std::list<AlphaSimplex2D>							Cycle;
+
+    public:
+									AlphaSimplex2D()					{}
+									AlphaSimplex2D(const Parent& p): 
+											Parent(p) 					{}
+									AlphaSimplex2D(const AlphaSimplex2D& s): 
+											Parent(s) 					{ attached_ = s.attached_; alpha_ = s.alpha_; }
+	    							AlphaSimplex2D(const ::Vertex& v);
+		
+								    AlphaSimplex2D(const Edge& e);
+								    AlphaSimplex2D(const Edge& e, const SimplexSet& simplices);
+		
+									AlphaSimplex2D(const Face& c);
+	    
+		RealType					value() const						{ return CGAL::to_double(alpha_); }
+		RealValue					alpha() const						{ return alpha_; }
+		bool						attached() const					{ return attached_; }
+		Cycle						boundary() const;
+
+		// Ordering
+		struct AlphaOrder
+		{ bool operator()(const AlphaSimplex2D& first, const AlphaSimplex2D& second) const; };
+		
+		std::ostream& 				operator<<(std::ostream& out) const;
+		
+	private:
+		RealValue 					alpha_;
+		bool 						attached_;
+};
+
+typedef 			std::vector<AlphaSimplex2D>								AlphaSimplex2DVector;
+void 				fill_alpha_order(const Delaunay& Dt, 
+									 AlphaSimplex2DVector& alpha_order);
+
+std::ostream& 		operator<<(std::ostream& out, const AlphaSimplex2D& s)	{ return s.operator<<(out); }
+
+#include "alphashapes2d.hpp"
+
+#endif // __ALPHASHAPES2D_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.hpp	Thu Oct 11 11:46:00 2007 -0400
@@ -0,0 +1,121 @@
+AlphaSimplex2D::	    
+AlphaSimplex2D(const ::Vertex& v): alpha_(0), attached_(false)
+{
+	for (int i = 0; i < 3; ++i)
+		if (v.face()->vertex(i)->point() == v.point())
+			Parent::add(v.face()->vertex(i));
+}
+
+AlphaSimplex2D::	    
+AlphaSimplex2D(const 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 Edge& e, const SimplexSet& simplices): 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);
+
+	VertexSet::const_iterator v = Parent::vertices().begin();
+	const Point& p1 = (*v++)->point();
+	const Point& p2 = (*v)->point();
+	
+	attached_ = false;
+	if (CGAL::side_of_bounded_circle(p1, p2, 
+									 f->vertex(e.second)->point()) == CGAL::ON_BOUNDED_SIDE)
+		attached_ = true;
+	else if (CGAL::side_of_bounded_circle(p1, p2,
+										  o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
+			attached_ = true;
+	else
+		alpha_ = squared_radius(p1, p2);
+
+	if (attached_)
+	{
+		SimplexSet::const_iterator f_iter = simplices.find(AlphaSimplex2D(*f));
+		SimplexSet::const_iterator o_iter = simplices.find(AlphaSimplex2D(*o));
+		if (f_iter == simplices.end())			// f is infinite
+			alpha_ = o_iter->alpha();
+		else if (o_iter == simplices.end())		// o is infinite
+			alpha_ = f_iter->alpha();
+		else
+			alpha_ = std::min(f_iter->alpha(), o_iter->alpha());
+	}
+}
+
+AlphaSimplex2D::	    
+AlphaSimplex2D(const Face& f): attached_(false)
+{
+	for (int i = 0; i < 3; ++i)
+		Parent::add(f.vertex(i));
+	VertexSet::const_iterator v = Parent::vertices().begin();
+	Point p1 = (*v++)->point();
+	Point p2 = (*v++)->point();
+	Point p3 = (*v)->point();
+	alpha_ = CGAL::squared_radius(p1, p2, p3);
+}
+
+AlphaSimplex2D::Cycle
+AlphaSimplex2D::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 
+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_alpha_order(const Delaunay& Dt, AlphaSimplex2DVector& alpha_order)
+{
+	// Compute all simplices with their alpha values and attachment information
+	AlphaSimplex2D::SimplexSet simplices;
+	for(Face_iterator cur = Dt.finite_faces_begin(); cur != Dt.finite_faces_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*cur));
+	std::cout << "Faces inserted" << std::endl;
+	for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*cur, simplices));
+	std::cout << "Edges inserted" << std::endl;
+	for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*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(), AlphaSimplex2D::AlphaOrder());
+}
+
--- a/include/topology/filtration.hpp	Thu Sep 20 08:02:49 2007 -0400
+++ b/include/topology/filtration.hpp	Thu Oct 11 11:46:00 2007 -0400
@@ -387,8 +387,8 @@
 
 	for (typename Simplex::Cycle::const_iterator cur = bdry.begin(); cur != bdry.end(); ++cur)
 	{
-		rLog(rlFiltrationTranspositions, "Appending in init_cycle_trail(): %s", 
-										 tostring(*cur).c_str());
+		rLog(rlFiltration, "Appending in init_cycle_trail(): %s", 
+						   tostring(*cur).c_str());
 		AssertMsg(get_index(*cur) != end(), "Non-existent simplex in the cycle");
 		j->cycle().append(get_index(*cur), get_consistency_cmp());
 	}