Made storing negative simplices in StaticPersistence parameter-based, so that DynamicPersistence works again
#include "alphashapes3d.h"
#include "../cohomology/wrappers.h"
#include <topology/cohomology-persistence.h>
#include <utilities/log.h>
#include <utilities/timer.h>
#include <iostream>
#include <fstream>
#include <boost/program_options.hpp>
#include <boost/progress.hpp>
#include <boost/foreach.hpp>
typedef boost::tuple<Dimension, RealType> BirthInfo;
typedef CohomologyPersistence<BirthInfo> Persistence;
typedef Persistence::SimplexIndex Index;
typedef Persistence::Death Death;
namespace po = boost::program_options;
int main(int argc, char** argv)
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/persistence") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/chain") );
#endif
std::string infilename, outfilename;
po::options_description hidden("Hidden options");
hidden.add_options()
("input-file", po::value<std::string>(&infilename), "Point set whose alpha shape filtration and persistence we want to compute")
("output-file", po::value<std::string>(&outfilename), "Where to write the collection of persistence diagrams");
po::positional_options_description pos;
pos.add("input-file", 1);
pos.add("output-file", 2);
po::options_description all; all.add(hidden);
po::variables_map vm;
po::store(po::command_line_parser(argc, argv).
options(all).positional(pos).run(), vm);
po::notify(vm);
if (!vm.count("input-file") || !vm.count("output-file"))
{
std::cout << "Usage: " << argv[0] << " input-file output-file" << std::endl;
// std::cout << hidden << std::endl;
return 1;
}
std::ofstream diagram_out(outfilename.c_str());
Timer total_timer; total_timer.start();
// Read in the point set and compute its Delaunay triangulation
std::ifstream in(infilename.c_str());
double x,y,z;
Delaunay3D Dt;
while(in)
{
in >> x >> y >> z;
Point p(x,y,z);
Dt.insert(p);
}
rInfo("Delaunay triangulation computed");
// Set up the alpha shape filtration
AlphaSimplex3DVector complex;
fill_complex(Dt, complex);
rInfo("Simplices: %d", complex.size());
std::sort(complex.begin(), complex.end(), AlphaSimplex3D::AlphaOrder());
Timer persistence_timer; persistence_timer.start();
std::map<AlphaSimplex3D, Index, AlphaSimplex3D::VertexComparison> complex_map;
Persistence p;
boost::progress_display show_progress(complex.size());
#ifdef COUNTERS
Counter::CounterType max_element_count = 0;
#endif
for(AlphaSimplex3DVector::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
{
const AlphaSimplex3D& s = *cur;
std::vector<Index> boundary;
for (AlphaSimplex3D::BoundaryIterator bcur = s.boundary_begin(); bcur != s.boundary_end(); ++bcur)
boundary.push_back(complex_map[*bcur]);
Index idx; Death d;
bool store = s.dimension() < 3;
boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(s.dimension(), s.value()), store);
// c[*cur] = idx;
if (store)
complex_map[s] = idx;
if (d && (s.value() - d->get<1>()) > 0)
{
AssertMsg(d->get<0>() == s.dimension() - 1, "Dimensions must match");
diagram_out << (s.dimension() - 1) << " " << d->get<1>() << " " << s.value() << std::endl;
}
++show_progress;
#ifdef COUNTERS
max_element_count = std::max(max_element_count, cCohomologyElementCount->count);
#endif
}
// output infinte persistence pairs
for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
persistence_timer.stop();
total_timer.stop();
persistence_timer.check("Persistence timer");
total_timer.check("Total timer");
#ifdef COUNTERS
std::cout << "Max element count: " << max_element_count << std::endl;
#endif
}