r7939: Initial working version of fftw3 conversion -- tested only with pjrec

[ctsim.git] / libctsim / imagefile.cpp

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index ef8b1e4bc7d262f4834aa549e2285b409a9a5568..ad13545d10bf918342bc12d28ebbaf3107e3592f 100644 (file)
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -725,7 +725,7 @@ ImageFile::fft (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+  fftw_complex* in = static_cast<fftw_complex*> (fftw_malloc (sizeof(fftw_complex) * m_nx * m_ny));
   
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
@@ -734,17 +734,17 @@ ImageFile::fft (ImageFile& result) const
   unsigned int iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      in[iArray].re = vReal[ix][iy];
+      in[iArray][0] = vReal[ix][iy];
       if (isComplex())
-        in[iArray].im = vImag[ix][iy];
+        in[iArray][1] = vImag[ix][iy];
       else
-        in[iArray].im = 0;
+        in[iArray][1] = 0;
       iArray++;
     }
   }
   
-  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);  
-  fftwnd_one (plan, in, NULL);
+  fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);  
+  fftw_execute (plan);
   
   ImageFileArray vRealResult = result.getArray();
   ImageFileArray vImagResult = result.getImaginaryArray();
@@ -752,13 +752,13 @@ ImageFile::fft (ImageFile& result) const
   unsigned int iScale = m_nx * m_ny;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      vRealResult[ix][iy] = in[iArray].re / iScale;
-      vImagResult[ix][iy] = in[iArray].im / iScale;
+      vRealResult[ix][iy] = in[iArray][0] / iScale;
+      vImagResult[ix][iy] = in[iArray][1] / iScale;
       iArray++;
     }
   }
   delete in;
-  fftwnd_destroy_plan (plan);
+  fftw_destroy_plan (plan);
   
   Fourier::shuffleFourierToNaturalOrder (result);
   
@@ -796,32 +796,31 @@ ImageFile::ifft (ImageFile& result) const
   
   Fourier::shuffleNaturalToFourierOrder (result);
   
-  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx * m_ny));
   
   unsigned int iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      in[iArray].re = vRealResult[ix][iy];
-      in[iArray].im = vImagResult[ix][iy];
+      in[iArray][0] = vRealResult[ix][iy];
+      in[iArray][1] = vImagResult[ix][iy];
       iArray++;
     }
   }
   
-  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   
-  fftwnd_one (plan, in, NULL);
+  fftw_execute (plan);
   
   iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      vRealResult[ix][iy] = in[iArray].re;
-      vImagResult[ix][iy] = in[iArray].im;
+      vRealResult[ix][iy] = in[iArray][0];
+      vImagResult[ix][iy] = in[iArray][1];
       iArray++;
     }
   }
-  fftwnd_destroy_plan (plan);
-  
-  delete in;
+  fftw_destroy_plan (plan);
+  fftw_free(in);
   
   return true;
 }
@@ -842,24 +841,24 @@ ImageFile::fftRows (ImageFile& result) const
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+  fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
 
-  fftw_complex* in = new fftw_complex [m_nx];
   std::complex<double>* pcRow = new std::complex<double> [m_nx];  
   for (unsigned int iy = 0; iy < m_ny; iy++) {
     unsigned int ix;
     for (ix = 0; ix < m_nx; ix++) {
-      in[ix].re = vReal[ix][iy];
+      in[ix][0] = vReal[ix][iy];
       if (isComplex())
-        in[ix].im = vImag[ix][iy];
+        in[ix][1] = vImag[ix][iy];
       else
-        in[ix].im = 0;
+        in[ix][1] = 0;
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (ix = 0; ix < m_nx; ix++)
-      pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
+      pcRow[ix] = std::complex<double>(in[ix][0], in[ix][1]);
     
     Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
     for (ix = 0; ix < m_nx; ix++) {
@@ -870,7 +869,7 @@ ImageFile::fftRows (ImageFile& result) const
   delete [] pcRow;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }     
@@ -888,12 +887,11 @@ ImageFile::ifftRows (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_nx];
-  
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
-  
-  fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+   
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+  fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   std::complex<double>* pcRow = new std::complex<double> [m_nx];
   
   unsigned int ix, iy;
@@ -909,21 +907,21 @@ ImageFile::ifftRows (ImageFile& result) const
     Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
     
     for (ix = 0; ix < m_nx; ix++) {
-      in[ix].re = pcRow[ix].real();
-      in[ix].im = pcRow[ix].imag();
+      in[ix][0] = pcRow[ix].real();
+      in[ix][1] = pcRow[ix].imag();
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (ix = 0; ix < m_nx; ix++) {
-      vReal[ix][iy] = in[ix].re;
-      vImag[ix][iy] = in[ix].im;
+      vReal[ix][iy] = in[ix][0];
+      vImag[ix][iy] = in[ix][1];
     }
   }
   delete [] pcRow;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }
@@ -944,24 +942,24 @@ ImageFile::fftCols (ImageFile& result) const
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+  fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
 
   std::complex<double>* pcCol = new std::complex<double> [m_ny];  
-  fftw_complex* in = new fftw_complex [m_ny];
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     unsigned int iy;
     for (iy = 0; iy < m_ny; iy++) {
-      in[iy].re = vReal[ix][iy];
+      in[iy][0] = vReal[ix][iy];
       if (isComplex())
-        in[iy].im = vImag[ix][iy];
+        in[iy][1] = vImag[ix][iy];
       else
-        in[iy].im = 0;
+        in[iy][1] = 0;
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (iy = 0; iy < m_ny; iy++)
-      pcCol[iy] = std::complex<double>(in[iy].re, in[iy].im);
+      pcCol[iy] = std::complex<double>(in[iy][0], in[iy][1]);
     
     Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
     for (iy = 0; iy < m_ny; iy++) {
@@ -972,7 +970,7 @@ ImageFile::fftCols (ImageFile& result) const
   delete [] pcCol;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }
@@ -990,12 +988,11 @@ ImageFile::ifftCols (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_ny];
-  
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+  fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   std::complex<double>* pcCol = new std::complex<double> [m_ny];
   
   unsigned int ix, iy;
@@ -1011,21 +1008,21 @@ ImageFile::ifftCols (ImageFile& result) const
     Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
     
     for (iy = 0; iy < m_ny; iy++) {
-      in[iy].re = pcCol[iy].real();
-      in[iy].im = pcCol[iy].imag();
+      in[iy][0] = pcCol[iy].real();
+      in[iy][1] = pcCol[iy].imag();
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (iy = 0; iy < m_ny; iy++) {
-      vReal[ix][iy] = in[iy].re;
-      vImag[ix][iy] = in[iy].im;
+      vReal[ix][iy] = in[iy][0];
+      vImag[ix][iy] = in[iy][1];
     }
   }
   delete [] pcCol;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }