regularData2D.h

Go to the documentation of this file.

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#ifndef BALL_DATATYPE_REGULARDATA2D_H
#define BALL_DATATYPE_REGULARDATA2D_H

#ifndef BALL_MATHS_VECTOR2_H
# include <BALL/MATHS/vector2.h>
#endif

#ifndef BALL_SYSTEM_FILE_H
# include <BALL/SYSTEM/file.h>
#endif

#ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
# include <BALL/SYSTEM/binaryFileAdaptor.h>
#endif

#include <iostream>
#include <fstream>
#include <iterator>
#include <algorithm>

namespace BALL 
{
  template <typename ValueType>
  class TRegularData2D 
  {
    public:

    BALL_CREATE(TRegularData2D<ValueType>)

    

    class IndexType
    {
      public:
      inline IndexType() : x(0), y(0) {}
      inline IndexType(Position p) : x(p), y(p) {}
      inline IndexType(Position p, Position q) : x(p), y(q) {}

      Position x;
      Position y;

    };

    typedef std::vector<ValueType> VectorType;
    typedef TVector2<float> CoordinateType;
    typedef typename std::vector<ValueType>::iterator Iterator;
    typedef typename std::vector<ValueType>::const_iterator ConstIterator;

    //  STL compatibility types
    //
    typedef ValueType value_type;
    typedef typename std::vector<ValueType>::iterator iterator;
    typedef typename std::vector<ValueType>::const_iterator const_iterator;
    typedef typename std::vector<ValueType>::reference reference;
    typedef typename std::vector<ValueType>::const_reference const_reference;
    typedef typename std::vector<ValueType>::pointer pointer;
    typedef typename std::vector<ValueType>::difference_type difference_type;
    typedef typename std::vector<ValueType>::size_type size_type;


    TRegularData2D();

    TRegularData2D(const TRegularData2D<ValueType>& data);  

    TRegularData2D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);

    /* Constructor.
     * This constructor takes the size of the grid (as the number of grid points per dimension)
     * as input and (optionally) the origin and the dimension (in grid <em>coordinates</em>).
     *  @throw  Exception::OutOfMemory if the memory for the grid could not be allocated
     */
    TRegularData2D(const IndexType& size, 
                   const CoordinateType& origin = CoordinateType(0.0), 
                   const CoordinateType& dimension = CoordinateType(1.0));

    virtual ~TRegularData2D();

    virtual void clear();


    TRegularData2D& operator = (const TRegularData2D<ValueType>& data);
      

    
    bool operator == (const TRegularData2D<ValueType>& data) const;

    BALL_INLINE bool operator != (const TRegularData2D<ValueType>& data) const { return !this->operator == (data); }

    BALL_INLINE bool empty() const { return data_.empty(); }

    bool isInside(const CoordinateType& x) const;

    BALL_INLINE ConstIterator begin() const { return data_.begin(); }
    BALL_INLINE ConstIterator end() const { return data_.end(); }
    BALL_INLINE Iterator begin() { return data_.begin(); }
    BALL_INLINE Iterator end() { return data_.end(); }

    // STL compatibility
    BALL_INLINE size_type size() const { return data_.size(); }
    BALL_INLINE size_type max_size() const { return data_.max_size(); }
    BALL_INLINE void swap(TRegularData2D<ValueType>& data) { std::swap(*this, data); }
    

    const ValueType& getData(const IndexType& index) const;

    ValueType& getData(const IndexType& index);

    const ValueType& getData(Position index) const;

    ValueType& getData(Position index);

    const ValueType& operator [] (const IndexType& index) const { return data_[index.x + size_.x * index.y]; }

    ValueType& operator [] (const IndexType& index) { return data_[index.x + size_.x * index.y]; }

    const ValueType& operator [] (Position index) const { return data_[index]; }

    ValueType& operator [] (Position index) { return data_[index]; }

    ValueType operator () (const CoordinateType& x) const;

    ValueType getInterpolatedValue(const CoordinateType& x) const;

    const ValueType& getClosestValue(const CoordinateType& x) const;

    ValueType& getClosestValue(const CoordinateType& x);

    IndexType getLowerIndex(const CoordinateType& v) const;

    IndexType getClosestIndex(const CoordinateType& v) const;

    inline const IndexType& getSize() const { return size_; }
        
    inline const CoordinateType& getOrigin() const { return origin_; }

    const CoordinateType& getSpacing() const {  return spacing_; }

    void setOrigin(const CoordinateType& origin);

    const CoordinateType& getDimension() const { return dimension_; }

    void setDimension(const CoordinateType& dimension) { dimension_ = dimension; }

    void resize(const IndexType& new_size);

    void rescale(const IndexType& new_size);

    CoordinateType getCoordinates(const IndexType& index) const;

    CoordinateType getCoordinates(Position index) const;

    void getEnclosingIndices
      (const CoordinateType& r, Position& ll, Position& lr, Position& ul, Position& ur) const;
                  
    void getEnclosingValues
      (const CoordinateType& r, ValueType& ll, ValueType& lr, ValueType& ul, ValueType& ur) const;
                  
    ValueType calculateMean() const;
    
    ValueType calculateSD() const;
    
    void binaryWrite(const String& filename) const;

    void binaryRead(const String& filename);
  
    
    protected:
      
    VectorType data_;

    CoordinateType  origin_;

    CoordinateType  dimension_;

    CoordinateType  spacing_;

    IndexType size_;

    typedef struct { ValueType bt[1024]; } BlockValueType;
  };

  typedef TRegularData2D<float> RegularData2D;

  // default constructor.
  template <class ValueType>
  TRegularData2D<ValueType>::TRegularData2D()
    : data_(),
      origin_(0.0),
      dimension_(0.0),
      spacing_(1.0),
      size_(0)
  {
  }

  // copy constructor
  template <class ValueType>
  TRegularData2D<ValueType>::TRegularData2D(const TRegularData2D<ValueType>& data)
    : data_(),
      origin_(data.origin_),
      dimension_(data.dimension_),
      spacing_(data.spacing_),
      size_(data.size_)
  {
    try
    {
      data_ = data.data_;
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, data.data_.size() * sizeof(ValueType));
    }
  }

  template <class ValueType>
  TRegularData2D<ValueType>::TRegularData2D
    (const typename TRegularData2D<ValueType>::IndexType& size,
     const typename TRegularData2D<ValueType>::CoordinateType& origin,
     const typename TRegularData2D<ValueType>::CoordinateType& dimension)
    : data_(),
      origin_(origin),
      dimension_(dimension),
      spacing_(0.0, 0.0),
      size_(size)
  {
    // Compute the grid spacing
    spacing_.x = dimension_.x / (double)(size_.x - 1);
    spacing_.y = dimension_.y / (double)(size_.y - 1);

    // Compute the number of grid points
    size_type number_of_points = size_.x * size_.y;
    try
    {
      data_.resize(number_of_points);
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_points * sizeof(ValueType));
    }
  }

  template <class ValueType>
  TRegularData2D<ValueType>::TRegularData2D
    (const typename TRegularData2D<ValueType>::CoordinateType& origin,
     const typename TRegularData2D<ValueType>::CoordinateType& dimension,
     const typename TRegularData2D<ValueType>::CoordinateType& spacing)
    : data_(),
      origin_(origin),
      dimension_(dimension),
      spacing_(spacing),
      size_(0)
  {
    // Compute the grid size
    size_.x = (Size)(dimension_.x / spacing_.x + 0.5) + 1;
    size_.y = (Size)(dimension_.y / spacing_.y + 0.5) + 1;
    
    // Compute the number of grid points
    size_type size = size_.x * size_.y;
    try
    {
      data_ .resize(size);
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
    }
    
    // Adjust the spacing -- dimension has precedence.
    spacing_.x = dimension_.x / (double)(size_.x - 1);
    spacing_.y = dimension_.y / (double)(size_.y - 1);
  }

  template <class ValueType>
  TRegularData2D<ValueType>::~TRegularData2D()
  {
  }

  // assignment operator
  template <typename ValueType>
  BALL_INLINE
  TRegularData2D<ValueType>& TRegularData2D<ValueType>::operator = 
    (const TRegularData2D<ValueType>& rhs)
  {
    // Avoid self assignment
    if (&rhs != this)
    {
      // Copy the coordinate-related attributes and
      // the size.
      origin_ = rhs.origin_;
      dimension_ = rhs.dimension_;
      spacing_ = rhs.spacing_;
      size_ = rhs.size_;

      // Copy the data itself and rethrow allocation exceptions.
      try
      {
        data_ = rhs.data_;
      }
      catch (std::bad_alloc&)
      {
        data_.resize(0);
        throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.data_.size() * sizeof(ValueType));
      }
    }

    return *this;
  }

  template <typename ValueType>
  void TRegularData2D<ValueType>::rescale(const typename TRegularData2D<ValueType>::IndexType& size)
  {
    // If the old size equals the new size, we're done.
    if ((size.x == size_.x) && (size_.y == size.y))
    {
      return;
    }
  
    // If the new grid is empty, this whole thing is quite easy.
    if ((size.x == 0) || (size.y == 0))
    {
      data_.resize(0);
      dimension_.set(0.0);
      return;
    }

    // Compute the new array size.
    size_type new_size = (size_type)(size.x * size.y);

    // Catch any bad_allocs thrown by vector::resize
    try
    {
      // Create a new temporary array.
      TRegularData2D<ValueType> old_data(*this);

      // Resize the data to its new size.
      data_.resize(new_size);
      spacing_.x = dimension_.x / (double)(size.x - 1);
      spacing_.y = dimension_.y / (double)(size.y - 1);
      
      // Walk over the new grid and copy the (interpolated) old stuff back.
      CoordinateType v;
      for (size_type i = 0; i < new_size; i++)
      { 
        Position x = i % size.x;
        Position y = i / size.x;
        v.x = origin_.x + x * spacing_.x;
        v.y = origin_.y + y * spacing_.y;
        data_[i] = old_data(v);
      }
      
      // Correct the grid dimension. Origin and spacing remain constant.
      size_ = size;
    }
    catch (std::bad_alloc&)
    {
      throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
    }
  }

  template <typename ValueType>
  void TRegularData2D<ValueType>::resize(const typename TRegularData2D<ValueType>::IndexType& size)
  {
    // If the old size equals the new size, we're done.
    if (size.x == size_.x && size_.y == size.y)
    {
      return;
    }
  
    // If the new grid is empty, this whole thing is quite easy.
    if ((size.x == 0) || (size.y == 0))
    {
      data_.resize(0);
      dimension_.set(0.0, 0.0);
      return;
    }

    // Compute the new array size.
    size_type new_size = (size_type)(size.x * size.y);

    // Catch any bad_allocs thrown by vector::resize
    try
    {
      // Create a new temporary array.
      std::vector<ValueType> old_data(data_);

      // Resize the data to its new size.
      data_.resize(new_size);
      
      // walk over the new grid and copy the old stuff back.
      static ValueType default_value = (ValueType)0;
      for (size_type i = 0; i < new_size; i++)
      { 
        size_type x = i % size.x;
        size_type y = i / size.x;
        if (x >= size_.x || y >= size_.y)
        {
          data_[i] = default_value;
        }
        else
        {
          data_[i] = old_data[x + y * size_.x];
        }
      }
      
      // Correct the grid dimension. Origin and spacing remain constant.
      dimension_.x *= (double)size.x / (double)size_.x;
      dimension_.y *= (double)size.y / (double)size_.y;
      size_ = size;
    }
    catch (std::bad_alloc&)
    {
      throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
    }
  }

  template <class ValueType>
  void TRegularData2D<ValueType>::setOrigin(const typename TRegularData2D<ValueType>::CoordinateType& origin)
  {
    origin_ = origin;
  }

  template <class ValueType> 
  BALL_INLINE
  bool TRegularData2D<ValueType>::isInside(const typename TRegularData2D<ValueType>::CoordinateType& r) const   
  {
    return ((r.x >= origin_.x) && (r.x <= (origin_.x + dimension_.x))
            && (r.y >= origin_.y) && (r.y <= (origin_.y + dimension_.y)));
  }

  template <class ValueType>
  BALL_INLINE 
  const ValueType& TRegularData2D<ValueType>::getData
    (const typename TRegularData2D<ValueType>::IndexType& index) const
  { 
    size_type pos = index.x + index.y * size_.x;
    if (pos >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }       
    return data_[pos];
  }

  template <class ValueType>
  BALL_INLINE 
  ValueType& TRegularData2D<ValueType>::getData(const typename TRegularData2D<ValueType>::IndexType& index)
  { 
    size_type pos = index.x + index.y * size_.x;
    if (pos >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }       
    return data_[pos];
  }

  template <class ValueType>
  BALL_INLINE 
  const ValueType& TRegularData2D<ValueType>::getData(Position index) const
  { 
    if (index >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }       
    return data_[index];
  }

  template <class ValueType>
  BALL_INLINE 
  ValueType& TRegularData2D<ValueType>::getData(Position index)
  { 
    if (index >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }       
    return data_[index];
  }

  template <class ValueType>
  BALL_INLINE 
  typename TRegularData2D<ValueType>::CoordinateType TRegularData2D<ValueType>::getCoordinates
    (const typename TRegularData2D<ValueType>::IndexType& index) const 
  {
    if ((index.x >= size_.x) || (index.y >= size_.y))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    CoordinateType r(origin_.x + index.x * spacing_.x,
                     origin_.y + index.y * spacing_.y);

    return r;
  }

  template <class ValueType> 
  BALL_INLINE 
  typename TRegularData2D<ValueType>::CoordinateType
  TRegularData2D<ValueType>::getCoordinates(Position position) const 
  {
    if (position >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    } 
    
    Position x = (Position)(position %  size_.x);
    Position y = (Position)(position / size_.x);

    return CoordinateType(origin_.x + (double)x * spacing_.x,
                          origin_.y + (double)y * spacing_.y);
  }

  template <typename ValueType>
  BALL_INLINE
  void TRegularData2D<ValueType>::getEnclosingIndices
    (const typename TRegularData2D<ValueType>::CoordinateType& r,
     Position& ll, Position& lr, Position& ul, Position& ur) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   

    // Calculate the grid indices of the lower left front corner
    // of the enclosing rectangle
    IndexType position;
    position.x = (Position)((r.x - origin_.x) / spacing_.x);
    position.y = (Position)((r.y - origin_.y) / spacing_.y);
    
    // Calculate the (linear) indices of the four rectangle corners
    ll = position.x + size_.x * position.y;
    lr = ll + 1;
    ul = ll + size_.x;
    ur = ul + 1;
  }

  template <typename ValueType>
  BALL_INLINE
  void TRegularData2D<ValueType>::getEnclosingValues
    (const typename TRegularData2D<ValueType>::CoordinateType& r,
     ValueType& ll, ValueType& lr, ValueType& ul, ValueType& ur) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    // compute the four grid indices forming the enclosing rectangle
    Position ll_id, lr_id, ul_id, ur_id;
    getEnclosingIndices(r, ll_id, lr_id, ul_id, ur_id);
        
    // Retrieve the grid values
    ll = data_[ll_id];
    lr = data_[lr_id];
    ul = data_[ul_id];
    ur = data_[ur_id];
  }

  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData2D<ValueType>::getInterpolatedValue
    (const typename TRegularData2D<ValueType>::CoordinateType& r) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    return this->operator () (r);
  }

  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData2D<ValueType>::operator ()
    (const typename TRegularData2D<ValueType>::CoordinateType& r) const
  {
    CoordinateType h(r - origin_);
    Position x = (Position)(h.x / spacing_.x);
    Position y = (Position)(h.y / spacing_.y);

    // correct for numerical inaccuracies
    if (x >= (size_.x - 1))
    {
      x = size_.x - 2;
    }
    if (y >= (size_.y - 1))
    {
      y = size_.y - 2;
    }

    Size l = x + size_.x * y;
    CoordinateType r_0(origin_.x + (double)x * spacing_.x,
                       origin_.y + (double)y * spacing_.y);

    double dx = 1.0 - ((r.x - r_0.x) / spacing_.x);
    double dy = 1.0 - ((r.y - r_0.y) / spacing_.y);

    return  data_[l] * dx * dy
          + data_[l + 1] * (1.0 - dx) * dy
          + data_[l + size_.x] * dx * (1.0 - dy)
          + data_[l + size_.x + 1] * (1.0 - dx) * (1.0 - dy);
  }

  template <typename ValueType>
  BALL_INLINE
  typename TRegularData2D<ValueType>::IndexType TRegularData2D<ValueType>::getLowerIndex
    (const typename TRegularData2D<ValueType>::CoordinateType& r) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    static IndexType position;
    position.x = (Position)((r.x - origin_.x) / spacing_.x);
    position.y = (Position)((r.y - origin_.y) / spacing_.y);
    
    return position;
  }

  template <typename ValueType>
  BALL_INLINE
  typename TRegularData2D<ValueType>::IndexType TRegularData2D<ValueType>::getClosestIndex
    (const typename TRegularData2D<ValueType>::CoordinateType& r) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    static IndexType position;
    position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
    position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
    
    return position;
  }

  template <typename ValueType>
  BALL_INLINE
  const ValueType& TRegularData2D<ValueType>::getClosestValue
    (const typename TRegularData2D<ValueType>::CoordinateType& r) const
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    static IndexType position;
    position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
    position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);

    return operator [] (position);
  }

  template <typename ValueType>
  BALL_INLINE
  ValueType& TRegularData2D<ValueType>::getClosestValue
    (const typename TRegularData2D<ValueType>::CoordinateType& r)
  {
    if (!isInside(r))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }   
    
    static IndexType position;
    position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
    position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);

    return operator [] (position);
  }

  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData2D<ValueType>::calculateMean() const
  {
    Position data_points  = (size_.x * size_.y);
    ValueType mean = 0;
    for (Position i = 0; i < data_points; i++)
    {
      mean += data_[i];
    }
    mean /= data_points;
    return mean;
  }
  
  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData2D<ValueType>::calculateSD() const
  {
    Position data_points  = (size_.x * size_.y);
    ValueType stddev = 0;
    ValueType mean = this->calculateMean();
    for (Position i = 0; i < data_points; i++)
    {
      stddev += (pow(data_(i)-mean,2));
    }
    stddev /= (data_points-1);
    stddev = sqrt(stddev);
    return stddev;
  }
    
  template <typename ValueType>
  void TRegularData2D<ValueType>::clear()
  {
    data_.resize(0);
    
    origin_.set(0.0);
    dimension_.set(0.0, 0.0);
    size_.x = 0;
    size_.y = 0;
    spacing_.set(1.0, 1.0);
  }

  template <typename ValueType> 
  bool TRegularData2D<ValueType>::operator == (const TRegularData2D<ValueType>& data) const
  {
    return ((origin_ == data.origin_)
            && (dimension_ == data.dimension_)
            && (size_.x == data.size_.x)
            && (size_.y == data.size_.y)
            && (data_ == data.data_));
  } 


  template <typename ValueType>
  std::ostream& operator << (std::ostream& os, const TRegularData2D<ValueType>& data)
  {
    // Write the grid origin, dimension, and number of grid points
    os << data.getOrigin().x << " " << data.getOrigin().y
      << std::endl
      << data.getOrigin().x + data.getDimension().x << " "
      << data.getOrigin().y + data.getDimension().y
      << std::endl
      << data.getSize().x - 1 << " " << data.getSize().y - 1
      << std::endl;

    // Write the array contents.
    std::copy(data.begin(), data.end(), std::ostream_iterator<ValueType>(os, "\n"));
    return os;
  }

  template <typename ValueType>
  std::istream& operator >> (std::istream& is, TRegularData2D<ValueType>& grid)
  {
    typename TRegularData2D<ValueType>::CoordinateType origin;
    typename TRegularData2D<ValueType>::CoordinateType dimension;
    typename TRegularData2D<ValueType>::IndexType size;

    is >> origin.x >> origin.y;
    is >> dimension.x >> dimension.y;
    is >> size.x >> size.y;
    
    dimension -= origin;
    size.x++;
    size.y++;

    grid.resize(size);
    grid.setOrigin(origin);
    grid.setDimension(dimension);

    std::copy(std::istream_iterator<ValueType>(is), 
              std::istream_iterator<ValueType>(), 
              grid.begin());
    //    std::copy_n(std::istream_iterator<ValueType>(is), grid.size(), grid.begin());
    
    return is;
  }

  template <typename ValueType>
  void TRegularData2D<ValueType>::binaryWrite(const String& filename) const
  {
    File outfile(filename.c_str(), std::ios::out|std::ios::binary);
    if (!outfile.isValid()) 
      throw Exception::FileNotFound(__FILE__, __LINE__, filename);

    // write all information we need to recreate the grid
    BinaryFileAdaptor<BlockValueType> adapt_block;
    BinaryFileAdaptor<ValueType>      adapt_single;
    BinaryFileAdaptor<float>          adapt_float;
    
    BinaryFileAdaptor<Size>           adapt_size;

    adapt_size.setData(data_.size());
    outfile << adapt_size;
    
    // NOTE: we do not use the binary file adaptor to write out the Vector2-variables here,
    //       the reason is a bit stupid: the data layout of Vector2 has a three-byte "hole"
    //       that is automatically padded by the compiler to ensure the correct alignment.
    //       valgrind is very unhappy when we write these uninitialized values out, even
    //       though they are entirely harmless. To prevent false-positives in the error checking,
    //       we write out the members of the vectors instead. This even saves some bytes in the
    //       output (not that it would matter...)
    adapt_float.setData(origin_.x);
    outfile << adapt_float;
    adapt_float.setData(origin_.y);
    outfile << adapt_float;

    adapt_float.setData(dimension_.x);
    outfile << adapt_float;
    adapt_float.setData(dimension_.y);
    outfile << adapt_float;

    adapt_float.setData(spacing_.x);
    outfile << adapt_float;
    adapt_float.setData(spacing_.y);
    outfile << adapt_float;

    BinaryFileAdaptor<IndexType> adapt_index;
    adapt_index.setData(size_);
    outfile << adapt_index;
  
    // we slide a window of BLOCK_SIZE over our data.
    const int BLOCK_SIZE = 1024;
    Index window_pos = 0;
    
    while (((int)data_.size() - (BLOCK_SIZE + window_pos)) >= 0)
    {
      adapt_block.setData(* (BlockValueType*)&(data_[window_pos]));
      outfile << adapt_block;
      window_pos += BLOCK_SIZE;
    }

    // Now we have to write the remaining data one by one.
    for (Size i = window_pos; i < data_.size(); i++)
    {
      adapt_single.setData(data_[i]);
      outfile << adapt_single;
    }

    // That's it. 
    outfile.close();
  }

  template <typename ValueType>
  void TRegularData2D<ValueType>::binaryRead(const String& filename)
  {
    File infile(filename, std::ios::in|std::ios::binary);
    if (!infile.isValid()) 
    {
      throw Exception::FileNotFound(__FILE__, __LINE__, filename);
    }
    
    BinaryFileAdaptor< BlockValueType > adapt_block;
    BinaryFileAdaptor< ValueType >      adapt_single;
    
    // read all information we need to recreate the grid
    BinaryFileAdaptor<Size>           adapt_size;
    BinaryFileAdaptor<float>          adapt_float;

    infile >> adapt_size;
    Size new_size = adapt_size.getData();
  
    infile >> adapt_float;
    origin_.x = adapt_float.getData();
    infile >> adapt_float;
    origin_.y = adapt_float.getData();

    infile >> adapt_float;
    dimension_.x = adapt_float.getData();
    infile >> adapt_float;
    dimension_.y = adapt_float.getData();

    infile >> adapt_float;
    spacing_.x = adapt_float.getData();
    infile >> adapt_float;
    spacing_.y = adapt_float.getData();

    BinaryFileAdaptor<IndexType> adapt_index;
    infile >> adapt_index;
    size_ =  adapt_index.getData();
  
    data_.resize(new_size);

    // we slide a window of size 1024 over our data
    Index window_pos = 0;
    
    while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
    {
      infile >> adapt_block;
      *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
      /*
      for (Size i=0; i<1024; i++)
      {
        data_[i+window_pos] = adapt_block.getData().bt[i];
      }
      */
      window_pos+=1024;
    }

    // now we have to read the remaining data one by one
    for (Size i=window_pos; i<data_.size(); i++)
    {
      infile >> adapt_single;
      data_[i] = adapt_single.getData();
    }

    // that's it. I hope...
    infile.close();
  } 
 } // namespace BALL

#endif // BALL_DATATYPE_TREGULARDATA2D_H