BALL: FFT3D.h Source File

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