Loading...
Searching...
No Matches
regularData3D.h
Go to the documentation of this file.
// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
#ifndef BALL_DATATYPE_REGULARDATA3D_H
#define BALL_DATATYPE_REGULARDATA3D_H
#ifndef BALL_MATHS_VECTOR3_H
# include <BALL/MATHS/vector3.h>
#endif
#ifndef BALL_SYSTEM_FILE_H
# include <BALL/SYSTEM/file.h>
#endif
#ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
# include <BALL/SYSTEM/binaryFileAdaptor.h>
#endif
#ifndef BALL_MATHS_COMMON_H
# include <BALL/MATHS/common.h>
#endif
#include <iostream>
#include <fstream>
#include <iterator>
#include <algorithm>
namespace BALL
{
template <typename ValueType>
class TRegularData3D
{
public:
BALL_CREATE(TRegularData3D<ValueType>)
class IndexType
{
public:
inline IndexType() : x(0), y(0), z(0) {}
inline IndexType(Position p) : x(p), y(p), z(p) {}
inline IndexType(Position p, Position q, Position r) : x(p), y(q), z(r) {}
Position x;
Position y;
Position z;
};
typedef std::vector<ValueType> VectorType;
typedef TVector3<float> CoordinateType;
typedef typename std::vector<ValueType>::iterator Iterator;
typedef typename std::vector<ValueType>::const_iterator ConstIterator;
// STL compatibility types
//
typedef ValueType value_type;
typedef typename std::vector<ValueType>::iterator iterator;
typedef typename std::vector<ValueType>::const_iterator const_iterator;
typedef typename std::vector<ValueType>::reference reference;
typedef typename std::vector<ValueType>::const_reference const_reference;
typedef typename std::vector<ValueType>::pointer pointer;
typedef typename std::vector<ValueType>::difference_type difference_type;
typedef typename std::vector<ValueType>::size_type size_type;
TRegularData3D();
TRegularData3D(const TRegularData3D<ValueType>& grid);
TRegularData3D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);
TRegularData3D(const CoordinateType& origin, const CoordinateType& x_axis,
const CoordinateType& y_axis, const CoordinateType& z_axis, const IndexType& size);
TRegularData3D
(const IndexType& size,
const CoordinateType& origin = CoordinateType(0.0),
const CoordinateType& dimension = CoordinateType(1.0));
virtual ~TRegularData3D();
virtual void clear();
TRegularData3D& operator = (const TRegularData3D<ValueType>& data);
bool operator == (const TRegularData3D<ValueType>& grid) const;
BALL_INLINE bool operator != (const TRegularData3D<ValueType>& grid) const { return !this->operator == (grid); }
BALL_INLINE bool empty() const { return data_.empty(); }
bool isInside(const CoordinateType& r) const;
BALL_INLINE bool isOrthogonal() const { return is_orthogonal_;}
BALL_INLINE ConstIterator begin() const { return data_.begin(); }
BALL_INLINE ConstIterator end() const { return data_.end(); }
BALL_INLINE Iterator begin() { return data_.begin(); }
BALL_INLINE Iterator end() { return data_.end(); }
// STL compatibility
BALL_INLINE size_type size() const { return data_.size(); }
BALL_INLINE size_type max_size() const { return data_.max_size(); }
BALL_INLINE void swap(TRegularData3D<ValueType>& grid) { std::swap(*this, grid); }
const vector<ValueType>& getData() const;
const ValueType& getData(const IndexType& index) const;
ValueType& getData(const IndexType& index);
const ValueType& getData(Position index) const;
ValueType& getData(Position index);
const ValueType& operator [] (const IndexType& index) const
{
return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y];
}
ValueType& operator [] (const IndexType& index)
{
return data_[index.x + size_.x * index.y + index.z * size_.x * size_.y];
}
const ValueType& operator [] (Position index) const { return data_[index]; }
ValueType& operator [] (Position index) { return data_[index]; }
ValueType operator () (const CoordinateType& x) const;
ValueType getInterpolatedValue(const CoordinateType& x) const;
const ValueType& getClosestValue(const CoordinateType& x) const;
ValueType& getClosestValue(const CoordinateType& x);
IndexType getClosestIndex(const CoordinateType& v) const;
IndexType getLowerIndex(const CoordinateType& v) const;
inline const IndexType& getSize() const { return size_; }
inline const CoordinateType& getOrigin() const { return origin_; }
const CoordinateType& getSpacing() const { return spacing_; }
void setOrigin(const CoordinateType& origin) { origin_ = origin; }
const CoordinateType& getDimension() const { return dimension_; }
void setDimension(const CoordinateType& dimension)
{
dimension_ = dimension;
if (is_orthogonal_)
{
spacing_.x = dimension_.x / (double)(size_.x - 1);
spacing_.y = dimension_.y / (double)(size_.y - 1);
spacing_.z = dimension_.z / (double)(size_.z - 1);
}
}
void resize(const IndexType& size);
void rescale(const IndexType& new_size);
CoordinateType getCoordinates(const IndexType& index) const;
CoordinateType getCoordinates(Position index) const;
void getEnclosingIndices
(const CoordinateType& r,
Position& llf, Position& rlf, Position& luf, Position& ruf,
Position& llb, Position& rlb, Position& lub, Position& rub) const;
void getEnclosingValues
(const CoordinateType& r,
ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const;
ValueType calculateMean() const;
ValueType calculateSD() const;
void binaryWrite(const String& filename) const;
void binaryWriteRaw(const String& filename) const;
void binaryRead(const String& filename);
protected:
const CoordinateType mapToCartesian_(CoordinateType r) const
{
Vector3 result;
r.x /= (size_.x - 1.);
r.y /= (size_.y - 1.);
r.z /= (size_.z - 1.);
result.x = mapping_[0] * r.x + mapping_[1] * r.y + mapping_[2] * r.z + origin_.x;
result.y = mapping_[3] * r.x + mapping_[4] * r.y + mapping_[5] * r.z + origin_.y;
result.z = mapping_[6] * r.x + mapping_[7] * r.y + mapping_[8] * r.z + origin_.z;
return result;
}
const CoordinateType mapInverse_(CoordinateType r) const
{
r -= origin_;
Vector3 result;
result.x = inverse_mapping_[0] * r.x + inverse_mapping_[1] * r.y + inverse_mapping_[2] * r.z;
result.y = inverse_mapping_[3] * r.x + inverse_mapping_[4] * r.y + inverse_mapping_[5] * r.z;
result.z = inverse_mapping_[6] * r.x + inverse_mapping_[7] * r.y + inverse_mapping_[8] * r.z;
result.x *= (size_.x - 1);
result.y *= (size_.y - 1);
result.z *= (size_.z - 1);
return result;
}
VectorType data_;
CoordinateType origin_;
CoordinateType dimension_;
CoordinateType spacing_;
IndexType size_;
typedef struct { ValueType bt[1024]; } BlockValueType;
bool is_orthogonal_;
std::vector<double> mapping_;
std::vector<double> inverse_mapping_;
};
typedef TRegularData3D<float> RegularData3D;
// default constructor.
template <class ValueType>
TRegularData3D<ValueType>::TRegularData3D()
: data_(0),
origin_(0.0),
dimension_(0.0),
spacing_(1.0),
size_(0, 0, 0),
is_orthogonal_(true)
{
}
// copy constructor
template <class ValueType>
TRegularData3D<ValueType>::TRegularData3D
(const TRegularData3D<ValueType>& data)
: data_(),
origin_(data.origin_),
dimension_(data.dimension_),
spacing_(data.spacing_),
size_(data.size_),
is_orthogonal_(data.is_orthogonal_),
mapping_(data.mapping_),
inverse_mapping_(data.inverse_mapping_)
{
try
{
data_ = data.data_;
}
catch (std::bad_alloc&)
{
data_.resize(0);
throw Exception::OutOfMemory(__FILE__, __LINE__, data.data_.size() * sizeof(ValueType));
}
}
template <class ValueType>
TRegularData3D<ValueType>::TRegularData3D
(const typename TRegularData3D<ValueType>::IndexType& size,
const typename TRegularData3D<ValueType>::CoordinateType& origin,
const typename TRegularData3D<ValueType>::CoordinateType& dimension)
: data_(),
origin_(origin),
dimension_(dimension),
spacing_(0.0),
size_(size),
is_orthogonal_(true)
{
// Compute the grid spacing
spacing_.x = dimension_.x / (double)(size_.x - 1);
spacing_.y = dimension_.y / (double)(size_.y - 1);
spacing_.z = dimension_.z / (double)(size_.z - 1);
// Compute the number of grid points
size_type number_of_points = size_.x * size_.y * size_.z;
try
{
data_.resize(number_of_points);
}
catch (std::bad_alloc&)
{
data_.resize(0);
throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_points * sizeof(ValueType));
}
}
template <class ValueType>
TRegularData3D<ValueType>::TRegularData3D
(const typename TRegularData3D<ValueType>::CoordinateType& origin,
const typename TRegularData3D<ValueType>::CoordinateType& dimension,
const typename TRegularData3D<ValueType>::CoordinateType& spacing)
: data_(),
origin_(origin),
dimension_(dimension),
spacing_(spacing),
size_(0),
is_orthogonal_(true)
{
// Compute the grid size
size_.x = (Size)(dimension_.x / spacing_.x + 0.5) + 1;
size_.y = (Size)(dimension_.y / spacing_.y + 0.5) + 1;
size_.z = (Size)(dimension_.z / spacing_.z + 0.5) + 1;
// Compute the number of grid points
size_type size = size_.x * size_.y * size_.z;
try
{
data_ .resize(size);
}
catch (std::bad_alloc&)
{
data_.resize(0);
throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
}
// Adjust the spacing -- dimension has precedence.
spacing_.x = dimension_.x / (double)(size_.x - 1);
spacing_.y = dimension_.y / (double)(size_.y - 1);
spacing_.z = dimension_.z / (double)(size_.z - 1);
}
template <class ValueType>
TRegularData3D<ValueType>::TRegularData3D
(const typename TRegularData3D<ValueType>::CoordinateType& origin,
const typename TRegularData3D<ValueType>::CoordinateType& x_axis,
const typename TRegularData3D<ValueType>::CoordinateType& y_axis,
const typename TRegularData3D<ValueType>::CoordinateType& z_axis,
const typename TRegularData3D<ValueType>::IndexType& new_size)
: data_(),
origin_(origin),
dimension_(x_axis+y_axis+z_axis),
spacing_(0.0),
size_(new_size),
is_orthogonal_(false)
{
// compute the spacing
spacing_.x = x_axis.getLength() / (new_size.x - 1.);
spacing_.y = y_axis.getLength() / (new_size.y - 1.);
spacing_.z = z_axis.getLength() / (new_size.z - 1.);
size_type size = size_.x * size_.y * size_.z;
try
{
data_.resize(size);
}
catch (std::bad_alloc&)
{
data_.resize(0);
throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
}
// prepare the mapping matrix and its inverse
mapping_.resize(9);
inverse_mapping_.resize(9);
mapping_[0] = x_axis.x; mapping_[1] = y_axis.x; mapping_[2] = z_axis.x;
mapping_[3] = x_axis.y; mapping_[4] = y_axis.y; mapping_[5] = z_axis.y;
mapping_[6] = x_axis.z; mapping_[7] = y_axis.z; mapping_[8] = z_axis.z;
// this is numerically instable and only works well for the "simple"
// cases. should be replaced by QR or an SVD
// just for readability
double a = mapping_[0]; double b = mapping_[1]; double c = mapping_[2];
double d = mapping_[3]; double e = mapping_[4]; double f = mapping_[5];
double g = mapping_[6]; double h = mapping_[7]; double i = mapping_[8];
double determinant = 1. / (a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g));
inverse_mapping_[0] = determinant * (e*i - f*h);
inverse_mapping_[1] = determinant * (b*i - c*h);
inverse_mapping_[2] = determinant * (b*f - c*e);
inverse_mapping_[3] = determinant * (f*g - d*i);
inverse_mapping_[4] = determinant * (a*i - c*g);
inverse_mapping_[5] = determinant * (c*d - a*f);
inverse_mapping_[6] = determinant * (d*h - e*g);
inverse_mapping_[7] = determinant * (b*g - a*h);
inverse_mapping_[8] = determinant * (a*e - b*d);
}
template <class ValueType>
TRegularData3D<ValueType>::~TRegularData3D()
{
}
template <typename ValueType>
BALL_INLINE
TRegularData3D<ValueType>& TRegularData3D<ValueType>::operator =
(const TRegularData3D<ValueType>& rhs)
{
// Avoid self-assignment.
if (&rhs != this)
{
// Copy the coordinate-related attributes and
// the size.
origin_ = rhs.origin_;
dimension_ = rhs.dimension_;
spacing_ = rhs.spacing_;
size_ = rhs.size_;
is_orthogonal_ = rhs.is_orthogonal_;
mapping_ = rhs.mapping_;
inverse_mapping_ = rhs.inverse_mapping_;
// Copy the data itself and rethrow allocation exceptions.
try
{
data_ = rhs.data_;
}
catch (std::bad_alloc&)
{
data_.resize(0);
throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.data_.size() * sizeof(ValueType));
}
}
return *this;
}
template <typename ValueType>
void TRegularData3D<ValueType>::resize(const typename TRegularData3D<ValueType>::IndexType& size)
{
if (is_orthogonal_)
{
// If the old size equals the new size, we're done.
if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
{
return;
}
// If the new grid is empty, this whole thing is quite easy.
if ((size.x == 0) || (size.y == 0) || (size.z == 0))
{
data_.resize(0);
dimension_.set(0.0, 0.0, 0.0);
return;
}
// Compute the new array size.
size_type new_size = (size_type)(size.x * size.y * size.z);
// Catch any bad_allocs thrown by vector::resize
try
{
// Create a new temporary array.
std::vector<ValueType> old_data(data_);
// Resize the data to its new size.
data_.resize(new_size);
// walk over the new grid and copy the old stuff back.
static ValueType default_value = (ValueType)0;
for (size_type i = 0; i < new_size; i++)
{
size_type x = i % size.x;
size_type y = (i % (size.x * size.y)) / size.x;
size_type z = i / (size.x * size.y);
if ((x >= size_.x) || (y >= size_.y) || (z >= size_.z))
{
data_[i] = default_value;
}
else
{
data_[i] = old_data[x + y * size_.x + z * size_.x * size_.y];
}
}
// Correct the grid dimension. Origin and spacing remain constant.
if ((size_.x != 0) && (size_.y != 0) && (size_.z != 0))
{
dimension_.x *= (double)size.x / (double)size_.x;
dimension_.y *= (double)size.y / (double)size_.y;
dimension_.z *= (double)size.z / (double)size_.z;
}
size_ = size;
}
catch (std::bad_alloc&)
{
throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
}
}
}
template <typename ValueType>
void TRegularData3D<ValueType>::rescale(const typename TRegularData3D<ValueType>::IndexType& size)
{
if (is_orthogonal_)
{
// If the old size equals the new size, we're done.
if ((size.x == size_.x) && (size_.y == size.y) && (size_.z == size.z))
{
return;
}
// If the new grid is empty, this whole thing is quite easy.
if ((size.x == 0) || (size.y == 0) || (size.z == 0))
{
data_.resize(0);
dimension_.set(0.0);
return;
}
// Compute the new array size.
size_type new_size = (size_type)(size.x * size.y * size.z);
// Catch any bad_allocs thrown by vector::resize
try
{
// Create a new temporary array.
TRegularData3D<ValueType> old_data(*this);
// Resize the data array to its new size.
data_.resize(new_size);
spacing_.x = dimension_.x / (double)(size.x - 1);
spacing_.y = dimension_.y / (double)(size.y - 1);
spacing_.z = dimension_.z / (double)(size.z - 1);
// Correct the grid size. Origin and dimension remain constant.
size_ = size;
// Walk over the new grid and copy the (interpolated) old stuff back.
for (size_type i = 0; i < data_.size(); i++)
{
try
{
data_[i] = old_data.getInterpolatedValue(getCoordinates(i));
}
catch (Exception::OutOfGrid&)
{
data_[i] = old_data.getClosestValue(getCoordinates(i));
}
}
}
catch (std::bad_alloc&)
{
throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * (Size)sizeof(ValueType));
}
}
}
template <class ValueType>
BALL_INLINE
bool TRegularData3D<ValueType>::isInside(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
if (is_orthogonal_)
{
return !((r.x > (origin_.x + dimension_.x ))
|| (r.y > (origin_.y + dimension_.y))
|| (r.z > (origin_.z + dimension_.z))
|| (r.x < origin_.x) || (r.y < origin_.y) || (r.z < origin_.z));
}
else
{
// compute A^-1 * pos and see whether the indices are part of the grid
CoordinateType ri = mapInverse_(r);
ri.x = Maths::round(ri.x);
ri.y = Maths::round(ri.y);
ri.z = Maths::round(ri.z);
return !( (ri.x < 0) || (ri.y < 0) || (ri.z < 0)
|| (ri.x >= size_.x) || (ri.y >= size_.y) || (ri.z >= size_.z) );
}
}
template <class ValueType>
BALL_INLINE
const ValueType& TRegularData3D<ValueType>::getData
(const typename TRegularData3D<ValueType>::IndexType& index) const
{
size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
if (pos >= data_.size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return data_[pos];
}
template <class ValueType>
BALL_INLINE
const vector<ValueType>& TRegularData3D<ValueType>::getData() const
{
return data_;
}
template <typename ValueType>
BALL_INLINE
ValueType& TRegularData3D<ValueType>::getData
(const typename TRegularData3D<ValueType>::IndexType& index)
{
size_type pos = index.x + index.y * size_.x + index.z * size_.x * size_.y;
if (pos >= data_.size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return data_[pos];
}
template <class ValueType>
BALL_INLINE
const ValueType& TRegularData3D<ValueType>::getData(Position index) const
{
if (index >= data_.size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return data_[index];
}
template <class ValueType>
BALL_INLINE
ValueType& TRegularData3D<ValueType>::getData(Position index)
{
if (index >= data_.size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return data_[index];
}
template <class ValueType>
BALL_INLINE
typename TRegularData3D<ValueType>::CoordinateType TRegularData3D<ValueType>::getCoordinates
(const typename TRegularData3D<ValueType>::IndexType& index) const
{
if ((index.x >= size_.x) || (index.y >= size_.y) || (index.z >= size_.z))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
if (is_orthogonal_)
{
CoordinateType r(origin_.x + index.x * spacing_.x,
origin_.y + index.y * spacing_.y,
origin_.z + index.z * spacing_.z);
return r;
}
else
{
CoordinateType r(index.x, index.y, index.z);
r = mapToCartesian_(r);
return r;
}
}
template <class ValueType>
BALL_INLINE
typename TRegularData3D<ValueType>::CoordinateType
TRegularData3D<ValueType>::getCoordinates(Position position) const
{
if (position >= data_.size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
Position x = (Position)(position % size_.x);
Position y = (Position)((position % (size_.x * size_.y))/ size_.x);
Position z = (Position)(position / (size_.x * size_.y));
if (is_orthogonal_)
{
return CoordinateType(origin_.x + (double)x * spacing_.x,
origin_.y + (double)y * spacing_.y,
origin_.z + (double)z * spacing_.z);
}
else
{
CoordinateType r(x,y,z);
r = mapToCartesian_(r);
return r;
}
}
template <typename ValueType>
BALL_INLINE
void TRegularData3D<ValueType>::getEnclosingIndices
(const typename TRegularData3D<ValueType>::CoordinateType& r,
Position& llf, Position& rlf, Position& luf, Position& ruf,
Position& llb, Position& rlb, Position& lub, Position& rub) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
// calculate the grid indices of the lower left front corner
// of the enclosing box
IndexType position;
if (is_orthogonal_)
{
position.x = (Position)((r.x - origin_.x) / spacing_.x);
position.y = (Position)((r.y - origin_.y) / spacing_.y);
position.z = (Position)((r.z - origin_.z) / spacing_.z);
}
else
{
CoordinateType pos = mapInverse_(r);
position.x = (Position) pos.x;
position.y = (Position) pos.y;
position.z = (Position) pos.z;
}
// calculate the (linear) indices of the eight
// box corners
llf = position.x + size_.x * position.y
+ size_.x * size_.y * position.z;
rlf = llf + 1;
luf = llf + size_.x;
ruf = luf + 1;
llb = llf + size_.x * size_.y;
rlb = llb + 1;
lub = llb + size_.x;
rub = lub + 1;
}
template <typename ValueType>
BALL_INLINE
void TRegularData3D<ValueType>::getEnclosingValues
(const typename TRegularData3D<ValueType>::CoordinateType& r,
ValueType& llf, ValueType& rlf, ValueType& luf, ValueType& ruf,
ValueType& llb, ValueType& rlb, ValueType& lub, ValueType& rub) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
// compute the eight grid indices forming the enclosing box
Position llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx;
getEnclosingIndices(r, llf_idx, rlf_idx, luf_idx, ruf_idx, llb_idx, rlb_idx, lub_idx, rub_idx);
// retrieve the grid values
llf = data_[llf_idx];
rlf = data_[rlf_idx];
luf = data_[luf_idx];
ruf = data_[ruf_idx];
llb = data_[llb_idx];
rlb = data_[rlb_idx];
lub = data_[lub_idx];
rub = data_[rub_idx];
}
template <typename ValueType>
BALL_INLINE
ValueType TRegularData3D<ValueType>::getInterpolatedValue
(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return operator () (r);
}
template <typename ValueType>
BALL_INLINE
ValueType TRegularData3D<ValueType>::operator ()
(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
Position x;
Position y;
Position z;
TVector3<double> r_0;
if (is_orthogonal_)
{
Vector3 h(r - origin_);
x = (Position)(h.x / spacing_.x);
y = (Position)(h.y / spacing_.y);
z = (Position)(h.z / spacing_.z);
while (x >= (size_.x - 1)) x--;
while (y >= (size_.y - 1)) y--;
while (z >= (size_.z - 1)) z--;
r_0.x = origin_.x + x*spacing_.x;
r_0.y = origin_.y + y*spacing_.y;
r_0.z = origin_.z + z*spacing_.z;
}
else
{
Vector3 pos = mapInverse_(r);
x = (Position) pos.x;
y = (Position) pos.y;
z = (Position) pos.z;
while (x >= (size_.x - 1)) x--;
while (y >= (size_.y - 1)) y--;
while (z >= (size_.z - 1)) z--;
// This can probably be done much faster...
Vector3 lower_pos(x,y,z);
lower_pos = mapToCartesian_(lower_pos);
r_0.x = lower_pos.x;
r_0.y = lower_pos.y;
r_0.z = lower_pos.z;
}
Position Nx = size_.x;
Position Nxy = size_.y * Nx;
Position l = x + Nx * y + Nxy * z;
double dx = 1. - ((double)(r.x - r_0.x) / spacing_.x);
double dy = 1. - ((double)(r.y - r_0.y) / spacing_.y);
double dz = 1. - ((double)(r.z - r_0.z) / spacing_.z);
return data_[l] * dx * dy * dz
+ data_[l + 1] * (1. - dx) * dy * dz
+ data_[l + Nx] * dx * (1. - dy) * dz
+ data_[l + Nx + 1] * (1. - dx) * (1. - dy) * dz
+ data_[l + Nxy] * dx * dy * (1. - dz)
+ data_[l + Nxy + 1] * (1. - dx) * dy * (1. - dz)
+ data_[l + Nxy + Nx] * dx * (1. - dy) * (1. - dz)
+ data_[l + Nxy + Nx + 1] * (1. - dx) * (1. - dy) * (1. - dz);
}
template <typename ValueType>
BALL_INLINE
typename TRegularData3D<ValueType>::IndexType TRegularData3D<ValueType>::getClosestIndex
(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
static IndexType position;
if (is_orthogonal_)
{
position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
}
else
{
Vector3 pos = mapInverse_(r);
position.x = (Position) Maths::round(pos.x);
position.y = (Position) Maths::round(pos.y);
position.z = (Position) Maths::round(pos.z);
}
return position;
}
template <typename ValueType>
BALL_INLINE
typename TRegularData3D<ValueType>::IndexType TRegularData3D<ValueType>::getLowerIndex
(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
static IndexType position;
if (is_orthogonal_)
{
position.x = (Position)((r.x - origin_.x) / spacing_.x);
position.y = (Position)((r.y - origin_.y) / spacing_.y);
position.z = (Position)((r.z - origin_.z) / spacing_.z);
}
else
{
Vector3 pos = mapInverse_(r);
position.x = (Position)pos.x;
position.y = (Position)pos.y;
position.z = (Position)pos.z;
}
return position;
}
template <typename ValueType>
BALL_INLINE
const ValueType& TRegularData3D<ValueType>::getClosestValue
(const typename TRegularData3D<ValueType>::CoordinateType& r) const
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
static IndexType position;
if (is_orthogonal_)
{
position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
}
else
{
static Vector3 pos = mapInverse_(r);
position.x = (Position) Maths::round(pos.x);
position.y = (Position) Maths::round(pos.y);
position.z = (Position) Maths::round(pos.z);
}
return operator [] (position);
}
template <typename ValueType>
BALL_INLINE
ValueType& TRegularData3D<ValueType>::getClosestValue
(const typename TRegularData3D<ValueType>::CoordinateType& r)
{
if (!isInside(r))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
static IndexType position;
if (is_orthogonal_)
{
position.x = (Position)((r.x - origin_.x) / spacing_.x + 0.5);
position.y = (Position)((r.y - origin_.y) / spacing_.y + 0.5);
position.z = (Position)((r.z - origin_.z) / spacing_.z + 0.5);
}
else
{
static Vector3 pos = mapInverse_(r);
position.x = (Position) Maths::round(pos.x);
position.y = (Position) Maths::round(pos.y);
position.z = (Position) Maths::round(pos.z);
}
return operator [] (position);
}
template <typename ValueType>
BALL_INLINE
ValueType TRegularData3D<ValueType>::calculateMean() const
{
Position data_points = (size_.x * size_.y * size_.z);
ValueType mean = 0;
for (Position i = 0; i < data_points; i++)
{
mean += data_[i];
}
mean /= data_points;
return mean;
}
template <typename ValueType>
BALL_INLINE
ValueType TRegularData3D<ValueType>::calculateSD() const
{
Position data_points = (size_.x * size_.y * size_.z);
ValueType stddev = 0;
ValueType mean = this->calculateMean();
for (Position i = 0; i < data_points; i++)
{
stddev += (pow(data_[i]-mean,2));
}
stddev /= (data_points-1);
stddev = sqrt(stddev);
return stddev;
}
template <typename ValueType>
void TRegularData3D<ValueType>::clear()
{
data_.resize(0);
origin_.set(0.0);
dimension_.set(0.0);
size_.x = 0;
size_.y = 0;
size_.z = 0;
spacing_.set(1.0);
is_orthogonal_ = true;
}
template <typename ValueType>
bool TRegularData3D<ValueType>::operator == (const TRegularData3D<ValueType>& grid) const
{
return ((origin_ == grid.origin_)
&& (dimension_ == grid.dimension_)
&& (size_.x == grid.size_.x)
&& (size_.y == grid.size_.y)
&& (size_.z == grid.size_.z)
&& (data_ == grid.data_)
&& (is_orthogonal_ == grid.is_orthogonal_));
}
template <typename ValueType>
void TRegularData3D<ValueType>::binaryWrite(const String& filename) const
{
File outfile(filename, std::ios::out|std::ios::binary);
if (!outfile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
BinaryFileAdaptor< BlockValueType > adapt_block;
BinaryFileAdaptor< ValueType > adapt_single;
// TODO: check for endiannes and swap bytes accordingly
// write all information we need to recreate the grid
BinaryFileAdaptor<CoordinateType> adapt_coordinate;
BinaryFileAdaptor<Size> adapt_size;
adapt_size.setData(data_.size());
outfile << adapt_size;
adapt_coordinate.setData(origin_);
outfile << adapt_coordinate;
adapt_coordinate.setData(dimension_);
outfile << adapt_coordinate;
adapt_coordinate.setData(spacing_);
outfile << adapt_coordinate;
BinaryFileAdaptor<IndexType> adapt_index;
adapt_index.setData(size_);
outfile << adapt_index;
// we slide a window of size 1024 over our data
Index window_pos = 0;
while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
{
adapt_block.setData(* (BlockValueType*)&(data_[window_pos]));
outfile << adapt_block;
window_pos+=1024;
}
// now we have to write the remaining data one by one
for (Size i=window_pos; i<data_.size(); i++)
{
adapt_single.setData(data_[i]);
outfile << adapt_single;
}
// that's it. I hope...
outfile.close();
}
template <typename ValueType>
void TRegularData3D<ValueType>::binaryRead(const String& filename)
{
File infile(filename, std::ios::in|std::ios::binary);
if (!infile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
BinaryFileAdaptor< BlockValueType > adapt_block;
BinaryFileAdaptor< ValueType > adapt_single;
// TODO: check for endiannes and swap bytes accordingly
// read all information we need to recreate the grid
BinaryFileAdaptor<CoordinateType> adapt_coordinate;
BinaryFileAdaptor<Size> adapt_size;
infile >> adapt_size;
Size new_size = adapt_size.getData();
infile >> adapt_coordinate;
origin_ = adapt_coordinate.getData();
infile >> adapt_coordinate;
dimension_ = adapt_coordinate.getData();
infile >> adapt_coordinate;
spacing_ = adapt_coordinate.getData();
BinaryFileAdaptor<IndexType> adapt_index;
infile >> adapt_index;
size_ = adapt_index.getData();
data_.resize(new_size);
// we slide a window of size 1024 over our data
Index window_pos = 0;
while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
{
infile >> adapt_block;
*(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
/*
for (Size i=0; i<1024; i++)
{
data_[i+window_pos] = adapt_block.getData().bt[i];
}
*/
window_pos+=1024;
}
// now we have to read the remaining data one by one
for (Size i=window_pos; i<data_.size(); i++)
{
infile >> adapt_single;
data_[i] = adapt_single.getData();
}
// that's it. I hope...
infile.close();
}
template <typename ValueType>
std::ostream& operator << (std::ostream& os, const TRegularData3D<ValueType>& grid)
{
// Write the grid origin, dimension, and number of grid points
os << grid.getOrigin().x << " " << grid.getOrigin().y << " " << grid.getOrigin().z
<< std::endl
<< grid.getOrigin().x + grid.getDimension().x << " "
<< grid.getOrigin().y + grid.getDimension().y << " "
<< grid.getOrigin().z + grid.getDimension().z
<< std::endl
<< grid.getSize().x - 1 << " " << grid.getSize().y - 1 << " " << grid.getSize().z - 1
<< std::endl;
// Write the array contents.
std::copy(grid.begin(), grid.end(), std::ostream_iterator<ValueType>(os, "\n"));
return os;
}
template <typename ValueType>
std::istream& operator >> (std::istream& is, TRegularData3D<ValueType>& grid)
{
typename TRegularData3D<ValueType>::CoordinateType origin;
typename TRegularData3D<ValueType>::CoordinateType dimension;
typename TRegularData3D<ValueType>::IndexType size;
is >> origin.x >> origin.y >> origin.z;
is >> dimension.x >> dimension.y >> dimension.z;
is >> size.x >> size.y >> size.z;
dimension -= origin;
size.x++;
size.y++;
size.z++;
grid.resize(size);
grid.setOrigin(origin);
grid.setDimension(dimension);
std::copy(std::istream_iterator<ValueType>(is), std::istream_iterator<ValueType>(), grid.begin());
return is;
}
} // namespace BALL
#endif // BALL_DATATYPE_REGULARDATA3D_H
common.h
binaryFileAdaptor.h
BALL::BinaryFileAdaptor
Definition binaryFileAdaptor.h:54
BALL::BinaryFileAdaptor::getData
const T & getData() const
Definition binaryFileAdaptor.h:134
BALL::BinaryFileAdaptor::setData
void setData(const T &data)
Definition binaryFileAdaptor.h:128
BALL::Exception::FileNotFound
Definition COMMON/exception.h:361
BALL::Exception::OutOfGrid
Definition COMMON/exception.h:351
BALL::Exception::OutOfMemory
Definition COMMON/exception.h:318
BALL::File
Definition file.h:136
BALL::File::isValid
bool isValid() const
BALL::File::close
void close()
BALL::ValueType
BALL::String
Definition string.h:57
BALL::IndexType
BALL::TRegularData3D::IndexType
The index type used to refer to a specific element in the grid (x-, y-, and z-index).
Definition regularData3D.h:57
BALL::TRegularData3D::IndexType::x
Position x
Definition regularData3D.h:64
BALL::TRegularData3D::IndexType::z
Position z
Definition regularData3D.h:68
BALL::TRegularData3D::IndexType::IndexType
IndexType(Position p, Position q, Position r)
Definition regularData3D.h:61
BALL::TRegularData3D::IndexType::IndexType
IndexType()
Definition regularData3D.h:59
BALL::TRegularData3D::IndexType::y
Position y
Definition regularData3D.h:66
BALL::TRegularData3D::IndexType::IndexType
IndexType(Position p)
Definition regularData3D.h:60
BALL::TRegularData3D
Definition regularData3D.h:46
BALL::TRegularData3D::getEnclosingValues
void getEnclosingValues(const CoordinateType &r, ValueType &llf, ValueType &rlf, ValueType &luf, ValueType &ruf, ValueType &llb, ValueType &rlb, ValueType &lub, ValueType &rub) const
Definition regularData3D.h:1017
BALL::TRegularData3D< float >::origin_
CoordinateType origin_
Definition regularData3D.h:491
BALL::TRegularData3D::begin
BALL_INLINE ConstIterator begin() const
Definition regularData3D.h:176
BALL::TRegularData3D::begin
BALL_INLINE Iterator begin()
Definition regularData3D.h:180
BALL::TRegularData3D::setOrigin
void setOrigin(const CoordinateType &origin)
Definition regularData3D.h:321
BALL::TRegularData3D< float >::inverse_mapping_
std::vector< double > inverse_mapping_
Definition regularData3D.h:510
BALL::TRegularData3D::VectorType
std::vector< ValueType > VectorType
The type containing an STL vector of the appropriate type.
Definition regularData3D.h:72
BALL::TRegularData3D::getSize
const IndexType & getSize() const
Definition regularData3D.h:306
BALL::TRegularData3D::isInside
bool isInside(const CoordinateType &r) const
Test if a given point is inside the grid.
Definition regularData3D.h:840
BALL::TRegularData3D::mapToCartesian_
const CoordinateType mapToCartesian_(CoordinateType r) const
map indices to Cartesian coordinates
Definition regularData3D.h:456
BALL::TRegularData3D::const_reference
std::vector< ValueType >::const_reference const_reference
Definition regularData3D.h:87
BALL::TRegularData3D::binaryWriteRaw
void binaryWriteRaw(const String &filename) const
BALL::TRegularData3D::size_type
std::vector< ValueType >::size_type size_type
Definition regularData3D.h:90
BALL::TRegularData3D::CoordinateType
TVector3< float > CoordinateType
The coordinate type.
Definition regularData3D.h:74
BALL::TRegularData3D::const_iterator
std::vector< ValueType >::const_iterator const_iterator
Definition regularData3D.h:85
BALL::TRegularData3D::operator!=
BALL_INLINE bool operator!=(const TRegularData3D< ValueType > &grid) const
Definition regularData3D.h:160
BALL::TRegularData3D::rescale
void rescale(const IndexType &new_size)
Definition regularData3D.h:782
BALL::TRegularData3D::getClosestValue
const ValueType & getClosestValue(const CoordinateType &x) const
BALL::TRegularData3D::empty
BALL_INLINE bool empty() const
Empty predicate.
Definition regularData3D.h:163
BALL::TRegularData3D::operator==
bool operator==(const TRegularData3D< ValueType > &grid) const
Definition regularData3D.h:1278
BALL::TRegularData3D::TRegularData3D
TRegularData3D()
Definition regularData3D.h:519
BALL::TRegularData3D::getData
ValueType & getData(Position index)
Definition regularData3D.h:908
BALL::TRegularData3D< float >::spacing_
CoordinateType spacing_
Definition regularData3D.h:497
BALL::TRegularData3D::calculateSD
ValueType calculateSD() const
Definition regularData3D.h:1248
BALL::TRegularData3D::TRegularData3D
TRegularData3D(const TRegularData3D< ValueType > &grid)
Definition regularData3D.h:532
BALL::TRegularData3D::resize
void resize(const IndexType &size)
Definition regularData3D.h:719
BALL::TRegularData3D::iterator
std::vector< ValueType >::iterator iterator
Definition regularData3D.h:84
BALL::TRegularData3D< float >::size_
IndexType size_
Definition regularData3D.h:500
BALL::TRegularData3D::mapInverse_
const CoordinateType mapInverse_(CoordinateType r) const
map Cartesian coordinates to indices (note: this does not yet convert the double values to Index)
Definition regularData3D.h:471
BALL::TRegularData3D< float >::dimension_
CoordinateType dimension_
Definition regularData3D.h:494
BALL::TRegularData3D::calculateMean
ValueType calculateMean() const
Definition regularData3D.h:1234
BALL::TRegularData3D::clear
virtual void clear()
Definition regularData3D.h:1263
BALL::TRegularData3D::TRegularData3D
TRegularData3D(const IndexType &size, const CoordinateType &origin=CoordinateType(0.0), const CoordinateType &dimension=CoordinateType(1.0))
BALL::TRegularData3D::operator=
TRegularData3D & operator=(const TRegularData3D< ValueType > &data)
Definition regularData3D.h:688
BALL::TRegularData3D::end
BALL_INLINE Iterator end()
Definition regularData3D.h:182
BALL::TRegularData3D::getClosestValue
ValueType & getClosestValue(const CoordinateType &x)
BALL::TRegularData3D< float >::is_orthogonal_
bool is_orthogonal_
Definition regularData3D.h:506
BALL::TRegularData3D::binaryRead
void binaryRead(const String &filename)
Definition regularData3D.h:1342
BALL::TRegularData3D::getDimension
const CoordinateType & getDimension() const
Definition regularData3D.h:328
BALL::TRegularData3D::getData
const vector< ValueType > & getData() const
Get the full data.
Definition regularData3D.h:877
BALL::TRegularData3D::ConstIterator
std::vector< ValueType >::const_iterator ConstIterator
A nonmutable iterator.
Definition regularData3D.h:78
BALL::TRegularData3D::max_size
BALL_INLINE size_type max_size() const
Definition regularData3D.h:191
BALL::TRegularData3D::getEnclosingIndices
void getEnclosingIndices(const CoordinateType &r, Position &llf, Position &rlf, Position &luf, Position &ruf, Position &llb, Position &rlb, Position &lub, Position &rub) const
Definition regularData3D.h:975
BALL::TRegularData3D::getClosestIndex
IndexType getClosestIndex(const CoordinateType &v) const
Definition regularData3D.h:1119
BALL::TRegularData3D::getLowerIndex
IndexType getLowerIndex(const CoordinateType &v) const
Definition regularData3D.h:1148
BALL::TRegularData3D::getCoordinates
CoordinateType getCoordinates(Position index) const
Definition regularData3D.h:946
BALL::TRegularData3D::Iterator
std::vector< ValueType >::iterator Iterator
A mutable iterator.
Definition regularData3D.h:76
BALL::TRegularData3D< float >::mapping_
std::vector< double > mapping_
Definition regularData3D.h:509
BALL::TRegularData3D::getData
const ValueType & getData(Position index) const
Definition regularData3D.h:897
BALL::TRegularData3D::difference_type
std::vector< ValueType >::difference_type difference_type
Definition regularData3D.h:89
BALL::TRegularData3D::end
BALL_INLINE ConstIterator end() const
Definition regularData3D.h:178
BALL::TRegularData3D::~TRegularData3D
virtual ~TRegularData3D()
Definition regularData3D.h:681
BALL::TRegularData3D::isOrthogonal
BALL_INLINE bool isOrthogonal() const
State whether grid is orthogonal or not.
Definition regularData3D.h:170
BALL::TRegularData3D::reference
std::vector< ValueType >::reference reference
Definition regularData3D.h:86
BALL::TRegularData3D::TRegularData3D
TRegularData3D(const CoordinateType &origin, const CoordinateType &dimension, const CoordinateType &spacing)
BALL::TRegularData3D::getOrigin
const CoordinateType & getOrigin() const
Definition regularData3D.h:312
BALL::TRegularData3D::operator[]
const ValueType & operator[](const IndexType &index) const
Definition regularData3D.h:226
BALL::TRegularData3D::getInterpolatedValue
ValueType getInterpolatedValue(const CoordinateType &x) const
Definition regularData3D.h:1044
BALL::TRegularData3D::getData
const ValueType & getData(const IndexType &index) const
BALL::TRegularData3D::getCoordinates
CoordinateType getCoordinates(const IndexType &index) const
BALL::TRegularData3D::pointer
std::vector< ValueType >::pointer pointer
Definition regularData3D.h:88
BALL::TRegularData3D::binaryWrite
void binaryWrite(const String &filename) const
Definition regularData3D.h:1290
BALL::TRegularData3D::value_type
ValueType value_type
Definition regularData3D.h:83
BALL::TRegularData3D::TRegularData3D
TRegularData3D(const CoordinateType &origin, const CoordinateType &x_axis, const CoordinateType &y_axis, const CoordinateType &z_axis, const IndexType &size)
BALL::TRegularData3D::getSpacing
const CoordinateType & getSpacing() const
Definition regularData3D.h:318
BALL::TRegularData3D::swap
BALL_INLINE void swap(TRegularData3D< ValueType > &grid)
Definition regularData3D.h:192
BALL::TRegularData3D::getData
ValueType & getData(const IndexType &index)
BALL::TRegularData3D::operator()
ValueType operator()(const CoordinateType &x) const
Definition regularData3D.h:1056
BALL::TRegularData3D::setDimension
void setDimension(const CoordinateType &dimension)
Definition regularData3D.h:338
BALL::TRegularData3D< float >::data_
VectorType data_
Definition regularData3D.h:488
BALL::TRegularData3D< float >::size
BALL_INLINE size_type size() const
Definition regularData3D.h:190
BALL::TVector3
Definition vector3.h:71
BALL::TVector3::x
T x
Definition vector3.h:489
BALL::TVector3::y
T y
Definition vector3.h:493
BALL::TVector3::z
T z
Definition vector3.h:497
double
BALL_CREATE
#define BALL_CREATE(name)
Definition create.h:62
BALL_INLINE
#define BALL_INLINE
Definition debug.h:15
file.h
BALL::operator<<
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
BALL::Constants::h
BALL_EXTERN_VARIABLE const double h
Definition constants.h:102
BALL::Constants::c
BALL_EXTERN_VARIABLE const double c
Definition constants.h:149
BALL::Maths::round
T round(const T &t)
Definition MATHS/common.h:127
BALL
Definition constants.h:13
BALL::Size
BALL_SIZE_TYPE Size
Definition COMMON/global.h:114
BALL::RegularData3D
TRegularData3D< float > RegularData3D
Definition regularData3D.h:515
BALL::Position
BALL_SIZE_TYPE Position
Definition COMMON/global.h:140
BALL::Vector3
TVector3< float > Vector3
Definition vector3.h:1084
BALL::operator>>
std::istream & operator>>(std::istream &is, TRegularData1D< ValueType > &grid)
Input operator.
Definition regularData1D.h:802
BALL::Index
BALL_INDEX_TYPE Index
Definition COMMON/global.h:105
BALL::TRegularData3D::BlockValueType
The block data type for reading and writing binary data.
Definition regularData3D.h:503
BALL::TRegularData3D::BlockValueType::bt
ValueType bt[1024]
Definition regularData3D.h:503
vector3.h