BALL  1.4.1
regularData2D.h
Go to the documentation of this file.
00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 
00005 #ifndef BALL_DATATYPE_REGULARDATA2D_H
00006 #define BALL_DATATYPE_REGULARDATA2D_H
00007 
00008 #ifndef BALL_MATHS_VECTOR2_H
00009 # include <BALL/MATHS/vector2.h>
00010 #endif
00011 
00012 #ifndef BALL_SYSTEM_FILE_H
00013 # include <BALL/SYSTEM/file.h>
00014 #endif
00015 
00016 #ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
00017 # include <BALL/SYSTEM/binaryFileAdaptor.h>
00018 #endif
00019 
00020 #include <iostream>
00021 #include <fstream>
00022 #include <iterator>
00023 #include <algorithm>
00024 
00025 namespace BALL 
00026 {
00036   template <typename ValueType>
00037   class TRegularData2D 
00038   {
00039     public:
00040 
00041     BALL_CREATE(TRegularData2D<ValueType>)
00042 
00043     
00046 
00048     class IndexType
00049     {
00050       public:
00051       inline IndexType() : x(0), y(0) {}
00052       inline IndexType(Position p) : x(p), y(p) {}
00053       inline IndexType(Position p, Position q) : x(p), y(q) {}
00054 
00056       Position x;
00058       Position y;
00059 
00060     };
00061 
00063     typedef std::vector<ValueType> VectorType;
00065     typedef TVector2<float> CoordinateType;
00067     typedef typename std::vector<ValueType>::iterator Iterator;
00069     typedef typename std::vector<ValueType>::const_iterator ConstIterator;
00071 
00072     //  STL compatibility types
00073     //
00074     typedef ValueType value_type;
00075     typedef typename std::vector<ValueType>::iterator iterator;
00076     typedef typename std::vector<ValueType>::const_iterator const_iterator;
00077     typedef typename std::vector<ValueType>::reference reference;
00078     typedef typename std::vector<ValueType>::const_reference const_reference;
00079     typedef typename std::vector<ValueType>::pointer pointer;
00080     typedef typename std::vector<ValueType>::difference_type difference_type;
00081     typedef typename std::vector<ValueType>::size_type size_type;
00082 
00086 
00090     TRegularData2D();
00091 
00095     TRegularData2D(const TRegularData2D<ValueType>& data);  
00096 
00103     TRegularData2D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);
00104 
00105     /* Constructor.
00106      * This constructor takes the size of the grid (as the number of grid points per dimension)
00107      * as input and (optionally) the origin and the dimension (in grid <em>coordinates</em>).
00108      *  @throw  Exception::OutOfMemory if the memory for the grid could not be allocated
00109      */
00110     TRegularData2D(const IndexType& size, 
00111                    const CoordinateType& origin = CoordinateType(0.0), 
00112                    const CoordinateType& dimension = CoordinateType(1.0));
00113 
00116     virtual ~TRegularData2D();
00117 
00121     virtual void clear();
00123 
00127 
00132     TRegularData2D& operator = (const TRegularData2D<ValueType>& data);
00134       
00135 
00139     
00144     bool operator == (const TRegularData2D<ValueType>& data) const;
00145 
00147     BALL_INLINE bool operator != (const TRegularData2D<ValueType>& data) const { return !this->operator == (data); }
00148 
00150     BALL_INLINE bool empty() const { return data_.empty(); }
00151 
00153     bool isInside(const CoordinateType& x) const;
00155 
00159 
00160     BALL_INLINE ConstIterator begin() const { return data_.begin(); }
00162     BALL_INLINE ConstIterator end() const { return data_.end(); }
00164     BALL_INLINE Iterator begin() { return data_.begin(); }
00166     BALL_INLINE Iterator end() { return data_.end(); }
00168 
00172     // STL compatibility
00173     BALL_INLINE size_type size() const { return data_.size(); }
00174     BALL_INLINE size_type max_size() const { return data_.max_size(); }
00175     BALL_INLINE void swap(TRegularData2D<ValueType>& data) { std::swap(*this, data); }
00176     
00177 
00182     const ValueType& getData(const IndexType& index) const;
00183 
00188     ValueType& getData(const IndexType& index);
00189 
00194     const ValueType& getData(Position index) const;
00195 
00200     ValueType& getData(Position index);
00201 
00206     const ValueType& operator [] (const IndexType& index) const { return data_[index.x + size_.x * index.y]; }
00207 
00212     ValueType& operator [] (const IndexType& index) { return data_[index.x + size_.x * index.y]; }
00213 
00218     const ValueType& operator [] (Position index) const { return data_[index]; }
00219 
00224     ValueType& operator [] (Position index) { return data_[index]; }
00225 
00234     ValueType operator () (const CoordinateType& x) const;
00235 
00242     ValueType getInterpolatedValue(const CoordinateType& x) const;
00243 
00250     const ValueType& getClosestValue(const CoordinateType& x) const;
00251 
00258     ValueType& getClosestValue(const CoordinateType& x);
00259 
00263     IndexType getLowerIndex(const CoordinateType& v) const;
00264 
00270     IndexType getClosestIndex(const CoordinateType& v) const;
00271 
00277     inline const IndexType& getSize() const { return size_; }
00278         
00283     inline const CoordinateType& getOrigin() const { return origin_; }
00284 
00289     const CoordinateType& getSpacing() const {  return spacing_; }
00290 
00293     void setOrigin(const CoordinateType& origin);
00294 
00300     const CoordinateType& getDimension() const { return dimension_; }
00301 
00307     void setDimension(const CoordinateType& dimension) { dimension_ = dimension; }
00308 
00321     void resize(const IndexType& new_size);
00322 
00332     void rescale(const IndexType& new_size);
00333 
00338     CoordinateType getCoordinates(const IndexType& index) const;
00339 
00344     CoordinateType getCoordinates(Position index) const;
00345 
00359     void getEnclosingIndices
00360       (const CoordinateType& r, Position& ll, Position& lr, Position& ul, Position& ur) const;
00361                   
00366     void getEnclosingValues
00367       (const CoordinateType& r, ValueType& ll, ValueType& lr, ValueType& ul, ValueType& ur) const;
00368                   
00372     ValueType calculateMean() const;
00373     
00377     ValueType calculateSD() const;
00378     
00382     void binaryWrite(const String& filename) const;
00383 
00387     void binaryRead(const String& filename);
00389   
00390     
00391     protected:
00392       
00394     VectorType data_;
00395 
00397     CoordinateType  origin_;
00398 
00400     CoordinateType  dimension_;
00401 
00403     CoordinateType  spacing_;
00404 
00406     IndexType size_;
00407 
00409     typedef struct { ValueType bt[1024]; } BlockValueType;
00410   };
00411 
00414   typedef TRegularData2D<float> RegularData2D;
00415 
00416   // default constructor.
00417   template <class ValueType>
00418   TRegularData2D<ValueType>::TRegularData2D()
00419     : data_(),
00420       origin_(0.0),
00421       dimension_(0.0),
00422       spacing_(1.0),
00423       size_(0)
00424   {
00425   }
00426 
00427   // copy constructor
00428   template <class ValueType>
00429   TRegularData2D<ValueType>::TRegularData2D(const TRegularData2D<ValueType>& data)
00430     : data_(),
00431       origin_(data.origin_),
00432       dimension_(data.dimension_),
00433       spacing_(data.spacing_),
00434       size_(data.size_)
00435   {
00436     try
00437     {
00438       data_ = data.data_;
00439     }
00440     catch (std::bad_alloc&)
00441     {
00442       data_.resize(0);
00443       throw Exception::OutOfMemory(__FILE__, __LINE__, data.data_.size() * sizeof(ValueType));
00444     }
00445   }
00446 
00447   template <class ValueType>
00448   TRegularData2D<ValueType>::TRegularData2D
00449     (const typename TRegularData2D<ValueType>::IndexType& size,
00450      const typename TRegularData2D<ValueType>::CoordinateType& origin,
00451      const typename TRegularData2D<ValueType>::CoordinateType& dimension)
00452     : data_(),
00453       origin_(origin),
00454       dimension_(dimension),
00455       spacing_(0.0, 0.0),
00456       size_(size)
00457   {
00458     // Compute the grid spacing
00459     spacing_.x = dimension_.x / (double)(size_.x - 1);
00460     spacing_.y = dimension_.y / (double)(size_.y - 1);
00461 
00462     // Compute the number of grid points
00463     size_type number_of_points = size_.x * size_.y;
00464     try
00465     {
00466       data_.resize(number_of_points);
00467     }
00468     catch (std::bad_alloc&)
00469     {
00470       data_.resize(0);
00471       throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_points * sizeof(ValueType));
00472     }
00473   }
00474 
00475   template <class ValueType>
00476   TRegularData2D<ValueType>::TRegularData2D
00477     (const typename TRegularData2D<ValueType>::CoordinateType& origin,
00478      const typename TRegularData2D<ValueType>::CoordinateType& dimension,
00479      const typename TRegularData2D<ValueType>::CoordinateType& spacing)
00480     : data_(),
00481       origin_(origin),
00482       dimension_(dimension),
00483       spacing_(spacing),
00484       size_(0)
00485   {
00486     // Compute the grid size
00487     size_.x = (Size)(dimension_.x / spacing_.x + 0.5) + 1;
00488     size_.y = (Size)(dimension_.y / spacing_.y + 0.5) + 1;
00489     
00490     // Compute the number of grid points
00491     size_type size = size_.x * size_.y;
00492     try
00493     {
00494       data_ .resize(size);
00495     }
00496     catch (std::bad_alloc&)
00497     {
00498       data_.resize(0);
00499       throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
00500     }
00501     
00502     // Adjust the spacing -- dimension has precedence.
00503     spacing_.x = dimension_.x / (double)(size_.x - 1);
00504     spacing_.y = dimension_.y / (double)(size_.y - 1);
00505   }
00506 
00507   template <class ValueType>
00508   TRegularData2D<ValueType>::~TRegularData2D()
00509   {
00510   }
00511 
00512   // assignment operator
00513   template <typename ValueType>
00514   BALL_INLINE
00515   TRegularData2D<ValueType>& TRegularData2D<ValueType>::operator = 
00516     (const TRegularData2D<ValueType>& rhs)
00517   {
00518     // Avoid self assignment
00519     if (&rhs != this)
00520     {
00521       // Copy the coordinate-related attributes and
00522       // the size.
00523       origin_ = rhs.origin_;
00524       dimension_ = rhs.dimension_;
00525       spacing_ = rhs.spacing_;
00526       size_ = rhs.size_;
00527 
00528       // Copy the data itself and rethrow allocation exceptions.
00529       try
00530       {
00531         data_ = rhs.data_;
00532       }
00533       catch (std::bad_alloc&)
00534       {
00535         data_.resize(0);
00536         throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.data_.size() * sizeof(ValueType));
00537       }
00538     }
00539 
00540     return *this;
00541   }
00542 
00543   template <typename ValueType>
00544   void TRegularData2D<ValueType>::rescale(const typename TRegularData2D<ValueType>::IndexType& size)
00545   {
00546     // If the old size equals the new size, we're done.
00547     if ((size.x == size_.x) && (size_.y == size.y))
00548     {
00549       return;
00550     }
00551   
00552     // If the new grid is empty, this whole thing is quite easy.
00553     if ((size.x == 0) || (size.y == 0))
00554     {
00555       data_.resize(0);
00556       dimension_.set(0.0);
00557       return;
00558     }
00559 
00560     // Compute the new array size.
00561     size_type new_size = (size_type)(size.x * size.y);
00562 
00563     // Catch any bad_allocs thrown by vector::resize
00564     try
00565     {
00566       // Create a new temporary array.
00567       TRegularData2D<ValueType> old_data(*this);
00568 
00569       // Resize the data to its new size.
00570       data_.resize(new_size);
00571       spacing_.x = dimension_.x / (double)(size.x - 1);
00572       spacing_.y = dimension_.y / (double)(size.y - 1);
00573       
00574       // Walk over the new grid and copy the (interpolated) old stuff back.
00575       CoordinateType v;
00576       for (size_type i = 0; i < new_size; i++)
00577       { 
00578         Position x = i % size.x;
00579         Position y = i / size.x;
00580         v.x = origin_.x + x * spacing_.x;
00581         v.y = origin_.y + y * spacing_.y;
00582         data_[i] = old_data(v);
00583       }
00584       
00585       // Correct the grid dimension. Origin and spacing remain constant.
00586       size_ = size;
00587     }
00588     catch (std::bad_alloc&)
00589     {
00590       throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
00591     }
00592   }
00593 
00594   template <typename ValueType>
00595   void TRegularData2D<ValueType>::resize(const typename TRegularData2D<ValueType>::IndexType& size)
00596   {
00597     // If the old size equals the new size, we're done.
00598     if (size.x == size_.x && size_.y == size.y)
00599     {
00600       return;
00601     }
00602   
00603     // If the new grid is empty, this whole thing is quite easy.
00604     if ((size.x == 0) || (size.y == 0))
00605     {
00606       data_.resize(0);
00607       dimension_.set(0.0, 0.0);
00608       return;
00609     }
00610 
00611     // Compute the new array size.
00612     size_type new_size = (size_type)(size.x * size.y);
00613 
00614     // Catch any bad_allocs thrown by vector::resize
00615     try
00616     {
00617       // Create a new temporary array.
00618       std::vector<ValueType> old_data(data_);
00619 
00620       // Resize the data to its new size.
00621       data_.resize(new_size);
00622       
00623       // walk over the new grid and copy the old stuff back.
00624       static ValueType default_value = (ValueType)0;
00625       for (size_type i = 0; i < new_size; i++)
00626       { 
00627         size_type x = i % size.x;
00628         size_type y = i / size.x;
00629         if (x >= size_.x || y >= size_.y)
00630         {
00631           data_[i] = default_value;
00632         }
00633         else
00634         {
00635           data_[i] = old_data[x + y * size_.x];
00636         }
00637       }
00638       
00639       // Correct the grid dimension. Origin and spacing remain constant.
00640       dimension_.x *= (double)size.x / (double)size_.x;
00641       dimension_.y *= (double)size.y / (double)size_.y;
00642       size_ = size;
00643     }
00644     catch (std::bad_alloc&)
00645     {
00646       throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
00647     }
00648   }
00649 
00650   template <class ValueType>
00651   void TRegularData2D<ValueType>::setOrigin(const typename TRegularData2D<ValueType>::CoordinateType& origin)
00652   {
00653     origin_ = origin;
00654   }
00655 
00656   template <class ValueType> 
00657   BALL_INLINE
00658   bool TRegularData2D<ValueType>::isInside(const typename TRegularData2D<ValueType>::CoordinateType& r) const   
00659   {
00660     return ((r.x >= origin_.x) && (r.x <= (origin_.x + dimension_.x))
00661             && (r.y >= origin_.y) && (r.y <= (origin_.y + dimension_.y)));
00662   }
00663 
00664   template <class ValueType>
00665   BALL_INLINE 
00666   const ValueType& TRegularData2D<ValueType>::getData
00667     (const typename TRegularData2D<ValueType>::IndexType& index) const
00668   { 
00669     size_type pos = index.x + index.y * size_.x;
00670     if (pos >= data_.size())
00671     {
00672       throw Exception::OutOfGrid(__FILE__, __LINE__);
00673     }       
00674     return data_[pos];
00675   }
00676 
00677   template <class ValueType>
00678   BALL_INLINE 
00679   ValueType& TRegularData2D<ValueType>::getData(const typename TRegularData2D<ValueType>::IndexType& index)
00680   { 
00681     size_type pos = index.x + index.y * size_.x;
00682     if (pos >= data_.size())
00683     {
00684       throw Exception::OutOfGrid(__FILE__, __LINE__);
00685     }       
00686     return data_[pos];
00687   }
00688 
00689   template <class ValueType>
00690   BALL_INLINE 
00691   const ValueType& TRegularData2D<ValueType>::getData(Position index) const
00692   { 
00693     if (index >= data_.size())
00694     {
00695       throw Exception::OutOfGrid(__FILE__, __LINE__);
00696     }       
00697     return data_[index];
00698   }
00699 
00700   template <class ValueType>
00701   BALL_INLINE 
00702   ValueType& TRegularData2D<ValueType>::getData(Position index)
00703   { 
00704     if (index >= data_.size())
00705     {
00706       throw Exception::OutOfGrid(__FILE__, __LINE__);
00707     }       
00708     return data_[index];
00709   }
00710 
00711   template <class ValueType>
00712   BALL_INLINE 
00713   typename TRegularData2D<ValueType>::CoordinateType TRegularData2D<ValueType>::getCoordinates
00714     (const typename TRegularData2D<ValueType>::IndexType& index) const 
00715   {
00716     if ((index.x >= size_.x) || (index.y >= size_.y))
00717     {
00718       throw Exception::OutOfGrid(__FILE__, __LINE__);
00719     }   
00720     
00721     CoordinateType r(origin_.x + index.x * spacing_.x,
00722                      origin_.y + index.y * spacing_.y);
00723 
00724     return r;
00725   }
00726 
00727   template <class ValueType> 
00728   BALL_INLINE 
00729   typename TRegularData2D<ValueType>::CoordinateType
00730   TRegularData2D<ValueType>::getCoordinates(Position position) const 
00731   {
00732     if (position >= data_.size())
00733     {
00734       throw Exception::OutOfGrid(__FILE__, __LINE__);
00735     } 
00736     
00737     Position x = (Position)(position %  size_.x);
00738     Position y = (Position)(position / size_.x);
00739 
00740     return CoordinateType(origin_.x + (double)x * spacing_.x,
00741                           origin_.y + (double)y * spacing_.y);
00742   }
00743 
00744   template <typename ValueType>
00745   BALL_INLINE
00746   void TRegularData2D<ValueType>::getEnclosingIndices
00747     (const typename TRegularData2D<ValueType>::CoordinateType& r,
00748      Position& ll, Position& lr, Position& ul, Position& ur) const
00749   {
00750     if (!isInside(r))
00751     {
00752       throw Exception::OutOfGrid(__FILE__, __LINE__);
00753     }   
00754 
00755     // Calculate the grid indices of the lower left front corner
00756     // of the enclosing rectangle
00757     IndexType position;
00758     position.x = (Position)((r.x - origin_.x) / spacing_.x);
00759     position.y = (Position)((r.y - origin_.y) / spacing_.y);
00760     
00761     // Calculate the (linear) indices of the four rectangle corners
00762     ll = position.x + size_.x * position.y;
00763     lr = ll + 1;
00764     ul = ll + size_.x;
00765     ur = ul + 1;
00766   }
00767 
00768   template <typename ValueType>
00769   BALL_INLINE
00770   void TRegularData2D<ValueType>::getEnclosingValues
00771     (const typename TRegularData2D<ValueType>::CoordinateType& r,
00772      ValueType& ll, ValueType& lr, ValueType& ul, ValueType& ur) const
00773   {
00774     if (!isInside(r))
00775     {
00776       throw Exception::OutOfGrid(__FILE__, __LINE__);
00777     }   
00778     
00779     // compute the four grid indices forming the enclosing rectangle
00780     Position ll_id, lr_id, ul_id, ur_id;
00781     getEnclosingIndices(r, ll_id, lr_id, ul_id, ur_id);
00782         
00783     // Retrieve the grid values
00784     ll = data_[ll_id];
00785     lr = data_[lr_id];
00786     ul = data_[ul_id];
00787     ur = data_[ur_id];
00788   }
00789 
00790   template <typename ValueType>
00791   BALL_INLINE
00792   ValueType TRegularData2D<ValueType>::getInterpolatedValue
00793     (const typename TRegularData2D<ValueType>::CoordinateType& r) const
00794   {
00795     if (!isInside(r))
00796     {
00797       throw Exception::OutOfGrid(__FILE__, __LINE__);
00798     }   
00799     
00800     return this->operator () (r);
00801   }
00802 
00803   template <typename ValueType>
00804   BALL_INLINE
00805   ValueType TRegularData2D<ValueType>::operator ()
00806     (const typename TRegularData2D<ValueType>::CoordinateType& r) const
00807   {
00808     CoordinateType h(r - origin_);
00809     Position x = (Position)(h.x / spacing_.x);
00810     Position y = (Position)(h.y / spacing_.y);
00811 
00812     // correct for numerical inaccuracies
00813     if (x >= (size_.x - 1))
00814     {
00815       x = size_.x - 2;
00816     }
00817     if (y >= (size_.y - 1))
00818     {
00819       y = size_.y - 2;
00820     }
00821 
00822     Size l = x + size_.x * y;
00823     CoordinateType r_0(origin_.x + (double)x * spacing_.x,
00824                        origin_.y + (double)y * spacing_.y);
00825 
00826     double dx = 1.0 - ((r.x - r_0.x) / spacing_.x);
00827     double dy = 1.0 - ((r.y - r_0.y) / spacing_.y);
00828 
00829     return  data_[l] * dx * dy
00830           + data_[l + 1] * (1.0 - dx) * dy
00831           + data_[l + size_.x] * dx * (1.0 - dy)
00832           + data_[l + size_.x + 1] * (1.0 - dx) * (1.0 - dy);
00833   }
00834 
00835   template <typename ValueType>
00836   BALL_INLINE
00837   typename TRegularData2D<ValueType>::IndexType TRegularData2D<ValueType>::getLowerIndex
00838     (const typename TRegularData2D<ValueType>::CoordinateType& r) const
00839   {
00840     if (!isInside(r))
00841     {
00842       throw Exception::OutOfGrid(__FILE__, __LINE__);
00843     }   
00844     
00845     static IndexType position;
00846     position.x = (Position)((r.x - origin_.x) / spacing_.x);
00847     position.y = (Position)((r.y - origin_.y) / spacing_.y);
00848     
00849     return position;
00850   }
00851 
00852   template <typename ValueType>
00853   BALL_INLINE
00854   typename TRegularData2D<ValueType>::IndexType TRegularData2D<ValueType>::getClosestIndex
00855     (const typename TRegularData2D<ValueType>::CoordinateType& r) const
00856   {
00857     if (!isInside(r))
00858     {
00859       throw Exception::OutOfGrid(__FILE__, __LINE__);
00860     }   
00861     
00862     static IndexType position;
00863     position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
00864     position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
00865     
00866     return position;
00867   }
00868 
00869   template <typename ValueType>
00870   BALL_INLINE
00871   const ValueType& TRegularData2D<ValueType>::getClosestValue
00872     (const typename TRegularData2D<ValueType>::CoordinateType& r) const
00873   {
00874     if (!isInside(r))
00875     {
00876       throw Exception::OutOfGrid(__FILE__, __LINE__);
00877     }   
00878     
00879     static IndexType position;
00880     position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
00881     position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
00882 
00883     return operator [] (position);
00884   }
00885 
00886   template <typename ValueType>
00887   BALL_INLINE
00888   ValueType& TRegularData2D<ValueType>::getClosestValue
00889     (const typename TRegularData2D<ValueType>::CoordinateType& r)
00890   {
00891     if (!isInside(r))
00892     {
00893       throw Exception::OutOfGrid(__FILE__, __LINE__);
00894     }   
00895     
00896     static IndexType position;
00897     position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
00898     position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
00899 
00900     return operator [] (position);
00901   }
00902 
00903   template <typename ValueType>
00904   BALL_INLINE
00905   ValueType TRegularData2D<ValueType>::calculateMean() const
00906   {
00907     Position data_points  = (size_.x * size_.y);
00908     ValueType mean = 0;
00909     for (Position i = 0; i < data_points; i++)
00910     {
00911       mean += data_[i];
00912     }
00913     mean /= data_points;
00914     return mean;
00915   }
00916   
00917   template <typename ValueType>
00918   BALL_INLINE
00919   ValueType TRegularData2D<ValueType>::calculateSD() const
00920   {
00921     Position data_points  = (size_.x * size_.y);
00922     ValueType stddev = 0;
00923     ValueType mean = this->calculateMean();
00924     for (Position i = 0; i < data_points; i++)
00925     {
00926       stddev += (pow(data_(i)-mean,2));
00927     }
00928     stddev /= (data_points-1);
00929     stddev = sqrt(stddev);
00930     return stddev;
00931   }
00932     
00933   template <typename ValueType>
00934   void TRegularData2D<ValueType>::clear()
00935   {
00936     data_.resize(0);
00937     
00938     origin_.set(0.0);
00939     dimension_.set(0.0, 0.0);
00940     size_.x = 0;
00941     size_.y = 0;
00942     spacing_.set(1.0, 1.0);
00943   }
00944 
00945   template <typename ValueType> 
00946   bool TRegularData2D<ValueType>::operator == (const TRegularData2D<ValueType>& data) const
00947   {
00948     return ((origin_ == data.origin_)
00949             && (dimension_ == data.dimension_)
00950             && (size_.x == data.size_.x)
00951             && (size_.y == data.size_.y)
00952             && (data_ == data.data_));
00953   } 
00954 
00955 
00958 
00959   template <typename ValueType>
00960   std::ostream& operator << (std::ostream& os, const TRegularData2D<ValueType>& data)
00961   {
00962     // Write the grid origin, dimension, and number of grid points
00963     os << data.getOrigin().x << " " << data.getOrigin().y
00964       << std::endl
00965       << data.getOrigin().x + data.getDimension().x << " "
00966       << data.getOrigin().y + data.getDimension().y
00967       << std::endl
00968       << data.getSize().x - 1 << " " << data.getSize().y - 1
00969       << std::endl;
00970 
00971     // Write the array contents.
00972     std::copy(data.begin(), data.end(), std::ostream_iterator<ValueType>(os, "\n"));
00973     return os;
00974   }
00975 
00977   template <typename ValueType>
00978   std::istream& operator >> (std::istream& is, TRegularData2D<ValueType>& grid)
00979   {
00980     typename TRegularData2D<ValueType>::CoordinateType origin;
00981     typename TRegularData2D<ValueType>::CoordinateType dimension;
00982     typename TRegularData2D<ValueType>::IndexType size;
00983 
00984     is >> origin.x >> origin.y;
00985     is >> dimension.x >> dimension.y;
00986     is >> size.x >> size.y;
00987     
00988     dimension -= origin;
00989     size.x++;
00990     size.y++;
00991 
00992     grid.resize(size);
00993     grid.setOrigin(origin);
00994     grid.setDimension(dimension);
00995 
00996     std::copy(std::istream_iterator<ValueType>(is), 
00997               std::istream_iterator<ValueType>(), 
00998               grid.begin());
00999     //    std::copy_n(std::istream_iterator<ValueType>(is), grid.size(), grid.begin());
01000     
01001     return is;
01002   }
01004 
01005   template <typename ValueType>
01006   void TRegularData2D<ValueType>::binaryWrite(const String& filename) const
01007   {
01008     File outfile(filename.c_str(), std::ios::out|std::ios::binary);
01009     if (!outfile.isValid()) 
01010       throw Exception::FileNotFound(__FILE__, __LINE__, filename);
01011 
01012     // write all information we need to recreate the grid
01013     BinaryFileAdaptor<BlockValueType> adapt_block;
01014     BinaryFileAdaptor<ValueType>      adapt_single;
01015     BinaryFileAdaptor<float>          adapt_float;
01016     
01017     BinaryFileAdaptor<Size>           adapt_size;
01018 
01019     adapt_size.setData(data_.size());
01020     outfile << adapt_size;
01021     
01022     // NOTE: we do not use the binary file adaptor to write out the Vector2-variables here,
01023     //       the reason is a bit stupid: the data layout of Vector2 has a three-byte "hole"
01024     //       that is automatically padded by the compiler to ensure the correct alignment.
01025     //       valgrind is very unhappy when we write these uninitialized values out, even
01026     //       though they are entirely harmless. To prevent false-positives in the error checking,
01027     //       we write out the members of the vectors instead. This even saves some bytes in the
01028     //       output (not that it would matter...)
01029     adapt_float.setData(origin_.x);
01030     outfile << adapt_float;
01031     adapt_float.setData(origin_.y);
01032     outfile << adapt_float;
01033 
01034     adapt_float.setData(dimension_.x);
01035     outfile << adapt_float;
01036     adapt_float.setData(dimension_.y);
01037     outfile << adapt_float;
01038 
01039     adapt_float.setData(spacing_.x);
01040     outfile << adapt_float;
01041     adapt_float.setData(spacing_.y);
01042     outfile << adapt_float;
01043 
01044     BinaryFileAdaptor<IndexType> adapt_index;
01045     adapt_index.setData(size_);
01046     outfile << adapt_index;
01047   
01048     // we slide a window of BLOCK_SIZE over our data.
01049     const int BLOCK_SIZE = 1024;
01050     Index window_pos = 0;
01051     
01052     while (((int)data_.size() - (BLOCK_SIZE + window_pos)) >= 0)
01053     {
01054       adapt_block.setData(* (BlockValueType*)&(data_[window_pos]));
01055       outfile << adapt_block;
01056       window_pos += BLOCK_SIZE;
01057     }
01058 
01059     // Now we have to write the remaining data one by one.
01060     for (Size i = window_pos; i < data_.size(); i++)
01061     {
01062       adapt_single.setData(data_[i]);
01063       outfile << adapt_single;
01064     }
01065 
01066     // That's it. 
01067     outfile.close();
01068   }
01069 
01070   template <typename ValueType>
01071   void TRegularData2D<ValueType>::binaryRead(const String& filename)
01072   {
01073     File infile(filename, std::ios::in|std::ios::binary);
01074     if (!infile.isValid()) 
01075     {
01076       throw Exception::FileNotFound(__FILE__, __LINE__, filename);
01077     }
01078     
01079     BinaryFileAdaptor< BlockValueType > adapt_block;
01080     BinaryFileAdaptor< ValueType >      adapt_single;
01081     
01082     // read all information we need to recreate the grid
01083     BinaryFileAdaptor<Size>           adapt_size;
01084     BinaryFileAdaptor<float>          adapt_float;
01085 
01086     infile >> adapt_size;
01087     Size new_size = adapt_size.getData();
01088   
01089     infile >> adapt_float;
01090     origin_.x = adapt_float.getData();
01091     infile >> adapt_float;
01092     origin_.y = adapt_float.getData();
01093 
01094     infile >> adapt_float;
01095     dimension_.x = adapt_float.getData();
01096     infile >> adapt_float;
01097     dimension_.y = adapt_float.getData();
01098 
01099     infile >> adapt_float;
01100     spacing_.x = adapt_float.getData();
01101     infile >> adapt_float;
01102     spacing_.y = adapt_float.getData();
01103 
01104     BinaryFileAdaptor<IndexType> adapt_index;
01105     infile >> adapt_index;
01106     size_ =  adapt_index.getData();
01107   
01108     data_.resize(new_size);
01109 
01110     // we slide a window of size 1024 over our data
01111     Index window_pos = 0;
01112     
01113     while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
01114     {
01115       infile >> adapt_block;
01116       *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
01117       /*
01118       for (Size i=0; i<1024; i++)
01119       {
01120         data_[i+window_pos] = adapt_block.getData().bt[i];
01121       }
01122       */
01123       window_pos+=1024;
01124     }
01125 
01126     // now we have to read the remaining data one by one
01127     for (Size i=window_pos; i<data_.size(); i++)
01128     {
01129       infile >> adapt_single;
01130       data_[i] = adapt_single.getData();
01131     }
01132 
01133     // that's it. I hope...
01134     infile.close();
01135   } 
01136  } // namespace BALL
01137 
01138 #endif // BALL_DATATYPE_TREGULARDATA2D_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines