Loading...
Searching...
No Matches
FFT3D.h
Go to the documentation of this file.
// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
#ifndef BALL_MATHS_TFFT3D_H
#define BALL_MATHS_TFFT3D_H
#ifndef BALL_COMMON_EXCEPTION_H
# include <BALL/COMMON/exception.h>
#endif
#ifndef BALL_DATATYPE_REGULARDATA3D_H
# include <BALL/DATATYPE/regularData3D.h>
#endif
//#ifndef BALL_MATHS_VECTOR2_H
//# include <BALL/MATHS/vector3.h>
//#endif
#include <BALL/MATHS/fftwCommon.h>
#include <cmath>
#include <complex>
#include <fftw3.h>
namespace BALL
{
template <typename ComplexTraits>
class TFFT3D
: public TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> >
{
public:
typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
typedef TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
BALL_CREATE(TFFT3D)
TFFT3D();
TFFT3D(const TFFT3D &data);
// AR: ldn is not any longer the binary logarithm but the absolute number of grid points
TFFT3D(Size ldnX, Size ldnY, Size ldnZ, double stepPhysX=1., double stepPhysY=1., double stepPhysZ=1., Vector3 origin=Vector3(0.,0.,0), bool inFourierSpace=false);
virtual ~TFFT3D();
const TFFT3D& operator = (const TFFT3D& fft_3d);
virtual void clear();
virtual void destroy();
bool operator == (const TFFT3D& fft3d) const;
// @name Accessors
void doFFT();
void doiFFT();
bool translate(const Vector3& trans_origin);
bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z);
double getPhysStepWidthX() const;
double getPhysStepWidthY() const;
double getPhysStepWidthZ() const;
double getFourierStepWidthX() const;
double getFourierStepWidthY() const;
double getFourierStepWidthZ() const;
double getPhysSpaceMinX() const;
double getPhysSpaceMinY() const
;
double getPhysSpaceMinZ() const;
double getPhysSpaceMaxX() const;
double getPhysSpaceMaxY() const;
double getPhysSpaceMaxZ() const;
double getFourierSpaceMinX() const;
double getFourierSpaceMinY() const;
double getFourierSpaceMinZ() const;
double getFourierSpaceMaxX() const;
double getFourierSpaceMaxY() const;
double getFourierSpaceMaxZ() const;
Size getMaxXIndex() const;
Size getMaxYIndex() const;
Size getMaxZIndex() const;
Size getNumberOfInverseTransforms() const;
Vector3 getGridCoordinates(Position position) const;
Complex getData(const Vector3& pos) const;
Complex getInterpolatedValue(const Vector3& pos) const;
void setData(const Vector3& pos, Complex val);
Complex& operator[](const Vector3& pos);
const Complex& operator[](const Vector3& pos) const;
Complex& operator[](const Position& pos)
{
return TRegularData3D<Complex>::operator [] (pos);
}
const Complex& operator[](const Position& pos) const
{
return TRegularData3D<Complex>::operator [] (pos);
}
// AR:
void setNumberOfFFTTransforms(Size num)
{
numPhysToFourier_ = num;
}
// AR:
void setNumberOfiFFTTransforms(Size num)
{
numFourierToPhys_ = num;
}
Complex phase(const Vector3& pos) const;
bool isInFourierSpace() const;
protected:
Size lengthX_, lengthY_, lengthZ_;
bool inFourierSpace_;
Size numPhysToFourier_;
Size numFourierToPhys_;
Vector3 origin_;
double stepPhysX_, stepPhysY_, stepPhysZ_;
double stepFourierX_, stepFourierY_, stepFourierZ_;
Vector3 minPhys_, maxPhys_;
Vector3 minFourier_, maxFourier_;
// AR: new version for FFTW3
typename ComplexTraits::FftwPlan planForward_;
typename ComplexTraits::FftwPlan planBackward_;
// AR: to control plan calculation with new fftw3
Size dataLength_;
Complex *dataAdress_;
bool planCalculated_;
};
typedef TFFT3D<BALL_FFTW_DEFAULT_TRAITS> FFT3D;
template <typename ComplexTraits>
const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator <<
(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from);
template <typename ComplexTraits>
const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from);
template <typename ComplexTraits>
TFFT3D<ComplexTraits>::TFFT3D()
: TRegularData3D<Complex>(),
dataLength_(0),
dataAdress_(0),
planCalculated_(false)
{
}
template <typename ComplexTraits>
bool TFFT3D<ComplexTraits>::operator == (const TFFT3D& fft3D) const
{
// AR: test whether data_.size() == fft3D.data_.size()
// instead of testing 3 lengths. Better for vector handling.
if (lengthX_ == fft3D.lengthX_ &&
lengthY_ == fft3D.lengthY_ &&
lengthZ_ == fft3D.lengthZ_ &&
origin_ == fft3D.origin_ &&
stepPhysX_ == fft3D.stepPhysX_ &&
stepPhysY_ == fft3D.stepPhysY_ &&
stepPhysZ_ == fft3D.stepPhysZ_ &&
stepFourierX_ == fft3D.stepFourierX_ &&
stepFourierY_ == fft3D.stepFourierY_ &&
stepFourierZ_ == fft3D.stepFourierZ_ &&
minPhys_ == fft3D.minPhys_ &&
maxPhys_ == fft3D.maxPhys_ &&
minFourier_ == fft3D.minFourier_ &&
maxFourier_ == fft3D.maxFourier_ &&
numPhysToFourier_ == fft3D.numPhysToFourier_ &&
numFourierToPhys_ == fft3D.numFourierToPhys_ &&
planCalculated_ == fft3D.planCalculated_)
{
Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_;
Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;
for (double posX=min.x; posX<=max.x; posX+=stepX)
{
for (double posY=min.y; posY<=max.y; posY+=stepY)
{
for (double posZ=min.z; posZ<=max.z; posZ+=stepZ)
{
if (getData(Vector3(posX,posY,posZ)) != fft3D.getData(Vector3(posX,posY,posZ)))
{
return false;
}
}
}
}
return true;
}
return false;
}
template <typename ComplexTraits>
bool TFFT3D<ComplexTraits>::translate(const Vector3& trans_origin)
{
Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_);
Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_);
Position internalOriginZ = (Position) Maths::rint(trans_origin.z*stepPhysZ_);
if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_) && (internalOriginZ <= lengthZ_))
{
origin_.x = trans_origin.x;
origin_.y = trans_origin.y;
origin_.z = trans_origin.z;
minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
return true;
}
else
{
return false;
}
}
template <typename ComplexTraits>
bool TFFT3D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
{
if ((new_width_x <= 0) || (new_width_y <= 0) || (new_width_z <= 0))
{
return false;
}
else
{
stepPhysX_ = new_width_x;
stepPhysY_ = new_width_y;
stepPhysZ_ = new_width_z;
stepFourierX_ = 2.*Constants::PI/(stepPhysX_*lengthX_);
stepFourierY_ = 2.*Constants::PI/(stepPhysY_*lengthY_);
stepFourierZ_ = 2.*Constants::PI/(stepPhysZ_*lengthZ_);
minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
return true;
}
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysStepWidthX() const
{
return stepPhysX_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysStepWidthY() const
{
return stepPhysY_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysStepWidthZ() const
{
return stepPhysZ_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierStepWidthX() const
{
return stepFourierX_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierStepWidthY() const
{
return stepFourierY_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierStepWidthZ() const
{
return stepFourierZ_;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMinX() const
{
return minPhys_.x;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMinY() const
{
return minPhys_.y;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMinZ() const
{
return minPhys_.z;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMaxX() const
{
return maxPhys_.x;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMaxY() const
{
return maxPhys_.y;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getPhysSpaceMaxZ() const
{
return maxPhys_.z;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMinX() const
{
return minFourier_.x;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMinY() const
{
return minFourier_.y;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMinZ() const
{
return minFourier_.z;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMaxX() const
{
return maxFourier_.x;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMaxY() const
{
return maxFourier_.y;
}
template <typename ComplexTraits>
double TFFT3D<ComplexTraits>::getFourierSpaceMaxZ() const
{
return maxFourier_.z;
}
template <typename ComplexTraits>
Size TFFT3D<ComplexTraits>::getMaxXIndex() const
{
return (lengthX_ - 1);
}
template <typename ComplexTraits>
Size TFFT3D<ComplexTraits>::getMaxYIndex() const
{
return (lengthY_ - 1);
}
template <typename ComplexTraits>
Size TFFT3D<ComplexTraits>::getMaxZIndex() const
{
return (lengthZ_ - 1);
}
template <typename ComplexTraits>
Size TFFT3D<ComplexTraits>::getNumberOfInverseTransforms() const
{
return numFourierToPhys_;
}
template <typename ComplexTraits>
Vector3 TFFT3D<ComplexTraits>::getGridCoordinates(Position position) const
{
if (!inFourierSpace_)
{
if (position >= ComplexVector::size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
Vector3 r;
Position x, y, z;
z = position % lengthZ_;
y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
x = position / (lengthY_ * lengthZ_);
r.set(-origin_.x + (float)x * stepPhysX_,
-origin_.y + (float)y * stepPhysY_,
-origin_.z + (float)z * stepPhysZ_);
return r;
}
else
{
if (position >= ComplexVector::size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
Vector3 r;
Index x, y, z;
z = position % lengthZ_;
y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
x = position / (lengthY_ * lengthZ_);
if (x>=lengthX_/2.)
{
x-=lengthX_;
}
if (y>=lengthY_/2.)
{
y-=lengthY_;
}
if (z>=lengthZ_/2.)
{
z-=lengthZ_;
}
r.set((float)x * stepFourierX_,
(float)y * stepFourierY_,
(float)z * stepFourierZ_);
return r;
}
}
template <typename ComplexTraits>
typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getData(const Vector3& pos) const
{
Complex result;
double normalization=1.;
if (!inFourierSpace_)
{
result = (*this)[pos];
normalization=1./((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
}
else
{
// AR:
//old: result = (*this)[pos];
result = (*this)[pos]*phase(pos);
//normalization=1./pow(sqrt(2.*Constants::PI),3)*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
normalization=1./pow(sqrt(2.*Constants::PI),3)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
//normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
}
result *= normalization;
return result;
}
template <typename ComplexTraits>
typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getInterpolatedValue(const Vector3& pos) const
{
Complex result;
Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_;
Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;
if ( (pos.x < min.x) || (pos.y < min.y) || (pos.z < min.z)
|| (pos.x > max.x) || (pos.y > max.y) || (pos.z > max.z) )
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
Vector3 h(pos.x - min.x, pos.y - min.y, pos.z - min.z);
double modX = fmod((double)h.x,stepX);
double modY = fmod((double)h.y,stepY);
double modZ = fmod((double)h.z,stepZ);
if (modX==0 && modY==0 && modZ==0) // we are on the grid
{
return getData(pos);
}
double beforeX = floor(h.x/stepX)*stepX+min.x;
double beforeY = floor(h.y/stepY)*stepY+min.y;
double beforeZ = floor(h.z/stepZ)*stepZ+min.z;
double afterX = ceil(h.x/stepX)*stepX+min.x;
double afterY = ceil(h.y/stepY)*stepY+min.y;
double afterZ = ceil(h.z/stepZ)*stepZ+min.z;
double tx = (pos.x - beforeX)/stepX;
double ty = (pos.y - beforeY)/stepY;
double tz = (pos.z - beforeZ)/stepZ;
result = getData(Vector3(beforeX,beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*(1.-tz));
result += getData(Vector3(afterX, beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)*(1.-tz));
result += getData(Vector3(beforeX,afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty *(1.-tz));
result += getData(Vector3(beforeX,beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)* tz );
result += getData(Vector3(afterX, afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx * ty *(1.-tz));
result += getData(Vector3(afterX, beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)* tz );
result += getData(Vector3(beforeX,afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty * tz );
result += getData(Vector3(afterX, afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx * ty * tz );
return result;
}
template <typename ComplexTraits>
void TFFT3D<ComplexTraits>::setData(const Vector3& pos, Complex val)
{
Complex dummy;
if (!inFourierSpace_)
{
dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)),
val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)));
(*this)[pos]=dummy;
}
else
{
val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((pow(sqrt(2*Constants::PI),3)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
/*val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*Constants::PI)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));*/
dummy = val;
(*this)[pos]=dummy;
}
}
template <typename ComplexTraits>
typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos)
{
Index internalPos;
if (!inFourierSpace_)
{
Index i, j, k;
i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
internalPos = (k + (j + i*lengthY_)*lengthZ_);
/*(Index) rint( (pos.z+origin_.z)/stepPhysZ_
+ ( (pos.y+origin_.y)/stepPhysY_
+ (pos.x+origin_.x)/stepPhysX_*lengthY_
) * lengthZ_
); */
}
else
{
Index i, j, k;
i = (Index) Maths::rint(pos.x/stepFourierX_);
j = (Index) Maths::rint(pos.y/stepFourierY_);
k = (Index) Maths::rint(pos.z/stepFourierZ_);
if (i<0)
{
i+=lengthX_;
}
if (j<0)
{
j+=lengthY_;
}
if (k<0)
{
k+=lengthZ_;
}
internalPos = (k + (j + i*lengthY_)*lengthZ_);
}
if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return TRegularData3D<Complex>::operator[]((Position)internalPos);
}
template <typename ComplexTraits>
const typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos) const
{
Index internalPos;
if (!inFourierSpace_)
{
Index i, j, k;
i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
internalPos = (k + (j + i*lengthY_)*lengthZ_);
/*(Index) rint( (pos.z+origin_.z)/stepPhysZ_
+ ( (pos.y+origin_.y)/stepPhysY_
+ (pos.x+origin_.x)/stepPhysX_*lengthY_
) * lengthZ_
); */
}
else
{
Index i, j, k;
i = (Index) Maths::rint(pos.x/stepFourierX_);
j = (Index) Maths::rint(pos.y/stepFourierY_);
k = (Index) Maths::rint(pos.z/stepFourierZ_);
if (i<0)
{
i+=lengthX_;
}
if (j<0)
{
j+=lengthY_;
}
if (k<0)
{
k+=lengthZ_;
}
internalPos = (k + (j + i*lengthY_)*lengthZ_);
}
if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return TRegularData3D<Complex>::operator[]((Position)internalPos);
}
/*Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos)
{
return operator [] (pos);
}
const Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos) const
{
return operator [] (pos);
}
*/
template <typename ComplexTraits>
typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::phase(const Vector3& pos) const
{
// AR: old version: -2.*Constants::PI...
double phase = 2.*Constants::PI*( (Maths::rint(pos.x/stepFourierX_))*(Maths::rint(origin_.x/stepPhysX_))
/lengthX_
+ (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
/lengthY_
+ (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
/lengthZ_ );
Complex result = Complex(cos(phase), sin(phase));
return result;
/*double phase = -2.*Constants::PI*( (rint(pos.x/stepFourierX_))*(rint(origin_.x/stepPhysX_))
/lengthX_
+ (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
/lengthY_
+ (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
/lengthZ_ );
Complex result = Complex(cos(phase), sin(phase));
return result;*/
}
template <typename ComplexTraits>
bool TFFT3D<ComplexTraits>::isInFourierSpace() const
{
return inFourierSpace_;
}
template <typename ComplexTraits>
const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator<<
(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from)
{
// first decide if the TFFT3D data is in Fourier space.
if (!from.isInFourierSpace())
{
// create a new grid
Size lengthX = from.getMaxXIndex()+1;
Size lengthY = from.getMaxYIndex()+1;
Size lengthZ = from.getMaxZIndex()+1;
TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
// and fill it
double normalization=1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
typename TFFT3D<ComplexTraits>::Complex dataIn;
typename TFFT3D<ComplexTraits>::Complex dataOut;
for (Position i = 0; i < from.size(); i++)
{
Position x, y, z;
z = i % lengthZ;
y = (i % (lengthY * lengthZ)) / lengthZ;
x = i / (lengthY * lengthZ);
dataIn = from[i];
dataOut = dataIn;
newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization;
}
to = newGrid;
return to;
}
else
{
// we are in Fourier space
// create a new grid
Size lengthX = from.getMaxXIndex()+1;
Size lengthY = from.getMaxYIndex()+1;
Size lengthZ = from.getMaxZIndex()+1;
//float stepPhysX = from.getPhysStepWidthX();
//float stepPhysY = from.getPhysStepWidthY();
//float stepPhysZ = from.getPhysStepWidthZ();
float stepFourierX = from.getFourierStepWidthX();
float stepFourierY = from.getFourierStepWidthY();
float stepFourierZ = from.getFourierStepWidthZ();
TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
Vector3(from.getFourierSpaceMinX(),
from.getFourierSpaceMinY(),
from.getFourierSpaceMinZ()),
Vector3(from.getFourierSpaceMaxX(),
from.getFourierSpaceMaxY(),
from.getFourierSpaceMaxZ()));
// and fill it
// AR: old double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
double normalization=1./pow(sqrt(2.*Constants::PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
Index x, y, z;
Vector3 r;
typename TFFT3D<ComplexTraits>::Complex dataIn;
typename TFFT3D<ComplexTraits>::Complex dataOut;
for (Position i = 0; i < from.size(); i++)
{
z = i % lengthZ;
y = (i % (lengthY * lengthZ)) / lengthZ;
x = i / (lengthY * lengthZ);
if (x>lengthX/2.)
{
x-=lengthX;
}
if (y>lengthY/2.)
{
y-=lengthY;
}
if (z>lengthZ/2.)
{
z-=lengthZ;
}
r.set((float)x * stepFourierX,
(float)y * stepFourierY,
(float)z * stepFourierZ);
dataIn = from[i];
dataOut = dataIn;
newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
}
to = newGrid;
return to;
}
}
template <typename ComplexTraits>
const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from)
{
// first decide if the TFFT3D data is in Fourier space.
if (!from.isInFourierSpace())
{
// create a new grid
Size lengthX = from.getMaxXIndex()+1;
Size lengthY = from.getMaxYIndex()+1;
Size lengthZ = from.getMaxZIndex()+1;
RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
// and fill it
double normalization = 1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
typename TFFT3D<ComplexTraits>::Complex dataIn;
typename TFFT3D<ComplexTraits>::Complex dataOut;
for (Position i = 0; i < from.size(); i++)
{
Position x, y, z;
z = i % lengthZ;
y = (i % (lengthY * lengthZ)) / lengthZ;
x = i / (lengthY * lengthZ);
dataIn = from[i];
dataOut = dataIn;
newGrid[x + (y + z*lengthY)*lengthZ] = dataOut.real()*normalization;
}
to = newGrid;
return to;
}
else
{
// we are in Fourier space
// create a new grid
Size lengthX = from.getMaxXIndex()+1;
Size lengthY = from.getMaxYIndex()+1;
Size lengthZ = from.getMaxZIndex()+1;
//float stepPhysX = from.getPhysStepWidthX();
//float stepPhysY = from.getPhysStepWidthY();
//float stepPhysZ = from.getPhysStepWidthZ();
float stepFourierX = from.getFourierStepWidthX();
float stepFourierY = from.getFourierStepWidthY();
float stepFourierZ = from.getFourierStepWidthZ();
RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getFourierSpaceMinX(), from.getFourierSpaceMinY(), from.getFourierSpaceMinZ()), Vector3(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY(), from.getFourierSpaceMaxZ()));
// and fill it
// AR: old version double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
double normalization=1./pow(sqrt(2.*Constants::PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
Index x, y, z;
signed int xp, yp, zp;
Vector3 r;
typename TFFT3D<ComplexTraits>::Complex dataIn;
typename TFFT3D<ComplexTraits>::Complex dataOut;
for (Position i = 0; i < from.size(); i++)
{
z = i % lengthZ;
y = (i % (lengthY * lengthZ)) / lengthZ;
x = i / (lengthY * lengthZ);
xp = x;
yp = y;
zp = z;
if (xp>=lengthX/2.)
{
xp-=(int)lengthX;
}
if (yp>=lengthY/2.)
{
yp-=(int)lengthY;
}
if (zp>=lengthZ/2.)
{
zp-=(int)lengthZ;
}
if (x>=lengthX/2.)
{
x-=(int)(lengthX/2.);
}
else
{
x+=(int)(lengthX/2.);
}
if (y>=lengthY/2.)
{
y-=(int)(lengthY/2.);
}
else
{
y+=(int)(lengthY/2.);
}
if (z>=lengthZ/2.)
{
z-=(int)(lengthZ/2.);
}
else
{
z+=(int)(lengthZ/2.);
}
r.set((float)xp * stepFourierX,
(float)yp * stepFourierY,
(float)zp * stepFourierZ);
dataIn = from[i];
dataOut = dataIn;
newGrid[x + (y + z*lengthY)*lengthZ] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
}
to = newGrid;
return to;
}
}
}
#endif // BALL_MATHS_TFFT3D_H
exception.h
BALL::Exception::OutOfGrid
Definition COMMON/exception.h:351
BALL::TFFT3D
Definition FFT3D.h:41
BALL::TFFT3D::Complex
std::complex< typename ComplexTraits::ComplexPrecision > Complex
Definition FFT3D.h:44
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysStepWidthX
double getPhysStepWidthX() const
BALL::TFFT3D::setNumberOfiFFTTransforms
void setNumberOfiFFTTransforms(Size num)
Definition FFT3D.h:299
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMinX
double getFourierSpaceMinX() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMaxZ
double getPhysSpaceMaxZ() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepPhysY_
double stepPhysY_
Definition FFT3D.h:328
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getGridCoordinates
Vector3 getGridCoordinates(Position position) const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::planCalculated_
bool planCalculated_
Definition FFT3D.h:341
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::destroy
virtual void destroy()
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMinY
double getPhysSpaceMinY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMaxY
double getPhysSpaceMaxY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::lengthY_
Size lengthY_
Definition FFT3D.h:323
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::minPhys_
Vector3 minPhys_
Definition FFT3D.h:330
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::dataLength_
Size dataLength_
Definition FFT3D.h:339
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::planForward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planForward_
Definition FFT3D.h:335
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMinY
double getFourierSpaceMinY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMaxZ
double getFourierSpaceMaxZ() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getMaxYIndex
Size getMaxYIndex() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepFourierY_
double stepFourierY_
Definition FFT3D.h:329
BALL::TFFT3D::setNumberOfFFTTransforms
void setNumberOfFFTTransforms(Size num)
Definition FFT3D.h:293
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepPhysZ_
double stepPhysZ_
Definition FFT3D.h:328
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMaxX
double getFourierSpaceMaxX() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::maxPhys_
Vector3 maxPhys_
Definition FFT3D.h:330
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::doFFT
void doFFT()
BALL::TFFT3D::operator[]
Complex & operator[](const Vector3 &pos)
Definition FFT3D.h:771
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMaxX
double getPhysSpaceMaxX() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::setData
void setData(const Vector3 &pos, Complex val)
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::doiFFT
void doiFFT()
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepPhysX_
double stepPhysX_
Definition FFT3D.h:328
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMaxY
double getFourierSpaceMaxY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::lengthX_
Size lengthX_
Definition FFT3D.h:323
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::lengthZ_
Size lengthZ_
Definition FFT3D.h:323
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierStepWidthX
double getFourierStepWidthX() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::origin_
Vector3 origin_
Definition FFT3D.h:327
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMinX
double getPhysSpaceMinX() const
BALL::TFFT3D::operator==
bool operator==(const TFFT3D &fft3d) const
Definition FFT3D.h:373
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getNumberOfInverseTransforms
Size getNumberOfInverseTransforms() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::planBackward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planBackward_
Definition FFT3D.h:336
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepFourierX_
double stepFourierX_
Definition FFT3D.h:329
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::inFourierSpace_
bool inFourierSpace_
Definition FFT3D.h:324
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::minFourier_
Vector3 minFourier_
Definition FFT3D.h:331
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMinZ
double getPhysSpaceMinZ() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysStepWidthZ
double getPhysStepWidthZ() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::dataAdress_
Complex * dataAdress_
Definition FFT3D.h:340
BALL::TFFT3D::isInFourierSpace
bool isInFourierSpace() const
Definition FFT3D.h:921
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::numFourierToPhys_
Size numFourierToPhys_
Definition FFT3D.h:326
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getMaxXIndex
Size getMaxXIndex() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::setPhysStepWidth
bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
BALL::TFFT3D::phase
Complex phase(const Vector3 &pos) const
Definition FFT3D.h:891
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::translate
bool translate(const Vector3 &trans_origin)
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierStepWidthZ
double getFourierStepWidthZ() const
BALL::TFFT3D::operator[]
const Complex & operator[](const Position &pos) const
Definition FFT3D.h:287
BALL::TFFT3D::ComplexVector
TRegularData3D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
Definition FFT3D.h:45
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::stepFourierZ_
double stepFourierZ_
Definition FFT3D.h:329
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierStepWidthY
double getFourierStepWidthY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::numPhysToFourier_
Size numPhysToFourier_
Definition FFT3D.h:325
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getMaxZIndex
Size getMaxZIndex() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::clear
virtual void clear()
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getData
Complex getData(const Vector3 &pos) const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::maxFourier_
Vector3 maxFourier_
Definition FFT3D.h:331
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getPhysStepWidthY
double getPhysStepWidthY() const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getInterpolatedValue
Complex getInterpolatedValue(const Vector3 &pos) const
BALL::TFFT3D::TFFT3D
TFFT3D()
Default constructor.
Definition FFT3D.h:364
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMinZ
double getFourierSpaceMinZ() const
BALL::IndexType
BALL::TRegularData3D
Definition regularData3D.h:46
BALL::TRegularData3D< std::complex< ComplexTraits::ComplexPrecision > >::TRegularData3D
TRegularData3D()
Definition regularData3D.h:519
BALL::TRegularData3D< std::complex< ComplexTraits::ComplexPrecision > >::getData
const vector< std::complex< ComplexTraits::ComplexPrecision > > & getData() const
Definition regularData3D.h:877
BALL::TRegularData3D::operator[]
const ValueType & operator[](const IndexType &index) const
Definition regularData3D.h:226
BALL::TRegularData3D< std::complex< typename ComplexTraits::ComplexPrecision > >::size
BALL_INLINE size_type size() const
Definition regularData3D.h:190
BALL::TVector3::x
T x
Definition vector3.h:489
BALL::TVector3::y
T y
Definition vector3.h:493
BALL::TVector3::set
void set(const T *ptr)
Definition vector3.h:615
BALL::TVector3::z
T z
Definition vector3.h:497
BALL_CREATE
#define BALL_CREATE(name)
Definition create.h:62
fftwCommon.h
BALL::operator<<
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
BALL::Constants::PI
BALL_EXTERN_VARIABLE const double PI
PI.
Definition constants.h:35
BALL::Maths::rint
double rint(double x)
round to integral value in floating-point format
Definition MATHS/common.h:327
BALL
Definition constants.h:13
BALL::Size
BALL_SIZE_TYPE Size
Definition COMMON/global.h:114
BALL::FFT3D
TFFT3D< BALL_FFTW_DEFAULT_TRAITS > FFT3D
Definition FFT3D.h:347
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::Complex
std::complex< BALL_COMPLEX_PRECISION > Complex
Definition complex.h:21
BALL::Index
BALL_INDEX_TYPE Index
Definition COMMON/global.h:105
regularData3D.h