BALL  1.4.1
regularData3D.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_REGULARDATA3D_H
00006 #define BALL_DATATYPE_REGULARDATA3D_H
00007 
00008 #ifndef BALL_MATHS_VECTOR3_H
00009 # include <BALL/MATHS/vector3.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 #ifndef BALL_MATHS_COMMON_H
00021 # include <BALL/MATHS/common.h>
00022 #endif
00023 
00024 #include <iostream>
00025 #include <fstream>
00026 #include <iterator>
00027 #include <algorithm>
00028 
00029 namespace BALL 
00030 {
00044   template <typename ValueType>
00045   class TRegularData3D 
00046   {
00047     public:
00048 
00049     BALL_CREATE(TRegularData3D<ValueType>)
00050 
00051     
00054 
00056     class IndexType
00057     {
00058       public:
00059       inline IndexType() : x(0), y(0), z(0) {}
00060       inline IndexType(Position p) : x(p), y(p), z(p) {}
00061       inline IndexType(Position p, Position q, Position r) : x(p), y(q), z(r) {}
00062 
00064       Position x;
00066       Position y;
00068       Position z;
00069     };
00070 
00072     typedef std::vector<ValueType> VectorType;
00074     typedef TVector3<float> CoordinateType;
00076     typedef typename std::vector<ValueType>::iterator Iterator;
00078     typedef typename std::vector<ValueType>::const_iterator ConstIterator;
00080 
00081     //  STL compatibility types
00082     //
00083     typedef ValueType value_type;
00084     typedef typename std::vector<ValueType>::iterator iterator;
00085     typedef typename std::vector<ValueType>::const_iterator const_iterator;
00086     typedef typename std::vector<ValueType>::reference reference;
00087     typedef typename std::vector<ValueType>::const_reference const_reference;
00088     typedef typename std::vector<ValueType>::pointer pointer;
00089     typedef typename std::vector<ValueType>::difference_type difference_type;
00090     typedef typename std::vector<ValueType>::size_type size_type;
00091 
00095 
00099     TRegularData3D(); 
00100   
00104     TRegularData3D(const TRegularData3D<ValueType>& grid);
00105 
00109     TRegularData3D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);
00110 
00114     TRegularData3D(const CoordinateType& origin, const CoordinateType& x_axis, 
00115                    const CoordinateType& y_axis, const CoordinateType& z_axis, const IndexType& size);
00116 
00120     TRegularData3D
00121       (const IndexType& size, 
00122        const CoordinateType& origin = CoordinateType(0.0),
00123        const CoordinateType& dimension = CoordinateType(1.0));
00124 
00127     virtual ~TRegularData3D();
00128 
00132     virtual void clear();
00134 
00138 
00143     TRegularData3D& operator = (const TRegularData3D<ValueType>& data);
00145 
00146 
00150 
00156     bool operator == (const TRegularData3D<ValueType>& grid) const;
00157 
00160     BALL_INLINE bool operator != (const TRegularData3D<ValueType>& grid) const {  return !this->operator == (grid); }
00161 
00163     BALL_INLINE bool empty() const { return data_.empty(); }
00164 
00166     bool isInside(const CoordinateType& r) const;
00168 
00170     BALL_INLINE bool isOrthogonal() const { return is_orthogonal_;}
00171 
00175 
00176     BALL_INLINE ConstIterator begin() const { return data_.begin(); }
00178     BALL_INLINE ConstIterator end() const { return data_.end(); }
00180     BALL_INLINE Iterator begin() { return data_.begin(); }
00182     BALL_INLINE Iterator end() { return data_.end(); }
00184 
00188 
00189     // STL compatibility
00190     BALL_INLINE size_type size() const { return data_.size(); }
00191     BALL_INLINE size_type max_size() const { return data_.max_size(); }
00192     BALL_INLINE void swap(TRegularData3D<ValueType>& grid) { std::swap(*this, grid); }
00193 
00194 
00196     const vector<ValueType>& getData() const;
00197 
00202     const ValueType& getData(const IndexType& index) const;
00203 
00208     ValueType& getData(const IndexType& index);
00209 
00214     const ValueType& getData(Position index) const;
00215 
00220     ValueType& getData(Position index);
00221 
00226     const ValueType& operator [] (const IndexType& index) const 
00227     {   
00228       return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y]; 
00229     }
00230 
00235     ValueType& operator [] (const IndexType& index)
00236     {   
00237       return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y]; 
00238     }
00239 
00244     const ValueType& operator [] (Position index) const { return data_[index]; }
00245 
00250     ValueType& operator [] (Position index) { return data_[index]; }
00251 
00262     ValueType operator () (const CoordinateType& x) const;
00263     
00270     ValueType getInterpolatedValue(const CoordinateType& x) const;
00271 
00278     const ValueType& getClosestValue(const CoordinateType& x) const;
00279 
00286     ValueType& getClosestValue(const CoordinateType& x);
00287 
00293     IndexType getClosestIndex(const CoordinateType& v) const;
00294 
00299     IndexType getLowerIndex(const CoordinateType& v) const;
00300 
00306     inline const IndexType& getSize() const { return size_; }
00307     
00312     inline const CoordinateType& getOrigin() const { return origin_; }
00313 
00318     const CoordinateType& getSpacing() const {  return spacing_; }
00321     void setOrigin(const CoordinateType& origin) { origin_ = origin; }
00322 
00328     const CoordinateType& getDimension() const { return dimension_; }
00329 
00338     void setDimension(const CoordinateType& dimension)  
00339     {
00340       dimension_ = dimension;
00341       if (is_orthogonal_)
00342       {
00343         spacing_.x = dimension_.x / (double)(size_.x - 1);
00344         spacing_.y = dimension_.y / (double)(size_.y - 1);
00345         spacing_.z = dimension_.z / (double)(size_.z - 1);
00346       }
00347     }
00348 
00364     void resize(const IndexType& size);
00365 
00378     void rescale(const IndexType& new_size);
00379 
00384     CoordinateType getCoordinates(const IndexType& index) const;
00385 
00390     CoordinateType getCoordinates(Position index) const;
00391 
00407     void getEnclosingIndices
00408       (const CoordinateType& r,
00409       Position& llf, Position& rlf, Position& luf, Position& ruf,
00410       Position& llb, Position& rlb, Position& lub, Position& rub) const;
00411 
00416     void getEnclosingValues
00417       (const CoordinateType& r,
00418        ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
00419        ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const;
00420     
00424     ValueType calculateMean() const;
00425     
00429     ValueType calculateSD() const;
00430     
00435     void binaryWrite(const String& filename) const;   
00436 
00444     void binaryWriteRaw(const String& filename) const;
00445 
00450     void binaryRead(const String& filename);
00452     
00453     protected:
00454       
00456     const CoordinateType mapToCartesian_(CoordinateType r) const
00457     {
00458       Vector3 result;
00459       r.x /= (size_.x - 1.);
00460       r.y /= (size_.y - 1.);
00461       r.z /= (size_.z - 1.);
00462 
00463       result.x = mapping_[0] * r.x + mapping_[1] * r.y + mapping_[2] * r.z + origin_.x;
00464       result.y = mapping_[3] * r.x + mapping_[4] * r.y + mapping_[5] * r.z + origin_.y;
00465       result.z = mapping_[6] * r.x + mapping_[7] * r.y + mapping_[8] * r.z + origin_.z;
00466 
00467       return result;
00468     }
00469     
00471     const CoordinateType mapInverse_(CoordinateType r) const
00472     {
00473       r -= origin_;
00474 
00475       Vector3 result;
00476       result.x = inverse_mapping_[0] * r.x + inverse_mapping_[1] * r.y + inverse_mapping_[2] * r.z;
00477       result.y = inverse_mapping_[3] * r.x + inverse_mapping_[4] * r.y + inverse_mapping_[5] * r.z;
00478       result.z = inverse_mapping_[6] * r.x + inverse_mapping_[7] * r.y + inverse_mapping_[8] * r.z;
00479 
00480       result.x *= (size_.x - 1);
00481       result.y *= (size_.y - 1);
00482       result.z *= (size_.z - 1);
00483 
00484       return result;
00485     }
00486     
00488     VectorType data_;
00489 
00491     CoordinateType  origin_;
00492 
00494     CoordinateType  dimension_;
00495 
00497     CoordinateType  spacing_;
00498 
00500     IndexType size_;
00501 
00503     typedef struct { ValueType bt[1024]; } BlockValueType;
00504 
00506     bool is_orthogonal_;
00507 
00509     std::vector<double> mapping_;
00510     std::vector<double> inverse_mapping_;
00511   };
00512 
00515   typedef TRegularData3D<float> RegularData3D;
00516 
00517   // default constructor.
00518   template <class ValueType>
00519   TRegularData3D<ValueType>::TRegularData3D()
00520     : data_(0),
00521       origin_(0.0),
00522       dimension_(0.0),
00523       spacing_(1.0),
00524       size_(0, 0, 0),
00525       is_orthogonal_(true)
00526   {
00527   }
00528 
00529   // copy constructor
00530   template <class ValueType>
00531   TRegularData3D<ValueType>::TRegularData3D
00532     (const TRegularData3D<ValueType>& data)
00533     : data_(),
00534       origin_(data.origin_),
00535       dimension_(data.dimension_),
00536       spacing_(data.spacing_),
00537       size_(data.size_),
00538       is_orthogonal_(data.is_orthogonal_),
00539       mapping_(data.mapping_),
00540       inverse_mapping_(data.inverse_mapping_)
00541   {
00542     try
00543     {
00544       data_ = data.data_;
00545     }
00546     catch (std::bad_alloc&)
00547     {
00548       data_.resize(0);
00549       throw Exception::OutOfMemory(__FILE__, __LINE__, data.data_.size() * sizeof(ValueType));
00550     }
00551   }
00552 
00553   template <class ValueType>
00554   TRegularData3D<ValueType>::TRegularData3D
00555     (const typename TRegularData3D<ValueType>::IndexType& size,
00556      const typename TRegularData3D<ValueType>::CoordinateType& origin,
00557      const typename TRegularData3D<ValueType>::CoordinateType& dimension)
00558     : data_(),
00559       origin_(origin),
00560       dimension_(dimension),
00561       spacing_(0.0),
00562       size_(size),
00563       is_orthogonal_(true)
00564   {
00565     // Compute the grid spacing
00566     spacing_.x = dimension_.x / (double)(size_.x - 1);
00567     spacing_.y = dimension_.y / (double)(size_.y - 1);
00568     spacing_.z = dimension_.z / (double)(size_.z - 1);
00569 
00570     // Compute the number of grid points
00571     size_type number_of_points = size_.x * size_.y * size_.z;
00572     try
00573     {
00574       data_.resize(number_of_points);
00575     }
00576     catch (std::bad_alloc&)
00577     {
00578       data_.resize(0);
00579       throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_points * sizeof(ValueType));
00580     }
00581   }
00582 
00583 
00584   template <class ValueType>
00585   TRegularData3D<ValueType>::TRegularData3D
00586     (const typename TRegularData3D<ValueType>::CoordinateType& origin,
00587      const typename TRegularData3D<ValueType>::CoordinateType& dimension,
00588      const typename TRegularData3D<ValueType>::CoordinateType& spacing)
00589     : data_(),
00590       origin_(origin),
00591       dimension_(dimension),
00592       spacing_(spacing),
00593       size_(0),
00594       is_orthogonal_(true)
00595   {
00596     // Compute the grid size
00597     size_.x = (Size)(dimension_.x / spacing_.x + 0.5) + 1;
00598     size_.y = (Size)(dimension_.y / spacing_.y + 0.5) + 1;
00599     size_.z = (Size)(dimension_.z / spacing_.z + 0.5) + 1;
00600 
00601     // Compute the number of grid points
00602     size_type size = size_.x * size_.y * size_.z;
00603     try
00604     {
00605       data_ .resize(size);
00606     }
00607     catch (std::bad_alloc&)
00608     {
00609       data_.resize(0);
00610       throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
00611     }
00612 
00613     // Adjust the spacing -- dimension has precedence.
00614     spacing_.x = dimension_.x / (double)(size_.x - 1);
00615     spacing_.y = dimension_.y / (double)(size_.y - 1);
00616     spacing_.z = dimension_.z / (double)(size_.z - 1);
00617   }
00618 
00619   template <class ValueType>
00620   TRegularData3D<ValueType>::TRegularData3D
00621     (const typename TRegularData3D<ValueType>::CoordinateType& origin, 
00622      const typename TRegularData3D<ValueType>::CoordinateType& x_axis,
00623      const typename TRegularData3D<ValueType>::CoordinateType& y_axis, 
00624      const typename TRegularData3D<ValueType>::CoordinateType& z_axis, 
00625      const typename TRegularData3D<ValueType>::IndexType& new_size)
00626     : data_(),
00627       origin_(origin),
00628       dimension_(x_axis+y_axis+z_axis),
00629       spacing_(0.0),
00630       size_(new_size),
00631       is_orthogonal_(false)
00632   {
00633     // compute the spacing
00634     spacing_.x = x_axis.getLength() / (new_size.x - 1.);
00635     spacing_.y = y_axis.getLength() / (new_size.y - 1.);
00636     spacing_.z = z_axis.getLength() / (new_size.z - 1.);  
00637 
00638     size_type size = size_.x * size_.y * size_.z;
00639     try
00640     {
00641       data_.resize(size);
00642     }
00643     catch (std::bad_alloc&)
00644     {
00645       data_.resize(0);
00646       throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
00647     }
00648 
00649     // prepare the mapping matrix and its inverse
00650     mapping_.resize(9);
00651     inverse_mapping_.resize(9);
00652 
00653     mapping_[0] = x_axis.x; mapping_[1] = y_axis.x; mapping_[2] = z_axis.x;
00654     mapping_[3] = x_axis.y; mapping_[4] = y_axis.y; mapping_[5] = z_axis.y;
00655     mapping_[6] = x_axis.z; mapping_[7] = y_axis.z; mapping_[8] = z_axis.z;
00656 
00657     // this is numerically instable and only works well for the "simple"
00658     // cases. should be replaced by QR or an SVD
00659 
00660     // just for readability
00661     double a = mapping_[0]; double b = mapping_[1]; double c = mapping_[2];
00662     double d = mapping_[3]; double e = mapping_[4]; double f = mapping_[5];
00663     double g = mapping_[6]; double h = mapping_[7]; double i = mapping_[8];
00664 
00665     double determinant = 1. / (a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g));
00666 
00667     inverse_mapping_[0] = determinant * (e*i - f*h);
00668     inverse_mapping_[1] = determinant * (b*i - c*h);
00669     inverse_mapping_[2] = determinant * (b*f - c*e);
00670 
00671     inverse_mapping_[3] = determinant * (f*g - d*i);
00672     inverse_mapping_[4] = determinant * (a*i - c*g);
00673     inverse_mapping_[5] = determinant * (c*d - a*f);
00674 
00675     inverse_mapping_[6] = determinant * (d*h - e*g);
00676     inverse_mapping_[7] = determinant * (b*g - a*h);
00677     inverse_mapping_[8] = determinant * (a*e - b*d);
00678   }   
00679 
00680   template <class ValueType>
00681   TRegularData3D<ValueType>::~TRegularData3D()
00682   {
00683   }
00684 
00685   template <typename ValueType>
00686   BALL_INLINE
00687   TRegularData3D<ValueType>& TRegularData3D<ValueType>::operator =
00688     (const TRegularData3D<ValueType>& rhs)
00689   {
00690     // Avoid self-assignment.
00691     if (&rhs != this)
00692     {
00693       // Copy the coordinate-related attributes and
00694       // the size.
00695       origin_ = rhs.origin_;
00696       dimension_ = rhs.dimension_;
00697       spacing_ = rhs.spacing_;
00698       size_ = rhs.size_;
00699       
00700       is_orthogonal_ = rhs.is_orthogonal_;
00701       mapping_ = rhs.mapping_;
00702       inverse_mapping_ = rhs.inverse_mapping_;
00703       // Copy the data itself and rethrow allocation exceptions.
00704       try
00705       {
00706         data_ = rhs.data_;
00707       }
00708       catch (std::bad_alloc&)
00709       {
00710         data_.resize(0);
00711         throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.data_.size() * sizeof(ValueType));
00712       }
00713     }
00714 
00715     return *this;
00716   }
00717 
00718   template <typename ValueType>
00719   void TRegularData3D<ValueType>::resize(const typename TRegularData3D<ValueType>::IndexType& size)
00720   {
00721     if (is_orthogonal_)
00722     {
00723       // If the old size equals the new size, we're done.
00724       if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
00725       {
00726         return;
00727       }
00728 
00729       // If the new grid is empty, this whole thing is quite easy.
00730       if ((size.x == 0) || (size.y == 0) || (size.z == 0))
00731       {
00732         data_.resize(0);
00733         dimension_.set(0.0, 0.0, 0.0);
00734         return;
00735       }
00736       // Compute the new array size.
00737       size_type new_size = (size_type)(size.x * size.y * size.z);
00738 
00739       // Catch any bad_allocs thrown by vector::resize
00740       try
00741       {
00742         // Create a new temporary array.
00743         std::vector<ValueType> old_data(data_);
00744 
00745         // Resize the data to its new size.
00746         data_.resize(new_size);
00747 
00748         // walk over the new grid and copy the old stuff back.
00749         static ValueType default_value = (ValueType)0;
00750         for (size_type i = 0; i < new_size; i++)
00751         {
00752           size_type x = i % size.x;
00753           size_type y = (i % (size.x * size.y)) / size.x;
00754           size_type z = i / (size.x * size.y);
00755           if ((x >= size_.x) || (y >= size_.y) || (z >= size_.z))
00756           {
00757             data_[i] = default_value;
00758           }
00759           else
00760           {
00761             data_[i] = old_data[x + y * size_.x + z * size_.x * size_.y];
00762           }
00763         }
00764 
00765         // Correct the grid dimension. Origin and spacing remain constant.
00766         if ((size_.x != 0) && (size_.y != 0) && (size_.z != 0))
00767         {
00768           dimension_.x *= (double)size.x / (double)size_.x;
00769           dimension_.y *= (double)size.y / (double)size_.y;
00770           dimension_.z *= (double)size.z / (double)size_.z;
00771         }
00772         size_ = size;
00773       }
00774       catch (std::bad_alloc&)
00775       {
00776         throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
00777       }
00778     }
00779   }
00780 
00781   template <typename ValueType>
00782   void TRegularData3D<ValueType>::rescale(const typename TRegularData3D<ValueType>::IndexType& size)
00783   {
00784     if (is_orthogonal_)
00785     {
00786       // If the old size equals the new size, we're done.
00787       if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
00788       {
00789         return;
00790       }
00791 
00792       // If the new grid is empty, this whole thing is quite easy.
00793       if ((size.x == 0) || (size.y == 0) || (size.z == 0))
00794       {
00795         data_.resize(0);
00796         dimension_.set(0.0);
00797         return;
00798       }
00799 
00800       // Compute the new array size.
00801       size_type new_size = (size_type)(size.x * size.y * size.z);
00802 
00803       // Catch any bad_allocs thrown by vector::resize
00804       try
00805       {
00806         // Create a new temporary array.
00807         TRegularData3D<ValueType> old_data(*this);
00808 
00809         // Resize the data array to its new size.
00810         data_.resize(new_size);
00811         spacing_.x = dimension_.x / (double)(size.x - 1);
00812         spacing_.y = dimension_.y / (double)(size.y - 1);
00813         spacing_.z = dimension_.z / (double)(size.z - 1);
00814 
00815         // Correct the grid size. Origin and dimension remain constant.
00816         size_ = size;
00817 
00818         // Walk over the new grid and copy the (interpolated) old stuff back.
00819         for (size_type i = 0; i < data_.size(); i++)
00820         {
00821           try
00822           {
00823             data_[i] = old_data.getInterpolatedValue(getCoordinates(i));
00824           }
00825           catch (Exception::OutOfGrid&)
00826           {
00827             data_[i] = old_data.getClosestValue(getCoordinates(i));
00828           }
00829         }
00830       }
00831       catch (std::bad_alloc&)
00832       {
00833         throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
00834       }
00835     }
00836   }
00837 
00838   template <class ValueType> 
00839   BALL_INLINE
00840   bool TRegularData3D<ValueType>::isInside(const typename TRegularData3D<ValueType>::CoordinateType& r) const   
00841   {
00842     if (is_orthogonal_)
00843     {
00844       return !((r.x > (origin_.x + dimension_.x )) 
00845           || (r.y > (origin_.y + dimension_.y))  
00846           || (r.z > (origin_.z + dimension_.z))
00847           || (r.x < origin_.x) || (r.y < origin_.y) || (r.z < origin_.z));
00848     } 
00849     else 
00850     {
00851       // compute A^-1 * pos and see whether the indices are part of the grid
00852       CoordinateType ri = mapInverse_(r);
00853       ri.x = Maths::round(ri.x);
00854       ri.y = Maths::round(ri.y);
00855       ri.z = Maths::round(ri.z);
00856       
00857       return !(   (ri.x < 0) || (ri.y < 0) || (ri.z < 0)
00858                || (ri.x >= size_.x) || (ri.y >= size_.y) || (ri.z >= size_.z) );
00859     }
00860   }
00861 
00862   template <class ValueType>
00863   BALL_INLINE
00864   const ValueType& TRegularData3D<ValueType>::getData
00865     (const typename TRegularData3D<ValueType>::IndexType& index) const
00866   {
00867     size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
00868     if (pos >= data_.size())
00869     {
00870       throw Exception::OutOfGrid(__FILE__, __LINE__);
00871     }
00872     return data_[pos];
00873   }
00874 
00875   template <class ValueType>
00876   BALL_INLINE
00877   const vector<ValueType>& TRegularData3D<ValueType>::getData() const
00878   {
00879     return data_;
00880   }
00881 
00882   template <typename ValueType>
00883   BALL_INLINE
00884   ValueType& TRegularData3D<ValueType>::getData
00885     (const typename TRegularData3D<ValueType>::IndexType& index)
00886   {
00887     size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
00888     if (pos >= data_.size())
00889     {
00890       throw Exception::OutOfGrid(__FILE__, __LINE__);
00891     }
00892     return data_[pos];
00893   }
00894 
00895   template <class ValueType>
00896   BALL_INLINE
00897   const ValueType& TRegularData3D<ValueType>::getData(Position index) const
00898   {
00899     if (index >= data_.size())
00900     {
00901       throw Exception::OutOfGrid(__FILE__, __LINE__);
00902     }
00903     return data_[index];
00904   }
00905 
00906   template <class ValueType>
00907   BALL_INLINE
00908   ValueType& TRegularData3D<ValueType>::getData(Position index)
00909   {
00910     if (index >= data_.size())
00911     {
00912       throw Exception::OutOfGrid(__FILE__, __LINE__);
00913     }
00914     return data_[index];
00915   }
00916 
00917   template <class ValueType>
00918   BALL_INLINE
00919   typename TRegularData3D<ValueType>::CoordinateType TRegularData3D<ValueType>::getCoordinates
00920     (const typename TRegularData3D<ValueType>::IndexType& index) const
00921   {
00922     if ((index.x >= size_.x) || (index.y >= size_.y) || (index.z >= size_.z))
00923     {
00924       throw Exception::OutOfGrid(__FILE__, __LINE__);
00925     }
00926 
00927     if (is_orthogonal_)
00928     {
00929       CoordinateType r(origin_.x + index.x * spacing_.x,
00930           origin_.y + index.y * spacing_.y,
00931           origin_.z + index.z * spacing_.z);
00932       return r;
00933     }
00934     else
00935     {
00936       CoordinateType r(index.x, index.y, index.z);
00937       r = mapToCartesian_(r);
00938 
00939       return r;
00940     }
00941   }
00942 
00943   template <class ValueType>
00944   BALL_INLINE
00945   typename TRegularData3D<ValueType>::CoordinateType
00946   TRegularData3D<ValueType>::getCoordinates(Position position) const
00947   {
00948     if (position >= data_.size())
00949     {
00950       throw Exception::OutOfGrid(__FILE__, __LINE__);
00951     }
00952 
00953     Position x = (Position)(position %  size_.x);
00954     Position y = (Position)((position % (size_.x * size_.y))/ size_.x);
00955     Position z = (Position)(position / (size_.x * size_.y));
00956 
00957     if (is_orthogonal_)
00958     {
00959       return CoordinateType(origin_.x + (double)x * spacing_.x,
00960                             origin_.y + (double)y * spacing_.y,
00961                             origin_.z + (double)z * spacing_.z);
00962     }
00963     else
00964     {
00965       CoordinateType r(x,y,z);
00966       r = mapToCartesian_(r);
00967 
00968       return r;
00969     }
00970   }
00971 
00972   template <typename ValueType>
00973   BALL_INLINE
00974   void TRegularData3D<ValueType>::getEnclosingIndices
00975     (const typename TRegularData3D<ValueType>::CoordinateType& r,
00976     Position& llf, Position& rlf, Position& luf, Position& ruf,
00977     Position& llb, Position& rlb, Position& lub, Position& rub) const
00978   {
00979     if (!isInside(r))
00980     {
00981       throw Exception::OutOfGrid(__FILE__, __LINE__);
00982     }   
00983 
00984     // calculate the grid indices of the lower left front corner
00985     // of the enclosing box
00986     IndexType position;
00987     if (is_orthogonal_)
00988     {
00989       position.x = (Position)((r.x - origin_.x) / spacing_.x);
00990       position.y = (Position)((r.y - origin_.y) / spacing_.y);
00991       position.z = (Position)((r.z - origin_.z) / spacing_.z);
00992     }
00993     else
00994     {
00995       CoordinateType pos = mapInverse_(r);
00996       position.x = (Position) pos.x;  
00997       position.y = (Position) pos.y;  
00998       position.z = (Position) pos.z;  
00999     } 
01000     
01001     // calculate the (linear) indices of the eight
01002     // box corners
01003     llf = position.x + size_.x * position.y 
01004         + size_.x * size_.y * position.z;
01005     rlf = llf + 1;
01006     luf = llf + size_.x;
01007     ruf = luf + 1;
01008     llb = llf + size_.x * size_.y;
01009     rlb = llb + 1;
01010     lub = llb + size_.x;
01011     rub = lub + 1;
01012   }
01013 
01014   template <typename ValueType>
01015   BALL_INLINE
01016   void TRegularData3D<ValueType>::getEnclosingValues
01017     (const typename TRegularData3D<ValueType>::CoordinateType& r,
01018     ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
01019     ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const
01020   {
01021     if (!isInside(r))
01022     {
01023       throw Exception::OutOfGrid(__FILE__, __LINE__);
01024     }   
01025     
01026     // compute the eight grid indices forming the enclosing box
01027     Position llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx;
01028     getEnclosingIndices(r, llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx);
01029         
01030     // retrieve the grid values
01031     llf = data_[llf_idx];
01032     rlf = data_[rlf_idx];
01033     luf = data_[luf_idx];
01034     ruf = data_[ruf_idx];
01035     llb = data_[llb_idx];
01036     rlb = data_[rlb_idx];
01037     lub = data_[lub_idx];
01038     rub = data_[rub_idx];
01039   }
01040 
01041   template <typename ValueType>
01042   BALL_INLINE
01043   ValueType TRegularData3D<ValueType>::getInterpolatedValue
01044     (const typename TRegularData3D<ValueType>::CoordinateType& r) const
01045   {
01046     if (!isInside(r))
01047     {
01048       throw Exception::OutOfGrid(__FILE__, __LINE__);
01049     }   
01050     return operator () (r);
01051   }
01052 
01053   template <typename ValueType>
01054   BALL_INLINE
01055   ValueType TRegularData3D<ValueType>::operator ()
01056     (const typename TRegularData3D<ValueType>::CoordinateType& r) const
01057   {
01058     Position x;
01059     Position y;   
01060     Position z;   
01061     TVector3<double> r_0;
01062     
01063     if (is_orthogonal_)
01064     { 
01065       Vector3 h(r - origin_);
01066 
01067       x = (Position)(h.x / spacing_.x);
01068       y = (Position)(h.y / spacing_.y);
01069       z = (Position)(h.z / spacing_.z);
01070 
01071       while (x >= (size_.x - 1)) x--;
01072       while (y >= (size_.y - 1)) y--;
01073       while (z >= (size_.z - 1)) z--;
01074 
01075       r_0.x = origin_.x + x*spacing_.x;
01076       r_0.y = origin_.y + y*spacing_.y;
01077       r_0.z = origin_.z + z*spacing_.z;
01078     }
01079     else
01080     {
01081       Vector3 pos = mapInverse_(r);
01082       x = (Position) pos.x;
01083       y = (Position) pos.y;
01084       z = (Position) pos.z;
01085 
01086       while (x >= (size_.x - 1)) x--;
01087       while (y >= (size_.y - 1)) y--;
01088       while (z >= (size_.z - 1)) z--;
01089 
01090       // This can probably be done much faster...
01091       Vector3 lower_pos(x,y,z);
01092       lower_pos = mapToCartesian_(lower_pos);
01093       r_0.x = lower_pos.x;
01094       r_0.y = lower_pos.y;
01095       r_0.z = lower_pos.z;
01096     }
01097     
01098     Position Nx = size_.x;
01099     Position Nxy = size_.y * Nx;
01100     Position l = x + Nx * y + Nxy * z;
01101 
01102     double dx = 1. - ((double)(r.x - r_0.x) / spacing_.x);
01103     double dy = 1. - ((double)(r.y - r_0.y) / spacing_.y);
01104     double dz = 1. - ((double)(r.z - r_0.z) / spacing_.z);
01105 
01106     return  data_[l] * dx * dy * dz
01107           + data_[l + 1] * (1. - dx) * dy * dz
01108           + data_[l + Nx] * dx * (1. - dy) * dz
01109           + data_[l + Nx + 1] * (1. - dx) * (1. - dy) * dz
01110           + data_[l + Nxy] * dx * dy * (1. - dz)
01111           + data_[l + Nxy + 1] * (1. - dx) * dy * (1. - dz)
01112           + data_[l + Nxy + Nx] * dx * (1. - dy) * (1. - dz)
01113           + data_[l + Nxy + Nx + 1] * (1. - dx) * (1. - dy) * (1. - dz);
01114   }
01115 
01116   template <typename ValueType>
01117   BALL_INLINE
01118   typename TRegularData3D<ValueType>::IndexType TRegularData3D<ValueType>::getClosestIndex
01119     (const typename TRegularData3D<ValueType>::CoordinateType& r) const
01120   {
01121     if (!isInside(r))
01122     {
01123       throw Exception::OutOfGrid(__FILE__, __LINE__);
01124     }
01125 
01126     static IndexType position;
01127 
01128     if (is_orthogonal_)
01129     {
01130       position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
01131       position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
01132       position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
01133     }
01134     else
01135     {
01136       Vector3 pos = mapInverse_(r);
01137       position.x = (Position) Maths::round(pos.x);
01138       position.y = (Position) Maths::round(pos.y);
01139       position.z = (Position) Maths::round(pos.z);
01140     }
01141       
01142     return position;
01143   }
01144 
01145   template <typename ValueType>
01146   BALL_INLINE
01147   typename TRegularData3D<ValueType>::IndexType  TRegularData3D<ValueType>::getLowerIndex
01148     (const typename TRegularData3D<ValueType>::CoordinateType& r) const
01149   {
01150     if (!isInside(r))
01151     {
01152       throw Exception::OutOfGrid(__FILE__, __LINE__);
01153     }
01154 
01155     static IndexType position;
01156     if (is_orthogonal_)
01157     {
01158       position.x = (Position)((r.x - origin_.x) / spacing_.x);
01159       position.y = (Position)((r.y - origin_.y) / spacing_.y);
01160       position.z = (Position)((r.z - origin_.z) / spacing_.z);
01161     }
01162     else
01163     {
01164       Vector3 pos = mapInverse_(r);
01165       position.x = (Position)pos.x;
01166       position.y = (Position)pos.y;
01167       position.z = (Position)pos.z;
01168     }
01169 
01170     return position;
01171   }
01172 
01173   template <typename ValueType>
01174   BALL_INLINE
01175   const ValueType& TRegularData3D<ValueType>::getClosestValue
01176     (const typename TRegularData3D<ValueType>::CoordinateType& r) const
01177   {
01178     if (!isInside(r))
01179     {
01180       throw Exception::OutOfGrid(__FILE__, __LINE__);
01181     }
01182 
01183     static IndexType position;
01184 
01185     if (is_orthogonal_)
01186     {
01187       position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
01188       position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
01189       position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
01190     }
01191     else
01192     {
01193       static Vector3 pos = mapInverse_(r);
01194       position.x = (Position) Maths::round(pos.x);
01195       position.y = (Position) Maths::round(pos.y);
01196       position.z = (Position) Maths::round(pos.z);
01197     }
01198 
01199     return operator [] (position);
01200   }
01201 
01202   template <typename ValueType>
01203   BALL_INLINE
01204   ValueType& TRegularData3D<ValueType>::getClosestValue
01205     (const typename TRegularData3D<ValueType>::CoordinateType& r)
01206   {
01207     if (!isInside(r))
01208     {
01209       throw Exception::OutOfGrid(__FILE__, __LINE__);
01210     }
01211 
01212     static IndexType position;
01213 
01214     if (is_orthogonal_)
01215     {
01216       position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
01217       position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
01218       position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
01219     }
01220     else
01221     {
01222       static Vector3 pos = mapInverse_(r);
01223       position.x = (Position) Maths::round(pos.x);
01224       position.y = (Position) Maths::round(pos.y);
01225       position.z = (Position) Maths::round(pos.z);
01226     }
01227 
01228 
01229     return operator [] (position);
01230   }
01231 
01232   template <typename ValueType>
01233   BALL_INLINE
01234   ValueType TRegularData3D<ValueType>::calculateMean() const
01235   {
01236     Position data_points  = (size_.x * size_.y * size_.z);
01237     ValueType mean = 0;
01238     for (Position i = 0; i < data_points; i++)
01239     {
01240       mean += data_[i];
01241     }
01242     mean /= data_points;
01243     return mean;
01244   }
01245 
01246   template <typename ValueType>
01247   BALL_INLINE
01248   ValueType TRegularData3D<ValueType>::calculateSD() const
01249   {
01250     Position data_points  = (size_.x * size_.y * size_.z);
01251     ValueType stddev = 0;
01252     ValueType mean = this->calculateMean();
01253     for (Position i = 0; i < data_points; i++)
01254     {
01255       stddev += (pow(data_[i]-mean,2));
01256     }
01257     stddev /= (data_points-1);
01258     stddev = sqrt(stddev);
01259     return stddev;
01260   }
01261   
01262   template <typename ValueType>
01263   void TRegularData3D<ValueType>::clear()
01264   {
01265     data_.resize(0);
01266 
01267     origin_.set(0.0);
01268     dimension_.set(0.0);
01269     size_.x = 0;
01270     size_.y = 0;
01271     size_.z = 0;
01272     spacing_.set(1.0);    
01273     is_orthogonal_ = true;
01274   }
01275 
01276   
01277   template <typename ValueType> 
01278   bool TRegularData3D<ValueType>::operator == (const TRegularData3D<ValueType>& grid) const 
01279   {
01280     return ((origin_ == grid.origin_)
01281             && (dimension_ == grid.dimension_)
01282             && (size_.x == grid.size_.x)
01283             && (size_.y == grid.size_.y)
01284             && (size_.z == grid.size_.z)
01285             && (data_ == grid.data_)
01286             && (is_orthogonal_ == grid.is_orthogonal_));
01287   } 
01288 
01289   template <typename ValueType>
01290   void TRegularData3D<ValueType>::binaryWrite(const String& filename) const
01291   {
01292     File outfile(filename, std::ios::out|std::ios::binary);
01293     if (!outfile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
01294     
01295     BinaryFileAdaptor< BlockValueType > adapt_block;
01296     BinaryFileAdaptor< ValueType >       adapt_single;
01297     
01298     // TODO: check for endiannes and swap bytes accordingly
01299 
01300     // write all information we need to recreate the grid
01301     BinaryFileAdaptor<CoordinateType> adapt_coordinate;
01302     BinaryFileAdaptor<Size>           adapt_size;
01303 
01304     adapt_size.setData(data_.size());
01305     outfile << adapt_size;
01306     
01307     adapt_coordinate.setData(origin_);
01308     outfile << adapt_coordinate;
01309 
01310     adapt_coordinate.setData(dimension_);
01311     outfile << adapt_coordinate;
01312 
01313     adapt_coordinate.setData(spacing_);
01314     outfile << adapt_coordinate;
01315 
01316     BinaryFileAdaptor<IndexType> adapt_index;
01317     adapt_index.setData(size_);
01318     outfile << adapt_index;
01319   
01320     // we slide a window of size 1024 over our data
01321     Index window_pos = 0;
01322     
01323     while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
01324     {
01325       adapt_block.setData(* (BlockValueType*)&(data_[window_pos]));
01326       outfile << adapt_block;
01327       window_pos+=1024;
01328     }
01329 
01330     // now we have to write the remaining data one by one
01331     for (Size i=window_pos; i<data_.size(); i++)
01332     {
01333       adapt_single.setData(data_[i]);
01334       outfile << adapt_single;
01335     }
01336 
01337     // that's it. I hope...
01338     outfile.close();
01339   }
01340 
01341   template <typename ValueType>
01342   void TRegularData3D<ValueType>::binaryRead(const String& filename)
01343   {
01344     File infile(filename, std::ios::in|std::ios::binary);
01345     if (!infile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
01346     
01347     BinaryFileAdaptor< BlockValueType > adapt_block;
01348     BinaryFileAdaptor< ValueType >      adapt_single;
01349     
01350     // TODO: check for endiannes and swap bytes accordingly
01351 
01352     // read all information we need to recreate the grid
01353     BinaryFileAdaptor<CoordinateType> adapt_coordinate;
01354     BinaryFileAdaptor<Size>           adapt_size;
01355 
01356     infile >> adapt_size;
01357     Size new_size = adapt_size.getData();
01358   
01359     infile >> adapt_coordinate;
01360     origin_ = adapt_coordinate.getData();
01361 
01362     infile >> adapt_coordinate;
01363     dimension_ = adapt_coordinate.getData();
01364 
01365     infile >> adapt_coordinate;
01366     spacing_ = adapt_coordinate.getData();
01367 
01368     BinaryFileAdaptor<IndexType> adapt_index;
01369     infile >> adapt_index;
01370     size_ =  adapt_index.getData();
01371   
01372     data_.resize(new_size);
01373 
01374     // we slide a window of size 1024 over our data
01375     Index window_pos = 0;
01376     
01377     while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
01378     {
01379       infile >> adapt_block;
01380       *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
01381       /*
01382       for (Size i=0; i<1024; i++)
01383       {
01384         data_[i+window_pos] = adapt_block.getData().bt[i];
01385       }
01386       */
01387       window_pos+=1024;
01388     }
01389 
01390     // now we have to read the remaining data one by one
01391     for (Size i=window_pos; i<data_.size(); i++)
01392     {
01393       infile >> adapt_single;
01394       data_[i] = adapt_single.getData();
01395     }
01396 
01397     // that's it. I hope...
01398     infile.close();
01399   }
01400 
01403 
01404   template <typename ValueType>
01405   std::ostream& operator << (std::ostream& os, const TRegularData3D<ValueType>& grid)
01406   {
01407     // Write the grid origin, dimension, and number of grid points
01408     os << grid.getOrigin().x << " " << grid.getOrigin().y << " " << grid.getOrigin().z
01409        << std::endl
01410        << grid.getOrigin().x + grid.getDimension().x << " " 
01411        << grid.getOrigin().y + grid.getDimension().y << " " 
01412        << grid.getOrigin().z + grid.getDimension().z
01413        << std::endl
01414        << grid.getSize().x - 1 << " " << grid.getSize().y - 1 << " " << grid.getSize().z - 1 
01415        << std::endl;
01416 
01417     // Write the array contents.
01418     std::copy(grid.begin(), grid.end(), std::ostream_iterator<ValueType>(os, "\n"));
01419 
01420     return os;
01421   }
01422 
01424   template <typename ValueType>
01425   std::istream& operator >> (std::istream& is, TRegularData3D<ValueType>& grid)
01426   {
01427     typename TRegularData3D<ValueType>::CoordinateType origin;
01428     typename TRegularData3D<ValueType>::CoordinateType dimension;
01429     typename TRegularData3D<ValueType>::IndexType size;
01430 
01431     is >> origin.x >> origin.y >> origin.z;
01432     is >> dimension.x >> dimension.y >> dimension.z;
01433     is >> size.x >> size.y >> size.z;
01434 
01435     dimension -= origin;
01436     size.x++;
01437     size.y++;
01438     size.z++;
01439 
01440     grid.resize(size);
01441     grid.setOrigin(origin);
01442     grid.setDimension(dimension);
01443     std::copy(std::istream_iterator<ValueType>(is), std::istream_iterator<ValueType>(), grid.begin());
01444 
01445     return is;
01446   }
01448   
01449  } // namespace BALL
01450 
01451 #endif // BALL_DATATYPE_REGULARDATA3D_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines