5 #ifndef BALL_MATHS_TFFT3D_H 6 #define BALL_MATHS_TFFT3D_H 8 #ifndef BALL_COMMON_EXCEPTION_H 13 #ifndef BALL_DATATYPE_REGULARDATA3D_H 38 template <
typename ComplexTraits>
40 :
public TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> >
44 typedef std::complex<typename ComplexTraits::ComplexPrecision>
Complex;
73 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);
129 bool setPhysStepWidth(
double new_width_x,
double new_width_y,
double new_width_z);
351 template <
typename ComplexTraits>
359 template <
typename ComplexTraits>
360 const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from);
363 template <
typename ComplexTraits>
372 template <
typename ComplexTraits>
403 for (
double posX=min.
x; posX<=max.
x; posX+=stepX)
405 for (
double posY=min.
y; posY<=max.
y; posY+=stepY)
407 for (
double posZ=min.
z; posZ<=max.
z; posZ+=stepZ)
423 template <
typename ComplexTraits>
449 template <
typename ComplexTraits>
452 if ((new_width_x <= 0) || (new_width_y <= 0) || (new_width_z <= 0))
474 template <
typename ComplexTraits>
480 template <
typename ComplexTraits>
486 template <
typename ComplexTraits>
492 template <
typename ComplexTraits>
498 template <
typename ComplexTraits>
504 template <
typename ComplexTraits>
510 template <
typename ComplexTraits>
516 template <
typename ComplexTraits>
522 template <
typename ComplexTraits>
528 template <
typename ComplexTraits>
534 template <
typename ComplexTraits>
540 template <
typename ComplexTraits>
546 template <
typename ComplexTraits>
553 template <
typename ComplexTraits>
559 template <
typename ComplexTraits>
565 template <
typename ComplexTraits>
571 template <
typename ComplexTraits>
577 template <
typename ComplexTraits>
583 template <
typename ComplexTraits>
589 template <
typename ComplexTraits>
595 template <
typename ComplexTraits>
601 template <
typename ComplexTraits>
607 template <
typename ComplexTraits>
667 template <
typename ComplexTraits>
671 double normalization=1.;
675 result = (*this)[pos];
682 result = (*this)[pos]*
phase(pos);
690 result *= normalization;
695 template <
typename ComplexTraits>
706 if ( (pos.
x < min.
x) || (pos.
y < min.
y) || (pos.
z < min.
z)
707 || (pos.
x > max.
x) || (pos.
y > max.
y) || (pos.
z > max.
z) )
713 double modX = fmod((
double)h.
x,stepX);
714 double modY = fmod((
double)h.
y,stepY);
715 double modZ = fmod((
double)h.
z,stepZ);
717 if (modX==0 && modY==0 && modZ==0)
722 double beforeX =
floor(h.
x/stepX)*stepX+min.
x;
723 double beforeY =
floor(h.
y/stepY)*stepY+min.
y;
724 double beforeZ =
floor(h.
z/stepZ)*stepZ+min.
z;
725 double afterX = ceil(h.
x/stepX)*stepX+min.
x;
726 double afterY = ceil(h.
y/stepY)*stepY+min.
y;
727 double afterZ = ceil(h.
z/stepZ)*stepZ+min.
z;
729 double tx = (pos.
x - beforeX)/stepX;
730 double ty = (pos.
y - beforeY)/stepY;
731 double tz = (pos.
z - beforeZ)/stepZ;
733 result =
getData(
Vector3(beforeX,beforeY,beforeZ))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*(1.-tz));
734 result +=
getData(
Vector3(afterX, beforeY,beforeZ))*(
typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)*(1.-tz));
735 result +=
getData(
Vector3(beforeX,afterY, beforeZ))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)* ty *(1.-tz));
736 result +=
getData(
Vector3(beforeX,beforeY,afterZ ))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)* tz );
737 result +=
getData(
Vector3(afterX, afterY, beforeZ))*(
typename ComplexTraits::ComplexPrecision)( tx * ty *(1.-tz));
738 result +=
getData(
Vector3(afterX, beforeY,afterZ ))*(
typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)* tz );
739 result +=
getData(
Vector3(beforeX,afterY, afterZ ))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)* ty * tz );
740 result +=
getData(
Vector3(afterX, afterY, afterZ ))*(
typename ComplexTraits::ComplexPrecision)( tx * ty * tz );
745 template <
typename ComplexTraits>
770 template <
typename ComplexTraits>
825 template <
typename ComplexTraits>
890 template <
typename ComplexTraits>
903 Complex result =
Complex(cos(phase), sin(phase));
920 template <
typename ComplexTraits>
926 template <
typename ComplexTraits>
931 if (!from.isInFourierSpace())
934 Size lengthX = from.getMaxXIndex()+1;
935 Size lengthY = from.getMaxYIndex()+1;
936 Size lengthZ = from.getMaxZIndex()+1;
939 Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
940 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
943 double normalization=1./(pow((
float)(lengthX*lengthY*lengthZ),(
int)from.getNumberOfInverseTransforms()));
947 for (
Position i = 0; i < from.size(); i++)
952 y = (i % (lengthY * lengthZ)) / lengthZ;
953 x = i / (lengthY * lengthZ);
958 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(
typename ComplexTraits::ComplexPrecision)normalization;
970 Size lengthX = from.getMaxXIndex()+1;
971 Size lengthY = from.getMaxYIndex()+1;
972 Size lengthZ = from.getMaxZIndex()+1;
976 float stepFourierX = from.getFourierStepWidthX();
977 float stepFourierY = from.getFourierStepWidthY();
978 float stepFourierZ = from.getFourierStepWidthZ();
983 Vector3(from.getFourierSpaceMinX(),
984 from.getFourierSpaceMinY(),
985 from.getFourierSpaceMinZ()),
986 Vector3(from.getFourierSpaceMaxX(),
987 from.getFourierSpaceMaxY(),
988 from.getFourierSpaceMaxZ()));
992 double normalization=1./pow(sqrt(2.*
Constants::PI),3)/(pow((
float)(lengthX*lengthY*lengthZ),(
int)from.getNumberOfInverseTransforms()));
1000 for (
Position i = 0; i < from.size(); i++)
1003 y = (i % (lengthY * lengthZ)) / lengthZ;
1004 x = i / (lengthY * lengthZ);
1021 r.
set((
float)x * stepFourierX,
1022 (
float)y * stepFourierY,
1023 (
float)z * stepFourierZ);
1028 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(
typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
1037 template <
typename ComplexTraits>
1038 const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from)
1041 if (!from.isInFourierSpace())
1044 Size lengthX = from.getMaxXIndex()+1;
1045 Size lengthY = from.getMaxYIndex()+1;
1046 Size lengthZ = from.getMaxZIndex()+1;
1048 RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ),
Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
1049 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
1052 double normalization = 1./(pow((
float)(lengthX*lengthY*lengthZ),(
int)from.getNumberOfInverseTransforms()));
1056 for (
Position i = 0; i < from.size(); i++)
1061 y = (i % (lengthY * lengthZ)) / lengthZ;
1062 x = i / (lengthY * lengthZ);
1067 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut.real()*normalization;
1079 Size lengthX = from.getMaxXIndex()+1;
1080 Size lengthY = from.getMaxYIndex()+1;
1081 Size lengthZ = from.getMaxZIndex()+1;
1085 float stepFourierX = from.getFourierStepWidthX();
1086 float stepFourierY = from.getFourierStepWidthY();
1087 float stepFourierZ = from.getFourierStepWidthZ();
1091 RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ),
Vector3(from.getFourierSpaceMinX(), from.getFourierSpaceMinY(), from.getFourierSpaceMinZ()),
Vector3(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY(), from.getFourierSpaceMaxZ()));
1095 double normalization=1./pow(sqrt(2.*
Constants::PI),3)/(pow((
float)(lengthX*lengthY*lengthZ),(
int)from.getNumberOfInverseTransforms()));
1098 signed int xp, yp, zp;
1103 for (
Position i = 0; i < from.size(); i++)
1106 y = (i % (lengthY * lengthZ)) / lengthZ;
1107 x = i / (lengthY * lengthZ);
1128 x-=(int)(lengthX/2.);
1132 x+=(int)(lengthX/2.);
1137 y-=(int)(lengthY/2.);
1141 y+=(int)(lengthY/2.);
1146 z-=(int)(lengthZ/2.);
1150 z+=(int)(lengthZ/2.);
1153 r.
set((
float)xp * stepFourierX,
1154 (
float)yp * stepFourierY,
1155 (
float)zp * stepFourierZ);
1160 newGrid[x + (y + z*lengthY)*lengthZ] = (dataOut*(
typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
1170 #endif // BALL_MATHS_TFFT3D_H const Complex & operator[](const Position &pos) const
double getFourierSpaceMinX() const
double getPhysStepWidthZ() const
T max(const T &a, const T &b)
double getFourierStepWidthX() const
TFFT3D< BALL_FFTW_DEFAULT_TRAITS > FFT3D
virtual ~TFFT3D()
Destructor.
double getFourierStepWidthZ() const
std::complex< typename ComplexTraits::ComplexPrecision > Complex
const TFFT3D & operator=(const TFFT3D &fft_3d)
Assignment operator.
double getPhysSpaceMinY() const
double getPhysStepWidthY() const
Size getMaxZIndex() const
const ValueType & operator[](const IndexType &index) const
ComplexTraits::FftwPlan planForward_
TRegularData3D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
double getFourierSpaceMaxY() const
bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
bool isInFourierSpace() const
Size getMaxYIndex() const
Complex & operator[](const Vector3 &pos)
BALL_INLINE size_type size() const
T min(const T &a, const T &b)
Complex & operator[](const Position &pos)
const vector< std::complex< ComplexTraits::ComplexPrecision > > & getData() const
Get the full data.
double getPhysSpaceMinZ() const
ComplexTraits::FftwPlan planBackward_
double getFourierStepWidthY() const
bool translate(const Vector3 &trans_origin)
BALL_EXTERN_VARIABLE const double PI
PI.
Complex getInterpolatedValue(const Vector3 &pos) const
TVector3< float > Vector3
double rint(double x)
round to integral value in floating-point format
double getFourierSpaceMinY() const
double getFourierSpaceMaxZ() const
Vector3 getGridCoordinates(Position position) const
double getFourierSpaceMaxX() const
BALL_EXTERN_VARIABLE const double k
TFFT3D()
Default constructor.
Complex getData(const Vector3 &pos) const
double getPhysSpaceMaxX() const
void setNumberOfFFTTransforms(Size num)
bool operator==(const TFFT3D &fft3d) const
void setData(const Vector3 &pos, Complex val)
BALL_EXTERN_VARIABLE const double h
double getPhysSpaceMaxY() const
void setNumberOfiFFTTransforms(Size num)
double getFourierSpaceMinZ() const
double getPhysStepWidthX() const
Size getMaxXIndex() const
double getPhysSpaceMinX() const
double getPhysSpaceMaxZ() const
#define BALL_CREATE(name)
Complex phase(const Vector3 &pos) const
Size getNumberOfInverseTransforms() const