BALL 1.5.0
regularData3D.h
Go to the documentation of this file.
1// -*- Mode: C++; tab-width: 2; -*-
2// vi: set ts=2:
3//
4
5#ifndef BALL_DATATYPE_REGULARDATA3D_H
6#define BALL_DATATYPE_REGULARDATA3D_H
7
8#ifndef BALL_MATHS_VECTOR3_H
9# include <BALL/MATHS/vector3.h>
10#endif
11
12#ifndef BALL_SYSTEM_FILE_H
13# include <BALL/SYSTEM/file.h>
14#endif
15
16#ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
18#endif
19
20#ifndef BALL_MATHS_COMMON_H
21# include <BALL/MATHS/common.h>
22#endif
23
24#include <iostream>
25#include <fstream>
26#include <iterator>
27#include <algorithm>
28
29namespace BALL
30{
44 template <typename ValueType>
46 {
47 public:
48
50
51
54
55
57 {
58 public:
59 inline IndexType() : x(0), y(0), z(0) {}
60 inline IndexType(Position p) : x(p), y(p), z(p) {}
61 inline IndexType(Position p, Position q, Position r) : x(p), y(q), z(r) {}
62
69 };
70
72 typedef std::vector<ValueType> VectorType;
76 typedef typename std::vector<ValueType>::iterator Iterator;
78 typedef typename std::vector<ValueType>::const_iterator ConstIterator;
80
81 // STL compatibility types
82 //
83 typedef ValueType value_type;
84 typedef typename std::vector<ValueType>::iterator iterator;
85 typedef typename std::vector<ValueType>::const_iterator const_iterator;
86 typedef typename std::vector<ValueType>::reference reference;
87 typedef typename std::vector<ValueType>::const_reference const_reference;
88 typedef typename std::vector<ValueType>::pointer pointer;
89 typedef typename std::vector<ValueType>::difference_type difference_type;
90 typedef typename std::vector<ValueType>::size_type size_type;
91
95
100
105
109 TRegularData3D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);
110
114 TRegularData3D(const CoordinateType& origin, const CoordinateType& x_axis,
115 const CoordinateType& y_axis, const CoordinateType& z_axis, const IndexType& size);
116
121 (const IndexType& size,
122 const CoordinateType& origin = CoordinateType(0.0),
123 const CoordinateType& dimension = CoordinateType(1.0));
124
128
132 virtual void clear();
134
138
145
146
150
157
160 BALL_INLINE bool operator != (const TRegularData3D<ValueType>& grid) const { return !this->operator == (grid); }
161
163 BALL_INLINE bool empty() const { return data_.empty(); }
164
166 bool isInside(const CoordinateType& r) const;
168
171
176 BALL_INLINE ConstIterator begin() const { return data_.begin(); }
178 BALL_INLINE ConstIterator end() const { return data_.end(); }
180 BALL_INLINE Iterator begin() { return data_.begin(); }
182 BALL_INLINE Iterator end() { return data_.end(); }
184
188
189 // STL compatibility
190 BALL_INLINE size_type size() const { return data_.size(); }
191 BALL_INLINE size_type max_size() const { return data_.max_size(); }
192 BALL_INLINE void swap(TRegularData3D<ValueType>& grid) { std::swap(*this, grid); }
193
194
196 const vector<ValueType>& getData() const;
197
202 const ValueType& getData(const IndexType& index) const;
203
208 ValueType& getData(const IndexType& index);
209
214 const ValueType& getData(Position index) const;
215
220 ValueType& getData(Position index);
221
226 const ValueType& operator [] (const IndexType& index) const
227 {
228 return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y];
229 }
230
235 ValueType& operator [] (const IndexType& index)
236 {
237 return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y];
238 }
239
244 const ValueType& operator [] (Position index) const { return data_[index]; }
245
250 ValueType& operator [] (Position index) { return data_[index]; }
251
262 ValueType operator () (const CoordinateType& x) const;
263
270 ValueType getInterpolatedValue(const CoordinateType& x) const;
271
278 const ValueType& getClosestValue(const CoordinateType& x) const;
279
286 ValueType& getClosestValue(const CoordinateType& x);
287
293 IndexType getClosestIndex(const CoordinateType& v) const;
294
299 IndexType getLowerIndex(const CoordinateType& v) const;
300
306 inline const IndexType& getSize() const { return size_; }
307
312 inline const CoordinateType& getOrigin() const { return origin_; }
313
318 const CoordinateType& getSpacing() const { return spacing_; }
321 void setOrigin(const CoordinateType& origin) { origin_ = origin; }
322
328 const CoordinateType& getDimension() const { return dimension_; }
329
338 void setDimension(const CoordinateType& dimension)
339 {
340 dimension_ = dimension;
341 if (is_orthogonal_)
342 {
343 spacing_.x = dimension_.x / (double)(size_.x - 1);
344 spacing_.y = dimension_.y / (double)(size_.y - 1);
345 spacing_.z = dimension_.z / (double)(size_.z - 1);
346 }
347 }
348
364 void resize(const IndexType& size);
365
378 void rescale(const IndexType& new_size);
379
385
391
408 (const CoordinateType& r,
409 Position& llf, Position& rlf, Position& luf, Position& ruf,
410 Position& llb, Position& rlb, Position& lub, Position& rub) const;
411
417 (const CoordinateType& r,
418 ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
419 ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const;
420
424 ValueType calculateMean() const;
425
429 ValueType calculateSD() const;
430
435 void binaryWrite(const String& filename) const;
436
444 void binaryWriteRaw(const String& filename) const;
445
450 void binaryRead(const String& filename);
452
453 protected:
454
457 {
458 Vector3 result;
459 r.x /= (size_.x - 1.);
460 r.y /= (size_.y - 1.);
461 r.z /= (size_.z - 1.);
462
463 result.x = mapping_[0] * r.x + mapping_[1] * r.y + mapping_[2] * r.z + origin_.x;
464 result.y = mapping_[3] * r.x + mapping_[4] * r.y + mapping_[5] * r.z + origin_.y;
465 result.z = mapping_[6] * r.x + mapping_[7] * r.y + mapping_[8] * r.z + origin_.z;
466
467 return result;
468 }
469
472 {
473 r -= origin_;
474
475 Vector3 result;
476 result.x = inverse_mapping_[0] * r.x + inverse_mapping_[1] * r.y + inverse_mapping_[2] * r.z;
477 result.y = inverse_mapping_[3] * r.x + inverse_mapping_[4] * r.y + inverse_mapping_[5] * r.z;
478 result.z = inverse_mapping_[6] * r.x + inverse_mapping_[7] * r.y + inverse_mapping_[8] * r.z;
479
480 result.x *= (size_.x - 1);
481 result.y *= (size_.y - 1);
482 result.z *= (size_.z - 1);
483
484 return result;
485 }
486
489
492
495
498
501
503 typedef struct { ValueType bt[1024]; } BlockValueType;
504
507
509 std::vector<double> mapping_;
510 std::vector<double> inverse_mapping_;
511 };
512
516
517 // default constructor.
518 template <class ValueType>
520 : data_(0),
521 origin_(0.0),
522 dimension_(0.0),
523 spacing_(1.0),
524 size_(0, 0, 0),
525 is_orthogonal_(true)
526 {
527 }
528
529 // copy constructor
530 template <class ValueType>
532 (const TRegularData3D<ValueType>& data)
533 : data_(),
534 origin_(data.origin_),
535 dimension_(data.dimension_),
536 spacing_(data.spacing_),
537 size_(data.size_),
538 is_orthogonal_(data.is_orthogonal_),
539 mapping_(data.mapping_),
540 inverse_mapping_(data.inverse_mapping_)
541 {
542 try
543 {
544 data_ = data.data_;
545 }
546 catch (std::bad_alloc&)
547 {
548 data_.resize(0);
549 throw Exception::OutOfMemory(__FILE__, __LINE__, data.data_.size() * sizeof(ValueType));
550 }
551 }
552
553 template <class ValueType>
555 (const typename TRegularData3D<ValueType>::IndexType& size,
556 const typename TRegularData3D<ValueType>::CoordinateType& origin,
557 const typename TRegularData3D<ValueType>::CoordinateType& dimension)
558 : data_(),
559 origin_(origin),
560 dimension_(dimension),
561 spacing_(0.0),
562 size_(size),
563 is_orthogonal_(true)
564 {
565 // Compute the grid spacing
566 spacing_.x = dimension_.x / (double)(size_.x - 1);
567 spacing_.y = dimension_.y / (double)(size_.y - 1);
568 spacing_.z = dimension_.z / (double)(size_.z - 1);
569
570 // Compute the number of grid points
571 size_type number_of_points = size_.x * size_.y * size_.z;
572 try
573 {
574 data_.resize(number_of_points);
575 }
576 catch (std::bad_alloc&)
577 {
578 data_.resize(0);
579 throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_points * sizeof(ValueType));
580 }
581 }
582
583
584 template <class ValueType>
586 (const typename TRegularData3D<ValueType>::CoordinateType& origin,
587 const typename TRegularData3D<ValueType>::CoordinateType& dimension,
588 const typename TRegularData3D<ValueType>::CoordinateType& spacing)
589 : data_(),
590 origin_(origin),
591 dimension_(dimension),
592 spacing_(spacing),
593 size_(0),
594 is_orthogonal_(true)
595 {
596 // Compute the grid size
597 size_.x = (Size)(dimension_.x / spacing_.x + 0.5) + 1;
598 size_.y = (Size)(dimension_.y / spacing_.y + 0.5) + 1;
599 size_.z = (Size)(dimension_.z / spacing_.z + 0.5) + 1;
600
601 // Compute the number of grid points
603 try
604 {
605 data_ .resize(size);
606 }
607 catch (std::bad_alloc&)
608 {
609 data_.resize(0);
610 throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
611 }
612
613 // Adjust the spacing -- dimension has precedence.
614 spacing_.x = dimension_.x / (double)(size_.x - 1);
615 spacing_.y = dimension_.y / (double)(size_.y - 1);
616 spacing_.z = dimension_.z / (double)(size_.z - 1);
617 }
618
619 template <class ValueType>
621 (const typename TRegularData3D<ValueType>::CoordinateType& origin,
622 const typename TRegularData3D<ValueType>::CoordinateType& x_axis,
623 const typename TRegularData3D<ValueType>::CoordinateType& y_axis,
624 const typename TRegularData3D<ValueType>::CoordinateType& z_axis,
625 const typename TRegularData3D<ValueType>::IndexType& new_size)
626 : data_(),
627 origin_(origin),
628 dimension_(x_axis+y_axis+z_axis),
629 spacing_(0.0),
630 size_(new_size),
631 is_orthogonal_(false)
632 {
633 // compute the spacing
634 spacing_.x = x_axis.getLength() / (new_size.x - 1.);
635 spacing_.y = y_axis.getLength() / (new_size.y - 1.);
636 spacing_.z = z_axis.getLength() / (new_size.z - 1.);
637
639 try
640 {
641 data_.resize(size);
642 }
643 catch (std::bad_alloc&)
644 {
645 data_.resize(0);
646 throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
647 }
648
649 // prepare the mapping matrix and its inverse
650 mapping_.resize(9);
651 inverse_mapping_.resize(9);
652
653 mapping_[0] = x_axis.x; mapping_[1] = y_axis.x; mapping_[2] = z_axis.x;
654 mapping_[3] = x_axis.y; mapping_[4] = y_axis.y; mapping_[5] = z_axis.y;
655 mapping_[6] = x_axis.z; mapping_[7] = y_axis.z; mapping_[8] = z_axis.z;
656
657 // this is numerically instable and only works well for the "simple"
658 // cases. should be replaced by QR or an SVD
659
660 // just for readability
661 double a = mapping_[0]; double b = mapping_[1]; double c = mapping_[2];
662 double d = mapping_[3]; double e = mapping_[4]; double f = mapping_[5];
663 double g = mapping_[6]; double h = mapping_[7]; double i = mapping_[8];
664
665 double determinant = 1. / (a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g));
666
667 inverse_mapping_[0] = determinant * (e*i - f*h);
668 inverse_mapping_[1] = determinant * (b*i - c*h);
669 inverse_mapping_[2] = determinant * (b*f - c*e);
670
671 inverse_mapping_[3] = determinant * (f*g - d*i);
672 inverse_mapping_[4] = determinant * (a*i - c*g);
673 inverse_mapping_[5] = determinant * (c*d - a*f);
674
675 inverse_mapping_[6] = determinant * (d*h - e*g);
676 inverse_mapping_[7] = determinant * (b*g - a*h);
677 inverse_mapping_[8] = determinant * (a*e - b*d);
678 }
679
680 template <class ValueType>
682 {
683 }
684
685 template <typename ValueType>
688 (const TRegularData3D<ValueType>& rhs)
689 {
690 // Avoid self-assignment.
691 if (&rhs != this)
692 {
693 // Copy the coordinate-related attributes and
694 // the size.
695 origin_ = rhs.origin_;
696 dimension_ = rhs.dimension_;
697 spacing_ = rhs.spacing_;
698 size_ = rhs.size_;
699
700 is_orthogonal_ = rhs.is_orthogonal_;
701 mapping_ = rhs.mapping_;
702 inverse_mapping_ = rhs.inverse_mapping_;
703 // Copy the data itself and rethrow allocation exceptions.
704 try
705 {
706 data_ = rhs.data_;
707 }
708 catch (std::bad_alloc&)
709 {
710 data_.resize(0);
711 throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.data_.size() * sizeof(ValueType));
712 }
713 }
714
715 return *this;
716 }
717
718 template <typename ValueType>
720 {
721 if (is_orthogonal_)
722 {
723 // If the old size equals the new size, we're done.
724 if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
725 {
726 return;
727 }
728
729 // If the new grid is empty, this whole thing is quite easy.
730 if ((size.x == 0) || (size.y == 0) || (size.z == 0))
731 {
732 data_.resize(0);
733 dimension_.set(0.0, 0.0, 0.0);
734 return;
735 }
736 // Compute the new array size.
737 size_type new_size = (size_type)(size.x * size.y * size.z);
738
739 // Catch any bad_allocs thrown by vector::resize
740 try
741 {
742 // Create a new temporary array.
743 std::vector<ValueType> old_data(data_);
744
745 // Resize the data to its new size.
746 data_.resize(new_size);
747
748 // walk over the new grid and copy the old stuff back.
749 static ValueType default_value = (ValueType)0;
750 for (size_type i = 0; i < new_size; i++)
751 {
752 size_type x = i % size.x;
753 size_type y = (i % (size.x * size.y)) / size.x;
754 size_type z = i / (size.x * size.y);
755 if ((x >= size_.x) || (y >= size_.y) || (z >= size_.z))
756 {
757 data_[i] = default_value;
758 }
759 else
760 {
761 data_[i] = old_data[x + y * size_.x + z * size_.x * size_.y];
762 }
763 }
764
765 // Correct the grid dimension. Origin and spacing remain constant.
766 if ((size_.x != 0) && (size_.y != 0) && (size_.z != 0))
767 {
768 dimension_.x *= (double)size.x / (double)size_.x;
769 dimension_.y *= (double)size.y / (double)size_.y;
770 dimension_.z *= (double)size.z / (double)size_.z;
771 }
772 size_ = size;
773 }
774 catch (std::bad_alloc&)
775 {
776 throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
777 }
778 }
779 }
780
781 template <typename ValueType>
783 {
784 if (is_orthogonal_)
785 {
786 // If the old size equals the new size, we're done.
787 if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
788 {
789 return;
790 }
791
792 // If the new grid is empty, this whole thing is quite easy.
793 if ((size.x == 0) || (size.y == 0) || (size.z == 0))
794 {
795 data_.resize(0);
796 dimension_.set(0.0);
797 return;
798 }
799
800 // Compute the new array size.
801 size_type new_size = (size_type)(size.x * size.y * size.z);
802
803 // Catch any bad_allocs thrown by vector::resize
804 try
805 {
806 // Create a new temporary array.
807 TRegularData3D<ValueType> old_data(*this);
808
809 // Resize the data array to its new size.
810 data_.resize(new_size);
811 spacing_.x = dimension_.x / (double)(size.x - 1);
812 spacing_.y = dimension_.y / (double)(size.y - 1);
813 spacing_.z = dimension_.z / (double)(size.z - 1);
814
815 // Correct the grid size. Origin and dimension remain constant.
816 size_ = size;
817
818 // Walk over the new grid and copy the (interpolated) old stuff back.
819 for (size_type i = 0; i < data_.size(); i++)
820 {
821 try
822 {
823 data_[i] = old_data.getInterpolatedValue(getCoordinates(i));
824 }
825 catch (Exception::OutOfGrid&)
826 {
827 data_[i] = old_data.getClosestValue(getCoordinates(i));
828 }
829 }
830 }
831 catch (std::bad_alloc&)
832 {
833 throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
834 }
835 }
836 }
837
838 template <class ValueType>
841 {
842 if (is_orthogonal_)
843 {
844 return !((r.x > (origin_.x + dimension_.x ))
845 || (r.y > (origin_.y + dimension_.y))
846 || (r.z > (origin_.z + dimension_.z))
847 || (r.x < origin_.x) || (r.y < origin_.y) || (r.z < origin_.z));
848 }
849 else
850 {
851 // compute A^-1 * pos and see whether the indices are part of the grid
852 CoordinateType ri = mapInverse_(r);
853 ri.x = Maths::round(ri.x);
854 ri.y = Maths::round(ri.y);
855 ri.z = Maths::round(ri.z);
856
857 return !( (ri.x < 0) || (ri.y < 0) || (ri.z < 0)
858 || (ri.x >= size_.x) || (ri.y >= size_.y) || (ri.z >= size_.z) );
859 }
860 }
861
862 template <class ValueType>
865 (const typename TRegularData3D<ValueType>::IndexType& index) const
866 {
867 size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
868 if (pos >= data_.size())
869 {
870 throw Exception::OutOfGrid(__FILE__, __LINE__);
871 }
872 return data_[pos];
873 }
874
875 template <class ValueType>
877 const vector<ValueType>& TRegularData3D<ValueType>::getData() const
878 {
879 return data_;
880 }
881
882 template <typename ValueType>
885 (const typename TRegularData3D<ValueType>::IndexType& index)
886 {
887 size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
888 if (pos >= data_.size())
889 {
890 throw Exception::OutOfGrid(__FILE__, __LINE__);
891 }
892 return data_[pos];
893 }
894
895 template <class ValueType>
897 const ValueType& TRegularData3D<ValueType>::getData(Position index) const
898 {
899 if (index >= data_.size())
900 {
901 throw Exception::OutOfGrid(__FILE__, __LINE__);
902 }
903 return data_[index];
904 }
905
906 template <class ValueType>
909 {
910 if (index >= data_.size())
911 {
912 throw Exception::OutOfGrid(__FILE__, __LINE__);
913 }
914 return data_[index];
915 }
916
917 template <class ValueType>
920 (const typename TRegularData3D<ValueType>::IndexType& index) const
921 {
922 if ((index.x >= size_.x) || (index.y >= size_.y) || (index.z >= size_.z))
923 {
924 throw Exception::OutOfGrid(__FILE__, __LINE__);
925 }
926
927 if (is_orthogonal_)
928 {
929 CoordinateType r(origin_.x + index.x * spacing_.x,
930 origin_.y + index.y * spacing_.y,
931 origin_.z + index.z * spacing_.z);
932 return r;
933 }
934 else
935 {
936 CoordinateType r(index.x, index.y, index.z);
937 r = mapToCartesian_(r);
938
939 return r;
940 }
941 }
942
943 template <class ValueType>
947 {
948 if (position >= data_.size())
949 {
950 throw Exception::OutOfGrid(__FILE__, __LINE__);
951 }
952
953 Position x = (Position)(position % size_.x);
954 Position y = (Position)((position % (size_.x * size_.y))/ size_.x);
955 Position z = (Position)(position / (size_.x * size_.y));
956
957 if (is_orthogonal_)
958 {
959 return CoordinateType(origin_.x + (double)x * spacing_.x,
960 origin_.y + (double)y * spacing_.y,
961 origin_.z + (double)z * spacing_.z);
962 }
963 else
964 {
965 CoordinateType r(x,y,z);
966 r = mapToCartesian_(r);
967
968 return r;
969 }
970 }
971
972 template <typename ValueType>
976 Position& llf, Position& rlf, Position& luf, Position& ruf,
977 Position& llb, Position& rlb, Position& lub, Position& rub) const
978 {
979 if (!isInside(r))
980 {
981 throw Exception::OutOfGrid(__FILE__, __LINE__);
982 }
983
984 // calculate the grid indices of the lower left front corner
985 // of the enclosing box
986 IndexType position;
987 if (is_orthogonal_)
988 {
989 position.x = (Position)((r.x - origin_.x) / spacing_.x);
990 position.y = (Position)((r.y - origin_.y) / spacing_.y);
991 position.z = (Position)((r.z - origin_.z) / spacing_.z);
992 }
993 else
994 {
995 CoordinateType pos = mapInverse_(r);
996 position.x = (Position) pos.x;
997 position.y = (Position) pos.y;
998 position.z = (Position) pos.z;
999 }
1000
1001 // calculate the (linear) indices of the eight
1002 // box corners
1003 llf = position.x + size_.x * position.y
1004 + size_.x * size_.y * position.z;
1005 rlf = llf + 1;
1006 luf = llf + size_.x;
1007 ruf = luf + 1;
1008 llb = llf + size_.x * size_.y;
1009 rlb = llb + 1;
1010 lub = llb + size_.x;
1011 rub = lub + 1;
1012 }
1013
1014 template <typename ValueType>
1018 ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
1019 ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const
1020 {
1021 if (!isInside(r))
1022 {
1023 throw Exception::OutOfGrid(__FILE__, __LINE__);
1024 }
1025
1026 // compute the eight grid indices forming the enclosing box
1027 Position llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx;
1028 getEnclosingIndices(r, llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx);
1029
1030 // retrieve the grid values
1031 llf = data_[llf_idx];
1032 rlf = data_[rlf_idx];
1033 luf = data_[luf_idx];
1034 ruf = data_[ruf_idx];
1035 llb = data_[llb_idx];
1036 rlb = data_[rlb_idx];
1037 lub = data_[lub_idx];
1038 rub = data_[rub_idx];
1039 }
1040
1041 template <typename ValueType>
1044 (const typename TRegularData3D<ValueType>::CoordinateType& r) const
1045 {
1046 if (!isInside(r))
1047 {
1048 throw Exception::OutOfGrid(__FILE__, __LINE__);
1049 }
1050 return operator () (r);
1051 }
1052
1053 template <typename ValueType>
1056 (const typename TRegularData3D<ValueType>::CoordinateType& r) const
1057 {
1058 Position x;
1059 Position y;
1060 Position z;
1061 TVector3<double> r_0;
1062
1063 if (is_orthogonal_)
1064 {
1065 Vector3 h(r - origin_);
1066
1067 x = (Position)(h.x / spacing_.x);
1068 y = (Position)(h.y / spacing_.y);
1069 z = (Position)(h.z / spacing_.z);
1070
1071 while (x >= (size_.x - 1)) x--;
1072 while (y >= (size_.y - 1)) y--;
1073 while (z >= (size_.z - 1)) z--;
1074
1075 r_0.x = origin_.x + x*spacing_.x;
1076 r_0.y = origin_.y + y*spacing_.y;
1077 r_0.z = origin_.z + z*spacing_.z;
1078 }
1079 else
1080 {
1081 Vector3 pos = mapInverse_(r);
1082 x = (Position) pos.x;
1083 y = (Position) pos.y;
1084 z = (Position) pos.z;
1085
1086 while (x >= (size_.x - 1)) x--;
1087 while (y >= (size_.y - 1)) y--;
1088 while (z >= (size_.z - 1)) z--;
1089
1090 // This can probably be done much faster...
1091 Vector3 lower_pos(x,y,z);
1092 lower_pos = mapToCartesian_(lower_pos);
1093 r_0.x = lower_pos.x;
1094 r_0.y = lower_pos.y;
1095 r_0.z = lower_pos.z;
1096 }
1097
1098 Position Nx = size_.x;
1099 Position Nxy = size_.y * Nx;
1100 Position l = x + Nx * y + Nxy * z;
1101
1102 double dx = 1. - ((double)(r.x - r_0.x) / spacing_.x);
1103 double dy = 1. - ((double)(r.y - r_0.y) / spacing_.y);
1104 double dz = 1. - ((double)(r.z - r_0.z) / spacing_.z);
1105
1106 return data_[l] * dx * dy * dz
1107 + data_[l + 1] * (1. - dx) * dy * dz
1108 + data_[l + Nx] * dx * (1. - dy) * dz
1109 + data_[l + Nx + 1] * (1. - dx) * (1. - dy) * dz
1110 + data_[l + Nxy] * dx * dy * (1. - dz)
1111 + data_[l + Nxy + 1] * (1. - dx) * dy * (1. - dz)
1112 + data_[l + Nxy + Nx] * dx * (1. - dy) * (1. - dz)
1113 + data_[l + Nxy + Nx + 1] * (1. - dx) * (1. - dy) * (1. - dz);
1114 }
1115
1116 template <typename ValueType>
1119 (const typename TRegularData3D<ValueType>::CoordinateType& r) const
1120 {
1121 if (!isInside(r))
1122 {
1123 throw Exception::OutOfGrid(__FILE__, __LINE__);
1124 }
1125
1126 static IndexType position;
1127
1128 if (is_orthogonal_)
1129 {
1130 position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
1131 position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
1132 position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
1133 }
1134 else
1135 {
1136 Vector3 pos = mapInverse_(r);
1137 position.x = (Position) Maths::round(pos.x);
1138 position.y = (Position) Maths::round(pos.y);
1139 position.z = (Position) Maths::round(pos.z);
1140 }
1141
1142 return position;
1143 }
1144
1145 template <typename ValueType>
1148 (const typename TRegularData3D<ValueType>::CoordinateType& r) const
1149 {
1150 if (!isInside(r))
1151 {
1152 throw Exception::OutOfGrid(__FILE__, __LINE__);
1153 }
1154
1155 static IndexType position;
1156 if (is_orthogonal_)
1157 {
1158 position.x = (Position)((r.x - origin_.x) / spacing_.x);
1159 position.y = (Position)((r.y - origin_.y) / spacing_.y);
1160 position.z = (Position)((r.z - origin_.z) / spacing_.z);
1161 }
1162 else
1163 {
1164 Vector3 pos = mapInverse_(r);
1165 position.x = (Position)pos.x;
1166 position.y = (Position)pos.y;
1167 position.z = (Position)pos.z;
1168 }
1169
1170 return position;
1171 }
1172
1173 template <typename ValueType>
1176 (const typename TRegularData3D<ValueType>::CoordinateType& r) const
1177 {
1178 if (!isInside(r))
1179 {
1180 throw Exception::OutOfGrid(__FILE__, __LINE__);
1181 }
1182
1183 static IndexType position;
1184
1185 if (is_orthogonal_)
1186 {
1187 position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
1188 position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
1189 position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
1190 }
1191 else
1192 {
1193 static Vector3 pos = mapInverse_(r);
1194 position.x = (Position) Maths::round(pos.x);
1195 position.y = (Position) Maths::round(pos.y);
1196 position.z = (Position) Maths::round(pos.z);
1197 }
1198
1199 return operator [] (position);
1200 }
1201
1202 template <typename ValueType>
1206 {
1207 if (!isInside(r))
1208 {
1209 throw Exception::OutOfGrid(__FILE__, __LINE__);
1210 }
1211
1212 static IndexType position;
1213
1214 if (is_orthogonal_)
1215 {
1216 position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
1217 position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
1218 position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
1219 }
1220 else
1221 {
1222 static Vector3 pos = mapInverse_(r);
1223 position.x = (Position) Maths::round(pos.x);
1224 position.y = (Position) Maths::round(pos.y);
1225 position.z = (Position) Maths::round(pos.z);
1226 }
1227
1228
1229 return operator [] (position);
1230 }
1231
1232 template <typename ValueType>
1235 {
1236 Position data_points = (size_.x * size_.y * size_.z);
1237 ValueType mean = 0;
1238 for (Position i = 0; i < data_points; i++)
1239 {
1240 mean += data_[i];
1241 }
1242 mean /= data_points;
1243 return mean;
1244 }
1245
1246 template <typename ValueType>
1249 {
1250 Position data_points = (size_.x * size_.y * size_.z);
1251 ValueType stddev = 0;
1252 ValueType mean = this->calculateMean();
1253 for (Position i = 0; i < data_points; i++)
1254 {
1255 stddev += (pow(data_[i]-mean,2));
1256 }
1257 stddev /= (data_points-1);
1258 stddev = sqrt(stddev);
1259 return stddev;
1260 }
1261
1262 template <typename ValueType>
1264 {
1265 data_.resize(0);
1266
1267 origin_.set(0.0);
1268 dimension_.set(0.0);
1269 size_.x = 0;
1270 size_.y = 0;
1271 size_.z = 0;
1272 spacing_.set(1.0);
1273 is_orthogonal_ = true;
1274 }
1275
1276
1277 template <typename ValueType>
1279 {
1280 return ((origin_ == grid.origin_)
1281 && (dimension_ == grid.dimension_)
1282 && (size_.x == grid.size_.x)
1283 && (size_.y == grid.size_.y)
1284 && (size_.z == grid.size_.z)
1285 && (data_ == grid.data_)
1286 && (is_orthogonal_ == grid.is_orthogonal_));
1287 }
1288
1289 template <typename ValueType>
1291 {
1292 File outfile(filename, std::ios::out|std::ios::binary);
1293 if (!outfile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
1294
1296 BinaryFileAdaptor< ValueType > adapt_single;
1297
1298 // TODO: check for endiannes and swap bytes accordingly
1299
1300 // write all information we need to recreate the grid
1301 BinaryFileAdaptor<CoordinateType> adapt_coordinate;
1302 BinaryFileAdaptor<Size> adapt_size;
1303
1304 adapt_size.setData(data_.size());
1305 outfile << adapt_size;
1306
1307 adapt_coordinate.setData(origin_);
1308 outfile << adapt_coordinate;
1309
1310 adapt_coordinate.setData(dimension_);
1311 outfile << adapt_coordinate;
1312
1313 adapt_coordinate.setData(spacing_);
1314 outfile << adapt_coordinate;
1315
1316 BinaryFileAdaptor<IndexType> adapt_index;
1317 adapt_index.setData(size_);
1318 outfile << adapt_index;
1319
1320 // we slide a window of size 1024 over our data
1321 Index window_pos = 0;
1322
1323 while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
1324 {
1325 adapt_block.setData(* (BlockValueType*)&(data_[window_pos]));
1326 outfile << adapt_block;
1327 window_pos+=1024;
1328 }
1329
1330 // now we have to write the remaining data one by one
1331 for (Size i=window_pos; i<data_.size(); i++)
1332 {
1333 adapt_single.setData(data_[i]);
1334 outfile << adapt_single;
1335 }
1336
1337 // that's it. I hope...
1338 outfile.close();
1339 }
1340
1341 template <typename ValueType>
1343 {
1344 File infile(filename, std::ios::in|std::ios::binary);
1345 if (!infile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
1346
1348 BinaryFileAdaptor< ValueType > adapt_single;
1349
1350 // TODO: check for endiannes and swap bytes accordingly
1351
1352 // read all information we need to recreate the grid
1353 BinaryFileAdaptor<CoordinateType> adapt_coordinate;
1354 BinaryFileAdaptor<Size> adapt_size;
1355
1356 infile >> adapt_size;
1357 Size new_size = adapt_size.getData();
1358
1359 infile >> adapt_coordinate;
1360 origin_ = adapt_coordinate.getData();
1361
1362 infile >> adapt_coordinate;
1363 dimension_ = adapt_coordinate.getData();
1364
1365 infile >> adapt_coordinate;
1366 spacing_ = adapt_coordinate.getData();
1367
1368 BinaryFileAdaptor<IndexType> adapt_index;
1369 infile >> adapt_index;
1370 size_ = adapt_index.getData();
1371
1372 data_.resize(new_size);
1373
1374 // we slide a window of size 1024 over our data
1375 Index window_pos = 0;
1376
1377 while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
1378 {
1379 infile >> adapt_block;
1380 *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
1381 /*
1382 for (Size i=0; i<1024; i++)
1383 {
1384 data_[i+window_pos] = adapt_block.getData().bt[i];
1385 }
1386 */
1387 window_pos+=1024;
1388 }
1389
1390 // now we have to read the remaining data one by one
1391 for (Size i=window_pos; i<data_.size(); i++)
1392 {
1393 infile >> adapt_single;
1394 data_[i] = adapt_single.getData();
1395 }
1396
1397 // that's it. I hope...
1398 infile.close();
1399 }
1400
1404 template <typename ValueType>
1405 std::ostream& operator << (std::ostream& os, const TRegularData3D<ValueType>& grid)
1406 {
1407 // Write the grid origin, dimension, and number of grid points
1408 os << grid.getOrigin().x << " " << grid.getOrigin().y << " " << grid.getOrigin().z
1409 << std::endl
1410 << grid.getOrigin().x + grid.getDimension().x << " "
1411 << grid.getOrigin().y + grid.getDimension().y << " "
1412 << grid.getOrigin().z + grid.getDimension().z
1413 << std::endl
1414 << grid.getSize().x - 1 << " " << grid.getSize().y - 1 << " " << grid.getSize().z - 1
1415 << std::endl;
1416
1417 // Write the array contents.
1418 std::copy(grid.begin(), grid.end(), std::ostream_iterator<ValueType>(os, "\n"));
1419
1420 return os;
1421 }
1422
1424 template <typename ValueType>
1425 std::istream& operator >> (std::istream& is, TRegularData3D<ValueType>& grid)
1426 {
1430
1431 is >> origin.x >> origin.y >> origin.z;
1432 is >> dimension.x >> dimension.y >> dimension.z;
1433 is >> size.x >> size.y >> size.z;
1434
1435 dimension -= origin;
1436 size.x++;
1437 size.y++;
1438 size.z++;
1439
1440 grid.resize(size);
1441 grid.setOrigin(origin);
1442 grid.setDimension(dimension);
1443 std::copy(std::istream_iterator<ValueType>(is), std::istream_iterator<ValueType>(), grid.begin());
1444
1445 return is;
1446 }
1448
1449 } // namespace BALL
1450
1451#endif // BALL_DATATYPE_REGULARDATA3D_H
#define BALL_INLINE
Definition: debug.h:15
#define BALL_CREATE(name)
Definition: create.h:62
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
Definition: constants.h:13
BALL_SIZE_TYPE Size
TRegularData3D< float > RegularData3D
BALL_SIZE_TYPE Position
TVector3< float > Vector3
Definition: vector3.h:1084
std::istream & operator>>(std::istream &is, TRegularData1D< ValueType > &grid)
Input operator.
BALL_EXTERN_VARIABLE const double h
Definition: constants.h:102
BALL_EXTERN_VARIABLE const double c
Definition: constants.h:149
std::ifstream infile
Definition: classTest.h:27
T round(const T &t)
Definition: MATHS/common.h:127
void getEnclosingValues(const CoordinateType &r, ValueType &llf, ValueType &rlf, ValueType &luf, ValueType &ruf, ValueType &llb, ValueType &rlb, ValueType &lub, ValueType &rub) const
CoordinateType origin_
Origin of the grid (offset)
BALL_INLINE ConstIterator begin() const
BALL_INLINE Iterator begin()
void setOrigin(const CoordinateType &origin)
std::vector< double > inverse_mapping_
std::vector< ValueType > VectorType
The type containing an STL vector of the appropriate type.
Definition: regularData3D.h:72
const IndexType & getSize() const
bool isInside(const CoordinateType &r) const
Test if a given point is inside the grid.
const CoordinateType mapToCartesian_(CoordinateType r) const
map indices to Cartesian coordinates
std::vector< ValueType >::const_reference const_reference
Definition: regularData3D.h:87
void binaryWriteRaw(const String &filename) const
std::vector< ValueType >::size_type size_type
Definition: regularData3D.h:90
TVector3< float > CoordinateType
The coordinate type.
Definition: regularData3D.h:74
std::vector< ValueType >::const_iterator const_iterator
Definition: regularData3D.h:85
BALL_INLINE bool operator!=(const TRegularData3D< ValueType > &grid) const
void rescale(const IndexType &new_size)
const ValueType & getClosestValue(const CoordinateType &x) const
BALL_INLINE bool empty() const
Empty predicate.
bool operator==(const TRegularData3D< ValueType > &grid) const
ValueType & getData(Position index)
CoordinateType spacing_
Grid spacing.
ValueType calculateSD() const
TRegularData3D(const TRegularData3D< ValueType > &grid)
void resize(const IndexType &size)
std::vector< ValueType >::iterator iterator
Definition: regularData3D.h:84
IndexType size_
The dimensions in grid points.
const CoordinateType mapInverse_(CoordinateType r) const
map Cartesian coordinates to indices (note: this does not yet convert the double values to Index)
CoordinateType dimension_
Dimension of the grid.
ValueType calculateMean() const
virtual void clear()
TRegularData3D(const IndexType &size, const CoordinateType &origin=CoordinateType(0.0), const CoordinateType &dimension=CoordinateType(1.0))
TRegularData3D & operator=(const TRegularData3D< ValueType > &data)
BALL_INLINE Iterator end()
ValueType & getClosestValue(const CoordinateType &x)
bool is_orthogonal_
A flag deciding whether the grid is orthogonal or not.
void binaryRead(const String &filename)
const CoordinateType & getDimension() const
const vector< ValueType > & getData() const
Get the full data.
std::vector< ValueType >::const_iterator ConstIterator
A nonmutable iterator.
Definition: regularData3D.h:78
BALL_INLINE size_type max_size() const
void getEnclosingIndices(const CoordinateType &r, Position &llf, Position &rlf, Position &luf, Position &ruf, Position &llb, Position &rlb, Position &lub, Position &rub) const
IndexType getClosestIndex(const CoordinateType &v) const
IndexType getLowerIndex(const CoordinateType &v) const
CoordinateType getCoordinates(Position index) const
std::vector< ValueType >::iterator Iterator
A mutable iterator.
Definition: regularData3D.h:76
std::vector< double > mapping_
mapping matrix and its inverse
const ValueType & getData(Position index) const
std::vector< ValueType >::difference_type difference_type
Definition: regularData3D.h:89
BALL_INLINE ConstIterator end() const
BALL_INLINE bool isOrthogonal() const
State whether grid is orthogonal or not.
std::vector< ValueType >::reference reference
Definition: regularData3D.h:86
TRegularData3D(const CoordinateType &origin, const CoordinateType &dimension, const CoordinateType &spacing)
const CoordinateType & getOrigin() const
const ValueType & operator[](const IndexType &index) const
ValueType getInterpolatedValue(const CoordinateType &x) const
const ValueType & getData(const IndexType &index) const
CoordinateType getCoordinates(const IndexType &index) const
std::vector< ValueType >::pointer pointer
Definition: regularData3D.h:88
void binaryWrite(const String &filename) const
TRegularData3D(const CoordinateType &origin, const CoordinateType &x_axis, const CoordinateType &y_axis, const CoordinateType &z_axis, const IndexType &size)
const CoordinateType & getSpacing() const
BALL_INLINE void swap(TRegularData3D< ValueType > &grid)
ValueType & getData(const IndexType &index)
ValueType operator()(const CoordinateType &x) const
void setDimension(const CoordinateType &dimension)
VectorType data_
The grid data.
BALL_INLINE size_type size() const
The index type used to refer to a specific element in the grid (x-, y-, and z-index)
Definition: regularData3D.h:57
IndexType(Position p, Position q, Position r)
Definition: regularData3D.h:61
The block data type for reading and writing binary data.
const T & getData() const
void setData(const T &data)
bool isValid() const
void close()