BALL
1.4.1
|
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