r384: Added first vesion of EZPlotDialog

[ctsim.git] / libctsim / imagefile.cpp

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index df0e51ec073981d312bdf265891a909d7125cc3a..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.33 2001/01/02 16:02:13 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
@@ -84,21 +84,29 @@ F64Image::F64Image (void)
 }
 
 void 
-ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
 {
-  ImageFileArray v = getArray();
-  SignalFilter filter (filterName, domainName, bw, filt_param);
-  
-  int iXCenter, iYCenter;
   if (isEven (m_nx))
     iXCenter = m_nx / 2;
   else
     iXCenter = (m_nx - 1) / 2;
+
   if (isEven (m_ny))
     iYCenter = m_ny / 2;
   else
     iYCenter = (m_ny - 1) / 2;
+}
+
+
+void 
+ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+{
+  ImageFileArray v = getArray();
+  SignalFilter filter (filterName, domainName, bw, filt_param);
   
+  unsigned int iXCenter, iYCenter;
+  getCenterCoordinates (iXCenter, iYCenter);
+
   for (unsigned int ix = 0; ix < m_nx; ix++)
     for (unsigned int iy = 0; iy < m_ny; iy++) {
       long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));
@@ -693,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;
         }
@@ -824,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
 
 
@@ -1297,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);