#pragma once #include #include "Bitmap_cubical_complex_base.h" using namespace std; namespace Gudhi { namespace Cubical_complex { //in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic boundary conditions). But, we set up the iterators and the procedures //to compute boundary and coboundary in the way that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all the cells on the //right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base template class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base { public: //constructors that take an extra parameter: Bitmap_cubical_complex_periodic_boundary_conditions_base(){}; Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseusStyleFile ); Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); //overwritten methods co compute boundary and coboundary virtual std::vector< size_t > get_boundary_of_a_cell( size_t cell )const; std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const; //inline unsigned get_dimension_of_a_cell( size_t cell )const; protected: std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed; void set_up_containers( const std::vector& sizes ) { unsigned multiplier = 1; for ( size_t i = 0 ; i != sizes.size() ; ++i ) { this->sizes.push_back(sizes[i]); this->multipliers.push_back(multiplier); if ( directions_in_which_periodic_b_cond_are_to_be_imposed[i] ) { multiplier *= 2*sizes[i]; } else { multiplier *= 2*sizes[i]+1; } } //std::reverse( this->sizes.begin() , this->sizes.end() ); this->data = std::vector(multiplier,std::numeric_limits::max()); this->total_number_of_cells = multiplier; } Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ); Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ); void construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ); }; template void Bitmap_cubical_complex_periodic_boundary_conditions_base::construct_complex_based_on_top_dimensional_cells( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ) { this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; this->set_up_containers( dimensions ); size_t i = 0; for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) { *it = topDimensionalCells[i]; ++i; } this->impose_lower_star_filtration(); } template Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ) { this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; this->set_up_containers( sizes ); } template Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseus_style_file ) { //for Perseus style files: bool dbg = false; ifstream inFiltration; inFiltration.open( perseus_style_file ); unsigned dimensionOfData; inFiltration >> dimensionOfData; this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensionOfData , false ); std::vector sizes; sizes.reserve( dimensionOfData ); for ( size_t i = 0 ; i != dimensionOfData ; ++i ) { int size_in_this_dimension; inFiltration >> size_in_this_dimension; if ( size_in_this_dimension < 0 ) { this->directions_in_which_periodic_b_cond_are_to_be_imposed[i] = true; } sizes.push_back( abs(size_in_this_dimension) ); } this->set_up_containers( sizes ); typename Bitmap_cubical_complex_periodic_boundary_conditions_base::Top_dimensional_cells_iterator it(*this); it = this->top_dimensional_cells_begin(); while ( !inFiltration.eof() ) { double filtrationLevel; inFiltration >> filtrationLevel; if ( inFiltration.eof() )break; if ( dbg ) { cerr << "Cell of an index : " << it.compute_index_in_bitmap() << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) << " get the value : " << filtrationLevel << endl; } *it = filtrationLevel; ++it; } inFiltration.close(); this->impose_lower_star_filtration(); /* char* filename = (char*)perseus_style_file; //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float"; ifstream file( filename , ios::binary | ios::ate ); unsigned realSizeOfFile = file.tellg(); file.close(); realSizeOfFile = realSizeOfFile/sizeof(T); unsigned w, h, d; w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) )); T* slice = new T[w*h*d]; if (slice == NULL) { cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n"; exit(EXIT_FAILURE); } FILE* fp; if ((fp=fopen( filename, "rb" )) == NULL ) { cerr << "Cannot open the file: " << filename << ". The program will now terminate \n"; exit(1); } fread( slice,4,w*h*d,fp ); fclose(fp); std::vector data(slice,slice+w*h*d); delete[] slice; std::vector< unsigned > sizes; sizes.push_back(w); sizes.push_back(w); sizes.push_back(w); std::vector< bool > directions; directions.push_back( true ); directions.push_back( true ); directions.push_back( true ); Bitmap_cubical_complex_periodic_boundary_conditions_base b( sizes, data, directions ); *this = b; */ } template Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector sizes ) { this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( sizes.size() , false ); this->set_up_containers( sizes ); } template Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells ) { std::vector directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector( dimensions.size() , false ); this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed ); } template Bitmap_cubical_complex_periodic_boundary_conditions_base::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector dimensions , std::vector topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed ) { this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed ); } //***********************Methods************************// template std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base::get_boundary_of_a_cell( size_t cell )const { bool dbg = false; if ( dbg ){cerr << "Computations of boundary of a cell : " << cell << endl;} std::vector< size_t > boundary_elements; size_t cell1 = cell; for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) { unsigned position = cell1/this->multipliers[i-1]; //this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. if ( position%2 == 1 ) { //if there are no periodic boundary conditions in this direction, we do not have to do anything. if ( !directions_in_which_periodic_b_cond_are_to_be_imposed[i-1] ) { //cerr << "A\n"; boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} } else { //in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. if ( position != 2*this->sizes[ i-1 ]-1 ) { //cerr << "B\n"; boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); boundary_elements.push_back( cell + this->multipliers[ i-1 ] ); if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";} } else { //cerr << "C\n"; boundary_elements.push_back( cell - this->multipliers[ i-1 ] ); boundary_elements.push_back( cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] ); if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] << " ";} } } } cell1 = cell1%this->multipliers[i-1]; } return boundary_elements; } template std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base::get_coboundary_of_a_cell( size_t cell )const { std::vector counter = this->compute_counter_for_given_cell( cell ); std::vector< size_t > coboundary_elements; size_t cell1 = cell; for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) { unsigned position = cell1/this->multipliers[i-1]; //if the cell has zero length in this direction, then it will have cbd in this direction. if ( position%2 == 0 ) { if ( !this->directions_in_which_periodic_b_cond_are_to_be_imposed[i-1] ) { //no periodic boundary conditions in this direction if ( (counter[i-1] != 0) && (cell > this->multipliers[i-1]) ) { coboundary_elements.push_back( cell - this->multipliers[i-1] ); } if ( (counter[i-1] != 2*this->sizes[i-1]) && (cell + this->multipliers[i-1] < this->data.size()) ) { coboundary_elements.push_back( cell + this->multipliers[i-1] ); } } else { //we want to have periodic boundary conditions in this direction if ( counter[i-1] != 0 ) { coboundary_elements.push_back( cell - this->multipliers[i-1] ); coboundary_elements.push_back( cell + this->multipliers[i-1] ); } else { //in this case counter[i-1] == 0. coboundary_elements.push_back( cell + this->multipliers[i-1] ); coboundary_elements.push_back( cell + (2*this->sizes[ i-1 ]-1)*this->multipliers[i-1] ); } } } cell1 = cell1%this->multipliers[i-1]; } return coboundary_elements; } }//Cubical_complex }//namespace Gudhi