Merge of bitmap feature from bitmap branch

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@1080 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 95d40c06d5ca0bffc9166edb04681bfd16ed116c

13 files changed, 3668 insertions, 0 deletions

diff --git a/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h
new file mode 100644
index 00000000..cde0b2fc
--- /dev/null
+++ b/src/Bitmap_cubical_complex/doc/Gudhi_Cubical_Complex_doc.h
@@ -0,0 +1,156 @@
+/*    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 DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_

+#define DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_

+

+namespace Gudhi {

+

+namespace Cubical_complex {

+

+/**  \defgroup cubical_complex Cubical complex

+ *

+ * \author   Pawel Dlotko

+ *

+ * @{

+ *

+

+ * Bitmap_cubical_complex is an example of a structured complex useful in computational mathematics (specially rigorous

+ * numerics) and image analysis. The presented implementation of cubical complexes is based on the following

+ * definition.

+ *

+ * An <em>elementary interval</em> is an interval of a form \f$ [n,n+1] \f$, or \f$[n,n]\f$, for \f$ n \in \mathcal{Z}

+ * \f$. The first one is called <em>non-degenerate</em>, while the second one is \a degenerate interval. A

+ * <em>boundary of a elementary interval</em> is a chain  \f$\partial [n,n+1] = [n+1,n+1]-[n,n] \f$ in case of

+ * non-degenerated elementary interval and \f$\partial [n,n] = 0 \f$ in case of degenerate elementary interval. An

+ * <em>elementary cube</em> \f$ C \f$ is a product of elementary intervals, \f$C=I_1 \times \ldots \times I_n\f$.

+ * <em>Embedding dimension</em> of a cube is n, the number of elementary intervals (degenerate or not) in the product.

+ * A <em>dimension of a cube</em> \f$C=I_1 \times ... \times I_n\f$ is the number of non degenerate elementary

+ * intervals in the product. A <em>boundary of a cube</em> \f$C=I_1 \times \ldots \times I_n\f$ is a chain obtained

+ * in the following way:

+ * \f[\partial C = (\partial I_1 \times \ldots \times I_n) + (I_1 \times \partial I_2 \times \ldots \times I_n) +

+ * \ldots + (I_1 \times I_2 \times \ldots \times \partial I_n).\f]

+ * A <em>cubical complex</em> \f$\mathcal{K}\f$ is a collection of cubes closed under operation of taking boundary

+ * (i.e. boundary of every cube from the collection is in the collection). A cube \f$C\f$ in cubical complex

+ * \f$\mathcal{K}\f$ is <em>maximal</em> if it is not in a boundary of any other cube in \f$\mathcal{K}\f$. A \a

+ * support of a cube \f$C\f$ is the set in \f$\mathbb{R}^n\f$ occupied by \f$C\f$ (\f$n\f$ is the embedding dimension

+ * of \f$C\f$).

+ *

+ * Cubes may be equipped with a filtration values in which case we have filtered cubical complex. All the cubical

+ * complexes considered in this implementation are filtered cubical complexes (although, the range of a filtration may

+ * be a set of two elements).

+ *

+ * For further details and theory of cubical complexes, please consult \cite kaczynski2004computational as well as the

+ * following paper \cite peikert2012topological .

+ *

+ * \section datastructure Data structure.

+ *

+ * The implementation of Cubical complex provides a representation of complexes that occupy a rectangular region in

+ * \f$\mathbb{R}^n\f$. This extra assumption allows for a memory efficient way of storing cubical complexes in a form

+ * of so called bitmaps. Let \f$R = [b_1,e_1] \times \ldots \times [b_n,e_n]\f$, for \f$b_1,...b_n,e_1,...,e_n \in

+ * \mathbb{Z}\f$, \f$b_i \leq d_i\f$ be the considered rectangular region and let \f$\mathcal{K}\f$ be a filtered

+ * cubical complex having the rectangle \f$R\f$ as its support. Note that the structure of the coordinate system gives

+ * a way a lexicographical ordering of cells of \f$\mathcal{K}\f$. This ordering is a base of the presented

+ * bitmap-based implementation. In this implementation, the whole cubical complex is stored as a vector of the values

+ * of filtration. This, together with dimension of \f$\mathcal{K}\f$ and the sizes of \f$\mathcal{K}\f$ in all

+ * directions, allows to determine, dimension, neighborhood, boundary and coboundary of every cube \f$C \in

+ * \mathcal{K}\f$.

+ *

+ * \image html "bitmapAllCubes.png" "Cubical complex.

+ *

+ * Note that the cubical complex in the figure above is, in a natural way, a product of one dimensional cubical

+ * complexes in \f$\mathbb{R}\f$. The number of all cubes in each direction is equal \f$2n+1\f$, where \f$n\f$ is the

+ * number of maximal cubes in the considered direction. Let us consider a cube at the position \f$k\f$ in the bitmap.

+ * Knowing the sizes of the bitmap, by a series of modulo operation, we can determine which elementary intervals are

+ * present in the product that gives the cube \f$C\f$. In a similar way, we can compute boundary and the coboundary of

+ * each cube. Further details can be found in the literature.

+ *

+ * \section inputformat Input Format.

+ *

+ * In the current implantation, filtration is given at the maximal cubes, and it is then extended by the lower star

+ * filtration to all cubes. There are a number of constructors that can be used to construct cubical complex by users

+ * who want to use the code directly. They can be found in the \a Bitmap_cubical_complex class.

+ * Currently one input from a text file is used. It uses a format used already in Perseus software

+ * (http://www.sas.upenn.edu/~vnanda/perseus/) by Vidit Nanda.

+ * Below we are providing a description of the format. The first line contains a number d begin the dimension of the

+ * bitmap (2 in the example below). Next d lines are the numbers of top dimensional cubes in each dimensions (3 and 3

+ * in the example below). Next, in lexicographical order, the filtration of top dimensional cubes is given (1 4 6 8

+ * 20 4 7 6 5 in the example below).

+ *

+ *

+ * \image html "exampleBitmap.png" "Example of a input data."

+ *

+ * The input file for the following complex is:

+ * \verbatim

+2

+3

+3

+1

+4

+6

+8

+20

+4

+7

+6

+5

+\endverbatim

+

+ * \section PeriodicBoundaryConditions Periodic boundary conditions

+ * Often one would like to impose periodic boundary conditions to the cubical complex. Let \f$ I_1\times ... \times

+ * I_n \f$ be a box that is decomposed with a cubical complex \f$ \mathcal{K} \f$. Imposing periodic boundary

+ * conditions in the direction i, means that the left and the right side of a complex \f$ \mathcal{K} \f$ are

+ * considered the same. In particular, if for a bitmap \f$ \mathcal{K} \f$ periodic boundary conditions are imposed

+ * in all directions, then complex \f$ \mathcal{K} \f$ became n-dimensional torus. One can use various constructors

+ * from the file Bitmap_cubical_complex_periodic_boundary_conditions_base.h to construct cubical complex with periodic

+ * boundary conditions. One can also use Perseus style input files. To indicate periodic boundary conditions in a

+ * given direction, then number of top dimensional cells in this direction have to be multiplied by -1. For instance:

+

+ *\verbatim

+2

+-3

+3

+1

+4

+6

+8

+20

+4

+7

+6

+5

+\endverbatim

+

+ * Indicate that we have imposed periodic boundary conditions in the direction x, but not in the direction y.

+

+ * \section BitmapExamples Examples

+ * End user programs are available in example/Bitmap_cubical_complex folder.

+

+ */

+/** @} */  // end defgroup cubical_complex

+

+}  // namespace Cubical_complex

+

+}  // namespace Gudhi

+

+#endif  // DOC_GUDHI_CUBICAL_COMPLEX_COMPLEX_H_

diff --git a/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png b/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png
new file mode 100644
index 00000000..77167b13
--- /dev/null
+++ b/src/Bitmap_cubical_complex/doc/bitmapAllCubes.png
diff --git a/src/Bitmap_cubical_complex/doc/exampleBitmap.png b/src/Bitmap_cubical_complex/doc/exampleBitmap.png
new file mode 100644
index 00000000..069c6eb2
--- /dev/null
+++ b/src/Bitmap_cubical_complex/doc/exampleBitmap.png
diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
new file mode 100644
index 00000000..7c7e8ac5
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
@@ -0,0 +1,72 @@
+/*    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 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/>.
+ */
+
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " <<
+      "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a" <<
+      "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let " <<
+      "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are " <<
+      "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See " <<
+      "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n" << std::endl;
+
+  int p = 2;
+  double min_persistence = 0;
+
+  if (argc != 2) {
+    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " <<
+        "the input. The program will now terminate.\n";
+    return 1;
+  }
+
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base;
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex;
+  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
+  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
+
+  Bitmap_cubical_complex b(argv[1]);
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(b);
+  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  std::stringstream ss;
+  ss << argv[1] << "_persistence";
+  std::ofstream out(ss.str().c_str());
+  pcoh.output_diagram(out);
+  out.close();
+
+  std::cout << "Result in file: " << ss.str().c_str() << "\n";
+
+  return 0;
+}
+
diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
new file mode 100644
index 00000000..2c5d7fd3
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
@@ -0,0 +1,74 @@
+/*    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 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/>.
+ */
+
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes persistent homology, by using " <<
+      "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in " <<
+      "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 " <<
+      "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers " <<
+      "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We " <<
+      "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for" <<
+      " example.\n" << std::endl;
+
+  int p = 2;
+  double min_persistence = 0;
+
+  if (argc != 2) {
+    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " <<
+        "the input. The program will now terminate.\n";
+    return 1;
+  }
+
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> Bitmap_base;
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex< Bitmap_base > Bitmap_cubical_complex;
+
+  Bitmap_cubical_complex b(argv[1]);
+
+  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
+  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(b, true);
+  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  std::stringstream ss;
+  ss << argv[1] << "_persistence";
+  std::ofstream out(ss.str().c_str());
+  pcoh.output_diagram(out);
+  out.close();
+
+  std::cout << "Result in file: " << ss.str().c_str() << "\n";
+
+  return 0;
+}
+
diff --git a/src/Bitmap_cubical_complex/example/CMakeLists.txt b/src/Bitmap_cubical_complex/example/CMakeLists.txt
new file mode 100644
index 00000000..8f9cfa80
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/CMakeLists.txt
@@ -0,0 +1,17 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBitmap)
+
+add_executable ( Bitmap_cubical_complex Bitmap_cubical_complex.cpp )
+target_link_libraries(Bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY})	
+add_test(Bitmap_cubical_complex_one_sphere ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt)
+add_test(Bitmap_cubical_complex_two_sphere ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt)
+
+add_executable ( Random_bitmap_cubical_complex Random_bitmap_cubical_complex.cpp )
+target_link_libraries(Random_bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY})	
+add_test(Random_bitmap_cubical_complex ${CMAKE_CURRENT_BINARY_DIR}/Random_bitmap_cubical_complex 2 100 100)
+
+add_executable ( Bitmap_cubical_complex_periodic_boundary_conditions Bitmap_cubical_complex_periodic_boundary_conditions.cpp )
+target_link_libraries(Bitmap_cubical_complex_periodic_boundary_conditions ${Boost_SYSTEM_LIBRARY})	
+add_test(Bitmap_cubical_complex_periodic_2d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt)
+add_test(Bitmap_cubical_complex_periodic_3d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt)
+
diff --git a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
new file mode 100644
index 00000000..416ad8f2
--- /dev/null
+++ b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp
@@ -0,0 +1,83 @@
+/*    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 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/>.
+ */
+
+
+// for persistence algorithm
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+int main(int argc, char** argv) {
+  srand(time(0));
+
+  std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " <<
+      "complexes. The first parameter of the program is the dimension D of the bitmap. The next D parameters are " <<
+      "number of top dimensional cubes in each dimension of the bitmap. The program will create random cubical " <<
+      "complex of that sizes and compute persistent homology of it." << std::endl;
+
+  int p = 2;
+  double min_persistence = 0;
+
+  if (argc < 3) {
+    std::cerr << "Wrong number of parameters, the program will now terminate\n";
+    return 1;
+  }
+
+  size_t dimensionOfBitmap = (size_t) atoi(argv[1]);
+  std::vector< unsigned > sizes;
+  size_t multipliers = 1;
+  for (size_t dim = 0; dim != dimensionOfBitmap; ++dim) {
+    unsigned sizeInThisDimension = (unsigned) atoi(argv[2 + dim]);
+    sizes.push_back(sizeInThisDimension);
+    multipliers *= sizeInThisDimension;
+  }
+
+  std::vector< double > data;
+  for (size_t i = 0; i != multipliers; ++i) {
+    data.push_back(rand() / static_cast<double>(RAND_MAX));
+  }
+
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base;
+  typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex;
+  Bitmap_cubical_complex b(sizes, data);
+
+  // Compute the persistence diagram of the complex
+  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
+  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
+  Persistent_cohomology pcoh(b);
+  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  std::stringstream ss;
+  ss << "randomComplex_persistence";
+  std::ofstream out(ss.str().c_str());
+  pcoh.output_diagram(out);
+  out.close();
+
+  return 0;
+}
+
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
new file mode 100644
index 00000000..7bf9617e
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -0,0 +1,592 @@
+/*    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_H_
+#define BITMAP_CUBICAL_COMPLEX_H_
+
+#include <gudhi/Bitmap_cubical_complex_base.h>
+#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>
+
+#ifdef GUDHI_USE_TBB
+#include <tbb/parallel_sort.h>
+#endif
+
+#include <limits>
+#include <utility>  // for pair<>
+#include <algorithm>  // for sort
+#include <vector>
+#include <numeric>  // for iota
+
+namespace Gudhi {
+
+namespace Cubical_complex {
+
+// global variable, was used just for debugging.
+const bool globalDbg = false;
+
+template <typename T> class is_before_in_filtration;
+
+/**
+* This is a Bitmap_cubical_complex class. It joints a functionalities of Bitmap_cubical_complex_base and Bitmap_cubical_complex_periodic_boundary_conditions_base classes into
+* Gudhi persistent homology engine. It is a template class that inherit from its template parameter. The template parameter is supposed to be either Bitmap_cubical_complex_base or Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+**/
+
+/**
+ *@class Bitmap_cubical_complex
+ *@brief Cubical complex represented as a bitmap.
+ *@ingroup cubical_complex
+ */
+template <typename T>
+class Bitmap_cubical_complex : public T {
+ public:
+  //*********************************************//
+  // Typedefs and typenames
+  //*********************************************//
+  typedef size_t Simplex_key;
+  typedef typename T::filtration_type Filtration_value;
+  typedef Simplex_key Simplex_handle;
+
+
+  //*********************************************//
+  // Constructors
+  //*********************************************//
+  // Over here we need to define various input types. I am proposing the following ones:
+  // Perseus style
+  // TODO(PD) H5 files?
+  // TODO(PD) binary files with little endiangs / big endians ?
+  // TODO(PD) constructor from a vector of elements of a type T. ?
+
+  /**
+   * Constructor form a Perseus-style file.
+   **/
+  Bitmap_cubical_complex(const char* perseus_style_file) :
+      T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) {
+    if (globalDbg) {
+      std::cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n";
+    }
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // 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->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * 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(const std::vector<unsigned>& dimensions,
+                         const std::vector<typename T::filtration_type>& top_dimensional_cells) :
+      T(dimensions, top_dimensional_cells),
+      key_associated_to_simplex(this->total_number_of_cells + 1) {
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // 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->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * 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::filtration_type
+   * with filtration on top dimensional cells. The last parameter of the constructor is a vector of bools of a length equal to the dimension of cubical complex.
+   * If the position i on this vector is true, then we impose periodic boundary conditions in this direction.
+   **/
+  Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
+                         const std::vector<typename T::filtration_type>& top_dimensional_cells,
+                         std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) :
+      T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed),
+      key_associated_to_simplex(this->total_number_of_cells + 1) {
+    for (size_t i = 0; i != this->total_number_of_cells; ++i) {
+      this->key_associated_to_simplex[i] = i;
+    }
+    // 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->initialize_simplex_associated_to_key();
+  }
+
+  /**
+   * Destructor of the Bitmap_cubical_complex class.
+   **/
+  virtual ~Bitmap_cubical_complex() {}
+
+  //*********************************************//
+  // Other 'easy' functions
+  //*********************************************//
+
+  /**
+   * Returns number of all cubes in the complex.
+   **/
+  size_t num_simplices()const {
+    return this->total_number_of_cells;
+  }
+
+  /**
+   * Returns a Simplex_handle to a cube that do not exist in this complex.
+   **/
+  static Simplex_handle null_simplex() {
+    if (globalDbg) {
+      std::cerr << "Simplex_handle null_simplex()\n";
+    }
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Returns dimension of the complex.
+   **/
+  inline size_t dimension()const {
+    return this->sizes.size();
+  }
+
+  /**
+   * Return dimension of a cell pointed by the Simplex_handle.
+   **/
+  inline unsigned dimension(Simplex_handle sh)const {
+    if (globalDbg) {
+      std::cerr << "unsigned dimension(const Simplex_handle& sh)\n";
+    }
+    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh);
+    return -1;
+  }
+
+  /**
+   * Return the filtration of a cell pointed by the Simplex_handle.
+   **/
+  typename T::filtration_type filtration(Simplex_handle sh) {
+    if (globalDbg) {
+      std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n";
+    }
+    // Returns the filtration value of a simplex.
+    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh];
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Return a key which is not a key of any cube in the considered data structure.
+   **/
+  static Simplex_key null_key() {
+    if (globalDbg) {
+      std::cerr << "Simplex_key null_key()\n";
+    }
+    return std::numeric_limits<Simplex_handle>::max();
+  }
+
+  /**
+   * Return the key of a cube pointed by the Simplex_handle.
+   **/
+  Simplex_key key(Simplex_handle sh)const {
+    if (globalDbg) {
+      std::cerr << "Simplex_key key(const Simplex_handle& sh)\n";
+    }
+    if (sh != std::numeric_limits<Simplex_handle>::max()) {
+      return this->key_associated_to_simplex[sh];
+    }
+    return this->null_key();
+  }
+
+  /**
+   * Return the Simplex_handle given the key of the cube.
+   **/
+  Simplex_handle simplex(Simplex_key key) {
+    if (globalDbg) {
+      std::cerr << "Simplex_handle simplex(Simplex_key key)\n";
+    }
+    if (key != std::numeric_limits<Simplex_handle>::max()) {
+      return this->simplex_associated_to_key[ key ];
+    }
+    return null_simplex();
+  }
+
+  /**
+   * Assign key to a cube pointed by the Simplex_handle
+   **/
+  void assign_key(Simplex_handle sh, Simplex_key key) {
+    if (globalDbg) {
+      std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
+    }
+    if (key == std::numeric_limits<Simplex_handle>::max()) return;
+    this->key_associated_to_simplex[sh] = key;
+    this->simplex_associated_to_key[key] = sh;
+  }
+
+  /**
+   * Function called from a constructor. It is needed for Filtration_simplex_iterator to work.
+   **/
+  void initialize_simplex_associated_to_key();
+
+  //*********************************************//
+  // Iterators
+  //*********************************************//
+
+  /**
+   * Boundary_simplex_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator;
+  typedef typename std::vector< Simplex_handle > Boundary_simplex_range;
+
+  /**
+   * 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), position(0) { }
+
+    Filtration_simplex_iterator operator++() {
+      if (globalDbg) {
+        std::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) {
+        std::cerr << "Filtration_simplex_iterator operator =\n";
+      }
+      this->b = rhs.b;
+      this->position = rhs.position;
+      return (*this);
+    }
+
+    bool operator==(const Filtration_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";
+      }
+      return ( this->position == rhs.position);
+    }
+
+    bool operator!=(const Filtration_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n";
+      }
+      return !(*this == rhs);
+    }
+
+    Simplex_handle operator*() {
+      if (globalDbg) {
+        std::cerr << "Simplex_handle operator*()\n";
+      }
+      return this->b->simplex_associated_to_key[ this->position ];
+    }
+
+    friend class Filtration_simplex_range;
+
+   private:
+    Bitmap_cubical_complex<T>* 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:
+    typedef Filtration_simplex_iterator const_iterator;
+    typedef Filtration_simplex_iterator iterator;
+
+    Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { }
+
+    Filtration_simplex_iterator begin() {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator begin() \n";
+      }
+      return Filtration_simplex_iterator(this->b);
+    }
+
+    Filtration_simplex_iterator end() {
+      if (globalDbg) {
+        std::cerr << "Filtration_simplex_iterator end()\n";
+      }
+      Filtration_simplex_iterator it(this->b);
+      it.position = this->b->simplex_associated_to_key.size();
+      return it;
+    }
+
+   private:
+    Bitmap_cubical_complex<T>* 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) {
+    return this->get_boundary_of_a_cell(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) {
+      std::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(PD) 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<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
+    std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh);
+    if (globalDbg) {
+      std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n";
+      std::cerr << "bdry.size() : " << bdry.size() << std::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 have less than two elements in the "
+            "boundary.");
+    return std::make_pair(bdry[0], 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) {
+        std::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), position(0), dimension(0) { }
+
+    Skeleton_simplex_iterator operator++() {
+      if (globalDbg) {
+        std::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) {
+        std::cerr << "Skeleton_simplex_iterator operator =\n";
+      }
+      this->b = rhs.b;
+      this->position = rhs.position;
+      this->dimension = rhs.dimension;
+      return (*this);
+    }
+
+    bool operator==(const Skeleton_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator ==\n";
+      }
+      return ( this->position == rhs.position);
+    }
+
+    bool operator!=(const Skeleton_simplex_iterator& rhs)const {
+      if (globalDbg) {
+        std::cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n";
+      }
+      return !(*this == rhs);
+    }
+
+    Simplex_handle operator*() {
+      if (globalDbg) {
+        std::cerr << "Simplex_handle operator*() \n";
+      }
+      return this->position;
+    }
+
+    friend class Skeleton_simplex_range;
+   private:
+    Bitmap_cubical_complex<T>* b;
+    size_t position;
+    unsigned 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:
+    typedef Skeleton_simplex_iterator const_iterator;
+    typedef Skeleton_simplex_iterator iterator;
+
+    Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) { }
+
+    Skeleton_simplex_iterator begin() {
+      if (globalDbg) {
+        std::cerr << "Skeleton_simplex_iterator begin()\n";
+      }
+      return Skeleton_simplex_iterator(this->b, this->dimension);
+    }
+
+    Skeleton_simplex_iterator end() {
+      if (globalDbg) {
+        std::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<T>* b;
+    unsigned dimension;
+  };
+
+  /**
+   * Function needed for compatibility with Gudhi. Not useful for other purposes.
+   **/
+  Skeleton_simplex_range skeleton_simplex_range(unsigned dimension) {
+    if (globalDbg) {
+      std::cerr << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n";
+    }
+    return Skeleton_simplex_range(this, dimension);
+  }
+
+  friend class is_before_in_filtration<T>;
+
+ protected:
+  std::vector< size_t > key_associated_to_simplex;
+  std::vector< size_t > simplex_associated_to_key;
+};  // Bitmap_cubical_complex
+
+template <typename T>
+void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() {
+  if (globalDbg) {
+    std::cerr << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n";
+  }
+  this->simplex_associated_to_key = std::vector<size_t>(this->data.size());
+  std::iota(std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0);
+#ifdef GUDHI_USE_TBB
+  tbb::parallel_sort(filtration_vect_, is_before_in_filtration(this));
+#else
+  std::sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(), is_before_in_filtration<T>(this));
+#endif
+
+  // we still need to deal here with a key_associated_to_simplex:
+  for ( size_t i = 0  ; i != simplex_associated_to_key.size() ; ++i ) {
+    this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i;
+  }
+}
+
+template <typename T>
+class is_before_in_filtration {
+ public:
+  explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC)
+      : CC_(CC) { }
+
+  bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1,
+                  const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const {
+    // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
+    typename T::filtration_type fil1 = CC_->data[sh1];
+    typename T::filtration_type fil2 = CC_->data[sh2];
+    if (fil1 != fil2) {
+      return fil1 < fil2;
+    }
+    // in this case they are on the same filtration level, so the dimension decide.
+    size_t dim1 = CC_->get_dimension_of_a_cell(sh1);
+    size_t dim2 = CC_->get_dimension_of_a_cell(sh2);
+    if (dim1 != dim2) {
+      return dim1 < dim2;
+    }
+    // in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply
+    // compare their positions in the bitmap:
+    return sh1 < sh2;
+  }
+
+ protected:
+  Bitmap_cubical_complex<T>* CC_;
+};
+
+}  // namespace Cubical_complex
+
+}  // namespace Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h
new file mode 100644
index 00000000..266ce051
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex/counter.h
@@ -0,0 +1,143 @@
+/*    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_COUNTER_H_
+#define BITMAP_CUBICAL_COMPLEX_COUNTER_H_
+
+#include <iostream>
+#include <vector>
+
+namespace Gudhi {
+
+namespace Cubical_complex {
+
+/**
+ * This is an implementation of a counter being a vector of integers. 
+ * The constructor of the class takes as an input two vectors W and V. 
+ * It assumes that W < V coordinatewise.
+ * If the initial counter W is not specified, it is assumed to be vector of zeros. 
+ * The class allows to iterate between W and V by using increment() function.
+ * The increment() function returns a bool value. 
+ * The current counter reach the end counter V if the value returned by the increment function is FALSE.
+ * This class is needed for the implementation of a bitmapCubicalComplex.
+ **/
+
+class counter {
+ public:
+  /**
+   * Constructor of a counter class. It takes only the parameter which is the end value of the counter. 
+   * The default beginning value is a vector of the same length as the endd, filled-in with zeros.
+   **/
+  counter(const std::vector<unsigned>& endd) : begin(endd.size(), 0), end(endd), current(endd.size(), 0) { }
+
+  /**
+   * Constructor of a counter class. It takes as the input beginn and end vector. 
+   * It assumes that begin vector is lexicographically below the end vector.
+   **/
+  counter(const std::vector< unsigned >& beginn, const std::vector< unsigned >& endd) : begin(beginn), end(endd), current(endd.size(), 0) {
+    if (beginn.size() != endd.size())
+      throw "In constructor of a counter, begin and end vectors do not have the same size. Program terminate";
+  }
+
+  /**
+   * Function to increment the counter. If the value returned by the function is true, 
+   * then the incrementation process was successful.
+   * If the value of the function is false, that means, that the counter have reached its end-value.
+   **/
+  bool increment() {
+    size_t i = 0;
+    while ((i != this->end.size()) && (this->current[i] == this->end[i])) {
+      ++i;
+    }
+
+    if (i == this->end.size())return false;
+    ++this->current[i];
+    for (size_t j = 0; j != i; ++j) {
+      this->current[j] = this->begin[j];
+    }
+    return true;
+  }
+
+  /**
+   * Function to check if we are at the end of counter.
+   **/
+  bool isFinal() {
+    for (size_t i = 0; i != this->current.size(); ++i) {
+      if (this->current[i] == this->end[i])return true;
+    }
+    return false;
+  }
+
+  /**
+   * Function required in the implementation of bitmapCubicalComplexWPeriodicBoundaryCondition. 
+   * Its aim is to find an counter corresponding to the element the following
+   * boundary element is identified with when periodic boundary conditions are imposed.
+   **/
+  std::vector< unsigned > find_opposite(const std::vector< bool >& directionsForPeriodicBCond) {
+    std::vector< unsigned > result;
+    for (size_t i = 0; i != this->current.size(); ++i) {
+      if ((this->current[i] == this->end[i]) && (directionsForPeriodicBCond[i] == true)) {
+        result.push_back(this->begin[i]);
+      } else {
+        result.push_back(this->current[i]);
+      }
+    }
+    return result;
+  }
+
+  /**
+   * Function checking at which positions the current value of a counter is the final value of the counter.
+   **/
+  std::vector< bool > directions_of_finals() {
+    std::vector< bool > result;
+    for (size_t i = 0; i != this->current.size(); ++i) {
+      if (this->current[i] == this->end[i]) {
+        result.push_back(true);
+      } else {
+        result.push_back(false);
+      }
+    }
+    return result;
+  }
+
+  /**
+   * Function to write counter to the stream.
+   **/
+  friend std::ostream& operator<<(std::ostream& out, const counter& c) {
+    // std::cerr << "c.current.size() : " << c.current.size() << endl;
+    for (size_t i = 0; i != c.current.size(); ++i) {
+      out << c.current[i] << " ";
+    }
+    return out;
+  }
+
+ private:
+  std::vector< unsigned > begin;
+  std::vector< unsigned > end;
+  std::vector< unsigned > current;
+};
+
+}  // namespace Cubical_complex
+
+}  // namespace Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_COUNTER_H_
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
new file mode 100644
index 00000000..7294da98
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -0,0 +1,819 @@
+/*    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/Bitmap_cubical_complex/counter.h>
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+#include <limits>
+#include <utility>  // for pair<>
+
+namespace Gudhi {
+
+namespace Cubical_complex {
+
+/**
+ * @class Bitmap_cubical_complex_base
+ * @brief Cubical complex represented as a bitmap, class with basic implementation.
+ * @ingroup 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:
+  typedef T filtration_type;
+
+  /**
+   *Default constructor
+   **/
+  Bitmap_cubical_complex_base() :
+      total_number_of_cells(0) { }
+  /**
+   * 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(const 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(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells);
+
+  /**
+   * Destructor of the Bitmap_cubical_complex_base class.
+   **/
+  virtual ~Bitmap_cubical_complex_base() { }
+
+  /**
+   * 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.
+   */
+  virtual 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.
+   **/
+  virtual 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.
+   * This allows reading and changing the value of filtration. Note that if the value of a filtration is changed, the
+   * code do not check if we have a filtration or not. i.e. it do not check if the value of a filtration of a cell is
+   * not smaller than the value of a filtration of its boundary and not greater than the value of its coboundary.
+   **/
+  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()const {
+    return this->data.size();
+  }
+
+  /**
+   * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the
+   * structure. This procedure is used for debugging purposes.
+   **/
+  template <typename K>
+  friend std::ostream& operator<<(std::ostream & os, const Bitmap_cubical_complex_base<K>& b);
+
+  /**
+   * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values
+   * equally distributed in the range of data.
+   * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute
+   * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to
+   * some number of bins. The function put_data_to_bins( size_t number_of_bins ) is designed for that purpose.
+   * The parameter of the function is the number of bins (distinct values) we want to have in the cubical complex.
+   **/
+  void put_data_to_bins(size_t number_of_bins);
+
+  /**
+   * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values
+   * equally distributed in the range of data.
+   * Sometimes if most of the cells have different birth-death times, the performance of the algorithms to compute
+   * persistence gets worst. When dealing with this type of data, one may want to put different values on cells to
+   * some number of bins. The function put_data_to_bins( T diameter_of_bin ) is designed for that purpose.
+   * The parameter of it is the diameter of each bin. Note that the bottleneck distance between the persistence
+   * diagram of the cubical complex before and after using such a function will be bounded by the parameter
+   * diameter_of_bin.
+   **/
+  void put_data_to_bins(T diameter_of_bin);
+
+  /**
+   * Functions to find min and max values of filtration.
+   **/
+  std::pair< T, T > min_max_filtration();
+
+  // ITERATORS
+
+  /**
+   * Iterator through all cells in the complex (in order they appear in the structure -- i.e.
+   * in lexicographical order).
+   **/
+  class All_cells_iterator : std::iterator< std::input_iterator_tag, T > {
+   public:
+    All_cells_iterator() {
+      this->counter = 0;
+    }
+
+    All_cells_iterator operator++() {
+      // first find first element of the counter that can be increased:
+      ++this->counter;
+      return *this;
+    }
+
+    All_cells_iterator operator++(int) {
+      All_cells_iterator result = *this;
+      ++(*this);
+      return result;
+    }
+
+    All_cells_iterator& operator=(const All_cells_iterator& rhs) {
+      this->counter = rhs.counter;
+      return *this;
+    }
+
+    bool operator==(const All_cells_iterator& rhs)const {
+      if (this->counter != rhs.counter)return false;
+      return true;
+    }
+
+    bool operator!=(const All_cells_iterator& rhs)const {
+      return !(*this == rhs);
+    }
+
+    /*
+     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
+     * an argument of the following functions:
+     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
+     * boundary and coboundary and dimension
+     * and in function get_cell_data to get a filtration of a cell.
+     */
+    size_t operator*() {
+      return this->counter;
+    }
+    friend class Bitmap_cubical_complex_base;
+   protected:
+    size_t counter;
+  };
+
+  /**
+   * Function returning a All_cells_iterator to the first cell of the bitmap.
+   **/
+  All_cells_iterator all_cells_iterator_begin() {
+    All_cells_iterator a;
+    return a;
+  }
+
+  /**
+   * Function returning a All_cells_iterator to the last cell of the bitmap.
+   **/
+  All_cells_iterator all_cells_iterator_end() {
+    All_cells_iterator a;
+    a.counter = this->data.size();
+    return a;
+  }
+
+  /**
+   * All_cells_range class provides ranges for All_cells_iterator
+   **/
+  class All_cells_range {
+   public:
+    All_cells_range(Bitmap_cubical_complex_base* b) : b(b) { }
+
+    All_cells_iterator begin() {
+      return b->all_cells_iterator_begin();
+    }
+
+    All_cells_iterator end() {
+      return b->all_cells_iterator_end();
+    }
+   private:
+    Bitmap_cubical_complex_base<T>* b;
+  };
+
+  All_cells_range all_cells_range() {
+    return All_cells_range(this);
+  }
+
+
+  /**
+   * Boundary_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< size_t >::const_iterator Boundary_iterator;
+  typedef typename std::vector< size_t > Boundary_range;
+
+  /**
+   * boundary_simplex_range creates an object of a Boundary_simplex_range class
+   * that provides ranges for the Boundary_simplex_iterator.
+   **/
+  Boundary_range boundary_range(size_t sh) {
+    return this->get_boundary_of_a_cell(sh);
+  }
+
+  /**
+   * Coboundary_range class provides ranges for boundary iterators.
+   **/
+  typedef typename std::vector< size_t >::const_iterator Coboundary_iterator;
+  typedef typename std::vector< size_t > Coboundary_range;
+
+  /**
+   * boundary_simplex_range creates an object of a Boundary_simplex_range class
+   * that provides ranges for the Boundary_simplex_iterator.
+   **/
+  Coboundary_range coboundary_range(size_t sh) {
+    return this->get_coboundary_of_a_cell(sh);
+  }
+
+  /**
+   * 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, T > {
+   public:
+    Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) {
+      this->counter = std::vector<size_t>(b.dimension());
+      // std::fill( this->counter.begin() , this->counter.end() , 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);
+    }
+
+    /*
+     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
+     * an argument of the following functions:
+     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
+     * boundary and coboundary and dimension
+     * and in function get_cell_data to get a filtration of a cell.
+     */
+    size_t operator*() {
+      return this->compute_index_in_bitmap();
+    }
+
+    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) {
+        std::cout << this->counter[i] << " ";
+      }
+    }
+    friend class Bitmap_cubical_complex_base;
+   protected:
+    std::vector< size_t > counter;
+    Bitmap_cubical_complex_base& b;
+  };
+
+  /**
+   * Function returning a Top_dimensional_cells_iterator to the first top dimensional cell of the bitmap.
+   **/
+  Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin() {
+    Top_dimensional_cells_iterator a(*this);
+    return a;
+  }
+
+  /**
+   * Function returning a Top_dimensional_cells_iterator to the last top dimensional cell of the bitmap.
+   **/
+  Top_dimensional_cells_iterator top_dimensional_cells_iterator_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;
+  }
+
+  /**
+   * Top_dimensional_cells_iterator_range class provides ranges for Top_dimensional_cells_iterator_range
+   **/
+  class Top_dimensional_cells_range {
+   public:
+    Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) { }
+
+    Top_dimensional_cells_iterator begin() {
+      return b->top_dimensional_cells_iterator_begin();
+    }
+
+    Top_dimensional_cells_iterator end() {
+      return b->top_dimensional_cells_iterator_end();
+    }
+   private:
+    Bitmap_cubical_complex_base<T>* b;
+  };
+
+  Top_dimensional_cells_range top_dimensional_cells_range() {
+    return Top_dimensional_cells_range(this);
+  }
+
+
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+
+  inline size_t number_cells()const {
+    return this->total_number_of_cells;
+  }
+
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+  //****************************************************************************************************************//
+
+ protected:
+  std::vector<unsigned> sizes;
+  std::vector<unsigned> multipliers;
+  std::vector<T> data;
+  size_t total_number_of_cells;
+
+  void set_up_containers(const 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;
+    }
+    this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max());
+    this->total_number_of_cells = multiplier;
+  }
+
+  size_t compute_position_in_bitmap(const 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;
+    counter.reserve(this->sizes.size());
+    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;
+  }
+  void read_perseus_style_file(const char* perseus_style_file);
+  void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions,
+                                                        const std::vector<T>& top_dimensional_cells);
+  Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions);
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions);
+  Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions,
+                              const std::vector<T>& top_dimensional_cells,
+                              std::vector<bool> directions);
+};
+
+template <typename T>
+void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) {
+  bool bdg = false;
+
+  std::pair< T, T > min_max = this->min_max_filtration();
+  T dx = (min_max.second - min_max.first) / (T) number_of_bins;
+
+  // now put the data into the appropriate bins:
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (bdg) {
+      std::cerr << "Before binning : " << this->data[i] << std::endl;
+    }
+    this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins;
+    if (bdg) {
+      std::cerr << "After binning : " << this->data[i] << std::endl;
+      getchar();
+    }
+  }
+}
+
+template <typename T>
+void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) {
+  bool bdg = false;
+  std::pair< T, T > min_max = this->min_max_filtration();
+
+  size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin;
+  // now put the data into the appropriate bins:
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (bdg) {
+      std::cerr << "Before binning : " << this->data[i] << std::endl;
+    }
+    this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins;
+    if (bdg) {
+      std::cerr << "After binning : " << this->data[i] << std::endl;
+      getchar();
+    }
+  }
+}
+
+template <typename T>
+std::pair< T, T > Bitmap_cubical_complex_base<T>::min_max_filtration() {
+  std::pair< T, T > min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min());
+  for (size_t i = 0; i != this->data.size(); ++i) {
+    if (this->data[i] < min_max.first)min_max.first = this->data[i];
+    if (this->data[i] > min_max.second)min_max.second = this->data[i];
+  }
+  return min_max;
+}
+
+template <typename K>
+std::ostream& operator<<(std::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
+(const std::vector<unsigned>& sizes) {
+  this->set_up_containers(sizes);
+}
+
+template <typename T>
+void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(const 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()) {
+    std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<size_t> sizes_in_following_directions"
+        << ", std::vector<T> 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."
+        << std::endl;
+    throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions,"
+           "std::vector<T> 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_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
+    this->get_cell_data(*it) = top_dimensional_cells[index];
+    ++index;
+  }
+  this->impose_lower_star_filtration();
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
+(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) {
+  this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
+}
+
+template <typename T>
+void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) {
+  bool dbg = false;
+  std::ifstream inFiltration;
+  inFiltration.open(perseus_style_file);
+  unsigned dimensionOfData;
+  inFiltration >> dimensionOfData;
+
+  if (dbg) {
+    std::cerr << "dimensionOfData : " << dimensionOfData << std::endl;
+    getchar();
+  }
+
+  std::vector<unsigned> sizes;
+  sizes.reserve(dimensionOfData);
+  for (size_t i = 0; i != dimensionOfData; ++i) {
+    unsigned size_in_this_dimension;
+    inFiltration >> size_in_this_dimension;
+    sizes.push_back(size_in_this_dimension);
+    if (dbg) {
+      std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
+    }
+  }
+  this->set_up_containers(sizes);
+
+  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
+  it = this->top_dimensional_cells_iterator_begin();
+
+  while (!inFiltration.eof()) {
+    T filtrationLevel;
+    inFiltration >> filtrationLevel;
+    if (dbg) {
+      std::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 << std::endl;
+    }
+    this->get_cell_data(*it) = filtrationLevel;
+    ++it;
+  }
+  inFiltration.close();
+  this->impose_lower_star_filtration();
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file,
+                                                            std::vector<bool> directions) {
+  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
+  // conditions.
+  // It ignores the last parameter of the function.
+  this->read_perseus_style_file(perseus_style_file);
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes,
+                                                            std::vector<bool> directions) {
+  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
+  // conditions.
+  // It ignores the last parameter of the function.
+  this->set_up_containers(sizes);
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions,
+                                                            const std::vector<T>& top_dimensional_cells,
+                                                            std::vector<bool> directions) {
+  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
+  // conditions.
+  // It ignores the last parameter of the function.
+  this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, top_dimensional_cells);
+}
+
+template <typename T>
+Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file) {
+  this->read_perseus_style_file(perseus_style_file);
+}
+
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const {
+  std::vector< size_t > boundary_elements;
+
+  // Speed traded of for memory. Check if it is better in practice.
+  boundary_elements.reserve(this->dimension()*2);
+
+  size_t cell1 = cell;
+  for (size_t i = this->multipliers.size(); i != 0; --i) {
+    unsigned position = cell1 / this->multipliers[i - 1];
+    if (position % 2 == 1) {
+      boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+      boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+    }
+    cell1 = cell1 % this->multipliers[i - 1];
+  }
+  return boundary_elements;
+}
+
+template <typename T>
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell)const {
+  std::vector<unsigned> 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 (position % 2 == 0) {
+      if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
+        coboundary_elements.push_back(cell - this->multipliers[i - 1]);
+      }
+      if (
+          (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
+        coboundary_elements.push_back(cell + this->multipliers[i - 1]);
+      }
+    }
+    cell1 = cell1 % this->multipliers[i - 1];
+  }
+  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) 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 / this->multipliers[i - 1];
+
+    if (dbg) {
+      std::cerr << "i-1 :" << i - 1 << std::endl;
+      std::cerr << "cell : " << cell << std::endl;
+      std::cerr << "position : " << position << std::endl;
+      std::cerr << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl;
+      getchar();
+    }
+
+    if (position % 2 == 1) {
+      if (dbg) std::cerr << "Nonzero length in this direction \n";
+      dimension++;
+    }
+    cell = cell % this->multipliers[i - 1];
+  }
+  return dimension;
+}
+
+template <typename T>
+inline 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);
+
+  size_t size_to_reserve = 1;
+  for (size_t i = 0; i != this->multipliers.size(); ++i) {
+    size_to_reserve *= (size_t) ((this->multipliers[i] - 1) / 2);
+  }
+
+  std::vector<size_t> indices_to_consider;
+  indices_to_consider.reserve(size_to_reserve);
+  // 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_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
+    indices_to_consider.push_back(it.compute_index_in_bitmap());
+  }
+
+  while (indices_to_consider.size()) {
+    if (dbg) {
+      std::cerr << "indices_to_consider in this iteration \n";
+      for (size_t i = 0; i != indices_to_consider.size(); ++i) {
+        std::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 (dbg) {
+          std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[ bd[boundaryIt] ]
+              << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[ indices_to_consider[i] ]
+              << std::endl;
+          getchar();
+        }
+        if (this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ]) {
+          this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];
+          if (dbg) {
+            std::cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : "
+                << this->data[ indices_to_consider[i] ] << std::endl;
+            getchar();
+          }
+        }
+        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;
+  }
+}
+
+}  // namespace Cubical_complex
+
+}  // namespace Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_BASE_H_
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
new file mode 100644
index 00000000..a446c0e8
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
@@ -0,0 +1,309 @@
+/*    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_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
+#define BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
+
+#include <gudhi/Bitmap_cubical_complex_base.h>
+
+#include <cmath>
+#include <limits>  // for numeric_limits<>
+#include <vector>
+
+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
+
+/**
+ * @class Bitmap_cubical_complex_periodic_boundary_conditions_base
+ * @brief Cubical complex with periodic boundary conditions represented as a bitmap.
+ * @ingroup cubical_complex
+ */
+/**
+ * This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are
+ * identical to the functions from Bitmap_cubical_complex_base.
+ * The ones that needed to be updated are the constructors and get_boundary_of_a_cell and get_coboundary_of_a_cell.
+ */
+template <typename T>
+class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T> {
+ public:
+  // constructors that take an extra parameter:
+
+  /**
+   * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base() { }
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following
+   * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If
+   * at i-th position of this vector there is true value, that means that periodic boundary conditions are to be
+   * imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction
+   * i.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes,
+                                                           const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus
+   * style file as an input. Please consult the documentation about the specification of the file.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseusStyleFile);
+  /**
+   * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following
+   * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of top dimensional
+   * cells (ordered lexicographically) and (3) vector of booleans. If at i-th position of this vector there is true
+   * value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the
+   * periodic boundary conditions will not be imposed in the direction i.
+   */
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
+                                                           const std::vector<T>& topDimensionalCells,
+                                                           const std::vector< bool >& directions_in_which_periodic_b_cond_are_to_be_imposed);
+
+  /**
+   * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   **/
+  virtual ~Bitmap_cubical_complex_periodic_boundary_conditions_base() {}
+
+  // overwritten methods co compute boundary and coboundary
+  /**
+   * A version of a function that return boundary of a given cell for an object of
+   * Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell) const;
+
+  /**
+   * A version of a function that return coboundary of a given cell for an object of
+   * Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   */
+  virtual std::vector< size_t > get_coboundary_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<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);
+
+      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<T>(multiplier, std::numeric_limits<T>::max());
+    this->total_number_of_cells = multiplier;
+  }
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes);
+  Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
+                                                           const std::vector<T>& topDimensionalCells);
+  void construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions,
+                                                        const std::vector<T>& topDimensionalCells,
+                                                        const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed);
+};
+
+template <typename T>
+void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions,
+                                                                                                                   const std::vector<T>& topDimensionalCells,
+                                                                                                                   const 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 (auto it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
+    this->get_cell_data(*it) = topDimensionalCells[i];
+    ++i;
+  }
+  this->impose_lower_star_filtration();
+}
+
+template <typename T>
+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes,
+                                                                                                                      const 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 <typename T>
+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseus_style_file) {
+  // for Perseus style files:
+  bool dbg = false;
+
+  std::ifstream inFiltration;
+  inFiltration.open(perseus_style_file);
+  unsigned dimensionOfData;
+  inFiltration >> dimensionOfData;
+
+  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensionOfData, false);
+
+  std::vector<unsigned> 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<T>::Top_dimensional_cells_iterator it(*this);
+  it = this->top_dimensional_cells_iterator_begin();
+
+  while (!inFiltration.eof()) {
+    double filtrationLevel;
+    inFiltration >> filtrationLevel;
+    if (inFiltration.eof())break;
+
+    if (dbg) {
+      std::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 << std::endl;
+    }
+    this->get_cell_data(*it) = filtrationLevel;
+    ++it;
+  }
+  inFiltration.close();
+  this->impose_lower_star_filtration();
+}
+
+template <typename T>
+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes) {
+  this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false);
+  this->set_up_containers(sizes);
+}
+
+template <typename T>
+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
+                                                                                                                      const std::vector<T>& topDimensionalCells) {
+  std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false);
+  this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells,
+                                                         directions_in_which_periodic_b_cond_are_to_be_imposed);
+}
+
+template <typename T>
+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::
+Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions,
+                                                         const std::vector<T>& topDimensionalCells,
+                                                         const 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 <typename T>
+std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell(size_t cell) const {
+  bool dbg = false;
+  if (dbg) {
+    std::cerr << "Computations of boundary of a cell : " << cell << std::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]) {
+        // std::cerr << "A\n";
+        boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+        boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+        if (dbg) {
+          std::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) {
+          // std::cerr << "B\n";
+          boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+          boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+          if (dbg) {
+            std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " ";
+          }
+        } else {
+          // std::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) {
+            std::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 <typename T>
+std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell(size_t cell) const {
+  std::vector<unsigned> 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;
+}
+
+}  // namespace Cubical_complex
+
+}  // namespace Gudhi
+
+#endif  // BITMAP_CUBICAL_COMPLEX_PERIODIC_BOUNDARY_CONDITIONS_BASE_H_
diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
new file mode 100644
index 00000000..a9162cee
--- /dev/null
+++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
@@ -0,0 +1,1378 @@
+/*    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 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/>.
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "cubical_complex"
+#include <boost/test/unit_test.hpp>
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+
+typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base;
+typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex;
+
+typedef Gudhi::Cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>
+Bitmap_cubical_complex_periodic_boundary_conditions_base;
+typedef Gudhi::Cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_periodic_boundary_conditions_base>
+Bitmap_cubical_complex_periodic_boundary_conditions;
+
+BOOST_AUTO_TEST_CASE(check_dimension) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  BOOST_CHECK(increasing.dimension() == 2);
+}
+
+BOOST_AUTO_TEST_CASE(topDimensionalCellsIterator_test) {
+  std::vector< double > expectedFiltrationValues1({0, 0, 0, 0, 100, 0, 0, 0, 0});
+
+  std::vector< double > expectedFiltrationValues2({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector< double > oneDimensionalCycle({0, 0, 0, 0, 100, 0, 0, 0, 0});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  Bitmap_cubical_complex hole(dimensions, oneDimensionalCycle);
+
+
+  int i = 0;
+  for (Bitmap_cubical_complex::Top_dimensional_cells_iterator
+       it = increasing.top_dimensional_cells_iterator_begin(); it != increasing.top_dimensional_cells_iterator_end(); ++it) {
+    BOOST_CHECK(increasing.get_cell_data(*it) == expectedFiltrationValues2[i]);
+    ++i;
+  }
+  i = 0;
+  for (Bitmap_cubical_complex::Top_dimensional_cells_iterator
+       it = hole.top_dimensional_cells_iterator_begin(); it != hole.top_dimensional_cells_iterator_end(); ++it) {
+    BOOST_CHECK(hole.get_cell_data(*it) == expectedFiltrationValues1[i]);
+    ++i;
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
+  std::vector<double> boundary0;
+  std::vector<double> boundary1;
+  boundary1.push_back(0);
+  boundary1.push_back(2);
+  std::vector<double> boundary2;
+  std::vector<double> boundary3;
+  boundary3.push_back(2);
+  boundary3.push_back(4);
+  std::vector<double> boundary4;
+  std::vector<double> boundary5;
+  boundary5.push_back(4);
+  boundary5.push_back(6);
+  std::vector<double> boundary6;
+  std::vector<double> boundary7;
+  boundary7.push_back(0);
+  boundary7.push_back(14);
+  std::vector<double> boundary8;
+  boundary8.push_back(1);
+  boundary8.push_back(15);
+  boundary8.push_back(7);
+  boundary8.push_back(9);
+  std::vector<double> boundary9;
+  boundary9.push_back(2);
+  boundary9.push_back(16);
+  std::vector<double> boundary10;
+  boundary10.push_back(3);
+  boundary10.push_back(17);
+  boundary10.push_back(9);
+  boundary10.push_back(11);
+  std::vector<double> boundary11;
+  boundary11.push_back(4);
+  boundary11.push_back(18);
+  std::vector<double> boundary12;
+  boundary12.push_back(5);
+  boundary12.push_back(19);
+  boundary12.push_back(11);
+  boundary12.push_back(13);
+  std::vector<double> boundary13;
+  boundary13.push_back(6);
+  boundary13.push_back(20);
+  std::vector<double> boundary14;
+  std::vector<double> boundary15;
+  boundary15.push_back(14);
+  boundary15.push_back(16);
+  std::vector<double> boundary16;
+  std::vector<double> boundary17;
+  boundary17.push_back(16);
+  boundary17.push_back(18);
+  std::vector<double> boundary18;
+  std::vector<double> boundary19;
+  boundary19.push_back(18);
+  boundary19.push_back(20);
+  std::vector<double> boundary20;
+  std::vector<double> boundary21;
+  boundary21.push_back(14);
+  boundary21.push_back(28);
+  std::vector<double> boundary22;
+  boundary22.push_back(15);
+  boundary22.push_back(29);
+  boundary22.push_back(21);
+  boundary22.push_back(23);
+  std::vector<double> boundary23;
+  boundary23.push_back(16);
+  boundary23.push_back(30);
+  std::vector<double> boundary24;
+  boundary24.push_back(17);
+  boundary24.push_back(31);
+  boundary24.push_back(23);
+  boundary24.push_back(25);
+  std::vector<double> boundary25;
+  boundary25.push_back(18);
+  boundary25.push_back(32);
+  std::vector<double> boundary26;
+  boundary26.push_back(19);
+  boundary26.push_back(33);
+  boundary26.push_back(25);
+  boundary26.push_back(27);
+  std::vector<double> boundary27;
+  boundary27.push_back(20);
+  boundary27.push_back(34);
+  std::vector<double> boundary28;
+  std::vector<double> boundary29;
+  boundary29.push_back(28);
+  boundary29.push_back(30);
+  std::vector<double> boundary30;
+  std::vector<double> boundary31;
+  boundary31.push_back(30);
+  boundary31.push_back(32);
+  std::vector<double> boundary32;
+  std::vector<double> boundary33;
+  boundary33.push_back(32);
+  boundary33.push_back(34);
+  std::vector<double> boundary34;
+  std::vector<double> boundary35;
+  boundary35.push_back(28);
+  boundary35.push_back(42);
+  std::vector<double> boundary36;
+  boundary36.push_back(29);
+  boundary36.push_back(43);
+  boundary36.push_back(35);
+  boundary36.push_back(37);
+  std::vector<double> boundary37;
+  boundary37.push_back(30);
+  boundary37.push_back(44);
+  std::vector<double> boundary38;
+  boundary38.push_back(31);
+  boundary38.push_back(45);
+  boundary38.push_back(37);
+  boundary38.push_back(39);
+  std::vector<double> boundary39;
+  boundary39.push_back(32);
+  boundary39.push_back(46);
+  std::vector<double> boundary40;
+  boundary40.push_back(33);
+  boundary40.push_back(47);
+  boundary40.push_back(39);
+  boundary40.push_back(41);
+  std::vector<double> boundary41;
+  boundary41.push_back(34);
+  boundary41.push_back(48);
+  std::vector<double> boundary42;
+  std::vector<double> boundary43;
+  boundary43.push_back(42);
+  boundary43.push_back(44);
+  std::vector<double> boundary44;
+  std::vector<double> boundary45;
+  boundary45.push_back(44);
+  boundary45.push_back(46);
+  std::vector<double> boundary46;
+  std::vector<double> boundary47;
+  boundary47.push_back(46);
+  boundary47.push_back(48);
+  std::vector<double> boundary48;
+  std::vector< std::vector<double> > boundaries;
+  boundaries.push_back(boundary0);
+  boundaries.push_back(boundary1);
+  boundaries.push_back(boundary2);
+  boundaries.push_back(boundary3);
+  boundaries.push_back(boundary4);
+  boundaries.push_back(boundary5);
+  boundaries.push_back(boundary6);
+  boundaries.push_back(boundary7);
+  boundaries.push_back(boundary8);
+  boundaries.push_back(boundary9);
+  boundaries.push_back(boundary10);
+  boundaries.push_back(boundary11);
+  boundaries.push_back(boundary12);
+  boundaries.push_back(boundary13);
+  boundaries.push_back(boundary14);
+  boundaries.push_back(boundary15);
+  boundaries.push_back(boundary16);
+  boundaries.push_back(boundary17);
+  boundaries.push_back(boundary18);
+  boundaries.push_back(boundary19);
+  boundaries.push_back(boundary20);
+  boundaries.push_back(boundary21);
+  boundaries.push_back(boundary22);
+  boundaries.push_back(boundary23);
+  boundaries.push_back(boundary24);
+  boundaries.push_back(boundary25);
+  boundaries.push_back(boundary26);
+  boundaries.push_back(boundary27);
+  boundaries.push_back(boundary28);
+  boundaries.push_back(boundary29);
+  boundaries.push_back(boundary30);
+  boundaries.push_back(boundary31);
+  boundaries.push_back(boundary32);
+  boundaries.push_back(boundary33);
+  boundaries.push_back(boundary34);
+  boundaries.push_back(boundary35);
+  boundaries.push_back(boundary36);
+  boundaries.push_back(boundary37);
+  boundaries.push_back(boundary38);
+  boundaries.push_back(boundary39);
+  boundaries.push_back(boundary40);
+  boundaries.push_back(boundary41);
+  boundaries.push_back(boundary42);
+  boundaries.push_back(boundary43);
+  boundaries.push_back(boundary44);
+  boundaries.push_back(boundary45);
+  boundaries.push_back(boundary46);
+  boundaries.push_back(boundary47);
+  boundaries.push_back(boundary48);
+
+
+
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+  for (size_t i = 0; i != increasing.size(); ++i) {
+    std::vector< size_t > bd = increasing.get_boundary_of_a_cell(i);
+    for (size_t j = 0; j != bd.size(); ++j) {
+      BOOST_CHECK(boundaries[i][j] == bd[j]);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_2) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+
+  std::vector<double> coboundaryElements;
+  coboundaryElements.push_back(7);
+  coboundaryElements.push_back(1);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(9);
+  coboundaryElements.push_back(1);
+  coboundaryElements.push_back(3);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(11);
+  coboundaryElements.push_back(3);
+  coboundaryElements.push_back(5);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(13);
+  coboundaryElements.push_back(5);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(7);
+  coboundaryElements.push_back(21);
+  coboundaryElements.push_back(15);
+  coboundaryElements.push_back(8);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(9);
+  coboundaryElements.push_back(23);
+  coboundaryElements.push_back(15);
+  coboundaryElements.push_back(17);
+  coboundaryElements.push_back(10);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(11);
+  coboundaryElements.push_back(25);
+  coboundaryElements.push_back(17);
+  coboundaryElements.push_back(19);
+  coboundaryElements.push_back(12);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(13);
+  coboundaryElements.push_back(27);
+  coboundaryElements.push_back(19);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(21);
+  coboundaryElements.push_back(35);
+  coboundaryElements.push_back(29);
+  coboundaryElements.push_back(22);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(23);
+  coboundaryElements.push_back(37);
+  coboundaryElements.push_back(29);
+  coboundaryElements.push_back(31);
+  coboundaryElements.push_back(24);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(25);
+  coboundaryElements.push_back(39);
+  coboundaryElements.push_back(31);
+  coboundaryElements.push_back(33);
+  coboundaryElements.push_back(26);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(27);
+  coboundaryElements.push_back(41);
+  coboundaryElements.push_back(33);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(35);
+  coboundaryElements.push_back(43);
+  coboundaryElements.push_back(36);
+  coboundaryElements.push_back(37);
+  coboundaryElements.push_back(43);
+  coboundaryElements.push_back(45);
+  coboundaryElements.push_back(38);
+  coboundaryElements.push_back(39);
+  coboundaryElements.push_back(45);
+  coboundaryElements.push_back(47);
+  coboundaryElements.push_back(40);
+  coboundaryElements.push_back(41);
+  coboundaryElements.push_back(47);
+  size_t number = 0;
+  for (size_t i = 0; i != increasing.size(); ++i) {
+    std::vector< size_t > bd = increasing.get_coboundary_of_a_cell(i);
+    for (size_t j = 0; j != bd.size(); ++j) {
+      BOOST_CHECK(coboundaryElements[number] == bd[j]);
+      ++number;
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(compute_boundary_test_3) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+  std::vector<unsigned> dim;
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(0);
+
+  for (size_t i = 0; i != increasing.size(); ++i) {
+    BOOST_CHECK(increasing.get_dimension_of_a_cell(i) == dim[i]);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) {
+  std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+  std::vector<unsigned> dimensions({3, 3});
+
+  Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells);
+
+  std::vector< unsigned > dim;
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+  dim.push_back(0);
+  dim.push_back(1);
+  dim.push_back(1);
+  dim.push_back(2);
+
+  std::vector<double> fil;
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(1);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(2);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(3);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(4);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(5);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(6);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(7);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(8);
+  fil.push_back(9);
+  fil.push_back(9);
+  fil.push_back(9);
+  fil.push_back(9);
+
+
+  Bitmap_cubical_complex::Filtration_simplex_range range = increasing.filtration_simplex_range();
+  size_t position = 0;
+  for (Bitmap_cubical_complex::Filtration_simplex_iterator it = range.begin(); it != range.end(); ++it) {
+    BOOST_CHECK(increasing.dimension(*it) == dim[position]);
+    BOOST_CHECK(increasing.filtration(*it) == fil[position]);
+    ++position;
+  }
+}
+
+BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) {
+  std::vector< double > filtration({0, 0, 0, 0});
+
+  std::vector<unsigned> dimensions({2, 2});
+
+  std::vector<bool> periodic_directions({true, true});
+
+  Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions);
+  BOOST_CHECK(cmplx.dimension() == 2);
+
+
+  std::vector<double> boundary0;
+  std::vector<double> boundary1;
+  boundary1.push_back(0);
+  boundary1.push_back(2);
+  std::vector<double> boundary2;
+  std::vector<double> boundary3;
+  boundary3.push_back(2);
+  boundary3.push_back(0);
+  std::vector<double> boundary4;
+  boundary4.push_back(0);
+  boundary4.push_back(8);
+  std::vector<double> boundary5;
+  boundary5.push_back(1);
+  boundary5.push_back(9);
+  boundary5.push_back(4);
+  boundary5.push_back(6);
+  std::vector<double> boundary6;
+  boundary6.push_back(2);
+  boundary6.push_back(10);
+  std::vector<double> boundary7;
+  boundary7.push_back(3);
+  boundary7.push_back(11);
+  boundary7.push_back(6);
+  boundary7.push_back(4);
+  std::vector<double> boundary8;
+  std::vector<double> boundary9;
+  boundary9.push_back(8);
+  boundary9.push_back(10);
+  std::vector<double> boundary10;
+  std::vector<double> boundary11;
+  boundary11.push_back(10);
+  boundary11.push_back(8);
+  std::vector<double> boundary12;
+  boundary12.push_back(8);
+  boundary12.push_back(0);
+  std::vector<double> boundary13;
+  boundary13.push_back(9);
+  boundary13.push_back(1);
+  boundary13.push_back(12);
+  boundary13.push_back(14);
+  std::vector<double> boundary14;
+  boundary14.push_back(10);
+  boundary14.push_back(2);
+  std::vector<double> boundary15;
+  boundary15.push_back(11);
+  boundary15.push_back(3);
+  boundary15.push_back(14);
+  boundary15.push_back(12);
+
+  std::vector< std::vector<double> > boundaries;
+  boundaries.push_back(boundary0);
+  boundaries.push_back(boundary1);
+  boundaries.push_back(boundary2);
+  boundaries.push_back(boundary3);
+  boundaries.push_back(boundary4);
+  boundaries.push_back(boundary5);
+  boundaries.push_back(boundary6);
+  boundaries.push_back(boundary7);
+  boundaries.push_back(boundary8);
+  boundaries.push_back(boundary9);
+  boundaries.push_back(boundary10);
+  boundaries.push_back(boundary11);
+  boundaries.push_back(boundary12);
+  boundaries.push_back(boundary13);
+  boundaries.push_back(boundary14);
+  boundaries.push_back(boundary15);
+
+  for (size_t i = 0; i != cmplx.size(); ++i) {
+    std::vector< size_t > bd = cmplx.get_boundary_of_a_cell(i);
+    for (size_t j = 0; j != bd.size(); ++j) {
+      BOOST_CHECK(boundaries[i][j] == bd[j]);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) {
+  std::vector< double > filtration({0, 0, 0, 0});
+
+  std::vector<unsigned> dimensions({2, 2});
+
+  std::vector<bool> periodic_directions({true, true});
+
+  Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions);
+  BOOST_CHECK(cmplx.dimension() == 2);
+
+
+  std::vector<double> coboundary0;
+  coboundary0.push_back(4);
+  coboundary0.push_back(12);
+  coboundary0.push_back(1);
+  coboundary0.push_back(3);
+  std::vector<double> coboundary1;
+  coboundary1.push_back(5);
+  coboundary1.push_back(13);
+  std::vector<double> coboundary2;
+  coboundary2.push_back(6);
+  coboundary2.push_back(14);
+  coboundary2.push_back(1);
+  coboundary2.push_back(3);
+  std::vector<double> coboundary3;
+  coboundary3.push_back(7);
+  coboundary3.push_back(15);
+  std::vector<double> coboundary4;
+  coboundary4.push_back(5);
+  coboundary4.push_back(7);
+  std::vector<double> coboundary5;
+  std::vector<double> coboundary6;
+  coboundary6.push_back(5);
+  coboundary6.push_back(7);
+  std::vector<double> coboundary7;
+  std::vector<double> coboundary8;
+  coboundary8.push_back(4);
+  coboundary8.push_back(12);
+  coboundary8.push_back(9);
+  coboundary8.push_back(11);
+  std::vector<double> coboundary9;
+  coboundary9.push_back(5);
+  coboundary9.push_back(13);
+  std::vector<double> coboundary10;
+  coboundary10.push_back(6);
+  coboundary10.push_back(14);
+  coboundary10.push_back(9);
+  coboundary10.push_back(11);
+  std::vector<double> coboundary11;
+  coboundary11.push_back(7);
+  coboundary11.push_back(15);
+  std::vector<double> coboundary12;
+  coboundary12.push_back(13);
+  coboundary12.push_back(15);
+  std::vector<double> coboundary13;
+  std::vector<double> coboundary14;
+  coboundary14.push_back(13);
+  coboundary14.push_back(15);
+  std::vector<double> coboundary15;
+
+  std::vector< std::vector<double> > coboundaries;
+  coboundaries.push_back(coboundary0);
+  coboundaries.push_back(coboundary1);
+  coboundaries.push_back(coboundary2);
+  coboundaries.push_back(coboundary3);
+  coboundaries.push_back(coboundary4);
+  coboundaries.push_back(coboundary5);
+  coboundaries.push_back(coboundary6);
+  coboundaries.push_back(coboundary7);
+  coboundaries.push_back(coboundary8);
+  coboundaries.push_back(coboundary9);
+  coboundaries.push_back(coboundary10);
+  coboundaries.push_back(coboundary11);
+  coboundaries.push_back(coboundary12);
+  coboundaries.push_back(coboundary13);
+  coboundaries.push_back(coboundary14);
+  coboundaries.push_back(coboundary15);
+
+  for (size_t i = 0; i != cmplx.size(); ++i) {
+    std::vector< size_t > cbd = cmplx.get_coboundary_of_a_cell(i);
+    for (size_t j = 0; j != cbd.size(); ++j) {
+      BOOST_CHECK(coboundaries[i][j] == cbd[j]);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) {
+  std::vector< double > filtrationOrg({0, 1, 2, 3});
+
+  std::vector<unsigned> dimensions({2, 2});
+
+  std::vector<bool> periodic_directions({true, true});
+
+  Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtrationOrg, periodic_directions);
+  BOOST_CHECK(cmplx.dimension() == 2);
+
+
+  std::vector<double> filtration;
+  filtration.push_back(0);  // 0
+  filtration.push_back(0);  // 1
+  filtration.push_back(0);  // 2
+  filtration.push_back(1);  // 3
+  filtration.push_back(0);  // 4
+  filtration.push_back(0);  // 5
+  filtration.push_back(0);  // 6
+  filtration.push_back(1);  // 7
+  filtration.push_back(0);  // 8
+  filtration.push_back(0);  // 9
+  filtration.push_back(0);  // 10
+  filtration.push_back(1);  // 11
+  filtration.push_back(2);  // 12
+  filtration.push_back(2);  // 13
+  filtration.push_back(2);  // 14
+  filtration.push_back(3);  // 15
+
+
+  for (size_t i = 0; i != cmplx.size(); ++i) {
+    BOOST_CHECK(filtration[i] == cmplx.get_cell_data(i));
+  }
+}
+BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check)
+{
+    std::vector< double > expected_filtration;
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+
+    std::vector<unsigned> expected_dimension;
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+
+    std::vector<size_t> expected_boundary;
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(1);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(5);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(3);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(9);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(21);
+    expected_boundary.push_back(15);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(23);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(19);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(24);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(24);
+
+
+    std::vector<size_t> expected_coboundary;
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(23);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(23);
+
+
+
+    std::vector< unsigned > sizes(2);
+    sizes[0] = 2;
+    sizes[1] = 2;
+
+    std::vector< double > data(4);
+    data[0] = 0;
+    data[1] = 1;
+    data[2] = 2;
+    data[3] = 3;
+
+    Bitmap_cubical_complex_base ba( sizes , data );
+    int i = 0;
+    int bd_it = 0;
+    int cbd_it = 0;
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {
+        BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) );
+        BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) );
+
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            ++bd_it;
+        }
+
+        Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd )
+        {
+             BOOST_CHECK( expected_coboundary[cbd_it] == *cbd );
+            ++cbd_it;
+        }
+        ++i;
+    }
+}
+
+
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check_2)
+{
+    std::vector< double > expected_filtration;
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+
+    std::vector<unsigned> expected_dimension;
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+
+    std::vector<size_t> expected_boundary;
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(1);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(5);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(3);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(9);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(21);
+    expected_boundary.push_back(15);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(23);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(19);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(24);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(24);
+
+
+    std::vector<size_t> expected_coboundary;
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(23);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(23);
+
+
+
+    std::vector< unsigned > sizes(2);
+    sizes[0] = 2;
+    sizes[1] = 2;
+
+    std::vector< double > data(4);
+    data[0] = 0;
+    data[1] = 1;
+    data[2] = 2;
+    data[3] = 3;
+
+    Bitmap_cubical_complex_base ba( sizes , data );
+    int i = 0;
+    int bd_it = 0;
+    int cbd_it = 0;
+
+    Bitmap_cubical_complex_base::All_cells_range range(&ba);
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it )
+    {
+        BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) );
+        BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) );
+
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            ++bd_it;
+        }
+
+        Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd )
+        {
+             BOOST_CHECK( expected_coboundary[cbd_it] == *cbd );
+            ++cbd_it;
+        }
+        ++i;
+    }
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check)
+{
+    std::vector< double > expected_filtration;
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+
+    std::vector<unsigned> expected_dimension;
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(2);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+    expected_dimension.push_back(1);
+    expected_dimension.push_back(0);
+
+    std::vector<size_t> expected_boundary;
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(0);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(1);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(5);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(2);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(3);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(7);
+    expected_boundary.push_back(9);
+    expected_boundary.push_back(4);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(10);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(11);
+    expected_boundary.push_back(21);
+    expected_boundary.push_back(15);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(12);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(13);
+    expected_boundary.push_back(23);
+    expected_boundary.push_back(17);
+    expected_boundary.push_back(19);
+    expected_boundary.push_back(14);
+    expected_boundary.push_back(24);
+    expected_boundary.push_back(20);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(22);
+    expected_boundary.push_back(24);
+
+
+    std::vector<size_t> expected_coboundary;
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(1);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(3);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(5);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(6);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(7);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(11);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(8);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(9);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(13);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(15);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(16);
+    expected_coboundary.push_back(17);
+    expected_coboundary.push_back(21);
+    expected_coboundary.push_back(23);
+    expected_coboundary.push_back(18);
+    expected_coboundary.push_back(19);
+    expected_coboundary.push_back(23);
+
+
+
+    std::vector< unsigned > sizes(2);
+    sizes[0] = 2;
+    sizes[1] = 2;
+
+    std::vector< double > data(4);
+    data[0] = 0;
+    data[1] = 1;
+    data[2] = 2;
+    data[3] = 3;
+
+    Bitmap_cubical_complex_base ba( sizes , data );
+    int i = 0;
+    int bd_it = 0;
+    int cbd_it = 0;
+
+    Bitmap_cubical_complex_base::All_cells_range range = ba.all_cells_range();
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it )
+    {
+        BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) );
+        BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) );
+
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            ++bd_it;
+        }
+
+        Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd )
+        {
+             BOOST_CHECK( expected_coboundary[cbd_it] == *cbd );
+            ++cbd_it;
+        }
+        ++i;
+    }
+}
+
+BOOST_AUTO_TEST_CASE(Top_dimensional_cells_iterator_range_check)
+{
+    std::vector< double > expected_filtration;
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(0);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(1);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(2);
+    expected_filtration.push_back(3);
+    expected_filtration.push_back(3);
+
+
+    std::vector< unsigned > sizes(2);
+    sizes[0] = 2;
+    sizes[1] = 2;
+
+    std::vector< double > data(4);
+    data[0] = 0;
+    data[1] = 1;
+    data[2] = 2;
+    data[3] = 3;
+
+    Bitmap_cubical_complex_base ba( sizes , data );
+    int i = 0;
+
+    Bitmap_cubical_complex_base::Top_dimensional_cells_range range = ba.top_dimensional_cells_range();
+    for ( Bitmap_cubical_complex_base::Top_dimensional_cells_iterator it = range.begin() ; it != range.end() ; ++it )
+    {
+        BOOST_CHECK( data[i] == ba.get_cell_data( *it ) );
+        BOOST_CHECK( ba.get_dimension_of_a_cell( *it ) == 2 );
+        ++i;
+    }
+}
+
diff --git a/src/Bitmap_cubical_complex/test/CMakeLists.txt b/src/Bitmap_cubical_complex/test/CMakeLists.txt
new file mode 100644
index 00000000..97c374e6
--- /dev/null
+++ b/src/Bitmap_cubical_complex/test/CMakeLists.txt
@@ -0,0 +1,25 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBitmapCCUT)
+
+if (GCOVR_PATH)
+  # for gcovr to make coverage reports - Corbera Jenkins plugin
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage")	
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage")
+endif()
+if (GPROF_PATH)
+  # for gprof to make coverage reports - Jenkins
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg")
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg")
+endif()
+
+add_executable ( BitmapCCUT Bitmap_test.cpp )
+target_link_libraries(BitmapCCUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(NAME BitmapCCUT
+  COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BitmapCCUT 
+    # XML format for Jenkins xUnit plugin 
+    --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BitmapCCUT.xml --log_level=test_suite --report_level=no)
+