BALL  1.4.1
FFT2D.h
Go to the documentation of this file.
00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 
00005 #ifndef BALL_MATHS_TFFT2D_H
00006 #define BALL_MATHS_TFFT2D_H
00007 
00008 #ifndef BALL_COMMON_EXCEPTION_H
00009 # include <BALL/COMMON/exception.h>
00010 #endif
00011 
00012 #ifndef BALL_DATATYPE_REGULARDATA2D_H
00013 # include <BALL/DATATYPE/regularData2D.h>
00014 #endif
00015 
00016 #ifndef BALL_MATHS_VECTOR2_H
00017 # include <BALL/MATHS/vector2.h>
00018 #endif
00019 
00020 #include <math.h>
00021 #include <complex>
00022 #include <fftw3.h>
00023 
00024 #include <BALL/MATHS/fftwCommon.h>
00025 
00026 
00027 namespace BALL
00028 {
00040   template <typename ComplexTraits>
00041   class TFFT2D 
00042     : public TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >
00043   {
00044     public:
00045     
00046       typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
00047       typedef TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
00048       typedef typename TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >::IndexType IndexType;
00049 
00050       BALL_CREATE(TFFT2D)
00051 
00052       
00055  
00056       
00057       TFFT2D();
00058 
00060       TFFT2D(const TFFT2D &data);
00061 
00071        // ldn is not any longer the binary logarithm but the absolute number of grid points
00072       TFFT2D(Size nX, Size nY, double stepPhysX=1., double stepPhysY=1., Vector2 origin=Vector2(0.,0.), bool inFourierSpace=false);
00073 
00075       virtual ~TFFT2D();
00076       
00078 
00082 
00084       const TFFT2D& operator = (const TFFT2D& fft_2d);
00085       
00088       virtual void clear();
00089       
00092       virtual void destroy();
00093 
00095 
00099 
00102       bool operator == (const TFFT2D& fft_2d) const;
00104       
00105       // @name Accessors
00106       
00109       void doFFT();
00110 
00113       void doiFFT();
00114 
00120       bool translate(const Vector2& trans_origin);
00121 
00127       bool setPhysStepWidth(double new_width_x, double new_width_y);
00128 
00131       double getPhysStepWidthX() const;
00132 
00135       double getPhysStepWidthY() const;
00136 
00139       double getFourierStepWidthX() const;
00140 
00143       double getFourierStepWidthY() const;
00144 
00147       double getPhysSpaceMinX() const;
00148 
00151       double getPhysSpaceMinY() const;
00152 
00155       double getPhysSpaceMaxX() const;
00156 
00159       double getPhysSpaceMaxY() const;
00160 
00163       double getFourierSpaceMinX() const;
00164 
00167       double getFourierSpaceMinY() const;
00168 
00171       double getFourierSpaceMaxX() const;
00172 
00175       double getFourierSpaceMaxY() const;
00176         
00182       Size getMaxXIndex() const;
00183 
00189       Size getMaxYIndex() const;
00190         
00194       Size getNumberOfInverseTransforms() const;
00195 
00198       Vector2 getGridCoordinates(Position position) const;
00199 
00206       Complex getData(const Vector2& pos) const;
00207 
00215       Complex getInterpolatedValue(const Vector2& pos) const;
00216 
00224       void setData(const Vector2& pos, Complex val);
00225 
00231       Complex& operator[](const Vector2& pos);
00232 
00238       const Complex& operator[](const Vector2& pos) const;
00239 
00244       const Complex& operator [] (const IndexType& index) const { return TRegularData2D<Complex>::operator [](index); }
00245 
00250       Complex& operator [] (const IndexType& index) {  return TRegularData2D<Complex>::operator [](index); }
00251       
00256       Complex& operator[](const Position& pos)
00257       {
00258         return TRegularData2D<Complex>::operator [] (pos);
00259       }
00260 
00265       const Complex& operator[](const Position& pos) const
00266       {
00267         return TRegularData2D<Complex>::operator [] (pos);
00268       }
00269       
00270       // 
00271       void setNumberOfFFTTransforms(Size num)
00272       {
00273         numPhysToFourier_ = num;
00274       }
00275       
00276       // 
00277       void setNumberOfiFFTTransforms(Size num)
00278       {
00279         numFourierToPhys_ = num;
00280       }
00281       
00286       Complex phase(const Vector2& pos) const;
00287         
00291       bool isInFourierSpace() const;
00292 
00293     protected:
00294       Size lengthX_, lengthY_;
00295       bool inFourierSpace_;
00296       Size numPhysToFourier_;
00297       Size numFourierToPhys_;
00298       Vector2 origin_;
00299       double stepPhysX_, stepPhysY_;
00300       double stepFourierX_, stepFourierY_;
00301       Vector2 minPhys_, maxPhys_;
00302       Vector2 minFourier_, maxFourier_;
00303       
00304       // new version for FFTW3
00305       typename ComplexTraits::FftwPlan planForward_;
00306       typename ComplexTraits::FftwPlan planBackward_;
00307 
00308       
00309       // to control plan calculation with new fftw3
00310       Size dataLength_;
00311       Complex *dataAdress_;
00312       bool planCalculated_;
00313       
00314   };
00315   
00318   typedef TFFT2D<BALL_FFTW_DEFAULT_TRAITS> FFT2D;
00319   
00322   template <typename ComplexTraits>
00323   const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator<< 
00324       (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from);
00325   
00330   template <typename ComplexTraits>
00331   const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from);
00332   
00333   template <typename ComplexTraits>
00334   TFFT2D<ComplexTraits>::TFFT2D()
00335     : TRegularData2D<Complex>(),
00336       dataLength_(0),
00337       dataAdress_(0),
00338       planCalculated_(false)
00339   {
00340   }
00341   
00342   template <typename ComplexTraits>
00343   bool TFFT2D<ComplexTraits>::operator == (const TFFT2D<ComplexTraits>& fft2D) const
00344   {
00345     // test whether data_.size() == fft2D.data_.size()
00346     // instead of testing 2 lengths. Better for vector handling.
00347     
00348     if (lengthX_ == fft2D.lengthX_ &&
00349         lengthY_ == fft2D.lengthY_ &&
00350         origin_ == fft2D.origin_ &&
00351         stepPhysX_ == fft2D.stepPhysX_ &&
00352         stepPhysY_ == fft2D.stepPhysY_ &&
00353         stepFourierX_ == fft2D.stepFourierX_ &&
00354         stepFourierY_ == fft2D.stepFourierY_ &&
00355         minPhys_ == fft2D.minPhys_ &&
00356         maxPhys_ == fft2D.maxPhys_ &&
00357         minFourier_ == fft2D.minFourier_ &&
00358         maxFourier_ == fft2D.maxFourier_ &&
00359         numPhysToFourier_ == fft2D.numPhysToFourier_ &&
00360         numFourierToPhys_ == fft2D.numFourierToPhys_)
00361     {
00362       Vector2 min  = inFourierSpace_ ?  minFourier_  :   minPhys_;
00363       Vector2 max  = inFourierSpace_ ?  maxFourier_  :   maxPhys_;
00364       double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
00365       double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;  
00366       
00367       for (double posX=min.x; posX<=max.x; posX+=stepX)
00368       {
00369         for (double posY=min.y; posY<=max.y; posY+=stepY)
00370         {
00371           if (getData(Vector2(posX,posY)) != fft2D.getData(Vector2(posX,posY)))
00372           {
00373             return false;
00374           }
00375         }
00376       }
00377       
00378       return true;
00379     }
00380   
00381     return false;
00382   }
00383   
00384   template <typename ComplexTraits>
00385   bool TFFT2D<ComplexTraits>::translate(const Vector2& trans_origin)
00386   {
00387     Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_);
00388     Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_);
00389     
00390     if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_))
00391     {
00392       origin_.x = trans_origin.x;
00393       origin_.y = trans_origin.y;
00394       
00395       minPhys_ = Vector2(-origin_.x,-origin_.y);
00396       maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
00397       minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
00398       maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
00399      
00400       return true;
00401     }
00402     else
00403     {
00404       return false;
00405     }
00406   }
00407 
00408   template <typename ComplexTraits>
00409   bool TFFT2D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y)
00410   {
00411     if ((new_width_x <= 0) || (new_width_y <= 0))
00412     {
00413       return false;
00414     }
00415     else
00416     {
00417       stepPhysX_ = new_width_x;
00418       stepPhysY_ = new_width_y;
00419       stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_);
00420       stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_);
00421 
00422       minPhys_ = Vector2(-origin_.x,-origin_.y);
00423       maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
00424       minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
00425       maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
00426   
00427       return true;
00428     }
00429   }
00430   
00431   template <typename ComplexTraits>
00432   double TFFT2D<ComplexTraits>::getPhysStepWidthX() const
00433   {
00434     return stepPhysX_;
00435   }
00436 
00437   template <typename ComplexTraits>
00438   double TFFT2D<ComplexTraits>::getPhysStepWidthY() const
00439   {
00440     return stepPhysY_;
00441   }
00442 
00443   template <typename ComplexTraits>
00444   double TFFT2D<ComplexTraits>::getFourierStepWidthX() const
00445   {
00446     return stepFourierX_;
00447   }
00448 
00449   template <typename ComplexTraits>
00450   double TFFT2D<ComplexTraits>::getFourierStepWidthY() const
00451   {
00452     return stepFourierY_;
00453   }
00454 
00455   template <typename ComplexTraits>
00456   double TFFT2D<ComplexTraits>::getPhysSpaceMinX() const
00457   {
00458     return minPhys_.x;
00459   }
00460 
00461   template <typename ComplexTraits>
00462   double TFFT2D<ComplexTraits>::getPhysSpaceMinY() const
00463   {
00464     return minPhys_.y;
00465   }
00466 
00467   template <typename ComplexTraits>
00468   double TFFT2D<ComplexTraits>::getPhysSpaceMaxX() const
00469   {
00470     return maxPhys_.x;
00471   }
00472 
00473   template <typename ComplexTraits>
00474   double TFFT2D<ComplexTraits>::getPhysSpaceMaxY() const
00475   {
00476     return maxPhys_.y;
00477   }
00478 
00479   template <typename ComplexTraits>
00480   double TFFT2D<ComplexTraits>::getFourierSpaceMinX() const
00481   {
00482     return minFourier_.x;
00483   }
00484 
00485   template <typename ComplexTraits>
00486   double TFFT2D<ComplexTraits>::getFourierSpaceMinY() const
00487   {
00488     return minFourier_.y;
00489   }
00490 
00491   template <typename ComplexTraits>
00492   double TFFT2D<ComplexTraits>::getFourierSpaceMaxX() const
00493   {
00494     return maxFourier_.x;
00495   }
00496 
00497   template <typename ComplexTraits>
00498   double TFFT2D<ComplexTraits>::getFourierSpaceMaxY() const
00499   {
00500     return maxFourier_.y;
00501   }
00502   
00503   template <typename ComplexTraits>
00504   Size TFFT2D<ComplexTraits>::getMaxXIndex() const
00505   {
00506     return (lengthX_ - 1);
00507   }
00508   
00509   template <typename ComplexTraits>
00510   Size TFFT2D<ComplexTraits>::getMaxYIndex() const
00511   {
00512     return (lengthY_ - 1);
00513   }
00514   
00515   template <typename ComplexTraits>
00516   Size TFFT2D<ComplexTraits>::getNumberOfInverseTransforms() const
00517   {
00518     return numFourierToPhys_;
00519   }
00520   
00521   // new for compatibility with FFT3D
00522   template <typename ComplexTraits>
00523   Vector2 TFFT2D<ComplexTraits>::getGridCoordinates(Position position) const
00524   {
00525     if (!inFourierSpace_)
00526     {
00527       if (position >= ComplexVector::size())
00528       {
00529         throw Exception::OutOfGrid(__FILE__, __LINE__);
00530       }
00531     
00532       Vector2 r;
00533       Position  x, y;
00534 
00535 
00536       // AR: ??????
00537       y = position % lengthY_;
00538       x = position / lengthY_;
00539 
00540       r.set(-origin_.x + (float)x * stepPhysX_,
00541             -origin_.y + (float)y * stepPhysY_);
00542 
00543       return r;
00544     }
00545     else
00546     {
00547       if (position >= ComplexVector::size())
00548       {
00549         throw Exception::OutOfGrid(__FILE__, __LINE__);
00550       }
00551     
00552       Vector2 r;
00553       Index x, y;
00554   
00555       // AR: ??????
00556       y = position % lengthY_;
00557       x = position / lengthY_;
00558 
00559       if (x>=lengthX_/2.)
00560       {
00561         x-=lengthX_;
00562       }
00563       
00564       if (y>=lengthY_/2.)
00565       {
00566         y-=lengthY_;
00567       }
00568 
00569       r.set((float)x * stepFourierX_,
00570             (float)y * stepFourierY_);
00571 
00572       return r;
00573     }
00574   }
00575   
00576   
00577   
00578   template <typename ComplexTraits>
00579   typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getData(const Vector2& pos) const
00580   {
00581     Complex result;
00582     double normalization=1.;
00583 
00584     if (!inFourierSpace_)
00585     {
00586       result = (*this)[pos];
00587       normalization=1./((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_));
00588     }
00589     else
00590     {
00591       result = (*this)[pos] * phase(pos);
00592       normalization=1./(2.*M_PI)*(stepPhysX_*stepPhysY_)/((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_));
00593     }
00594 
00595     result *= normalization;
00596     
00597     return result;
00598   }
00599 
00600   template <typename ComplexTraits>
00601   typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getInterpolatedValue(const Vector2& pos) const
00602   {
00603     Complex result;
00604     
00605     Vector2 min  = inFourierSpace_ ? minFourier_   :   minPhys_;
00606     Vector2 max  = inFourierSpace_ ? maxFourier_   :   maxPhys_;
00607     double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
00608     double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
00609     
00610     if (    (pos.x < min.x) || (pos.y < min.y)
00611          || (pos.x > max.x) || (pos.y > max.y)  )
00612     {
00613       throw Exception::OutOfGrid(__FILE__, __LINE__);
00614     }
00615 
00616     Vector2 h(pos.x - min.x, pos.y - min.y);
00617     double modX = fmod((double)h.x,stepX);
00618     double modY = fmod((double)h.y,stepY);
00619 
00620     if (modX==0 && modY ==0) // we are on the grid
00621     {
00622       return getData(pos);
00623     }
00624 
00625     double beforeX = floor(h.x/stepX)*stepX+ min.x;
00626     double beforeY = floor(h.y/stepY)*stepY+ min.y;
00627     double afterX  =  ceil(h.x/stepX)*stepX+ min.x;
00628     double afterY  =  ceil(h.y/stepY)*stepY+ min.y;
00629       
00630     double tx = (pos.x - beforeX)/stepX;
00631     double ty = (pos.y - beforeY)/stepY;
00632 
00633     result  = getData(Vector2(beforeX,beforeY))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty));
00634     result += getData(Vector2(afterX, beforeY))*(typename ComplexTraits::ComplexPrecision)(    tx *(1.-ty));
00635     result += getData(Vector2(beforeX,afterY ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*    ty );
00636     result += getData(Vector2(afterX, afterY ))*(typename ComplexTraits::ComplexPrecision)(    tx *    ty );
00637 
00638     return result;
00639   }
00640 
00641   template <typename ComplexTraits>
00642   void TFFT2D<ComplexTraits>::setData(const Vector2& pos, Complex val)
00643   {
00644     Complex dummy;
00645   
00646     if (!inFourierSpace_)
00647     {
00648       dummy = Complex(val.real()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)),
00649                         val.imag()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)));
00650   
00651       (*this)[pos]=dummy;
00652     }
00653     else
00654     {
00655       val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((2*M_PI/(stepPhysX_*stepPhysY_)))
00656                      *(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_),(int)numFourierToPhys_);
00657       
00658       dummy = val;
00659     
00660       (*this)[pos]=dummy;
00661     }
00662   }
00663 
00664   template <typename ComplexTraits>
00665   typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos)
00666   {
00667     Index internalPos;
00668 
00669     if (!inFourierSpace_)
00670     {
00671       Index i, j;
00672       
00673       i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
00674       j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
00675 
00676       internalPos = j + i*lengthY_;
00677     }
00678     else
00679     {
00680       Index i, j;
00681 
00682       i = (Index) Maths::rint(pos.x/stepFourierX_);
00683       j = (Index) Maths::rint(pos.y/stepFourierY_);
00684 
00685       if (i<0)
00686       {
00687         i+=lengthX_;
00688       }
00689 
00690       if (j<0)
00691       {
00692         j+=lengthY_;
00693       }
00694 
00695       internalPos = (j + i*lengthY_);
00696     }
00697 
00698     if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_)))
00699     {
00700       throw Exception::OutOfGrid(__FILE__, __LINE__);
00701     }
00702     
00703     return operator [] (internalPos);
00704   }
00705 
00706   template <typename ComplexTraits>
00707   const typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos) const
00708   {
00709     Index internalPos;
00710 
00711     if (!inFourierSpace_)
00712     {
00713       Index i, j;
00714       
00715       i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
00716       j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
00717 
00718       internalPos = j + i*lengthY_;
00719     }
00720     else
00721     {
00722       Index i, j;
00723 
00724       i = (Index) Maths::rint(pos.x/stepFourierX_);
00725       j = (Index) Maths::rint(pos.y/stepFourierY_);
00726 
00727       if (i<0)
00728       {
00729         i+=lengthX_;
00730       }
00731 
00732       if (j<0)
00733       {
00734         j+=lengthY_;
00735       }
00736 
00737       internalPos = (j + i*lengthY_);
00738     }
00739 
00740     if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_)))
00741     {
00742       throw Exception::OutOfGrid(__FILE__, __LINE__);
00743     }
00744     
00745     return operator [] (internalPos);
00746   }
00747   
00748   template <typename ComplexTraits>
00749   typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::phase(const Vector2& pos) const
00750   {
00751     double phase = 2.*M_PI*(  Maths::rint(pos.x/stepFourierX_)*Maths::rint(origin_.x/stepPhysX_)
00752                               /lengthX_
00753                             + Maths::rint(pos.y/stepFourierY_)*Maths::rint(origin_.y/stepPhysY_)
00754                               /lengthY_ );
00755 
00756     Complex result = Complex(cos(phase), sin(phase));
00757             
00758     return result;
00759   }
00760   
00761   template <typename ComplexTraits>
00762   bool TFFT2D<ComplexTraits>::isInFourierSpace() const
00763   {
00764     return inFourierSpace_;
00765   }
00766   
00767   template <typename ComplexTraits>
00768   const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator << 
00769     (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from)
00770   {
00771     // first decide if the FFT3D data is in Fourier space.
00772     if (!from.isInFourierSpace())
00773     {
00774       // create a new grid
00775       Size lengthX = from.getMaxXIndex()+1;
00776       Size lengthY = from.getMaxYIndex()+1;
00777       
00778       TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY),
00779                                       Vector2(from.getPhysSpaceMinX(), from.getPhysSpaceMinY()),
00780                                       Vector2(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY()));
00781 
00782       // and fill it
00783       double normalization=1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
00784       typename TFFT2D<ComplexTraits>::Complex dataIn;
00785       typename TFFT2D<ComplexTraits>::Complex dataOut;
00786       
00787       for (Position i = 0; i < from.size(); i++)
00788       {
00789         Position x, y;
00790 
00791         y =  i % lengthY;
00792         x =  i / lengthY;
00793 
00794         dataIn  = from[i];
00795         dataOut = dataIn;
00796         
00797         newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization;
00798       }
00799 
00800       to = newGrid;
00801 
00802       return to;
00803     }
00804     else
00805     {
00806       // we are in Fourier space
00807       
00808       // create a new grid
00809       Size lengthX = from.getMaxXIndex()+1;
00810       Size lengthY = from.getMaxYIndex()+1;
00811       //float stepPhysX = from.getPhysStepWidthX();
00812       //float stepPhysY = from.getPhysStepWidthY();
00813       float stepFourierX = from.getFourierStepWidthX();
00814       float stepFourierY = from.getFourierStepWidthY();
00815 
00816 
00817       
00818       TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY),
00819                                       Vector2(from.getFourierSpaceMinX(), 
00820                                               from.getFourierSpaceMinY()),
00821                                       Vector2(from.getFourierSpaceMaxX(),
00822                                               from.getFourierSpaceMaxY()));
00823 
00824       // and fill it
00825       // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
00826       double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
00827       
00828       
00829       Index x, y;
00830       Vector2 r;
00831       typename TFFT2D<ComplexTraits>::Complex dataIn;
00832       typename TFFT2D<ComplexTraits>::Complex dataOut;
00833   
00834       for (Position i = 0; i < from.size(); i++)
00835       {
00836         y =  i % lengthY;
00837         x =  i / lengthY;
00838 
00839         if (x>lengthX/2.)
00840         {
00841           x-=lengthX;
00842         }
00843 
00844         if (y>lengthY/2.)
00845         {
00846           y-=lengthY;
00847         }
00848 
00849         r.set((float)x * stepFourierX,
00850               (float)y * stepFourierY);
00851 
00852         dataIn = from[i];
00853         dataOut = dataIn;
00854         
00855         newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
00856       }
00857 
00858       to = newGrid;
00859 
00860       return to;
00861     }
00862   }
00863   
00864   template <typename ComplexTraits>
00865   const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from)
00866   {
00867     // first decide if the FFT2D data is in Fourier space.
00868     if (!from.isInFourierSpace())
00869     {
00870       // create a new grid
00871       Size lengthX = from.getMaxXIndex()+1;
00872       Size lengthY = from.getMaxYIndex()+1;
00873       
00874       RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY),
00875                             Vector2(from.getPhysSpaceMinX(), 
00876                                     from.getPhysSpaceMinY()),
00877                             Vector2(from.getPhysSpaceMaxX(),
00878                                     from.getPhysSpaceMaxY()));
00879 
00880       // and fill it
00881       double normalization = 1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
00882       typename TFFT2D<ComplexTraits>::Complex dataIn;
00883       typename TFFT2D<ComplexTraits>::Complex dataOut;
00884 
00885       typename TFFT2D<ComplexTraits>::IndexType current_index;
00886       typename RegularData2D::IndexType regdat_index;
00887       for (current_index.x = 0; current_index.x < lengthX; current_index.x++)
00888       {
00889         for (current_index.y = 0; current_index.y < lengthY; current_index.y++)
00890         {
00891           regdat_index.x = current_index.x;
00892           regdat_index.y = current_index.y;
00893 
00894           dataIn  = from[current_index];
00895           dataOut = dataIn;
00896         
00897           newGrid[regdat_index] = dataOut.real()*normalization;
00898         }
00899       }
00900 
00901       to = newGrid;
00902 
00903       return to;
00904     }
00905     else
00906     {
00907       // we are in Fourier space
00908       
00909       // create a new grid
00910       Size lengthX = from.getMaxXIndex()+1;
00911       Size lengthY = from.getMaxYIndex()+1;
00912       //float stepPhysX = from.getPhysStepWidthX();
00913       //float stepPhysY = from.getPhysStepWidthY();
00914       float stepFourierX = from.getFourierStepWidthX();
00915       float stepFourierY = from.getFourierStepWidthY();
00916   
00917       RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY), Vector2(from.getFourierSpaceMinX(), from.getFourierSpaceMinY()), Vector2(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY()));
00918 
00919       // and fill it
00920       // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
00921       double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
00922       
00923       Index x, y;
00924       signed int xp, yp;
00925       Vector2 r;
00926       typename TFFT2D<ComplexTraits>::Complex dataIn;
00927       typename TFFT2D<ComplexTraits>::Complex dataOut;
00928   
00929       RegularData2D::IndexType current_index;
00930       for (Position i = 0; i < from.size(); i++)
00931       {
00932         y =  i % lengthY;
00933         x =  i / lengthY;
00934         
00935         xp = x;
00936         yp = y;
00937 
00938         if (xp>=lengthX/2.)
00939         {
00940           xp-=(int)lengthX;
00941         }
00942         if (yp>=lengthY/2.)
00943         {
00944           yp-=(int)lengthY;
00945         }
00946 
00947         if (x>=lengthX/2.)
00948         {
00949           x-=(int)(lengthX/2.);
00950         }
00951         else
00952         {
00953           x+=(int)(lengthX/2.);
00954         }
00955 
00956         if (y>=lengthY/2.)
00957         {
00958           y-=(int)(lengthY/2.);
00959         }
00960         else
00961         {
00962           y+=(int)(lengthY/2.);
00963         }
00964 
00965 
00966         r.set((float)xp * stepFourierX,
00967               (float)yp * stepFourierY);
00968 
00969         dataIn = from[i];
00970         dataOut = dataIn;
00971 
00972         current_index.x = x;
00973         current_index.y = y;
00974 
00975         newGrid[current_index] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
00976       }
00977 
00978       to = newGrid;
00979 
00980       return to;
00981     }
00982   }
00983   
00984   
00985   
00986 }
00987 
00988 #endif // BALL_MATHS_TFFT2D_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines