examples/cech-complex/cech-complex.cpp
 author Dmitriy Morozov Wed, 23 May 2007 10:48:04 -0400 changeset 18 a4d403087fbc child 20 7bf6aa6b0ab6 permissions -rw-r--r--
```
#include <sys.h>
#include <debug.h>

#include <simplex.h>
#include <filtration.h>
#include "Miniball_dynamic_d.h"
#include <algorithm>

#include <iostream>
#include <fstream>

typedef 		std::vector<Point>										PointContainer;
typedef			unsigned int 											PointIndex;
typedef			SimplexWithValue<PointIndex>							Simplex;
typedef 		std::vector<Simplex> 									SimplexVector;
typedef 		Filtration<Simplex>										CechFiltration;

class DimensionValueComparison
{
public:
bool	operator()(const Simplex& s1, const Simplex& s2) const
{
if (s1.dimension() == s2.dimension())
return s1.get_value() < s2.get_value();
else
return s1.dimension() < s2.dimension();
}
};

int choose(int n, int k)
{
if (k > n/2) k = n-k;
int numerator = 1, denominator = 1;
for (int i = 0; i < k; ++i)
{ numerator *= (n-i); denominator *= (i+1); }
return numerator/denominator;
}

void add_simplices(SimplexVector& sv, int d, const PointContainer& points)
{
PointIndex indices[d+1];
for (int i = 0; i < d+1; ++i)
indices[i] = d - i;

while(indices[d] < points.size() - d)
{
Miniball mb(points[indices[0]].dim());
Simplex s;
for (int i = 0; i < d+1; ++i)
{
mb.check_in(points[indices[i]]);
}
mb.build();
sv.push_back(s);

for (int i = 0; i < d+1; ++i)
{
++indices[i];
if (indices[i] < points.size() - i)
{
for (int j = i-1; j >= 0; --j)
indices[j] = indices[j+1] + 1;
break;
}
}
}
}

int main(int argc, char** argv)
{
// Read in the point set and compute its Delaunay triangulation
std::istream& in = std::cin;
int ambient_d, homology_d;
in >> ambient_d >> homology_d;

std::cout << "Ambient dimension: " << ambient_d << std::endl;
std::cout << "Will compute PD up to dimension: " << homology_d << std::endl;

PointContainer points;
while(in)
{
Point p(ambient_d);
for (int i = 0; i < ambient_d; ++i)
in >> p[i];
points.push_back(p);
}
std::cout << "Points read: " << points.size() << std::endl;

// Compute Cech values
CechFiltration cf;
{										// resource acquisition is initialization
SimplexVector	sv;
int num_simplices = 0;
for (int i = 0; i <= homology_d + 1; ++i)
num_simplices += choose(points.size(), i+1);
sv.reserve(num_simplices);
std::cout << "Reserved SimplexVector of size: " << num_simplices << std::endl;

for (int i = 0; i <= homology_d + 1; ++i)
std::cout << "Size of SimplexVector: " << sv.size() << std::endl;

std::sort(sv.begin(), sv.end(), DimensionValueComparison());

for (SimplexVector::const_iterator cur = sv.begin(); cur != sv.end(); ++cur)
cf.append(*cur);
}

// Compute persistence
cf.fill_simplex_index_map();
cf.pair_simplices(cf.begin(), cf.end());
std::cout << "Simplices paired" << std::endl;

for (CechFiltration::Index i = cf.begin(); i != cf.end(); ++i)
if (i->is_paired())
{
if (i->sign())
std::cout << i->dimension() << " " << i->get_value() << " " << i->pair()->get_value() << std::endl;
} //else std::cout << i->value() << std::endl;

}

```