BALL: FFT3D.h Source File

00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 // $Id: FFT3D.h,v 1.17.16.1 2007/03/25 21:23:44 oliver Exp $
00005 //
00006 
00007 #ifndef BALL_MATHS_TFFT3D_H
00008 #define BALL_MATHS_TFFT3D_H
00009 
00010 #ifndef BALL_COMMON_EXCEPTION_H
00011 # include <BALL/COMMON/exception.h>
00012 #endif
00013 
00014 
00015 #ifndef BALL_DATATYPE_REGULARDATA3D_H
00016 # include <BALL/DATATYPE/regularData3D.h>
00017 #endif
00018 
00019 //#ifndef BALL_MATHS_VECTOR2_H
00020 //# include <BALL/MATHS/vector3.h>
00021 //#endif
00022 
00023 #include <BALL/MATHS/fftwCommon.h>
00024 #include <math.h>
00025 #include <complex>
00026 #include <fftw3.h>
00027 
00028 namespace BALL
00029 {
00040   template <typename ComplexTraits>
00041   class TFFT3D 
00042     : public TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> >
00043   {
00044     public:
00045     
00046       typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
00047       typedef TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
00048 
00049       BALL_CREATE(TFFT3D)
00050 
00051       
00054     
00055       
00056       TFFT3D();
00057     
00059       TFFT3D(const TFFT3D &data);
00060 
00074        // AR: ldn is not any longer the binary logarithm but the absolute number of grid points
00075       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);
00076 
00078       virtual ~TFFT3D();
00079           
00081     
00085     
00087       const TFFT3D& operator = (const TFFT3D& fft_3d);
00088     
00091       virtual void clear();
00092       
00095       virtual void destroy();
00096 
00098 
00102 
00105       bool operator == (const TFFT3D& fft3d) const;
00107       
00108       // @name Accessors
00109     
00111 
00113       void doFFT();
00114 
00117       void doiFFT();
00118 
00124       bool translate(const Vector3& trans_origin);
00125 
00131       bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z);
00132 
00135       double getPhysStepWidthX() const;
00136 
00139       double getPhysStepWidthY() const;
00140 
00143       double getPhysStepWidthZ() const;
00144 
00147       double getFourierStepWidthX() const;
00148 
00151       double getFourierStepWidthY() const;
00152 
00155       double getFourierStepWidthZ() const;
00156 
00159       double getPhysSpaceMinX() const;
00160 
00163       double getPhysSpaceMinY() const
00164 ;
00165 
00168       double getPhysSpaceMinZ() const;
00169 
00172       double getPhysSpaceMaxX() const;
00173 
00176       double getPhysSpaceMaxY() const;
00177 
00180       double getPhysSpaceMaxZ() const;
00181 
00184       double getFourierSpaceMinX() const;
00185 
00188       double getFourierSpaceMinY() const;
00189 
00192       double getFourierSpaceMinZ() const;
00193 
00196       double getFourierSpaceMaxX() const;
00197 
00200       double getFourierSpaceMaxY() const;
00201 
00204       double getFourierSpaceMaxZ() const;
00205 
00211       Size getMaxXIndex() const;
00212 
00218       Size getMaxYIndex() const;
00219 
00225       Size getMaxZIndex() const;
00226 
00230       Size getNumberOfInverseTransforms() const;
00231       
00234       Vector3 getGridCoordinates(Position position) const;
00235       
00240       Complex getData(const Vector3& pos) const
00241         throw(Exception::OutOfGrid);
00242 
00248       Complex getInterpolatedValue(const Vector3& pos) const
00249         throw(Exception::OutOfGrid);
00250 
00256       void setData(const Vector3& pos, Complex val)
00257         throw(Exception::OutOfGrid);
00258 
00262       Complex& operator[](const Vector3& pos)
00263         throw(Exception::OutOfGrid);
00264 
00268       const Complex& operator[](const Vector3& pos) const
00269         throw(Exception::OutOfGrid);
00270 
00273       Complex& operator[](const Position& pos)
00274         throw(Exception::OutOfGrid)
00275       {
00276         return TRegularData3D<Complex>::operator [] (pos);
00277       }
00278 
00281       const Complex& operator[](const Position& pos) const
00282         throw(Exception::OutOfGrid)
00283       {
00284         return TRegularData3D<Complex>::operator [] (pos);
00285       }
00286       
00287       // AR:
00288       void setNumberOfFFTTransforms(Size num)
00289       {
00290         numPhysToFourier_ = num;
00291       }
00292       
00293       // AR:
00294       void setNumberOfiFFTTransforms(Size num)
00295       {
00296         numFourierToPhys_ = num;
00297       }
00298       
00299         
00304       Complex phase(const Vector3& pos) const;
00306       
00310 
00314       bool isInFourierSpace() const;
00316       
00317     protected:
00318       Size lengthX_, lengthY_, lengthZ_;
00319       bool inFourierSpace_;
00320       Size numPhysToFourier_;
00321       Size numFourierToPhys_;
00322       Vector3 origin_;
00323       double stepPhysX_, stepPhysY_, stepPhysZ_;
00324       double stepFourierX_, stepFourierY_, stepFourierZ_;
00325       Vector3 minPhys_, maxPhys_;
00326       Vector3 minFourier_, maxFourier_;
00327       
00328       
00329       // AR: new version for FFTW3
00330       typename ComplexTraits::FftwPlan planForward_;
00331       typename ComplexTraits::FftwPlan planBackward_;
00332       
00333       // AR: to control plan calculation with new fftw3
00334       Size dataLength_;
00335       Complex *dataAdress_;
00336       bool planCalculated_;
00337       
00338   };
00339   
00342   typedef TFFT3D<BALL_FFTW_DEFAULT_TRAITS> FFT3D;
00343 
00346   template <typename ComplexTraits>
00347   const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator << 
00348     (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from);
00349   
00354   template <typename ComplexTraits>
00355   const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from);
00356     
00357     
00358   template <typename ComplexTraits>
00359   TFFT3D<ComplexTraits>::TFFT3D()
00360     : TRegularData3D<Complex>(),
00361       dataLength_(0),
00362       dataAdress_(0),
00363       planCalculated_(false)
00364   {
00365   }
00366   
00367   template <typename ComplexTraits>
00368   bool TFFT3D<ComplexTraits>::operator == (const TFFT3D& fft3D) const
00369   {
00370     
00371     // AR: test whether data_.size() == fft3D.data_.size()
00372     //     instead of testing 3 lengths. Better for vector handling.
00373     
00374     if (lengthX_ == fft3D.lengthX_ &&
00375         lengthY_ == fft3D.lengthY_ &&
00376         lengthZ_ == fft3D.lengthZ_ &&
00377         origin_ == fft3D.origin_ &&
00378         stepPhysX_ == fft3D.stepPhysX_ &&
00379         stepPhysY_ == fft3D.stepPhysY_ &&
00380         stepPhysZ_ == fft3D.stepPhysZ_ &&
00381         stepFourierX_ == fft3D.stepFourierX_ &&
00382         stepFourierY_ == fft3D.stepFourierY_ &&
00383         stepFourierZ_ == fft3D.stepFourierZ_ &&
00384         minPhys_ == fft3D.minPhys_ &&
00385         maxPhys_ == fft3D.maxPhys_ &&
00386         minFourier_ == fft3D.minFourier_ &&
00387         maxFourier_ == fft3D.maxFourier_ &&
00388         numPhysToFourier_ == fft3D.numPhysToFourier_ &&
00389         numFourierToPhys_ == fft3D.numFourierToPhys_ &&
00390         planCalculated_ == fft3D.planCalculated_)
00391     {
00392       Vector3 min  = inFourierSpace_ ?  minFourier_  :   minPhys_;
00393       Vector3 max  = inFourierSpace_ ?  maxFourier_  :   maxPhys_;
00394       double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
00395       double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;  
00396       double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;  
00397       
00398       for (double posX=min.x; posX<=max.x; posX+=stepX)
00399       {
00400         for (double posY=min.y; posY<=max.y; posY+=stepY)
00401         {
00402           for (double posZ=min.z; posZ<=max.z; posZ+=stepZ)
00403           {
00404             if (getData(Vector3(posX,posY,posZ)) != fft3D.getData(Vector3(posX,posY,posZ)))
00405             {
00406               return false;
00407             }
00408           }
00409         }
00410       }
00411       
00412       return true;
00413     }
00414   
00415     return false;
00416   }
00417   
00418   template <typename ComplexTraits>
00419   bool TFFT3D<ComplexTraits>::translate(const Vector3& trans_origin)
00420   {
00421     Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_);
00422     Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_);
00423     Position internalOriginZ = (Position) Maths::rint(trans_origin.z*stepPhysZ_);
00424     
00425     if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_) && (internalOriginZ <= lengthZ_))
00426     {
00427       origin_.x = trans_origin.x;
00428       origin_.y = trans_origin.y;
00429       origin_.z = trans_origin.z;
00430       
00431       minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
00432       maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
00433       minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
00434       maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
00435    
00436       return true;
00437     }
00438     else
00439     {
00440       return false;
00441     }
00442   }
00443 
00444   template <typename ComplexTraits>
00445   bool TFFT3D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
00446   {
00447     if ((new_width_x <= 0) || (new_width_y <= 0) || (new_width_z <= 0))
00448     {
00449       return false;
00450     }
00451     else
00452     {
00453       stepPhysX_ = new_width_x;
00454       stepPhysY_ = new_width_y;
00455       stepPhysZ_ = new_width_z;
00456       stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_);
00457       stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_);
00458       stepFourierZ_ = 2.*M_PI/(stepPhysZ_*lengthZ_);
00459 
00460       minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
00461       maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
00462       minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
00463       maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
00464   
00465       return true;
00466     }
00467   }
00468   
00469   template <typename ComplexTraits>
00470   double TFFT3D<ComplexTraits>::getPhysStepWidthX() const
00471   {
00472     return stepPhysX_;
00473   }
00474 
00475   template <typename ComplexTraits>
00476   double TFFT3D<ComplexTraits>::getPhysStepWidthY() const
00477   {
00478     return stepPhysY_;
00479   }
00480 
00481   template <typename ComplexTraits>
00482   double TFFT3D<ComplexTraits>::getPhysStepWidthZ() const
00483   {
00484     return stepPhysZ_;
00485   }
00486   
00487   template <typename ComplexTraits>
00488   double TFFT3D<ComplexTraits>::getFourierStepWidthX() const
00489   {
00490     return stepFourierX_;
00491   }
00492 
00493   template <typename ComplexTraits>
00494   double TFFT3D<ComplexTraits>::getFourierStepWidthY() const
00495   {
00496     return stepFourierY_;
00497   }
00498 
00499   template <typename ComplexTraits>
00500   double TFFT3D<ComplexTraits>::getFourierStepWidthZ() const
00501   {
00502     return stepFourierZ_;
00503   }
00504 
00505   template <typename ComplexTraits>
00506   double TFFT3D<ComplexTraits>::getPhysSpaceMinX() const
00507   {
00508     return minPhys_.x;
00509   }
00510 
00511   template <typename ComplexTraits>
00512   double TFFT3D<ComplexTraits>::getPhysSpaceMinY() const
00513   {
00514     return minPhys_.y;
00515   }
00516 
00517   template <typename ComplexTraits>
00518   double TFFT3D<ComplexTraits>::getPhysSpaceMinZ() const
00519   {
00520     return minPhys_.z;
00521   }
00522 
00523   template <typename ComplexTraits>
00524   double TFFT3D<ComplexTraits>::getPhysSpaceMaxX() const
00525   {
00526     return maxPhys_.x;
00527   }
00528 
00529   template <typename ComplexTraits>
00530   double TFFT3D<ComplexTraits>::getPhysSpaceMaxY() const
00531   {
00532     return maxPhys_.y;
00533   }
00534 
00535   template <typename ComplexTraits>
00536   double TFFT3D<ComplexTraits>::getPhysSpaceMaxZ() const
00537   {
00538     return maxPhys_.z;
00539   }
00540   
00541   template <typename ComplexTraits>
00542   double TFFT3D<ComplexTraits>::getFourierSpaceMinX() const
00543 
00544   {
00545     return minFourier_.x;
00546   }
00547 
00548   template <typename ComplexTraits>
00549   double TFFT3D<ComplexTraits>::getFourierSpaceMinY() const
00550   {
00551     return minFourier_.y;
00552   }
00553 
00554   template <typename ComplexTraits>
00555   double TFFT3D<ComplexTraits>::getFourierSpaceMinZ() const
00556   {
00557     return minFourier_.z;
00558   }
00559 
00560   template <typename ComplexTraits>
00561   double TFFT3D<ComplexTraits>::getFourierSpaceMaxX() const
00562   {
00563     return maxFourier_.x;
00564   }
00565 
00566   template <typename ComplexTraits>
00567   double TFFT3D<ComplexTraits>::getFourierSpaceMaxY() const
00568   {
00569     return maxFourier_.y;
00570   }
00571 
00572   template <typename ComplexTraits>
00573   double TFFT3D<ComplexTraits>::getFourierSpaceMaxZ() const
00574   {
00575     return maxFourier_.z;
00576   }
00577 
00578   template <typename ComplexTraits>
00579   Size TFFT3D<ComplexTraits>::getMaxXIndex() const
00580   {
00581     return (lengthX_ - 1);
00582   }
00583   
00584   template <typename ComplexTraits>
00585   Size TFFT3D<ComplexTraits>::getMaxYIndex() const
00586   {
00587     return (lengthY_ - 1);
00588   }
00589   
00590   template <typename ComplexTraits>
00591   Size TFFT3D<ComplexTraits>::getMaxZIndex() const
00592   {
00593     return (lengthZ_ - 1);
00594   }
00595   
00596   template <typename ComplexTraits>
00597   Size TFFT3D<ComplexTraits>::getNumberOfInverseTransforms() const
00598   {
00599     return numFourierToPhys_;
00600   }
00601 
00602   template <typename ComplexTraits>
00603   Vector3 TFFT3D<ComplexTraits>::getGridCoordinates(Position position) const
00604   {
00605     if (!inFourierSpace_)
00606     {
00607       if (position >= ComplexVector::size())
00608       {
00609         throw Exception::OutOfGrid(__FILE__, __LINE__);
00610       }
00611     
00612       Vector3 r;
00613       Position  x, y, z;
00614 
00615       z = position % lengthZ_;
00616       y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
00617       x =  position / (lengthY_ * lengthZ_);
00618 
00619       r.set(-origin_.x + (float)x * stepPhysX_,
00620             -origin_.y + (float)y * stepPhysY_,
00621             -origin_.z + (float)z * stepPhysZ_);
00622 
00623       return r;
00624     }
00625     else
00626     {
00627       if (position >= ComplexVector::size())
00628       {
00629         throw Exception::OutOfGrid(__FILE__, __LINE__);
00630       }
00631     
00632       Vector3 r;
00633       Index x, y, z;
00634   
00635       z = position % lengthZ_;
00636       y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
00637       x =  position / (lengthY_ * lengthZ_);
00638 
00639       if (x>=lengthX_/2.)
00640       {
00641         x-=lengthX_;
00642       }
00643       
00644       if (y>=lengthY_/2.)
00645       {
00646         y-=lengthY_;
00647       }
00648 
00649       if (z>=lengthZ_/2.)
00650       {
00651         z-=lengthZ_;
00652       }
00653 
00654       r.set((float)x * stepFourierX_,
00655             (float)y * stepFourierY_,
00656             (float)z * stepFourierZ_);
00657 
00658       return r;
00659     }
00660   }
00661   
00662   template <typename ComplexTraits>
00663   typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getData(const Vector3& pos) const
00664     throw(Exception::OutOfGrid)
00665   {
00666     Complex result;
00667     double normalization=1.;
00668 
00669     if (!inFourierSpace_)
00670     {
00671       result = (*this)[pos];
00672       normalization=1./((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
00673     }
00674     else
00675     {
00676       // AR:
00677       //old: result = (*this)[pos];
00678       result = (*this)[pos]*phase(pos);
00679       
00680       //normalization=1./pow(sqrt(2.*M_PI),3)*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
00681       
00682       normalization=1./pow(sqrt(2.*M_PI),3)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
00683       //normalization=1./(sqrt(2.*M_PI))*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
00684     }
00685 
00686     result *= normalization;
00687     
00688     return result;
00689   }
00690 
00691   template <typename ComplexTraits>
00692   typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getInterpolatedValue(const Vector3& pos) const
00693     throw(Exception::OutOfGrid)
00694   {
00695     Complex result;
00696     
00697     Vector3 min  = inFourierSpace_ ? minFourier_   :   minPhys_;
00698     Vector3 max  = inFourierSpace_ ? maxFourier_   :   maxPhys_;
00699     double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
00700     double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
00701     double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;
00702     
00703     if (    (pos.x < min.x) || (pos.y < min.y) || (pos.z < min.z)
00704          || (pos.x > max.x) || (pos.y > max.y) || (pos.z > max.z) )
00705     {
00706       throw Exception::OutOfGrid(__FILE__, __LINE__);
00707     }
00708 
00709     Vector3 h(pos.x - min.x, pos.y - min.y, pos.z - min.z);
00710     double modX = fmod((double)h.x,stepX);
00711     double modY = fmod((double)h.y,stepY);
00712     double modZ = fmod((double)h.z,stepZ);
00713 
00714     if (modX==0 && modY==0 && modZ==0) // we are on the grid
00715     {
00716       return getData(pos);
00717     }
00718 
00719     double beforeX = floor(h.x/stepX)*stepX+min.x;
00720     double beforeY = floor(h.y/stepY)*stepY+min.y;
00721     double beforeZ = floor(h.z/stepZ)*stepZ+min.z;
00722     double afterX  =  ceil(h.x/stepX)*stepX+min.x;
00723     double afterY  =  ceil(h.y/stepY)*stepY+min.y;
00724     double afterZ  =  ceil(h.z/stepZ)*stepZ+min.z;
00725       
00726     double tx = (pos.x - beforeX)/stepX;
00727     double ty = (pos.y - beforeY)/stepY;
00728     double tz = (pos.z - beforeZ)/stepZ;
00729 
00730     result  = getData(Vector3(beforeX,beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*(1.-tz));
00731     result += getData(Vector3(afterX, beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)(    tx *(1.-ty)*(1.-tz));
00732     result += getData(Vector3(beforeX,afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*    ty *(1.-tz));
00733     result += getData(Vector3(beforeX,beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*    tz );
00734     result += getData(Vector3(afterX, afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)(    tx *    ty *(1.-tz));
00735     result += getData(Vector3(afterX, beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)(    tx *(1.-ty)*    tz );
00736     result += getData(Vector3(beforeX,afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*    ty *    tz );
00737     result += getData(Vector3(afterX, afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)(    tx *    ty *    tz );
00738 
00739     return result;
00740   }
00741 
00742   template <typename ComplexTraits>
00743   void TFFT3D<ComplexTraits>::setData(const Vector3& pos, Complex val)
00744     throw(Exception::OutOfGrid)
00745   {
00746     Complex dummy;
00747   
00748     if (!inFourierSpace_)
00749     {
00750       dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)),
00751                         val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)));
00752   
00753       (*this)[pos]=dummy;
00754     }
00755     else
00756     {
00757       val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((pow(sqrt(2*M_PI),3)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
00758                      *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
00759       /*val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*M_PI)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
00760                      *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));*/
00761       
00762       dummy = val;
00763       
00764       (*this)[pos]=dummy;
00765     }
00766   }
00767 
00768   template <typename ComplexTraits>
00769   typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos)
00770     throw(Exception::OutOfGrid)
00771   {
00772     Index internalPos;
00773 
00774     if (!inFourierSpace_)
00775     {
00776       Index i, j, k;
00777 
00778       i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
00779       j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
00780       k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
00781       
00782       internalPos = (k + (j + i*lengthY_)*lengthZ_);
00783       
00784       /*(Index) rint(       (pos.z+origin_.z)/stepPhysZ_
00785                                   + (   (pos.y+origin_.y)/stepPhysY_
00786                                       + (pos.x+origin_.x)/stepPhysX_*lengthY_ 
00787                                     ) * lengthZ_
00788                                 ); */
00789     }
00790     else
00791     {
00792       Index i, j, k;
00793 
00794       i = (Index) Maths::rint(pos.x/stepFourierX_);
00795       j = (Index) Maths::rint(pos.y/stepFourierY_);
00796       k = (Index) Maths::rint(pos.z/stepFourierZ_);
00797 
00798       if (i<0)
00799       {
00800         i+=lengthX_;
00801       }
00802 
00803       if (j<0)
00804       {
00805         j+=lengthY_;
00806       }
00807 
00808       if (k<0)
00809       {
00810         k+=lengthZ_;
00811       }
00812 
00813       internalPos = (k + (j + i*lengthY_)*lengthZ_);
00814     }
00815 
00816     if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
00817     {
00818       throw Exception::OutOfGrid(__FILE__, __LINE__);
00819     }
00820     
00821     return TRegularData3D<Complex>::operator[]((Position)internalPos);
00822   }
00823 
00824   template <typename ComplexTraits>
00825   const typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos) const
00826     throw(Exception::OutOfGrid)
00827   {
00828     Index internalPos;
00829 
00830     if (!inFourierSpace_)
00831     {
00832       Index i, j, k;
00833 
00834       i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
00835       j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
00836       k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
00837       
00838       internalPos = (k + (j + i*lengthY_)*lengthZ_);
00839       
00840       /*(Index) rint(       (pos.z+origin_.z)/stepPhysZ_
00841                                   + (   (pos.y+origin_.y)/stepPhysY_
00842                                       + (pos.x+origin_.x)/stepPhysX_*lengthY_ 
00843                                     ) * lengthZ_
00844                                 ); */
00845     }
00846     else
00847     {
00848       Index i, j, k;
00849 
00850       i = (Index) Maths::rint(pos.x/stepFourierX_);
00851       j = (Index) Maths::rint(pos.y/stepFourierY_);
00852       k = (Index) Maths::rint(pos.z/stepFourierZ_);
00853 
00854       if (i<0)
00855       {
00856         i+=lengthX_;
00857       }
00858 
00859       if (j<0)
00860       {
00861         j+=lengthY_;
00862       }
00863 
00864       if (k<0)
00865       {
00866         k+=lengthZ_;
00867       }
00868 
00869       internalPos = (k + (j + i*lengthY_)*lengthZ_);
00870     }
00871 
00872     if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
00873     {
00874       throw Exception::OutOfGrid(__FILE__, __LINE__);
00875     }
00876     
00877     return TRegularData3D<Complex>::operator[]((Position)internalPos);
00878   }
00879 
00880   /*Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos)
00881     throw(Exception::OutOfGrid)
00882   {
00883     return operator [] (pos);
00884   }
00885 
00886   const Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos) const
00887     throw(Exception::OutOfGrid)
00888   {
00889     return operator [] (pos);
00890   }
00891 */
00892   template <typename ComplexTraits>
00893   typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::phase(const Vector3& pos) const
00894   {
00895     
00896     // AR: old version: -2.*M_PI...
00897     double phase = 2.*M_PI*(  (Maths::rint(pos.x/stepFourierX_))*(Maths::rint(origin_.x/stepPhysX_))
00898                               /lengthX_
00899                             + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
00900                               /lengthY_
00901                             + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
00902                               /lengthZ_ );
00903   
00904 
00905     Complex result = Complex(cos(phase), sin(phase));
00906             
00907     return result;
00908     
00909     /*double phase = -2.*M_PI*(  (rint(pos.x/stepFourierX_))*(rint(origin_.x/stepPhysX_))
00910                               /lengthX_
00911                             + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
00912                               /lengthY_
00913                             + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
00914                               /lengthZ_ );
00915   
00916 
00917     Complex result = Complex(cos(phase), sin(phase));
00918             
00919     return result;*/
00920   }
00921 
00922   template <typename ComplexTraits>
00923   bool TFFT3D<ComplexTraits>::isInFourierSpace() const
00924   {
00925     return inFourierSpace_;
00926   }
00927   
00928   template <typename ComplexTraits>
00929   const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator<< 
00930     (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from)
00931   {
00932     // first decide if the TFFT3D data is in Fourier space.
00933     if (!from.isInFourierSpace())
00934     {
00935       // create a new grid
00936       Size lengthX = from.getMaxXIndex()+1;
00937       Size lengthY = from.getMaxYIndex()+1;
00938       Size lengthZ = from.getMaxZIndex()+1;
00939       
00940       TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
00941                                       Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
00942                                       Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
00943 
00944       // and fill it
00945       double normalization=1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
00946       typename TFFT3D<ComplexTraits>::Complex dataIn;
00947       typename TFFT3D<ComplexTraits>::Complex dataOut;
00948       
00949       for (Position i = 0; i < from.size(); i++)
00950       {
00951         Position x, y, z;
00952 
00953         z =  i % lengthZ;
00954         y = (i % (lengthY * lengthZ)) / lengthZ;
00955         x =  i / (lengthY * lengthZ);
00956 
00957         dataIn  = from[i];
00958         dataOut = dataIn;
00959         
00960         newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization;
00961       }
00962 
00963       to = newGrid;
00964 
00965       return to;
00966     }
00967     else
00968     {
00969       // we are in Fourier space
00970       
00971       // create a new grid
00972       Size lengthX = from.getMaxXIndex()+1;
00973       Size lengthY = from.getMaxYIndex()+1;
00974       Size lengthZ = from.getMaxZIndex()+1;
00975       //float stepPhysX = from.getPhysStepWidthX();
00976       //float stepPhysY = from.getPhysStepWidthY();
00977       //float stepPhysZ = from.getPhysStepWidthZ();
00978       float stepFourierX = from.getFourierStepWidthX();
00979       float stepFourierY = from.getFourierStepWidthY();
00980       float stepFourierZ = from.getFourierStepWidthZ();
00981 
00982 
00983       
00984       TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
00985                                       Vector3(from.getFourierSpaceMinX(), 
00986                                               from.getFourierSpaceMinY(),
00987                                               from.getFourierSpaceMinZ()),
00988                                       Vector3(from.getFourierSpaceMaxX(),
00989                                               from.getFourierSpaceMaxY(),
00990                                               from.getFourierSpaceMaxZ()));
00991 
00992       // and fill it
00993       // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
00994       double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
00995       
00996       
00997       Index x, y, z;
00998       Vector3 r;
00999       typename TFFT3D<ComplexTraits>::Complex dataIn;
01000       typename TFFT3D<ComplexTraits>::Complex dataOut;
01001   
01002       for (Position i = 0; i < from.size(); i++)
01003       {
01004         z =  i % lengthZ;
01005         y = (i % (lengthY * lengthZ)) / lengthZ;
01006         x =  i / (lengthY * lengthZ);
01007 
01008         if (x>lengthX/2.)
01009         {
01010           x-=lengthX;
01011         }
01012 
01013         if (y>lengthY/2.)
01014         {
01015           y-=lengthY;
01016         }
01017 
01018         if (z>lengthZ/2.)
01019         {
01020           z-=lengthZ;
01021         }
01022 
01023         r.set((float)x * stepFourierX,
01024               (float)y * stepFourierY,
01025               (float)z * stepFourierZ);
01026 
01027         dataIn = from[i];
01028         dataOut = dataIn;
01029         
01030         newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
01031       }
01032 
01033       to = newGrid;
01034 
01035       return to;
01036     }
01037   }
01038   
01039   template <typename ComplexTraits>
01040   const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from)
01041   {
01042     // first decide if the TFFT3D data is in Fourier space.
01043     if (!from.isInFourierSpace())
01044     {
01045       // create a new grid
01046       Size lengthX = from.getMaxXIndex()+1;
01047       Size lengthY = from.getMaxYIndex()+1;
01048       Size lengthZ = from.getMaxZIndex()+1;
01049       
01050       RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
01051 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
01052 
01053       // and fill it
01054       double normalization = 1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
01055       typename TFFT3D<ComplexTraits>::Complex dataIn;
01056       typename TFFT3D<ComplexTraits>::Complex dataOut;
01057       
01058       for (Position i = 0; i < from.size(); i++)
01059       {
01060         Position x, y, z;
01061 
01062         z =  i % lengthZ;
01063         y = (i % (lengthY * lengthZ)) / lengthZ;
01064         x =  i / (lengthY * lengthZ);
01065 
01066         dataIn  = from[i];
01067         dataOut = dataIn;
01068         
01069         newGrid[x + (y + z*lengthY)*lengthZ] = dataOut.real()*normalization;
01070       }
01071 
01072       to = newGrid;
01073 
01074       return to;
01075     }
01076     else
01077     {
01078       // we are in Fourier space
01079       
01080       // create a new grid
01081       Size lengthX = from.getMaxXIndex()+1;
01082       Size lengthY = from.getMaxYIndex()+1;
01083       Size lengthZ = from.getMaxZIndex()+1;
01084       //float stepPhysX = from.getPhysStepWidthX();
01085       //float stepPhysY = from.getPhysStepWidthY();
01086       //float stepPhysZ = from.getPhysStepWidthZ();
01087       float stepFourierX = from.getFourierStepWidthX();
01088       float stepFourierY = from.getFourierStepWidthY();
01089       float stepFourierZ = from.getFourierStepWidthZ();
01090 
01091 
01092       
01093       RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getFourierSpaceMinX(), from.getFourierSpaceMinY(), from.getFourierSpaceMinZ()), Vector3(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY(), from.getFourierSpaceMaxZ()));
01094 
01095       // and fill it
01096       // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
01097       double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
01098       
01099       Index x, y, z;
01100       signed int xp, yp, zp;
01101       Vector3 r;
01102       typename TFFT3D<ComplexTraits>::Complex dataIn;
01103       typename TFFT3D<ComplexTraits>::Complex dataOut;
01104   
01105       for (Position i = 0; i < from.size(); i++)
01106       {
01107         z =  i % lengthZ;
01108         y = (i % (lengthY * lengthZ)) / lengthZ;
01109         x =  i / (lengthY * lengthZ);
01110         
01111         xp = x;
01112         yp = y;
01113         zp = z;
01114 
01115         if (xp>=lengthX/2.)
01116         {
01117           xp-=(int)lengthX;
01118         }
01119         if (yp>=lengthY/2.)
01120         {
01121           yp-=(int)lengthY;
01122         }
01123         if (zp>=lengthZ/2.)
01124         {
01125           zp-=(int)lengthZ;
01126         }
01127 
01128         if (x>=lengthX/2.)
01129         {
01130           x-=(int)(lengthX/2.);
01131         }
01132         else
01133         {
01134           x+=(int)(lengthX/2.);
01135         }
01136 
01137         if (y>=lengthY/2.)
01138         {
01139           y-=(int)(lengthY/2.);
01140         }
01141         else
01142         {
01143           y+=(int)(lengthY/2.);
01144         }
01145 
01146         if (z>=lengthZ/2.)
01147         {
01148           z-=(int)(lengthZ/2.);
01149         }
01150         else
01151         {
01152           z+=(int)(lengthZ/2.);
01153         }
01154 
01155         r.set((float)xp * stepFourierX,
01156               (float)yp * stepFourierY,
01157               (float)zp * stepFourierZ);
01158 
01159         dataIn = from[i];
01160         dataOut = dataIn;
01161 
01162         newGrid[x + (y + z*lengthY)*lengthZ] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
01163       }
01164 
01165       to = newGrid;
01166 
01167       return to;
01168     }
01169   } 
01170 }
01171 
01172 #endif // BALL_MATHS_TFFT3D_H