r384: Added first vesion of EZPlotDialog

[ctsim.git] / libctsim / imagefile.cpp

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index 946e4db75e17ee766a924939f2237b02a13ee1c7..b4564322c6f32745c75734a64efd559f04f3332d 100644 (file)
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.34 2001/01/04 21:28:41 kevin Exp $
+**  $Id: imagefile.cpp,v 1.36 2001/01/12 16:41:56 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -701,14 +701,16 @@ ImageFile::scaleImage (ImageFile& result) const
             vResultImag[ix][iy] = 0;
         }
       } else {
-        vResult[ix][iy] = vReal[scaleNX][scaleNY] + 
-          dXFrac * (vReal[scaleNX+1][scaleNY] - vReal[scaleNX][scaleNY]) + 
-          dYFrac * (vReal[scaleNX][scaleNY+1] - vReal[scaleNX][scaleNY]);
+        vResult[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vReal[scaleNX][scaleNY] + 
+          dXFrac * (1 - dYFrac) * vReal[scaleNX+1][scaleNY] + 
+          dYFrac * (1 - dXFrac) * vReal[scaleNX][scaleNY+1] +
+          dXFrac * dYFrac * vReal[scaleNX+1][scaleNY+1];
         if (result.isComplex()) {
           if (isComplex())
-            vResultImag[ix][iy] = vImag[scaleNX][scaleNY] + 
-            dXFrac * (vImag[scaleNX+1][scaleNY] - vImag[scaleNX][scaleNY]) + 
-            dYFrac * (vImag[scaleNX][scaleNY+1] - vImag[scaleNX][scaleNY]);
+            vResultImag[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vImag[scaleNX][scaleNY] + 
+          dXFrac * (1 - dYFrac) * vImag[scaleNX+1][scaleNY] + 
+          dYFrac * (1 - dXFrac) * vImag[scaleNX][scaleNY+1] +
+          dXFrac * dYFrac * vImag[scaleNX+1][scaleNY+1];
           else
             vResultImag[ix][iy] = 0;
         }
@@ -832,6 +834,123 @@ ImageFile::ifft (ImageFile& result) const
         
         return true;
 }
+
+bool
+ImageFile::fftRows (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+    return false;
+  }
+  
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.convertRealToComplex ())
+      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_FORWARD, FFTW_IN_PLACE);
+  std::complex<double>* pcRow = new std::complex<double> [m_nx];
+
+  unsigned int ix, iy;
+  unsigned int iArray = 0;
+  for (iy = 0; iy < m_ny; iy++) {
+    for (ix = 0; ix < m_nx; ix++) {
+      in[ix].re = vReal[ix][iy];
+      if (isComplex())
+        in[ix].im = vImag[ix][iy];
+      else
+        in[ix].im = 0;
+    }
+
+    fftw_one (plan, in, NULL);
+
+    for (ix = 0; ix < m_nx; ix++)
+      pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
+
+    Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
+    for (ix = 0; ix < m_nx; ix++) {
+      vReal[ix][iy] = pcRow[ix].real();
+      vImag[ix][iy] = pcRow[ix].imag();
+    }
+  }
+  delete [] pcRow;
+   
+   fftw_destroy_plan (plan);
+   delete in;
+
+   return true;
+}     
+
+bool
+ImageFile::ifftRows (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+    return false;
+  }
+  
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.convertRealToComplex ())
+      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);
+  std::complex<double>* pcRow = new std::complex<double> [m_nx];
+
+  unsigned int ix, iy;
+  unsigned int iArray = 0;
+  for (iy = 0; iy < m_ny; iy++) {
+    for (ix = 0; ix < m_nx; ix++) {
+      double dImag = 0;
+      if (isComplex())
+        dImag = vImag[ix][iy];
+      pcRow[ix] = std::complex<double> (vReal[ix][iy], dImag);
+    }
+
+    Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
+
+    for (ix = 0; ix < m_nx; ix++) {
+      in[ix].re = pcRow[ix].real();
+      in[ix].im = pcRow[ix].imag();
+    }
+
+    fftw_one (plan, in, NULL);
+
+    for (ix = 0; ix < m_nx; ix++) {
+      vReal[ix][iy] = in[ix].re;
+      vImag[ix][iy] = in[ix].im;
+    }
+  }
+  delete [] pcRow;
+  
+   fftw_destroy_plan (plan);
+   delete in;
+
+   return true;
+}
+
+bool
+ImageFile::fftCols (ImageFile& result) const
+{
+  return false;
+}
+
+bool
+ImageFile::ifftCols (ImageFile& result) const
+{
+  return false;
+}
+
 #endif // HAVE_FFTW
 
 
@@ -1305,7 +1424,7 @@ ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, dou
   
   FILE *out;
   if ((out = fopen (outfile,"w")) == NULL) {
-    sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
+    sys_error(ERR_SEVERE, "Error opening output file %s for writing", outfile);
     return false;
   }
   gdImageGif(gif,out);