/* This file is part of the Gudhi Library. The Gudhi library
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
* Author(s): Pawel Dlotko
*
* Copyright (C) 2015 INRIA Sophia-Saclay (France)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
#ifndef BITMAP_CUBICAL_COMPLEX_H_
#define BITMAP_CUBICAL_COMPLEX_H_
#include
#include
//global variable, was used just for debugging.
bool globalDbg = false;
template
class Bitmap_cubical_complex : public Bitmap_cubical_complex_base {
public:
//*********************************************************************************************************************************//
//Typedefs and typenames
//*********************************************************************************************************************************//
friend class Simplex_handle;
typedef size_t Simplex_key;
typedef T Filtration_value;
//*********************************************************************************************************************************//
//Simplex handle class
//*********************************************************************************************************************************//
/**
* Handle of a cell, required for compatibility with the function to compute persistence in Gudhi. Elements of this class are: the pointer to the bitmap B in which the considered cell is
* together with a position of this cell in B. Given this data, one can get all the information about the considered cell.
**/
class Simplex_handle {
public:
Simplex_handle() {
if (globalDbg) {
cerr << "Simplex_handle()\n";
}
this->b = 0;
this->position = 0;
}
Simplex_handle(Bitmap_cubical_complex* b) {
if (globalDbg) {
cerr << "Simplex_handle(Bitmap_cubical_complex* b)\n";
}
this->b = b;
this->position = 0;
}
Simplex_handle(const Simplex_handle& org) : b(org.b) {
if (globalDbg) {
cerr << "Simplex_handle( const Simplex_handle& org )\n";
}
this->position = org.position;
}
Simplex_handle& operator=(const Simplex_handle& rhs) {
if (globalDbg) {
cerr << "Simplex_handle operator = \n";
}
this->position = rhs.position;
this->b = rhs.b;
return *this;
}
Simplex_handle(Bitmap_cubical_complex* b, Simplex_key position) {
if (globalDbg) {
cerr << "Simplex_handle(Bitmap_cubical_complex* b , Simplex_key position)\n";
cerr << "Position : " << position << endl;
}
this->b = b;
this->position = position;
}
friend class Bitmap_cubical_complex;
private:
Bitmap_cubical_complex* b;
Simplex_key position; //Assumption -- this field always keep the REAL position of simplex in the bitmap, no matter what keys have been.
//to deal with the keys, the class Bitmap_cubical_complex have extra vectors: keyAssociatedToSimplex and simplexAssociatedToKey
//that allow to move between actual cell and the key assigned to it.
};
//*********************************************************************************************************************************//
//Constructors
//*********************************************************************************************************************************//
//Over here we need to definie various input types. I am proposing the following ones:
//Perseus style
//H5 files? TODO
//binary files with little endiangs / big endians? TODO
//constructor from a vector of elements of a type T. TODO
/**
* Constructor form a Perseus-style file.
**/
Bitmap_cubical_complex(char* perseusStyleFile) : Bitmap_cubical_complex_base(perseusStyleFile) {
if (globalDbg) {
cerr << "Bitmap_cubical_complex( char* perseusStyleFile )\n";
}
std::vector< size_t > keyAssociatedToSimplex(this->totalNumberOfCells + 1);
std::vector< size_t > simplexAssociatedToKey(this->totalNumberOfCells + 1);
for (size_t i = 0; i != this->totalNumberOfCells; ++i) {
keyAssociatedToSimplex[i] = simplexAssociatedToKey[i] = i;
}
this->keyAssociatedToSimplex = keyAssociatedToSimplex;
this->simplexAssociatedToKey = simplexAssociatedToKey;
//we initialize this only once, in each constructor, when the bitmap is constructed. If the user decide to change some elements of the bitmap, then this procedure need
//to be called again.
this->initializeElementsOrderedAccordingToFiltration();
}
/**
* Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells in the following directions and vector of element of a type T
* with filtration on top dimensional cells.
**/
Bitmap_cubical_complex(std::vector dimensions, std::vector topDimensionalCells) : Bitmap_cubical_complex_base(dimensions, topDimensionalCells) {
std::vector< size_t > keyAssociatedToSimplex(this->totalNumberOfCells + 1);
std::vector< size_t > simplexAssociatedToKey(this->totalNumberOfCells + 1);
for (size_t i = 0; i != this->totalNumberOfCells; ++i) {
keyAssociatedToSimplex[i] = simplexAssociatedToKey[i] = i;
}
this->keyAssociatedToSimplex = keyAssociatedToSimplex;
this->simplexAssociatedToKey = simplexAssociatedToKey;
//we initialize this only once, in each constructor, when the bitmap is constructed. If the user decide to change some elements of the bitmap, then this procedure need
//to be called again.
this->initializeElementsOrderedAccordingToFiltration();
}
//*********************************************************************************************************************************//
//Other 'easy' functions
//*********************************************************************************************************************************//
/**
* Returns number of all cubes in the complex.
**/
size_t num_simplices()const {
return this->totalNumberOfCells;
}
/**
* Returns a Simplex_handle to a cube that do not exist in this complex.
**/
Simplex_handle null_simplex() {
return Simplex_handle(this, this->data.size());
}
/**
* Returns dimension of the complex.
**/
size_t dimension() {
return this->sizes.size();
}
/**
* Return dimension of a cell pointed by the Simplex_handle.
**/
size_t dimension(const Simplex_handle& sh) {
if (globalDbg) {
cerr << "int dimension(const Simplex_handle& sh)\n";
}
if (sh.position != this->data.size()) return sh.b->get_dimension_of_a_cell(sh.position);
return std::numeric_limits::max();
}
/**
* Return the filtration of a cell pointed by the Simplex_handle.
**/
T filtration(const Simplex_handle& sh) {
if (globalDbg) {
cerr << "T filtration(const Simplex_handle& sh)\n";
}
//Returns the filtration value of a simplex.
if (sh.position != this->data.size()) return sh.b->data[ sh.position ];
return INT_MAX;
}
/**
* Return a key which is not a key of any cube in the considered data structure.
**/
Simplex_key null_key() {
if (globalDbg) {
cerr << "Simplex_key null_key()\n";
}
return this->data.size();
}
/**
* Return the key of a cube pointed by the Simplex_handle.
**/
Simplex_key key(const Simplex_handle& sh) {
if (globalDbg) {
cerr << "Simplex_key key(const Simplex_handle& sh)\n";
}
return sh.b->keyAssociatedToSimplex[ sh.position ];
}
/**
* Return the Simplex_handle given the key of the cube.
**/
Simplex_handle simplex(Simplex_key key) {
if (globalDbg) {
cerr << "Simplex_handle simplex(Simplex_key key)\n";
}
return Simplex_handle(this, this->simplexAssociatedToKey[ key ]);
}
/**
* Assign key to a cube pointed by the Simplex_handle
**/
void assign_key(Simplex_handle& sh, Simplex_key key) {
if (globalDbg) {
cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
}
this->keyAssociatedToSimplex[sh.position] = key;
this->simplexAssociatedToKey[key] = sh.position;
}
/**
* Function called from a constructor. It is needed for Filtration_simplex_iterator to work.
**/
void initializeElementsOrderedAccordingToFiltration();
//*********************************************************************************************************************************//
//Iterators
//*********************************************************************************************************************************//
/**
* Boundary_simplex_iterator class allows iteration on boundary of each cube.
**/
class Boundary_simplex_range;
class Boundary_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {
//Iterator on the simplices belonging to the boundary of a simplex.
//value_type must be 'Simplex_handle'.
public:
Boundary_simplex_iterator(Simplex_handle& sh) : sh(sh) {
if (globalDbg) {
cerr << "Boundary_simplex_iterator( Simplex_handle& sh )\n";
}
this->position = 0;
this->boundaryElements = this->sh.b->get_boundary_of_a_cell(this->sh.position);
}
Boundary_simplex_iterator operator++() {
if (globalDbg) {
cerr << "Boundary_simplex_iterator operator++()\n";
}
++this->position;
return *this;
}
Boundary_simplex_iterator operator++(int) {
Boundary_simplex_iterator result = *this;
++(*this);
return result;
}
Boundary_simplex_iterator operator=(const Boundary_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "Boundary_simplex_iterator operator =\n";
}
this->sh = rhs.sh;
this->boundaryElements.clear();
this->boundaryElementsinsert(this->boundaryElements.end(), rhs.boundaryElements.begin(), rhs.boundaryElements.end());
}
bool operator==(const Boundary_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator ==\n";
}
if (this->position == rhs.position) {
if (this->boundaryElements.size() != rhs.boundaryElements.size())return false;
for (size_t i = 0; i != this->boundaryElements.size(); ++i) {
if (this->boundaryElements[i] != rhs.boundaryElements[i])return false;
}
return true;
}
return false;
}
bool operator!=(const Boundary_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator != \n";
}
return !(*this == rhs);
}
Simplex_handle operator*() {
if (globalDbg) {
cerr << "Simplex_handle operator*\n";
}
return Simplex_handle(this->sh.b, this->boundaryElements[this->position]);
}
friend class Boundary_simplex_range;
private:
Simplex_handle sh;
std::vector< size_t > boundaryElements;
size_t position;
};
/**
* Boundary_simplex_range class provides ranges for boundary iterators.
**/
class Boundary_simplex_range {
//Range giving access to the simplices in the boundary of a simplex.
//.begin() and .end() return type Boundary_simplex_iterator.
public:
Boundary_simplex_range(const Simplex_handle& sh) : sh(sh) { };
Boundary_simplex_iterator begin() {
if (globalDbg) {
cerr << "Boundary_simplex_iterator begin\n";
}
Boundary_simplex_iterator it(this->sh);
return it;
}
Boundary_simplex_iterator end() {
if (globalDbg) {
cerr << "Boundary_simplex_iterator end()\n";
}
Boundary_simplex_iterator it(this->sh);
it.position = it.boundaryElements.size();
return it;
}
private:
Simplex_handle sh;
};
/**
* Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration. Secondary criteria for filtration are:
* (1) Dimension of a cube (lower dimensional comes first).
* (2) Position in the data structure (the ones that are earlies in the data structure comes first).
**/
class Filtration_simplex_range;
class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {
//Iterator over all simplices of the complex in the order of the indexing scheme.
//'value_type' must be 'Simplex_handle'.
public:
Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) { };
Filtration_simplex_iterator() : b(NULL) { };
Filtration_simplex_iterator operator++() {
if (globalDbg) {
cerr << "Filtration_simplex_iterator operator++\n";
}
++this->position;
return (*this);
}
Filtration_simplex_iterator operator++(int) {
Filtration_simplex_iterator result = *this;
++(*this);
return result;
}
Filtration_simplex_iterator operator=(const Filtration_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "Filtration_simplex_iterator operator =\n";
}
this->b = rhs.b;
this->position = rhs.position;
}
bool operator==(const Filtration_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";
}
if (this->position == rhs.position) {
return true;
}
return false;
}
bool operator!=(const Filtration_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n";
}
return !(*this == rhs);
}
Simplex_handle operator*() {
if (globalDbg) {
cerr << "Simplex_handle operator*()\n";
}
return Simplex_handle(this->b, this->b->elementsOrderedAccordingToFiltration[ this->position ]);
}
friend class Filtration_simplex_range;
private:
Bitmap_cubical_complex* b;
size_t position;
};
/**
* Filtration_simplex_range provides the ranges for Filtration_simplex_iterator.
**/
class Filtration_simplex_range {
//Range over the simplices of the complex in the order of the filtration.
//.begin() and .end() return type Filtration_simplex_iterator.
public:
Filtration_simplex_range(Bitmap_cubical_complex* b) : b(b) { };
Filtration_simplex_iterator begin() {
if (globalDbg) {
cerr << "Filtration_simplex_iterator begin() \n";
}
return Filtration_simplex_iterator(this->b);
}
Filtration_simplex_iterator end() {
if (globalDbg) {
cerr << "Filtration_simplex_iterator end()\n";
}
Filtration_simplex_iterator it(this->b);
it.position = this->b->elementsOrderedAccordingToFiltration.size();
return it;
}
private:
Bitmap_cubical_complex* b;
};
//*********************************************************************************************************************************//
//Methods to access iterators from the container:
/**
* boundary_simplex_range creates an object of a Boundary_simplex_range class that provides ranges for the Boundary_simplex_iterator.
**/
Boundary_simplex_range boundary_simplex_range(Simplex_handle& sh) {
if (globalDbg) {
cerr << "Boundary_simplex_range boundary_simplex_range(Simplex_handle& sh)\n";
}
//Returns a range giving access to all simplices of the boundary of a simplex, i.e. the set of codimension 1 subsimplices of the Simplex.
return Boundary_simplex_range(sh);
}
/**
* filtration_simplex_range creates an object of a Filtration_simplex_range class that provides ranges for the Filtration_simplex_iterator.
**/
Filtration_simplex_range filtration_simplex_range() {
if (globalDbg) {
cerr << "Filtration_simplex_range filtration_simplex_range()\n";
}
//Returns a range over the simplices of the complex in the order of the filtration
return Filtration_simplex_range(this);
}
//*********************************************************************************************************************************//
//*********************************************************************************************************************************//
//Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are there.
//TODO -- the file IndexingTag.h in the Gudhi library contains an empty structure, so I understand that this is something that was planned (for simplicial maps?)
//but was never finished. The only idea I have here is to use the same empty structure from IndexingTag.h file, but only if the compiler needs it. If the compiler
//do not need it, then I would rather not add here elements which I do not understand.
//typedef Indexing_tag
/**
* Function needed for compatibility with Gudhi. Not useful for other purposes.
**/
std::pair endpoints(Simplex_handle sh) {
std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh.position);
if (globalDbg) {
cerr << "std::pair endpoints( Simplex_handle sh )\n";
cerr << "bdry.size() : " << bdry.size() << endl;
}
//this method returns two first elements from the boundary of sh.
if (bdry.size() < 2)throw ("Error in endpoints in Bitmap_cubical_complex class. The cell for which this method was called have less than two elements in the boundary.");
return std::make_pair(Simplex_handle(this, bdry[0]), Simplex_handle(this, bdry[1]));
}
/**
* Class needed for compatibility with Gudhi. Not useful for other purposes.
**/
class Skeleton_simplex_range;
class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > {
//Iterator over all simplices of the complex in the order of the indexing scheme.
//'value_type' must be 'Simplex_handle'.
public:
Skeleton_simplex_iterator(Bitmap_cubical_complex* b, size_t d) : b(b), dimension(d) {
if (globalDbg) {
cerr << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , size_t d )\n";
}
//find the position of the first simplex of a dimension d
this->position = 0;
while ((this->position != b->data.size()) && (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) {
++this->position;
}
};
Skeleton_simplex_iterator() : b(NULL), dimension(0) { };
Skeleton_simplex_iterator operator++() {
if (globalDbg) {
cerr << "Skeleton_simplex_iterator operator++()\n";
}
//increment the position as long as you did not get to the next element of the dimension dimension.
++this->position;
while ((this->position != this->b->data.size()) && (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) {
++this->position;
}
return (*this);
}
Skeleton_simplex_iterator operator++(int) {
Skeleton_simplex_iterator result = *this;
++(*this);
return result;
}
Skeleton_simplex_iterator operator=(const Skeleton_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "Skeleton_simplex_iterator operator =\n";
}
this->b = rhs.b;
this->position = rhs.position;
}
bool operator==(const Skeleton_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator ==\n";
}
if (this->position == rhs.position) {
return true;
}
return false;
}
bool operator!=(const Skeleton_simplex_iterator& rhs) {
if (globalDbg) {
cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n";
}
return !(*this == rhs);
}
Simplex_handle operator*() {
if (globalDbg) {
cerr << "Simplex_handle operator*() \n";
}
return Simplex_handle(this->b, this->position);
}
friend class Skeleton_simplex_range;
private:
Bitmap_cubical_complex* b;
size_t position;
int dimension;
};
/**
* Class needed for compatibility with Gudhi. Not useful for other purposes.
**/
class Skeleton_simplex_range {
//Range over the simplices of the complex in the order of the filtration.
//.begin() and .end() return type Filtration_simplex_iterator.
public:
Skeleton_simplex_range(Bitmap_cubical_complex* b, int dimension) : b(b), dimension(dimension) { };
Skeleton_simplex_iterator begin() {
if (globalDbg) {
cerr << "Skeleton_simplex_iterator begin()\n";
}
return Skeleton_simplex_iterator(this->b, this->dimension);
}
Skeleton_simplex_iterator end() {
if (globalDbg) {
cerr << "Skeleton_simplex_iterator end()\n";
}
Skeleton_simplex_iterator it(this->b, this->dimension);
it.position = this->b->data.size();
return it;
}
private:
Bitmap_cubical_complex* b;
int dimension;
};
/**
* Function needed for compatibility with Gudhi. Not useful for other purposes.
**/
Skeleton_simplex_range skeleton_simplex_range(int dimension) {
if (globalDbg) {
cerr << "Skeleton_simplex_range skeleton_simplex_range( int dimension )\n";
}
return Skeleton_simplex_range(this, dimension);
}
//*********************************************************************************************************************************//
//functions used for debugging:
/**
* Function used for debugging purposes.
**/
void printKeyAssociatedToSimplex() {
for (size_t i = 0; i != this->data.size(); ++i) {
cerr << i << " -> " << this->simplexAssociatedToKey[i] << endl;
}
}
/**
* Function used for debugging purposes.
**/
size_t printRealPosition(const Simplex_handle& sh) {
return sh.position;
}
private:
std::vector< size_t > keyAssociatedToSimplex;
std::vector< size_t > simplexAssociatedToKey;
std::vector< size_t > elementsOrderedAccordingToFiltration; //needed by Filtration_simplex_iterator. If this iterator is not used, this field is not initialized.
}; //Bitmap_cubical_complex
template
bool compareElementsForElementsOrderedAccordingToFiltration(const std::pair< size_t, std::pair< T, char > >& f, const std::pair< size_t, std::pair< T, char > >& s) {
if (globalDbg) {
cerr << "ompareElementsForElementsOrderedAccordingToFiltration\n";
}
if (f.second.first < s.second.first) {
return true;
} else {
if (f.second.first > s.second.first) {
return false;
} else {
//in this case f.second.first == s.second.first, and we use dimension to compare:
if (f.second.second < s.second.second) {
return true;
} else {
if (f.second.second > s.second.second) {
return false;
} else {
//in this case, both the filtration value and the dimensions for those cells are the same. Since it may be nice to have a stable sorting procedure, in this case, we compare positions in the bitmap:
return ( f.first < s.first);
}
}
}
}
}
template
void Bitmap_cubical_complex::initializeElementsOrderedAccordingToFiltration() {
if (globalDbg) {
cerr << "void Bitmap_cubical_complex::initializeElementsOrderedAccordingToFiltration() \n";
}
//( position , (filtration , dimension) )
std::vector< std::pair< size_t, std::pair< T, char > > > dataOfElementsFromBitmap(this->data.size());
for (size_t i = 0; i != this->data.size(); ++i) {
//TODO -- this can be optimized by having a counter here. We do not need to re-compute the dimension for every cell from scratch
dataOfElementsFromBitmap[i] = std::make_pair(i, std::make_pair(this->data[i], this->get_dimension_of_a_cell(i)));
}
std::sort(dataOfElementsFromBitmap.begin(), dataOfElementsFromBitmap.end(), compareElementsForElementsOrderedAccordingToFiltration);
std::vector< size_t > elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap(this->data.size());
for (size_t i = 0; i != dataOfElementsFromBitmap.size(); ++i) {
elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap[i] = dataOfElementsFromBitmap[i].first;
}
this->elementsOrderedAccordingToFiltration = elementsOfBitmapOrderedAccordingToFiltrationThenAccordingToDimensionThenAccordingToPositionInBitmap;
}
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
//****************************************************************************************************************//
#endif // BITMAP_CUBICAL_COMPLEX_H_