5#ifndef BALL_MATHS_TFFT2D_H
6#define BALL_MATHS_TFFT2D_H
8#ifndef BALL_COMMON_EXCEPTION_H
12#ifndef BALL_DATATYPE_REGULARDATA2D_H
16#ifndef BALL_MATHS_VECTOR2_H
40 template <
typename ComplexTraits>
42 :
public TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >
46 typedef std::complex<typename ComplexTraits::ComplexPrecision>
Complex;
322 template <
typename ComplexTraits>
330 template <
typename ComplexTraits>
333 template <
typename ComplexTraits>
338 planCalculated_(false)
342 template <
typename ComplexTraits>
362 Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_;
363 Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
364 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
365 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
367 for (
double posX=min.x; posX<=max.x; posX+=stepX)
369 for (
double posY=min.y; posY<=max.y; posY+=stepY)
384 template <
typename ComplexTraits>
390 if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_))
392 origin_.x = trans_origin.
x;
393 origin_.y = trans_origin.
y;
395 minPhys_ =
Vector2(-origin_.x,-origin_.y);
396 maxPhys_ =
Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
397 minFourier_ =
Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
398 maxFourier_ =
Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
408 template <
typename ComplexTraits>
411 if ((new_width_x <= 0) || (new_width_y <= 0))
417 stepPhysX_ = new_width_x;
418 stepPhysY_ = new_width_y;
422 minPhys_ =
Vector2(-origin_.x,-origin_.y);
423 maxPhys_ =
Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
424 minFourier_ =
Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
425 maxFourier_ =
Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
431 template <
typename ComplexTraits>
437 template <
typename ComplexTraits>
443 template <
typename ComplexTraits>
446 return stepFourierX_;
449 template <
typename ComplexTraits>
452 return stepFourierY_;
455 template <
typename ComplexTraits>
461 template <
typename ComplexTraits>
467 template <
typename ComplexTraits>
473 template <
typename ComplexTraits>
479 template <
typename ComplexTraits>
482 return minFourier_.x;
485 template <
typename ComplexTraits>
488 return minFourier_.y;
491 template <
typename ComplexTraits>
494 return maxFourier_.x;
497 template <
typename ComplexTraits>
500 return maxFourier_.y;
503 template <
typename ComplexTraits>
506 return (lengthX_ - 1);
509 template <
typename ComplexTraits>
512 return (lengthY_ - 1);
515 template <
typename ComplexTraits>
518 return numFourierToPhys_;
522 template <
typename ComplexTraits>
525 if (!inFourierSpace_)
527 if (position >= ComplexVector::size())
537 y = position % lengthY_;
538 x = position / lengthY_;
540 r.
set(-origin_.x + (
float)x * stepPhysX_,
541 -origin_.y + (
float)y * stepPhysY_);
547 if (position >= ComplexVector::size())
556 y = position % lengthY_;
557 x = position / lengthY_;
569 r.
set((
float)x * stepFourierX_,
570 (
float)y * stepFourierY_);
578 template <
typename ComplexTraits>
582 double normalization=1.;
584 if (!inFourierSpace_)
586 result = (*this)[pos];
587 normalization=1./((
float)pow((
float)(lengthX_*lengthY_),(
int)numFourierToPhys_));
591 result = (*this)[pos] * phase(pos);
592 normalization=1./(2.*
Constants::PI)*(stepPhysX_*stepPhysY_)/((
float)pow((
float)(lengthX_*lengthY_),(
int)numFourierToPhys_));
595 result *= normalization;
600 template <
typename ComplexTraits>
605 Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_;
606 Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
607 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
608 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
610 if ( (pos.
x < min.x) || (pos.
y < min.y)
611 || (pos.
x > max.x) || (pos.
y > max.y) )
617 double modX = fmod((
double)h.x,stepX);
618 double modY = fmod((
double)h.y,stepY);
620 if (modX==0 && modY ==0)
625 double beforeX = floor(h.x/stepX)*stepX+ min.x;
626 double beforeY = floor(h.y/stepY)*stepY+ min.y;
627 double afterX = ceil(h.x/stepX)*stepX+ min.x;
628 double afterY = ceil(h.y/stepY)*stepY+ min.y;
630 double tx = (pos.
x - beforeX)/stepX;
631 double ty = (pos.
y - beforeY)/stepY;
633 result = getData(
Vector2(beforeX,beforeY))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty));
634 result += getData(
Vector2(afterX, beforeY))*(
typename ComplexTraits::ComplexPrecision)( tx *(1.-ty));
635 result += getData(
Vector2(beforeX,afterY ))*(
typename ComplexTraits::ComplexPrecision)((1.-tx)* ty );
636 result += getData(
Vector2(afterX, afterY ))*(
typename ComplexTraits::ComplexPrecision)( tx * ty );
641 template <
typename ComplexTraits>
646 if (!inFourierSpace_)
648 dummy =
Complex(val.real()*((
float)pow((
float)(lengthX_*lengthY_),(
int)numFourierToPhys_)),
649 val.imag()*((
float)pow((
float)(lengthX_*lengthY_),(
int)numFourierToPhys_)));
655 val*=phase(pos)*(
typename ComplexTraits::ComplexPrecision)((2*
Constants::PI/(stepPhysX_*stepPhysY_)))
656 *(
typename ComplexTraits::ComplexPrecision)pow((
typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_),(
int)numFourierToPhys_);
664 template <
typename ComplexTraits>
669 if (!inFourierSpace_)
676 internalPos = j + i*lengthY_;
695 internalPos = (j + i*lengthY_);
698 if ((internalPos < 0) || (internalPos>=(
Index) (lengthX_*lengthY_)))
703 return operator [] (internalPos);
706 template <
typename ComplexTraits>
711 if (!inFourierSpace_)
718 internalPos = j + i*lengthY_;
737 internalPos = (j + i*lengthY_);
740 if ((internalPos < 0) || (internalPos>=(
Index) (lengthX_*lengthY_)))
745 return operator [] (internalPos);
748 template <
typename ComplexTraits>
761 template <
typename ComplexTraits>
764 return inFourierSpace_;
767 template <
typename ComplexTraits>
772 if (!from.isInFourierSpace())
775 Size lengthX = from.getMaxXIndex()+1;
776 Size lengthY = from.getMaxYIndex()+1;
779 Vector2(from.getPhysSpaceMinX(), from.getPhysSpaceMinY()),
780 Vector2(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY()));
783 double normalization=1./(pow((
float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
787 for (
Position i = 0; i < from.size(); i++)
797 newGrid[x + y*lengthY] = dataOut*(
typename ComplexTraits::ComplexPrecision)normalization;
809 Size lengthX = from.getMaxXIndex()+1;
810 Size lengthY = from.getMaxYIndex()+1;
813 float stepFourierX = from.getFourierStepWidthX();
814 float stepFourierY = from.getFourierStepWidthY();
819 Vector2(from.getFourierSpaceMinX(),
820 from.getFourierSpaceMinY()),
821 Vector2(from.getFourierSpaceMaxX(),
822 from.getFourierSpaceMaxY()));
826 double normalization=1./(2.*
Constants::PI)/(pow((
float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
834 for (
Position i = 0; i < from.size(); i++)
849 r.
set((
float)x * stepFourierX,
850 (
float)y * stepFourierY);
855 newGrid[x + y*lengthY] = dataOut*(
typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
864 template <
typename ComplexTraits>
874 RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY),
886 typename RegularData2D::IndexType regdat_index;
887 for (current_index.
x = 0; current_index.
x < lengthX; current_index.
x++)
889 for (current_index.
y = 0; current_index.
y < lengthY; current_index.
y++)
891 regdat_index.
x = current_index.
x;
892 regdat_index.y = current_index.
y;
894 dataIn = from[current_index];
897 newGrid[regdat_index] = dataOut.real()*normalization;
929 RegularData2D::IndexType current_index;
949 x-=(int)(lengthX/2.);
953 x+=(int)(lengthX/2.);
958 y-=(int)(lengthY/2.);
962 y+=(int)(lengthY/2.);
966 r.
set((
float)xp * stepFourierX,
967 (
float)yp * stepFourierY);
975 newGrid[current_index] = (dataOut*(
typename ComplexTraits::ComplexPrecision)normalization*from.
phase(r)).real();
#define BALL_CREATE(name)
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
TVector2< float > Vector2
TFFT2D< BALL_FFTW_DEFAULT_TRAITS > FFT2D
std::complex< BALL_COMPLEX_PRECISION > Complex
BALL_EXTERN_VARIABLE const double PI
PI.
double rint(double x)
round to integral value in floating-point format
BALL_INLINE size_type size() const
const ValueType & operator[](const IndexType &index) const
The index type used to refer to a specific element in the grid (x-, and y-index)
const TFFT2D & operator=(const TFFT2D &fft_2d)
Assignment operator.
Size getMaxXIndex() const
void setNumberOfFFTTransforms(Size num)
double getFourierSpaceMinX() const
double getPhysSpaceMinY() const
Complex getData(const Vector2 &pos) const
TFFT2D(Size nX, Size nY, double stepPhysX=1., double stepPhysY=1., Vector2 origin=Vector2(0., 0.), bool inFourierSpace=false)
double getPhysStepWidthY() const
virtual ~TFFT2D()
Destructor.
double getPhysStepWidthX() const
void setData(const Vector2 &pos, Complex val)
double getPhysSpaceMaxY() const
const Complex & operator[](const Position &pos) const
double getFourierSpaceMaxY() const
TFFT2D(const TFFT2D &data)
Copy constructor.
double getPhysSpaceMinX() const
bool setPhysStepWidth(double new_width_x, double new_width_y)
bool isInFourierSpace() const
Complex & operator[](const Vector2 &pos)
Complex & operator[](const Position &pos)
void setNumberOfiFFTTransforms(Size num)
TRegularData2D< std::complex< typenameComplexTraits::ComplexPrecision > >::IndexType IndexType
Complex phase(const Vector2 &pos) const
Size getNumberOfInverseTransforms() const
Complex getInterpolatedValue(const Vector2 &pos) const
std::complex< typename ComplexTraits::ComplexPrecision > Complex
bool translate(const Vector2 &trans_origin)
Size getMaxYIndex() const
ComplexTraits::FftwPlan planForward_
double getFourierSpaceMaxX() const
double getFourierSpaceMinY() const
ComplexTraits::FftwPlan planBackward_
TFFT2D()
Default constructor.
double getFourierStepWidthX() const
double getFourierStepWidthY() const
TRegularData2D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
bool operator==(const TFFT2D &fft_2d) const
Vector2 getGridCoordinates(Position position) const
double getPhysSpaceMaxX() const