diff options
author | pdlotko <pdlotko@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-11-27 14:16:44 +0000 |
---|---|---|
committer | pdlotko <pdlotko@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-11-27 14:16:44 +0000 |
commit | 39ffc7af64543f0b2ce17487771afbd010a97349 (patch) | |
tree | 3f2736153e927cf6bc23f8c8bf0a0f868bf2cc6c /src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h | |
parent | 606115b399b010b8ebfbf13f6a8ca2b1c1bf17cd (diff) |
Adding changes to the cubcial complex class according to Marc and Vincent's comments.
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bitmap@930 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: cb9134a5c18fae3e2e63bbaf8329dd168d08d3b3
Diffstat (limited to 'src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h')
-rw-r--r-- | src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h | 1352 |
1 files changed, 759 insertions, 593 deletions
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index d9c91832..2c2bd481 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -1,593 +1,759 @@ -/* 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 <http://www.gnu.org/licenses/>.
- */
-
-#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_
-#define BITMAP_CUBICAL_COMPLEX_BASE_H_
-
-#include <gudhi/counter.h>
-
-#include <iostream>
-#include <string>
-#include <vector>
-#include <cstdlib>
-#include <climits>
-#include <fstream>
-#include <algorithm>
-#include <iterator>
-
-/**
- * This is a class implementing a basic bitmap data structure to store cubical complexes. It implements only the most basic subroutines.
- * The idea of the bitmap is the following. Our aim is to have a memory efficient data structure to store d-dimensional cubical complex C being a cubical decomposition
- * of a rectangular region of a space. This is achieved by storing C as a vector of bits (this is where the name 'bitmap' came from). Each cell is represented by a single
- * bit (in case of black and white bitmaps, or by a single element of a type T (here T is a filtration type of a bitmap, typically a double). All the informations needed for homology and
- * persistent homology computations (like dimension of a cell, boundary and coboundary elements of a cell, are then obtained from the position of the element in C.
- */
-template <typename T>
-class Bitmap_cubical_complex_base {
- public:
- /**
- * There are a few constructors of a Bitmap_cubical_complex_base class. First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal to the number of elements in the
- * input vector and size in the i-th dimension equal the number in the position i-of the input vector.
- */
- Bitmap_cubical_complex_base(std::vector<unsigned> sizes_);
- /**
- * The second constructor takes as a input a Perseus style file. For more details, please consult the documentations of Perseus software as well as examples attached to this
- * implementation.
- **/
- Bitmap_cubical_complex_base(char* perseusStyleFile_);
- /**
- * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one) together with vector of filtration values of top dimensional cells.
- **/
- Bitmap_cubical_complex_base(std::vector<unsigned> dimensions_, std::vector<T> topDimensionalCells_);
-
- /**
- * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell and get_cell_data are the basic functions that compute boundary / coboundary / dimension and the filtration
- * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a non-negative integer, indicating a position of a cube in the data structure.
- * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to the positions of (co)boundary element of the input cell.
- */
- inline std::vector< size_t > get_boundary_of_a_cell(size_t cell_);
- /**
- * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell and get_cell_data are the basic functions that compute boundary / coboundary / dimension and the filtration
- * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a non-negative integer, indicating a position of a cube in the data structure.
- * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to the positions of (co)boundary element of the input cell.
- **/
- inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell_);
- /**
- * In the case of get_dimension_of_a_cell function, the output is a non-negative integer indicating the dimension of a cell.
- **/
- inline unsigned get_dimension_of_a_cell(size_t cell_);
- /**
- * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position.
- **/
- inline T& get_cell_data(size_t cell_);
-
-
- /**
- * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. Then, there are a few ways one can pick the filtration of lower dimensional
- * cells. The most typical one is by so called lower star filtration. This function is always called by any constructor which takes the top dimensional cells. If you use such a constructor,
- * then there is no need to call this function. Call it only if you are putting the filtration of the cells by your own (for instance by using topDimensionalCellsIterator).
- **/
- void impose_lower_star_filtration(); // assume that top dimensional cells are already set.
-
- /**
- * Returns dimension of a complex.
- **/
- inline unsigned dimension() {
- return sizes.size();
- }
-
- /**
- * Returns number of all cubes in the data structure.
- **/
- inline unsigned size_of_bitmap() {
- return this->data.size();
- }
-
- /**
- * Writing to stream operator.
- **/
- template <typename K>
- friend std::ostream& operator<<(std::ostream & os_, const Bitmap_cubical_complex_base<K>& b_);
-
- // ITERATORS
-
- /**
- * Iterator through all cells in the complex (in order they appear in the structure -- i.e. in lexicographical order).
- **/
- typedef typename std::vector< T >::iterator all_cells_iterator;
-
- all_cells_iterator all_cells_begin()const {
- return this->data.begin();
- }
-
- all_cells_iterator all_cells_end()const {
- return this->data.end();
- }
-
-
- typedef typename std::vector< T >::const_iterator all_cells_const_iterator;
-
- all_cells_const_iterator all_cells_const_begin()const {
- return this->data.begin();
- }
-
- all_cells_const_iterator all_cells_const_end()const {
- return this->data.end();
- }
-
- /**
- * Iterator through top dimensional cells of the complex. The cells appear in order they are stored in the structure (i.e. in lexicographical order)
- **/
- class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, double > {
- public:
- Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b_) : b(b_) {
- for (size_t i = 0; i != b_.dimension(); ++i) {
- this->counter.push_back(0);
- }
- }
-
- Top_dimensional_cells_iterator operator++() {
- // first find first element of the counter that can be increased:
- size_t dim = 0;
- while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1))++dim;
-
- if (dim != this->b.dimension()) {
- ++this->counter[dim];
- for (size_t i = 0; i != dim; ++i) {
- this->counter[i] = 0;
- }
- } else {
- ++this->counter[0];
- }
- return *this;
- }
-
- Top_dimensional_cells_iterator operator++(int) {
- Top_dimensional_cells_iterator result = *this;
- ++(*this);
- return result;
- }
-
- Top_dimensional_cells_iterator operator=(const Top_dimensional_cells_iterator& rhs_) {
- this->counter = rhs_.counter;
- this->b = rhs_.b;
- return *this;
- }
-
- bool operator==(const Top_dimensional_cells_iterator& rhs_) {
- if (&this->b != &rhs_.b)return false;
- if (this->counter.size() != rhs_.counter.size())return false;
- for (size_t i = 0; i != this->counter.size(); ++i) {
- if (this->counter[i] != rhs_.counter[i])return false;
- }
- return true;
- }
-
- bool operator!=(const Top_dimensional_cells_iterator& rhs_) {
- return !(*this == rhs_);
- }
-
- T& operator*() {
- // given the counter, compute the index in the array and return this element.
- unsigned index = 0;
- for (size_t i = 0; i != this->counter.size(); ++i) {
- index += (2 * this->counter[i] + 1) * this->b.multipliers[i];
- }
- return this->b.data[index];
- }
-
- size_t computeIndexInBitmap() {
- size_t index = 0;
- for (size_t i = 0; i != this->counter.size(); ++i) {
- index += (2 * this->counter[i] + 1) * this->b.multipliers[i];
- }
- return index;
- }
-
- void printCounter() {
- for (size_t i = 0; i != this->counter.size(); ++i) {
- std::cout << this->counter[i] << " ";
- }
- }
- friend class Bitmap_cubical_complex_base;
- protected:
- std::vector< unsigned > counter;
- Bitmap_cubical_complex_base& b;
- };
-
- Top_dimensional_cells_iterator top_dimensional_cells_begin() {
- Top_dimensional_cells_iterator a(*this);
- return a;
- }
-
- Top_dimensional_cells_iterator top_dimensional_cells_end() {
- Top_dimensional_cells_iterator a(*this);
- for (size_t i = 0; i != this->dimension(); ++i) {
- a.counter[i] = this->sizes[i] - 1;
- }
- a.counter[0]++;
- return a;
- }
-
-
- //****************************************************************************************************************//
- //****************************************************************************************************************//
- //****************************************************************************************************************//
- //****************************************************************************************************************//
-
-
- //****************************************************************************************************************//
- //****************************************************************************************************************//
- //****************************************************************************************************************//
- //****************************************************************************************************************//
-
- protected:
- std::vector<unsigned> sizes;
- std::vector<unsigned> multipliers;
- std::vector<T> data;
- size_t totalNumberOfCells;
-
- void set_up_containers(std::vector<unsigned> sizes_) {
- unsigned multiplier = 1;
- for (size_t i = 0; i != sizes_.size(); ++i) {
- this->sizes.push_back(sizes_[i]);
- this->multipliers.push_back(multiplier);
- // multiplier *= 2*(sizes[i]+1)+1;
- multiplier *= 2 * sizes_[i] + 1;
- }
- // std::reverse( this->sizes.begin() , this->sizes.end() );
- std::vector<T> data(multiplier);
- std::fill(data.begin(), data.end(), INT_MAX);
- this->totalNumberOfCells = multiplier;
- this->data = data;
- }
-
- size_t compute_position_in_bitmap(std::vector< int > counter_) {
- size_t position = 0;
- for (size_t i = 0; i != this->multipliers.size(); ++i) {
- position += this->multipliers[i] * counter_[i];
- }
- return position;
- }
-
- std::vector<unsigned> compute_counter_for_given_cell(size_t cell_) {
- std::vector<unsigned> counter;
- for (size_t dim = this->sizes.size(); dim != 0; --dim) {
- counter.push_back(cell_ / this->multipliers[dim - 1]);
- cell_ = cell_ % this->multipliers[dim - 1];
- }
- std::reverse(counter.begin(), counter.end());
- return counter;
- }
-
- std::vector< size_t > generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions(std::vector< bool > directionsForPeriodicBCond_);
-};
-
-template <typename K>
-std::ostream& operator<<(std::ostream & out_, const Bitmap_cubical_complex_base<K>& b_) {
- // for ( typename bitmap<K>::all_cells_const_iterator it = b.all_cells_const_begin() ;
- // it != b.all_cells_const_end() ; ++it )
- for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator it = b_.all_cells_const_begin();
- it != b_.all_cells_const_end(); ++it) {
- out_ << *it << " ";
- }
- return out_;
-}
-
-template <typename T>
-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(std::vector<unsigned> sizes_) {
- this->set_up_containers(sizes_);
-}
-
-template <typename T>
-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(std::vector<unsigned> sizesInFollowingDirections_,
- std::vector<T> topDimensionalCells_) {
- this->set_up_containers(sizesInFollowingDirections_);
-
- size_t numberOfTopDimensionalElements = 1;
- for (size_t i = 0; i != sizesInFollowingDirections_.size(); ++i) {
- numberOfTopDimensionalElements *= sizesInFollowingDirections_[i];
- }
- if (numberOfTopDimensionalElements != topDimensionalCells_.size()) {
- std::cerr << "Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizesInFollowingDirections_ , "
- "std::vector<float> topDimensionalCells_ ). Number of top dimensional elements that follow from "
- "sizesInFollowingDirections vector is different than the size of topDimensionalCells vector." << std::endl;
- throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizesInFollowingDirections_ , "
- "std::vector<float> topDimensionalCells_ ). Number of top dimensional elements that follow from "
- "sizesInFollowingDirections vector is different than the size of topDimensionalCells vector.");
- }
-
- Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
- size_t index = 0;
- for (it = this->top_dimensional_cells_begin(); it != this->top_dimensional_cells_end(); ++it) {
- (*it) = topDimensionalCells_[index];
- ++index;
- }
- this->impose_lower_star_filtration();
-}
-
-template <typename T>
-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(char* perseusStyleFile_) {
- bool dbg = false;
- std::ifstream inFiltration, inIds;
- inFiltration.open(perseusStyleFile_);
- unsigned dimensionOfData;
- inFiltration >> dimensionOfData;
-
- if (dbg) {
- std::cerr << "dimensionOfData : " << dimensionOfData << std::endl;
- }
-
- std::vector<unsigned> sizes;
- for (size_t i = 0; i != dimensionOfData; ++i) {
- int sizeInThisDimension;
- inFiltration >> sizeInThisDimension;
- sizeInThisDimension = abs(sizeInThisDimension);
- sizes.push_back(sizeInThisDimension);
- if (dbg) {
- std::cerr << "sizeInThisDimension : " << sizeInThisDimension << std::endl;
- }
- }
- this->set_up_containers(sizes);
-
- Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
- it = this->top_dimensional_cells_begin();
-
- // TODO(Pawel Dlotko): Over here we also need to read id's of cell and put them to bitmapElement structure!
- while (!inFiltration.eof()) {
- double filtrationLevel;
- inFiltration >> filtrationLevel;
- if (dbg) {
- std::cerr << "Cell of an index : " << it.computeIndexInBitmap() << " and dimension: " <<
- this->get_dimension_of_a_cell(it.computeIndexInBitmap()) << " get the value : " <<
- filtrationLevel << std::endl;
- }
- *it = filtrationLevel;
- ++it;
- }
- inFiltration.close();
- this->impose_lower_star_filtration();
-}
-
-template <typename T>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell_) {
- bool bdg = false;
- // First of all, we need to take the list of coordinates in which the cell has nonzero length.
- // We do it by using modified version to compute dimension of a cell:
- std::vector< unsigned > dimensionsInWhichCellHasNonzeroLength;
- unsigned dimension = 0;
- size_t cell1 = cell_;
- for (size_t i = this->multipliers.size(); i != 0; --i) {
- unsigned position = cell1 / multipliers[i - 1];
- if (position % 2 == 1) {
- dimensionsInWhichCellHasNonzeroLength.push_back(i - 1);
- dimension++;
- }
- cell1 = cell1 % multipliers[i - 1];
- }
-
- if (bdg) {
- std::cerr << "dimensionsInWhichCellHasNonzeroLength : \n";
- for (size_t i = 0; i != dimensionsInWhichCellHasNonzeroLength.size(); ++i) {
- std::cerr << dimensionsInWhichCellHasNonzeroLength[i] << std::endl;
- }
- getchar();
- }
-
- std::vector< size_t > boundaryElements;
- if (dimensionsInWhichCellHasNonzeroLength.size() == 0)return boundaryElements;
- for (size_t i = 0; i != dimensionsInWhichCellHasNonzeroLength.size(); ++i) {
- boundaryElements.push_back(cell_ - multipliers[ dimensionsInWhichCellHasNonzeroLength[i] ]);
- boundaryElements.push_back(cell_ + multipliers[ dimensionsInWhichCellHasNonzeroLength[i] ]);
-
- if (bdg) std::cerr << "multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<
- multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;
- if (bdg) std::cerr << "cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<
- cell_ - multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;
- if (bdg) std::cerr << "cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] : " <<
- cell_ + multipliers[dimensionsInWhichCellHasNonzeroLength[i]] << std::endl;
- }
- return boundaryElements;
-}
-
-template <typename T>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell_) {
- bool bdg = false;
- // First of all, we need to take the list of coordinates in which the cell has nonzero length.
- // We do it by using modified version to compute dimension of a cell:
- std::vector< unsigned > dimensionsInWhichCellHasZeroLength;
- unsigned dimension = 0;
- size_t cell1 = cell_;
- for (size_t i = this->multipliers.size(); i != 0; --i) {
- unsigned position = cell1 / multipliers[i - 1];
- if (position % 2 == 0) {
- dimensionsInWhichCellHasZeroLength.push_back(i - 1);
- dimension++;
- }
- cell1 = cell1 % multipliers[i - 1];
- }
-
- std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell_);
- // reverse(counter.begin() , counter.end());
-
- if (bdg) {
- std::cerr << "dimensionsInWhichCellHasZeroLength : \n";
- for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {
- std::cerr << dimensionsInWhichCellHasZeroLength[i] << std::endl;
- }
- std::cerr << "\n counter : " << std::endl;
- for (size_t i = 0; i != counter.size(); ++i) {
- std::cerr << counter[i] << std::endl;
- }
- getchar();
- }
-
- std::vector< size_t > coboundaryElements;
- if (dimensionsInWhichCellHasZeroLength.size() == 0)return coboundaryElements;
- for (size_t i = 0; i != dimensionsInWhichCellHasZeroLength.size(); ++i) {
- if (bdg) {
- std::cerr << "Dimension : " << i << std::endl;
- if (counter[dimensionsInWhichCellHasZeroLength[i]] == 0) {
- std::cerr << "In dimension : " << i <<
- " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";
- }
- if (counter[dimensionsInWhichCellHasZeroLength[i]] == 2 * this->sizes[dimensionsInWhichCellHasZeroLength[i]]) {
- std::cerr << "In dimension : " << i <<
- " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n";
- }
- }
-
-
- if ((cell_ > multipliers[dimensionsInWhichCellHasZeroLength[i]]) &&
- (counter[dimensionsInWhichCellHasZeroLength[i]] != 0)) {
- // if ( counter[dimensionsInWhichCellHasZeroLength[i]] != 0 )
- if (bdg)std::cerr << "Subtracting : " << cell_ - multipliers[dimensionsInWhichCellHasZeroLength[i]] << std::endl;
- coboundaryElements.push_back(cell_ - multipliers[dimensionsInWhichCellHasZeroLength[i]]);
- }
- if ((cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]] < this->data.size()) &&
- (counter[dimensionsInWhichCellHasZeroLength[i]] != 2 * this->sizes[dimensionsInWhichCellHasZeroLength[i]])) {
- // if ( counter[dimensionsInWhichCellHasZeroLength[i]] != 2*this->sizes[dimensionsInWhichCellHasZeroLength[i]] )
- coboundaryElements.push_back(cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]]);
- if (bdg)std::cerr << "Adding : " << cell_ + multipliers[dimensionsInWhichCellHasZeroLength[i]] << std::endl;
- }
- }
- return coboundaryElements;
-}
-
-template <typename T>
-unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell_) {
- bool dbg = false;
- if (dbg)std::cerr << "\n\n\n Computing position o a cell of an index : " << cell_ << std::endl;
- unsigned dimension = 0;
- for (size_t i = this->multipliers.size(); i != 0; --i) {
- unsigned position = cell_ / multipliers[i - 1];
-
- if (dbg)std::cerr << "i-1 :" << i - 1 << std::endl;
- if (dbg)std::cerr << "cell_ : " << cell_ << std::endl;
- if (dbg)std::cerr << "position : " << position << std::endl;
- if (dbg)std::cerr << "multipliers[" << i - 1 << "] = " << multipliers[i - 1] << std::endl;
- if (dbg)getchar();
-
- if (position % 2 == 1) {
- if (dbg)std::cerr << "Nonzero length in this direction \n";
- dimension++;
- }
- cell_ = cell_ % multipliers[i - 1];
- }
- return dimension;
-}
-
-template <typename T>
-T& Bitmap_cubical_complex_base<T>::get_cell_data(size_t cell_) {
- return this->data[cell_];
-}
-
-template <typename T>
-void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() {
- bool dbg = false;
-
- // this vector will be used to check which elements have already been taken care of in imposing lower star filtration:
- std::vector<bool> isThisCellConsidered(this->data.size(), false);
-
- std::vector<size_t> indicesToConsider;
- // we assume here that we already have a filtration on the top dimensional cells
- // and we have to extend it to lower ones.
- typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
- for (it = this->top_dimensional_cells_begin(); it != this->top_dimensional_cells_end(); ++it) {
- indicesToConsider.push_back(it.computeIndexInBitmap());
- }
-
- while (indicesToConsider.size()) {
- if (dbg) {
- std::cerr << "indicesToConsider in this iteration \n";
- for (size_t i = 0; i != indicesToConsider.size(); ++i) {
- std::cout << indicesToConsider[i] << " ";
- }
- getchar();
- }
- std::vector<size_t> newIndicesToConsider;
- for (size_t i = 0; i != indicesToConsider.size(); ++i) {
- std::vector<size_t> bd = this->get_boundary_of_a_cell(indicesToConsider[i]);
- for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) {
- if (this->data[ bd[boundaryIt] ] > this->data[ indicesToConsider[i] ]) {
- this->data[ bd[boundaryIt] ] = this->data[ indicesToConsider[i] ];
- }
- if (isThisCellConsidered[ bd[boundaryIt] ] == false) {
- newIndicesToConsider.push_back(bd[boundaryIt]);
- isThisCellConsidered[ bd[boundaryIt] ] = true;
- }
- }
- }
- indicesToConsider.swap(newIndicesToConsider);
- }
-}
-
-template <typename T>
-std::vector< size_t >
-Bitmap_cubical_complex_base<T>::generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions(std::vector< bool > directionsForPeriodicBCond_) {
- bool dbg = false;
- if (this->sizes.size() != directionsForPeriodicBCond_.size())
- throw ("directionsForPeriodicBCond_ vector size is different from the size of the bitmap. "
- "The program will now terminate \n");
-
- std::vector<int> sizes(this->sizes.size());
- for (size_t i = 0; i != this->sizes.size(); ++i)sizes[i] = 2 * this->sizes[i];
-
- counter c(sizes);
-
- std::vector< size_t > result;
-
- for (size_t i = 0; i != this->data.size(); ++i) {
- size_t position;
- if (!c.isFinal()) {
- position = i;
- // result.push_back( i );
- } else {
- std::vector< bool > finals = c.directionsOfFinals();
- bool jumpInPosition = false;
- for (size_t dir = 0; dir != finals.size(); ++dir) {
- if (finals[dir] == false)continue;
- if (directionsForPeriodicBCond_[dir]) {
- jumpInPosition = true;
- }
- }
- if (jumpInPosition == true) {
- // in this case this guy is final, so we need to find 'the opposite one'
- position = compute_position_in_bitmap(c.findOpposite(directionsForPeriodicBCond_));
- } else {
- position = i;
- }
- }
- result.push_back(position);
- if (dbg) {
- std::cerr << " position : " << position << std::endl;
- std::cerr << c << std::endl;
- getchar();
- }
-
- c.increment();
- }
-
- return result;
-}
-
-#endif // BITMAP_CUBICAL_COMPLEX_BASE_H_
+/* 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 <http://www.gnu.org/licenses/>. + */ + +#pragma once + +#include <iostream> +#include <string> +#include <vector> +#include <string> +#include <fstream> +#include <algorithm> +#include <iterator> +#include <limits> +#include "counter.h" + + +using namespace std; + +namespace Gudhi +{ + +namespace Cubical_complex +{ + + + +/** + * This is a class implementing a basic bitmap data structure to store cubical complexes. + * It implements only the most basic subroutines. + * The idea of the bitmap is the following. Our aim is to have a memory efficient + * data structure to store d-dimensional cubical complex + * C being a cubical decomposition + * of a rectangular region of a space. This is achieved by storing C as a + * vector of bits (this is where the name 'bitmap' came from). + * Each cell is represented by a single + * bit (in case of black and white bitmaps, or by a single element of a type T + * (here T is a filtration type of a bitmap, typically a double). + * All the informations needed for homology and + * persistent homology computations (like dimension of a cell, boundary and + * coboundary elements of a cell, are then obtained from the + * position of the element in C. + * The default filtration used in this implementation is the lower star filtration. + */ +template <typename T> +class Bitmap_cubical_complex_base +{ +public: + /** + * There are a few constructors of a Bitmap_cubical_complex_base class. + * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal + * the number of elements in the + * input vector and size in the i-th dimension equal the number in the position i-of the input vector. + */ + Bitmap_cubical_complex_base( std::vector<unsigned>& sizes ); + /** + * The second constructor takes as a input a Perseus style file. For more details, + * please consult the documentations of + * Perseus software as well as examples attached to this + * implementation. + **/ + Bitmap_cubical_complex_base( const char* perseus_style_file ); + /** + * The last constructor of a Bitmap_cubical_complex_base class accepts vector of dimensions (as the first one) + * together with vector of filtration values of top dimensional cells. + **/ + Bitmap_cubical_complex_base( std::vector<unsigned>& dimensions , const std::vector<T>& top_dimensional_cells ); + + /** + * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell + * and get_cell_data are the basic + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. The input parameter of all of those function is a + * non-negative integer, indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to + * the positions of (co)boundary element of the input cell. + */ + inline std::vector< size_t > get_boundary_of_a_cell( size_t cell )const; + /** + * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell, + * get_dimension_of_a_cell and get_cell_data are the basic + * functions that compute boundary / coboundary / dimension and the filtration + * value form a position of a cell in the structure of a bitmap. + * The input parameter of all of those function is a non-negative integer, + * indicating a position of a cube in the data structure. + * In the case of functions that compute (co)boundary, the output is a vector if + * non-negative integers pointing to the + * positions of (co)boundary element of the input cell. + **/ + inline std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const; + /** + * In the case of get_dimension_of_a_cell function, the output is a non-negative integer + * indicating the dimension of a cell. + **/ + inline unsigned get_dimension_of_a_cell( size_t cell )const; + /** + * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position. + **/ + inline T& get_cell_data( size_t cell ); + + + /** + * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. + * Then, there are a few ways one can pick the filtration of lower dimensional + * cells. The most typical one is by so called lower star filtration. This function is always called by any + * constructor which takes the top dimensional cells. If you use such a constructor, + * then there is no need to call this function. Call it only if you are putting the filtration + * of the cells by your own (for instance by using Top_dimensional_cells_iterator). + **/ + void impose_lower_star_filtration();//assume that top dimensional cells are already set. + + /** + * Returns dimension of a complex. + **/ + inline unsigned dimension()const{ return sizes.size(); } + + /** + * Returns number of all cubes in the data structure. + **/ + inline unsigned size_of_bitmap()const + { + return this->data.size(); + } + + /** + * Writing to stream operator. + **/ + template <typename K> + friend ostream& operator << ( ostream & os , const Bitmap_cubical_complex_base<K>& b ); + + //ITERATORS + + /** + * Iterator through all cells in the complex (in order they appear in the structure -- i.e. + * in lexicographical order). + **/ + typedef typename std::vector< T >::iterator all_cells_iterator; + all_cells_iterator all_cells_begin()const + { + return this->data.begin(); + } + all_cells_iterator all_cells_end()const + { + return this->data.end(); + } + + + typedef typename std::vector< T >::const_iterator all_cells_const_iterator; + all_cells_const_iterator all_cells_const_begin()const + { + return this->data.begin(); + } + all_cells_const_iterator all_cells_const_end()const + { + return this->data.end(); + } + + /** + * Iterator through top dimensional cells of the complex. The cells appear in order they are stored + * in the structure (i.e. in lexicographical order) + **/ + class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, double > + { + public: + Top_dimensional_cells_iterator( Bitmap_cubical_complex_base& b ):b(b) + { + for ( size_t i = 0 ; i != b.dimension() ; ++i ) + { + this->counter.push_back(0); + } + } + Top_dimensional_cells_iterator operator++() + { + //first find first element of the counter that can be increased: + size_t dim = 0; + while ( ( dim != this->b.dimension() ) && ( this->counter[dim] == this->b.sizes[dim]-1 ) )++dim; + + if ( dim != this->b.dimension() ) + { + ++this->counter[dim]; + for ( size_t i = 0 ; i != dim ; ++i ) + { + this->counter[i] = 0; + } + } + else + { + ++this->counter[0]; + } + return *this; + } + Top_dimensional_cells_iterator operator++(int) + { + Top_dimensional_cells_iterator result = *this; + ++(*this); + return result; + } + Top_dimensional_cells_iterator operator =( const Top_dimensional_cells_iterator& rhs ) + { + this->counter = rhs.counter; + this->b = rhs.b; + return *this; + } + bool operator == ( const Top_dimensional_cells_iterator& rhs )const + { + if ( &this->b != &rhs.b )return false; + if ( this->counter.size() != rhs.counter.size() )return false; + for ( size_t i = 0 ; i != this->counter.size() ; ++i ) + { + if ( this->counter[i] != rhs.counter[i] )return false; + } + return true; + } + bool operator != ( const Top_dimensional_cells_iterator& rhs )const + { + return !(*this == rhs); + } + + T& operator*() + { + //given the counter, compute the index in the array and return this element. + unsigned index = 0; + for ( size_t i = 0 ; i != this->counter.size() ; ++i ) + { + index += (2*this->counter[i]+1)*this->b.multipliers[i]; + } + return this->b.data[index]; + } + + size_t compute_index_in_bitmap()const + { + size_t index = 0; + for ( size_t i = 0 ; i != this->counter.size() ; ++i ) + { + index += (2*this->counter[i]+1)*this->b.multipliers[i]; + } + return index; + } + + void print_counter()const + { + for ( size_t i = 0 ; i != this->counter.size() ; ++i ) + { + cout << this->counter[i] << " "; + } + } + friend class Bitmap_cubical_complex_base; + protected: + std::vector< unsigned > counter; + Bitmap_cubical_complex_base& b; + }; + Top_dimensional_cells_iterator top_dimensional_cells_begin() + { + Top_dimensional_cells_iterator a(*this); + return a; + } + Top_dimensional_cells_iterator top_dimensional_cells_end() + { + Top_dimensional_cells_iterator a(*this); + for ( size_t i = 0 ; i != this->dimension() ; ++i ) + { + a.counter[i] = this->sizes[i]-1; + } + a.counter[0]++; + return a; + } + + +//****************************************************************************************************************// +//****************************************************************************************************************// +//****************************************************************************************************************// +//****************************************************************************************************************// + + +//****************************************************************************************************************// +//****************************************************************************************************************// +//****************************************************************************************************************// +//****************************************************************************************************************// + +protected: + std::vector<unsigned> sizes; + std::vector<unsigned> multipliers; + std::vector<T> data; + size_t total_number_of_cells; + void set_up_containers( std::vector<unsigned>& sizes ) + { + unsigned multiplier = 1; + for ( size_t i = 0 ; i != sizes.size() ; ++i ) + { + this->sizes.push_back(sizes[i]); + this->multipliers.push_back(multiplier); + //multiplier *= 2*(sizes[i]+1)+1; + multiplier *= 2*sizes[i]+1; + } + //std::reverse( this->sizes.begin() , this->sizes.end() ); + std::vector<T> data(multiplier); + std::fill( data.begin() , data.end() , std::numeric_limits<int>::max() ); + this->total_number_of_cells = multiplier; + this->data = data; + } + + size_t compute_position_in_bitmap( std::vector< unsigned >& counter ) + { + size_t position = 0; + for ( size_t i = 0 ; i != this->multipliers.size() ; ++i ) + { + position += this->multipliers[i]*counter[i]; + } + return position; + } + + std::vector<unsigned> compute_counter_for_given_cell( size_t cell )const + { + std::vector<unsigned> counter; + for ( size_t dim = this->sizes.size() ; dim != 0 ; --dim ) + { + counter.push_back(cell/this->multipliers[dim-1]); + cell = cell%this->multipliers[dim-1]; + } + std::reverse( counter.begin() , counter.end() ); + return counter; + } + + std::vector< size_t > + generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions + ( std::vector< bool >& directions_for_periodic_b_cond ); +}; + + + + +template <typename K> +ostream& operator << ( ostream & out , const Bitmap_cubical_complex_base<K>& b ) +{ + for ( typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator + it = b.all_cells_const_begin() ; it != b.all_cells_const_end() ; ++it ) + { + out << *it << " "; + } + return out; +} + + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +( std::vector<unsigned>& sizes ) +{ + this->set_up_containers( sizes ); +} + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base +( std::vector<unsigned>& sizes_in_following_directions , const std::vector<T>& top_dimensional_cells ) +{ + this->set_up_containers( sizes_in_following_directions ); + + size_t number_of_top_dimensional_elements = 1; + for ( size_t i = 0 ; i != sizes_in_following_directions.size() ; ++i ) + { + number_of_top_dimensional_elements *= sizes_in_following_directions[i]; + } + if ( number_of_top_dimensional_elements != top_dimensional_cells.size() ) + { + cerr << + "Error in constructor\ + Bitmap_cubical_complex_base\ + ( std::vector<size_t> sizes_in_following_directions , std::vector<float> top_dimensional_cells ).\ + Number of top dimensional elements that follow from sizes_in_following_directions vector is different\ + than the size of top_dimensional_cells vector." << endl; + throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions,\ + std::vector<float> top_dimensional_cells )\ + . Number of top dimensional elements that follow from sizes_in_following_directions vector is different than the\ + size of top_dimensional_cells vector."); + } + + Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + size_t index = 0; + for ( it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) + { + (*it) = top_dimensional_cells[index]; + ++index; + } + this->impose_lower_star_filtration(); +} + + +template <typename T> +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus_style_file ) +{ + bool dbg = false; + ifstream inFiltration, inIds; + inFiltration.open( perseus_style_file ); + unsigned dimensionOfData; + inFiltration >> dimensionOfData; + + if (dbg){cerr << "dimensionOfData : " << dimensionOfData << endl;} + + std::vector<unsigned> sizes; + for ( size_t i = 0 ; i != dimensionOfData ; ++i ) + { + unsigned size_in_this_dimension; + inFiltration >> size_in_this_dimension; + size_in_this_dimension = abs(size_in_this_dimension); + sizes.push_back( size_in_this_dimension ); + if (dbg){cerr << "size_in_this_dimension : " << size_in_this_dimension << endl;} + } + this->set_up_containers( sizes ); + + Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + it = this->top_dimensional_cells_begin(); + + //TODO -- over here we also need to read id's of cell and put them to bitmapElement structure! + while ( !inFiltration.eof() ) + { + double filtrationLevel; + inFiltration >> filtrationLevel; + 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(); +} + + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell( size_t cell )const +{ + bool bdg = false; + //first of all, we need to take the list of coordinates in which the cell has nonzero length. + //We do it by using modified version to compute dimension of a cell: + std::vector< unsigned > dimensions_in_which_cell_has_nonzero_length; + unsigned dimension = 0; + size_t cell1 = cell; + for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) + { + unsigned position = cell1/multipliers[i-1]; + if ( position%2 == 1 ) + { + dimensions_in_which_cell_has_nonzero_length.push_back(i-1); + dimension++; + } + cell1 = cell1%multipliers[i-1]; + } + + if (bdg) + { + cerr << "dimensions_in_which_cell_has_nonzero_length : \n"; + for ( size_t i = 0 ; i != dimensions_in_which_cell_has_nonzero_length.size() ; ++i ) + { + cerr << dimensions_in_which_cell_has_nonzero_length[i] << endl; + } + getchar(); + } + + std::vector< size_t > boundary_elements; + if ( dimensions_in_which_cell_has_nonzero_length.size() == 0 )return boundary_elements; + for ( size_t i = 0 ; i != dimensions_in_which_cell_has_nonzero_length.size() ; ++i ) + { + boundary_elements.push_back( cell - multipliers[ dimensions_in_which_cell_has_nonzero_length[i] ] ); + boundary_elements.push_back( cell + multipliers[ dimensions_in_which_cell_has_nonzero_length[i] ] ); + + if (bdg) cerr << "multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " + << multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl; + if (bdg) cerr << "cell - multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " + << cell - multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl; + if (bdg) cerr << "cell + multipliers[dimensions_in_which_cell_has_nonzero_length[i]] : " + << cell + multipliers[dimensions_in_which_cell_has_nonzero_length[i]] << endl; + } + return boundary_elements; +} + + + + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell( size_t cell )const +{ + bool bdg = false; + //first of all, we need to take the list of coordinates in which the cell has nonzero length. + //We do it by using modified version to compute dimension of a cell: + std::vector< unsigned > dimensions_in_which_cell_has_zero_length; + unsigned dimension = 0; + size_t cell1 = cell; + for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) + { + unsigned position = cell1/multipliers[i-1]; + if ( position%2 == 0 ) + { + dimensions_in_which_cell_has_zero_length.push_back(i-1); + dimension++; + } + cell1 = cell1%multipliers[i-1]; + } + + std::vector<unsigned> counter = this->compute_counter_for_given_cell( cell ); + //reverse(counter.begin() , counter.end()); + + if (bdg) + { + cerr << "dimensions_in_which_cell_has_zero_length : \n"; + for ( size_t i = 0 ; i != dimensions_in_which_cell_has_zero_length.size() ; ++i ) + { + cerr << dimensions_in_which_cell_has_zero_length[i] << endl; + } + cerr << "\n counter : " << endl; + for ( size_t i = 0 ; i != counter.size() ; ++i ) + { + cerr << counter[i] << endl; + } + getchar(); + } + + std::vector< size_t > coboundary_elements; + if ( dimensions_in_which_cell_has_zero_length.size() == 0 )return coboundary_elements; + for ( size_t i = 0 ; i != dimensions_in_which_cell_has_zero_length.size() ; ++i ) + { + if ( bdg ) + { + cerr << "Dimension : " << i << endl; + if (counter[dimensions_in_which_cell_has_zero_length[i]] == 0) + { + cerr << "In dimension : " << i + << " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n"; + } + if ( counter[dimensions_in_which_cell_has_zero_length[i]] + == + 2*this->sizes[dimensions_in_which_cell_has_zero_length[i]] ) + { + cerr << "In dimension : " << i + << " we cannot substract, since we will jump out of a Bitmap_cubical_complex_base \n"; + } + } + + + if ( (cell > multipliers[dimensions_in_which_cell_has_zero_length[i]]) + && (counter[dimensions_in_which_cell_has_zero_length[i]] != 0) ) + //if ( counter[dimensions_in_which_cell_has_zero_length[i]] != 0 ) + { + if ( bdg ) + { + cerr << "Subtracting : " << cell - multipliers[dimensions_in_which_cell_has_zero_length[i]] << endl; + } + coboundary_elements.push_back( cell - multipliers[dimensions_in_which_cell_has_zero_length[i]] ); + } + if ( + (cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] < this->data.size()) && + (counter[dimensions_in_which_cell_has_zero_length[i]] + != + 2*this->sizes[dimensions_in_which_cell_has_zero_length[i]]) + ) + //if ( counter[dimensions_in_which_cell_has_zero_length[i]] != + //2*this->sizes[dimensions_in_which_cell_has_zero_length[i]] ) + { + coboundary_elements.push_back( cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] ); + if ( bdg )cerr << "Adding : " << cell + multipliers[dimensions_in_which_cell_has_zero_length[i]] << endl; + } + } + return coboundary_elements; +} + + + + + + +template <typename T> +unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell( size_t cell )const +{ + bool dbg = false; + if (dbg)cerr << "\n\n\n Computing position o a cell of an index : " << cell << endl; + unsigned dimension = 0; + for ( size_t i = this->multipliers.size() ; i != 0 ; --i ) + { + unsigned position = cell/multipliers[i-1]; + + if (dbg)cerr << "i-1 :" << i-1 << endl; + if (dbg)cerr << "cell : " << cell << endl; + if (dbg)cerr << "position : " << position << endl; + if (dbg)cerr << "multipliers["<<i-1<<"] = " << multipliers[i-1] << endl; + if (dbg)getchar(); + + if ( position%2 == 1 ) + { + if (dbg)cerr << "Nonzero length in this direction \n"; + dimension++; + } + cell = cell%multipliers[i-1]; + } + return dimension; +} + +template <typename T> +T& Bitmap_cubical_complex_base<T>::get_cell_data( size_t cell ) +{ + return this->data[cell]; +} + + +template <typename T> +void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() +{ + bool dbg = false; + + //this vector will be used to check which elements have already been taken care of + //in imposing lower star filtration: + std::vector<bool> is_this_cell_considered( this->data.size() , false ); + + std::vector<size_t> indices_to_consider; + //we assume here that we already have a filtration on the top dimensional cells and + //we have to extend it to lower ones. + typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); + for ( it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it ) + { + indices_to_consider.push_back( it.compute_index_in_bitmap() ); + } + + while ( indices_to_consider.size() ) + { + if ( dbg ) + { + cerr << "indices_to_consider in this iteration \n"; + for ( size_t i = 0 ; i != indices_to_consider.size() ; ++i ) + { + cout << indices_to_consider[i] << " "; + } + getchar(); + } + std::vector<size_t> new_indices_to_consider; + for ( size_t i = 0 ; i != indices_to_consider.size() ; ++i ) + { + std::vector<size_t> bd = this->get_boundary_of_a_cell( indices_to_consider[i] ); + for ( size_t boundaryIt = 0 ; boundaryIt != bd.size() ; ++boundaryIt ) + { + if ( this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ] ) + { + this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ]; + } + if ( is_this_cell_considered[ bd[boundaryIt] ] == false ) + { + new_indices_to_consider.push_back( bd[boundaryIt] ); + is_this_cell_considered[ bd[boundaryIt] ] = true; + } + } + } + indices_to_consider.swap(new_indices_to_consider); + } +} + + +template <typename T> +bool compareFirstElementsOfTuples( const std::pair< std::pair< T , size_t > , char >& first , + const std::pair< std::pair< T , size_t > , char >& second ) +{ + if ( first.first.first < second.first.first ) + { + return true; + } + else + { + if ( first.first.first > second.first.first ) + { + return false; + } + //in this case first.first.first == second.first.first, so we need to compare dimensions + return first.second < second.second; + } +} + + + +template <typename T> +std::vector< size_t > Bitmap_cubical_complex_base<T>:: +generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions +( std::vector< bool >& directions_for_periodic_b_cond ) +{ + bool dbg = false; + if ( this->sizes.size() != directions_for_periodic_b_cond.size() ) + throw "directions_for_periodic_b_cond vector size is different from the size of the bitmap. Program terminate \n"; + + std::vector<unsigned> sizes( this->sizes.size() ); + for ( size_t i = 0 ; i != this->sizes.size() ; ++i )sizes[i] = 2*this->sizes[i]; + + counter c( sizes ); + + std::vector< size_t > result; + + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + size_t position; + if ( !c.isFinal() ) + { + position = i; + //result.push_back( i ); + } + else + { + std::vector< bool > finals = c.directions_of_finals(); + bool jump_in_position = false; + for ( size_t dir = 0 ; dir != finals.size() ; ++dir ) + { + if ( finals[dir] == false )continue; + if ( directions_for_periodic_b_cond[dir] ) + { + jump_in_position = true; + } + } + if ( jump_in_position == true ) + { + //in this case this guy is final, so we need to find 'the opposite one' + position = compute_position_in_bitmap( c.find_opposite( directions_for_periodic_b_cond ) ); + } + else + { + position = i; + } + } + result.push_back( position ); + if ( dbg ) + { + cerr << " position : " << position << endl; + cerr << c << endl; + getchar(); + } + + c.increment(); + } + + return result; +} + +} + +}
\ No newline at end of file |