ImLib3D: FFT.cpp Source File

00001 /* ImLib3D
00002  * Copyright (c) 2001, ULP-IPB Strasbourg.
00003  *
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2 of the License, or (at
00007  * your option) any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017  */
00020 #include<ImLib3D/Image3D.hpp>
00021 
00022 #include<ImLib3D/FFT.hpp>
00023 #include<ImLib3D/Display.hpp>
00024 
00025 
00026 
00027 #include<fftw3.h>
00028 
00029 
00030 // void 
00031 // AssertCompatibleComplexType()
00032 // {
00033 //  complex_fftw std(111,222);
00034 //  fftw_complex wtyp;
00035 //  wtyp.re=111;
00036 //  wtyp.im=222;
00037 //  if(sizeof(wtyp)!=sizeof(std)){ThrowError("AssertCompatibleComplexType: incompatible sizes");}
00038 //  if(memcmp((void*)&std,(void*)&wtyp,sizeof(std)))
00039 //  {ThrowError("AssertCompatibleComplexType: incompatible storage");}
00040 // }
00041 
00042 void
00043 IP3D::FFT(const Image3Dcomplexf& src,Image3Dcomplexf &res0, bool resultFullSize)
00044 {
00045     ImageProcessorArgs<Image3Dcomplexf> args(ImArgs(src),res0);
00046     args.SameSize();
00047     Image3Dcomplexf &res=args.GetDest();
00048 
00049     
00050     double t0=DTime();
00051     fftwf_complex *in =(fftwf_complex *)((void *)src.GetData());
00052     fftwf_complex *out=(fftwf_complex *)((void *)res.GetData());
00053     fftwf_plan fftwPlan=fftwf_plan_dft_3d(src.Depth(),src.Height(),src.Width(),
00054                                           in,out,
00055                                           FFTW_FORWARD,FFTW_ESTIMATE);
00056     fftwf_execute(fftwPlan);
00057     fftwf_destroy_plan(fftwPlan);
00058     mprintf("fft complex:float elapsed time:%f\n",DTime(t0));
00059 }
00060 
00061 void
00062 IP3D::FFT(const Image3Dcomplexd& src,Image3Dcomplexd &res0, bool resultFullSize)
00063 {
00064     ImageProcessorArgs<Image3Dcomplexd> args(ImArgs(src),res0);
00065     args.SameSize();
00066     Image3Dcomplexd &res=args.GetDest();
00067 
00068     
00069     double t0=DTime();
00070     fftw_complex *in =(fftw_complex *)((void *)src.GetData());
00071     fftw_complex *out=(fftw_complex *)((void *)res.GetData());
00072     fftw_plan fftwPlan=fftw_plan_dft_3d(src.Depth(),src.Height(),src.Width(),
00073                                           in,out,
00074                                           FFTW_FORWARD,FFTW_ESTIMATE);
00075     fftw_execute(fftwPlan);
00076     fftw_destroy_plan(fftwPlan);
00077     mprintf("fft complex:double elapsed time:%f\n",DTime(t0));
00078 }
00079 
00080 
00081 void
00082 IP3D::FFT(const Image3Df& src, 
00083                                     Image3Dcomplexf &res0,
00084                                     bool resultFullSize)
00085 {
00086     Size3D halfSize(src.Width()/2+1,src.Height(),src.Depth());
00087     // if user wants fullsized result, we need a temp image for half sized
00088     Image3Dcomplexf *tres;
00089     if(resultFullSize){tres=new Image3Dcomplexf(halfSize);}
00090     else{res0.Resize(halfSize);tres=&res0;}
00091     Image3Dcomplexf &res=*tres;
00092 
00093     //FFTOperators::AssertCompatibleComplexType();
00094     double t0=DTime();
00095     float         *in =(float         *)         src.GetData();
00096     fftwf_complex *out=(fftwf_complex *)((void *)res.GetData());
00097     fftwf_plan fftwPlan=fftwf_plan_dft_r2c_3d(src.Depth(),
00098                                               src.Height(),
00099                                               src.Width(),
00100                                               in,out,
00101                                               FFTW_ESTIMATE);
00102     
00103     fftwf_execute(fftwPlan);
00104     mprintf("FFTRealToComplex:: (float) fft elapsed time:%f\n",DTime(t0));
00105 
00106     // if user asked for a full sized result, we must create result from temp
00107     if(resultFullSize)
00108     {
00109         res0.Resize(src.Size());
00110         Image3Dcomplexf::iteratorXYZ p;
00111         for(p=res.begin();p!=res.end();++p)
00112         {
00113             res0(p.Pos())=complexf((*p).real(),(*p).imag());
00114             if(p.x>0)
00115             {
00116                 res0(src.Width()-p.x,p.y,p.z)=complexf((*p).real(),(*p).imag());
00117             }
00118         }
00119         delete tres;
00120     }
00121 }
00122 
00123 
00124 void
00125 IP3D::FFT(const Image3Dd& src, Image3Dcomplexd &res0,bool resultFullSize)
00126 {
00127     Size3D halfSize(src.Width()/2+1,src.Height(),src.Depth());
00128     // if user wants fullsized result, we need a temp image for half sized
00129     Image3Dcomplexd *tres;
00130     if(resultFullSize){tres=new Image3Dcomplexd(halfSize);}
00131     else{res0.Resize(halfSize);tres=&res0;}
00132     Image3Dcomplexd &res=*tres;
00133 
00134     //FFTOperators::AssertCompatibleComplexType();
00135     double t0=DTime();
00136     double       *in =(double       *)         src.GetData();
00137     fftw_complex *out=(fftw_complex *)((void *)res.GetData());
00138     fftw_plan fftwPlan=fftw_plan_dft_r2c_3d(src.Depth(),
00139                                               src.Height(),
00140                                               src.Width(),
00141                                               in,out,
00142                                               FFTW_ESTIMATE);
00143     
00144     fftw_execute(fftwPlan);
00145     mprintf("FFTRealToComplex:: (float) fft elapsed time:%f\n",DTime(t0));
00146 
00147     // if user asked for a full sized result, we must create result from temp
00148     if(resultFullSize)
00149     {
00150         res0.Resize(src.Size());
00151         Image3Dcomplexd::iteratorXYZ p;
00152         for(p=res.begin();p!=res.end();++p)
00153         {
00154             res0(p.Pos())=complexf((*p).real(),(*p).imag());
00155             if(p.x>0)
00156             {
00157                 res0(src.Width()-p.x,p.y,p.z)=complexf((*p).real(),(*p).imag());
00158             }
00159         }
00160         delete tres;
00161     }
00162 }
00163 
00164 
00165 void
00166 IP3D::FFTInverse(const Image3Dcomplexf& src,Image3Dcomplexf &res0,
00167                  int resultWidth,bool overwriteSrc)
00168 {
00169     ImageProcessorArgs<Image3Dcomplexf> args(ImArgs(src),res0);
00170     args.SameSize();
00171     Image3Dcomplexf &res=args.GetDest();
00172 
00173     double t0=DTime();
00174     fftwf_complex *in =(fftwf_complex *)((void *)src.GetData());
00175     fftwf_complex *out=(fftwf_complex *)((void *)res.GetData());
00176     fftwf_plan fftwPlan=fftwf_plan_dft_3d(src.Depth(),src.Height(),src.Width(),
00177                                           in,out,
00178                                           FFTW_BACKWARD,FFTW_ESTIMATE);
00179     fftwf_execute(fftwPlan);
00180     fftwf_destroy_plan(fftwPlan);
00181     mprintf("FFTInverse:: complex<float> -> complex<float> elapsed time:%f\n",
00182             DTime(t0));
00183 }
00184 
00185 
00186 void
00187 IP3D::FFTInverse(const Image3Dcomplexd& src,Image3Dcomplexd &res0,
00188                  int resultWidth,bool overwriteSrc)
00189 {
00190     ImageProcessorArgs<Image3Dcomplexd> args(ImArgs(src),res0);
00191     args.SameSize();
00192     Image3Dcomplexd &res=args.GetDest();
00193 
00194     double t0=DTime();
00195     fftw_complex *in =(fftw_complex *)((void *)src.GetData());
00196     fftw_complex *out=(fftw_complex *)((void *)res.GetData());
00197     fftw_plan fftwPlan=fftw_plan_dft_3d(src.Depth(),src.Height(),src.Width(),
00198                                           in,out,
00199                                           FFTW_BACKWARD,FFTW_ESTIMATE);
00200     fftw_execute(fftwPlan);
00201     fftw_destroy_plan(fftwPlan);
00202     mprintf("FFTInverse:: complex<double> -> complex<double> elapsed time:%f\n",
00203             DTime(t0));
00204 }
00205 
00206 
00207 void 
00208 IP3D::FFTInverse(const Image3Dcomplexf& src0,Image3Df& res,
00209                  int resultWidth,bool overwriteSrc)
00210 {   
00211     // user actually wants the real part of the inverse complex transform
00212     if(resultWidth==src0.Width() || resultWidth==0)
00213     {
00214         Image3Dcomplexf tmp;
00215         IP3D::FFTInverse(src0,tmp);
00216         tmp.ComplexReal(res);
00217         return;
00218     }
00219 
00220     // ok. now we know this is a half-sided transform
00221     // the fftw routine actually overwrites the source
00222     // src: create tmp image if we're not allowed to overwrite
00223     // WARNING: this is a bit hacky... since fftw actually
00224     // modifies its source argument :-(
00225     Image3Dcomplexf *psrc=(Image3Dcomplexf *)&src0;
00226     if(!overwriteSrc){psrc=new Image3Dcomplexf(src0);}
00227     Image3Dcomplexf& src=*psrc;
00228 
00229     //  check if widths are compatible
00230     if(src.Width()!=resultWidth/2+1)
00231     {
00232         ThrowError("FFTInverse:: (ComplexToReal) src.Width incompatible with resultWidth argument. src.Width:%d resultWidth:%d",src.Width(),resultWidth);}
00233     res.Resize(Size3D(resultWidth,src.Height(),src.Depth()));   
00234 
00235     double t0=DTime();
00236     fftwf_complex *in =(fftwf_complex *)((void *)src.GetData());
00237     float         *out=(float         *)((void *)res.GetData());
00238     fftwf_plan fftwPlan=fftwf_plan_dft_c2r_3d(res.Depth(),
00239                                               res.Height(),
00240                                               res.Width(),
00241                                               in,out,
00242                                               FFTW_ESTIMATE);
00243     fftwf_execute(fftwPlan);
00244     fftwf_destroy_plan(fftwPlan);
00245     mprintf("FFTInverse:: ComplexToReal float, elapsed time:%f\n",DTime(t0));
00246 
00247     if(!overwriteSrc){delete psrc;}
00248 }
00249 
00250 void 
00251 IP3D::FFTInverse(const Image3Dcomplexd& src0,Image3Dd& res,
00252                  int resultWidth,bool overwriteSrc)
00253 {   
00254     // user actually wants the real part of the inverse complex transform
00255     if(resultWidth==src0.Width() || resultWidth==0)
00256     {
00257         Image3Dcomplexd tmp;
00258         IP3D::FFTInverse(src0,tmp);
00259         res.Resize(src0.Size());
00260         Image3Dcomplexd::const_iteratorFast pSrc;
00261         Image3Dd::iteratorFast pRes;
00262         for (pSrc = src0.begin(),pRes=res.begin(); pSrc != src0.end(); 
00263              pSrc++,pRes++)
00264         {
00265             *pRes=(*pSrc).real();
00266         }
00267         return;
00268     }
00269 
00270     // ok. now we know this is a half-sided transform
00271     // the fftw routine actually overwrites the source
00272     // src: create tmp image if we're not allowed to overwrite
00273     // WARNING: this is a bit hacky... since fftw actually
00274     // modifies its source argument :-(
00275     Image3Dcomplexd *psrc=(Image3Dcomplexd *)&src0;
00276     if(!overwriteSrc){psrc=new Image3Dcomplexd(src0);}
00277     Image3Dcomplexd& src=*psrc;
00278 
00279     //  check if widths are compatible
00280     if(src.Width()!=resultWidth/2+1)
00281     {
00282         ThrowError("FFTInverse:: (ComplexToReal) src.Width incompatible with resultWidth argument. src.Width:%d resultWidth:%d",src.Width(),resultWidth);}
00283     res.Resize(Size3D(resultWidth,src.Height(),src.Depth()));   
00284 
00285     double t0=DTime();
00286     fftw_complex *in =(fftw_complex *)((void *)src.GetData());
00287     double       *out=(double       *)((void *)res.GetData());
00288     fftw_plan fftwPlan=fftw_plan_dft_c2r_3d(res.Depth(),
00289                                               res.Height(),
00290                                               res.Width(),
00291                                               in,out,
00292                                               FFTW_ESTIMATE);
00293     fftw_execute(fftwPlan);
00294     fftw_destroy_plan(fftwPlan);
00295     mprintf("FFTInverse:: ComplexToReal double, elapsed time:%f\n",DTime(t0));
00296 
00297     if(!overwriteSrc){delete psrc;}
00298 }