Adding documentation of the class.

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

Former-commit-id: 6eec3b7216880475ed667660275c2314bf0feb0c

3 files changed, 427 insertions, 140 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
index b324d272..63edcadd 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -23,7 +23,8 @@
 
 #pragma once
 #include <limits>
-#include "Bitmap_cubical_complex_base.h"
+#include "Bitmap_cubical_complex_base.h"

+#include "Bitmap_cubical_complex_periodic_boundary_conditions_base.h"
 
 
 
@@ -38,15 +39,15 @@ const bool globalDbg = false;
 

 template <typename T> class is_before_in_filtration;
 
-template <typename T = double>
-class Bitmap_cubical_complex : public Bitmap_cubical_complex_base<T>
+template <typename T>
+class Bitmap_cubical_complex : public T
 {
 public:
 //*********************************************//
 //Typedefs and typenames
 //*********************************************//
     typedef size_t Simplex_key;
-    typedef T Filtration_value;

+    typedef typename T::filtration_type Filtration_value;

     typedef Simplex_key Simplex_handle;

 
 
@@ -63,7 +64,7 @@ public:
     * Constructor form a Perseus-style file.
     **/
     Bitmap_cubical_complex( const char* perseus_style_file ):
-    Bitmap_cubical_complex_base<T>(perseus_style_file),key_associated_to_simplex(this->total_number_of_cells+1)
+    T(perseus_style_file),key_associated_to_simplex(this->total_number_of_cells+1)
     {
         if ( globalDbg ){cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n";}
         for ( size_t i = 0 ; i != this->total_number_of_cells ; ++i )
@@ -82,8 +83,28 @@ public:
     * 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<T>& top_dimensional_cells ):
-    Bitmap_cubical_complex_base<T>(dimensions,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 )
@@ -138,9 +159,9 @@ public:
     /**
     * Return the filtration of a cell pointed by the Simplex_handle.
     **/
-    T filtration(Simplex_handle sh)
+    typename T::filtration_type filtration(Simplex_handle sh)
     {
-        if ( globalDbg ){cerr << "T filtration(const Simplex_handle& sh)\n";}
+        if ( globalDbg ){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();
@@ -493,8 +514,8 @@ public:
     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.
-      T fil1 = CC_->data[sh1];
-      T fil2 = CC_->data[sh2];
+      typename T::filtration_type fil1 = CC_->data[sh1];
+      typename T::filtration_type fil2 = CC_->data[sh2];
       if ( fil1 != fil2 )
       {
 	   return fil1 < fil2;
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index 54d60325..600f250d 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -63,7 +63,12 @@ namespace Cubical_complex
 template <typename T>
 class Bitmap_cubical_complex_base
 {
-public:
+public:

+    typedef T filtration_type;

+    /**

+    *Default constructor

+    **/

+    Bitmap_cubical_complex_base(){}
     /**
     * 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
@@ -93,7 +98,7 @@ public:
      * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to
      * the positions of (co)boundary element of the input cell.
     */
-    inline std::vector< size_t > get_boundary_of_a_cell( size_t cell )const;
+    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
@@ -286,7 +291,12 @@ public:
 //****************************************************************************************************************//
 //****************************************************************************************************************//
 //****************************************************************************************************************//
-
+

+

+inline size_t number_cells()const

+{

+    return this->total_number_of_cells;

+}
 
 //****************************************************************************************************************//
 //****************************************************************************************************************//
@@ -305,14 +315,10 @@ protected:
         {
             this->sizes.push_back(sizes[i]);
             this->multipliers.push_back(multiplier);
-            //multiplier *= 2*(sizes[i]+1)+1;
             multiplier *= 2*sizes[i]+1;
         }
-        //std::reverse( this->sizes.begin() , this->sizes.end() );
-        std::vector<T> data(multiplier);
-        std::fill( data.begin() , data.end() , std::numeric_limits<int>::max() );
+        this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max());
         this->total_number_of_cells = multiplier;
-        this->data = data;
     }
 
     size_t compute_position_in_bitmap( const std::vector< unsigned >& counter )
@@ -336,11 +342,12 @@ protected:
         }
         std::reverse( counter.begin() , counter.end() );
         return counter;
-    }
-
-    std::vector< size_t >
-    generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions
-    ( const std::vector< bool >& directions_for_periodic_b_cond );
+    }

+    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 );
 };
 
 
@@ -364,11 +371,10 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
 {
     this->set_up_containers( sizes );
 }
-
-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 )
-{
+

+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;
@@ -397,13 +403,19 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base
         (*it) = top_dimensional_cells[index];
         ++index;
     }
-    this->impose_lower_star_filtration();
-}
-
+    this->impose_lower_star_filtration();

+}

 
 template <typename T>
-Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus_style_file )
+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;
     ifstream inFiltration, inIds;
     inFiltration.open( perseus_style_file );
@@ -418,7 +430,7 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus
     {
         unsigned size_in_this_dimension;
         inFiltration >> size_in_this_dimension;
-        size_in_this_dimension = abs(size_in_this_dimension);
+        size_in_this_dimension = size_in_this_dimension;
         sizes.push_back( size_in_this_dimension );
         if (dbg){cerr << "size_in_this_dimension : " << size_in_this_dimension << endl;}
     }
@@ -427,7 +439,6 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus
     Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
     it = this->top_dimensional_cells_begin();
 
-    //TODO -- over here we also need to read id's of cell and put them to bitmapElement structure!
     while ( !inFiltration.eof() )
     {
         double filtrationLevel;
@@ -444,7 +455,37 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base( const char* perseus
         ++it;
     }
     inFiltration.close();
-    this->impose_lower_star_filtration();
+    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 bundary 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 bundary 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 bundary 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 );
 }
 
 
@@ -468,60 +509,39 @@ std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell( si
         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
 {
-    //first of all, we need to take the list of coordinates in which the cell has nonzero length.
-    //We do it by using modified version to compute dimension of a cell:
-    std::vector< unsigned > dimensions_in_which_cell_has_zero_length;

-

-    //Speed traded of for memory. Check if it is better in practice.

-    dimensions_in_which_cell_has_zero_length.reserve( this->dimension()*2 );

-
-    unsigned dimension = 0;
+    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 )
         {
-            dimensions_in_which_cell_has_zero_length.push_back(i-1);
-            dimension++;
+            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];
     }
-
-    std::vector<unsigned> counter = this->compute_counter_for_given_cell( cell );
-    std::vector< size_t > coboundary_elements;

+    return coboundary_elements;
+}

 

-    //Speed traded of for memory. Check if it is better in practice.

-    coboundary_elements.reserve ( dimensions_in_which_cell_has_zero_length.size()*2 );

 
-    if ( dimensions_in_which_cell_has_zero_length.size() == 0 )return coboundary_elements;
-    for ( size_t i = 0 ; i != dimensions_in_which_cell_has_zero_length.size() ; ++i )
-    {
-        if ( (cell > this->multipliers[dimensions_in_which_cell_has_zero_length[i]])
-              && (counter[dimensions_in_which_cell_has_zero_length[i]] != 0) )
-        {
-            coboundary_elements.push_back( cell - this->multipliers[dimensions_in_which_cell_has_zero_length[i]] );
-        }
-        if (
-           (cell + this->multipliers[dimensions_in_which_cell_has_zero_length[i]] < this->data.size()) &&
-           (counter[dimensions_in_which_cell_has_zero_length[i]]
-           !=
-           2*this->sizes[dimensions_in_which_cell_has_zero_length[i]])
-           )
-        {
-            coboundary_elements.push_back( cell + this->multipliers[dimensions_in_which_cell_has_zero_length[i]] );
-        }
-    }
-    return coboundary_elements;
-}
 
 
 
@@ -555,7 +575,7 @@ unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell( size_t cell )c
 }
 
 template <typename T>
-T& Bitmap_cubical_complex_base<T>::get_cell_data( size_t cell )
+inline T& Bitmap_cubical_complex_base<T>::get_cell_data( size_t cell )
 {
     return this->data[cell];
 }
@@ -602,10 +622,21 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration()
         {
             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 )

+                {

+                    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] ] << endl;

+                    getchar();

+

+                }
                 if ( this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ] )
                 {
-                    this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];
+                    this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ];

+                    if ( dbg )

+                    {

+                        cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " << this->data[ indices_to_consider[i] ] << endl;

+                        getchar();

+                    }
                 }
                 if ( is_this_cell_considered[ bd[boundaryIt] ] == false )
                 {
@@ -640,69 +671,6 @@ bool compareFirstElementsOfTuples( const std::pair< std::pair< T , size_t > , ch
 
 
 
-template <typename T>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::
-generate_vector_of_shifts_for_bitmaps_with_periodic_boundary_conditions
-( const std::vector< bool >& directions_for_periodic_b_cond )
-{
-    bool dbg = false;
-    if ( this->sizes.size() != directions_for_periodic_b_cond.size() )
-    throw "directions_for_periodic_b_cond vector size is different from the size of the bitmap. Program terminate \n";
-
-    std::vector<unsigned> sizes;

-    sizes.reserve( this->sizes.size() );
-    for ( size_t i = 0 ; i != this->sizes.size() ; ++i )

-    {

-        sizes.push_back(2*this->sizes[i]);

-    }
-
-    counter c( sizes );
-
-    std::vector< size_t > result;

-    result.reserve( this->data.size() );
-
-    for ( size_t i = 0 ; i != this->data.size() ; ++i )
-    {
-        size_t position;
-        if ( !c.isFinal() )
-        {
-            position = i;
-        }
-        else
-        {
-            std::vector< bool > finals = c.directions_of_finals();
-            bool jump_in_position = false;
-            for ( size_t dir = 0 ; dir != finals.size() ; ++dir )
-            {
-                if ( finals[dir] == false )continue;
-                if ( directions_for_periodic_b_cond[dir] )
-                {
-                    jump_in_position = true;
-                }
-            }
-            if ( jump_in_position == true )
-            {
-                //in this case this guy is final, so we need to find 'the opposite one'
-                position = compute_position_in_bitmap(  c.find_opposite( directions_for_periodic_b_cond ) );
-            }
-            else
-            {
-                position = i;
-            }
-        }
-        result.push_back( position );
-        if ( dbg )
-        {
-            cerr << " position : " << position << endl;
-            cerr << c << endl;
-            getchar();
-        }
-
-        c.increment();
-    }
-
-    return result;
-}
 
 }
 
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..28911ea8
--- /dev/null
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
@@ -0,0 +1,298 @@
+#pragma once

+#include <cmath>

+#include "Bitmap_cubical_complex_base.h"

+

+using namespace std;
+
+namespace Gudhi
+{
+
+namespace Cubical_complex
+{

+

+//in this class, we are storing all the elements which are in normal bitmap (i.e. the bitmap without the periodic boundary conditions). But, we set up the iterators and the procedures

+//to compute boundary and coboundary in the way that it is all right. We assume here that all the cells that are on the left / bottom and so on remains, while all the cells on the

+//right / top are not in the Bitmap_cubical_complex_periodic_boundary_conditions_base

+

+template <typename T>

+class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_cubical_complex_base<T>

+{

+public:

+    //constructors that take an extra parameter:

+    Bitmap_cubical_complex_periodic_boundary_conditions_base(){};

+    Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed );

+    Bitmap_cubical_complex_periodic_boundary_conditions_base( const char* perseusStyleFile );

+    Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed );

+

+    //overwritten methods co compute boundary and coboundary

+    virtual std::vector< size_t > get_boundary_of_a_cell( size_t cell )const;

+    std::vector< size_t > get_coboundary_of_a_cell( size_t cell )const;

+    //inline unsigned get_dimension_of_a_cell( size_t cell )const;

+

+protected:

+    std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed;

+    void set_up_containers( const std::vector<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( std::vector<unsigned> sizes );

+    Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells );

+    void construct_complex_based_on_top_dimensional_cells( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , 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( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed )

+{

+    this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;

+    this->set_up_containers( dimensions );

+

+    size_t i = 0;

+    for ( typename Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Top_dimensional_cells_iterator it = this->top_dimensional_cells_begin() ; it != this->top_dimensional_cells_end() ; ++it )

+    {

+        *it = topDimensionalCells[i];

+        ++i;

+    }

+    this->impose_lower_star_filtration();

+}

+

+template <typename T>

+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( std::vector<unsigned> sizes , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed )

+{

+    this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed;

+    this->set_up_containers( sizes );

+}

+

+template <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;

+

+    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_begin();
+
+    while ( !inFiltration.eof() )
+    {
+        double filtrationLevel;
+        inFiltration >> filtrationLevel;

+        if ( inFiltration.eof() )break;

+
+        if ( dbg )
+        {
+            cerr << "Cell of an index : "
+                 << it.compute_index_in_bitmap()
+                 << " and dimension: "
+                 << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
+                 << " get the value : " << filtrationLevel << endl;
+        }
+        *it = filtrationLevel;
+        ++it;
+    }
+    inFiltration.close();
+    this->impose_lower_star_filtration();

+

+    /*

+     char* filename = (char*)perseus_style_file;

+     //char* filename = "combustionWithPeriodicBoundaryConditions/v0/tV0_000000.float";

+     ifstream file( filename , ios::binary | ios::ate );

+     unsigned realSizeOfFile = file.tellg();

+     file.close();

+     realSizeOfFile = realSizeOfFile/sizeof(T);

+

+     unsigned w, h, d;

+

+     w = h = d = ceil(pow( realSizeOfFile , (double)(1/(double)3) ));

+

+     T* slice = new T[w*h*d];

+     if (slice == NULL)

+     {

+        cerr << "Allocation error, cannot allocate " << w*h*d*sizeof(T) << " bytes to store the data from the file. The program will now terminate \n";

+        exit(EXIT_FAILURE);

+     }

+

+     FILE* fp;

+     if ((fp=fopen( filename, "rb" )) == NULL )

+     {

+         cerr << "Cannot open the file: " << filename << ". The program will now terminate \n";

+         exit(1);

+     }

+     fread( slice,4,w*h*d,fp );

+     fclose(fp);

+     std::vector<T> data(slice,slice+w*h*d);

+     delete[] slice;

+     std::vector< unsigned > sizes;

+     sizes.push_back(w);

+     sizes.push_back(w);

+     sizes.push_back(w);

+

+     std::vector< bool > directions;

+     directions.push_back( true );

+     directions.push_back( true );

+     directions.push_back( true );

+     Bitmap_cubical_complex_periodic_boundary_conditions_base<T> b( sizes, data, directions );

+     *this = b;

+     */

+}

+

+template <typename T>

+Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( 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( std::vector<unsigned> dimensions , 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( std::vector<unsigned> dimensions , std::vector<T> topDimensionalCells , std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed )

+{

+    this->construct_complex_based_on_top_dimensional_cells( dimensions , topDimensionalCells , directions_in_which_periodic_b_cond_are_to_be_imposed );

+}

+

+//***********************Methods************************//

+

+template <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 ){cerr << "Computations of boundary of a cell : " << cell << endl;}

+

+    std::vector< size_t > boundary_elements;
+    size_t cell1 = cell;
+    for ( size_t i = this->multipliers.size() ; i != 0 ; --i )
+    {
+        unsigned position = cell1/this->multipliers[i-1];

+        //this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.

+
+        if ( position%2 == 1 )
+        {

+            //if there are no periodic boundary conditions in this direction, we do not have to do anything.

+            if ( !directions_in_which_periodic_b_cond_are_to_be_imposed[i-1] )

+            {

+                //cerr << "A\n";

+                boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

+                boundary_elements.push_back( cell + this->multipliers[ i-1 ] );

+                if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";}

+            }

+            else

+            {

+                //in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.

+                if ( position != 2*this->sizes[ i-1 ]-1 )

+                {

+                    //cerr << "B\n";

+                    boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

+                    boundary_elements.push_back( cell + this->multipliers[ i-1 ] );

+                    if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell + this->multipliers[ i-1 ] << " ";}

+                }

+                else

+                {

+                    //cerr << "C\n";

+                    boundary_elements.push_back( cell - this->multipliers[ i-1 ] );

+                    boundary_elements.push_back( cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] );

+                    if (dbg){cerr << cell - this->multipliers[ i-1 ] << " " << cell - (2*this->sizes[ i-1 ]-1)*this->multipliers[ i-1 ] << " ";}

+                }

+            }
+        }
+        cell1 = cell1%this->multipliers[i-1];
+    }
+    return boundary_elements;

+}

+

+template <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;
+}

+

+

+

+}//Cubical_complex

+}//namespace Gudhi