r640: no message

[ctsim.git] / libctsim / imagefile.cpp

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index 7c9e91e267b42159474d8b3bbd67abceb330df9c..1c8a063db7c05af76bde6ff2db25df9457bb557c 100644 (file)
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.40 2001/03/07 16:34:47 kevin Exp $
+**  $Id: imagefile.cpp,v 1.41 2001/03/18 18:08:25 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
@@ -36,18 +36,20 @@ const double ImageFile::s_dBlueGrayscaleFactor = 0.114;
 
 
 const int ImageFile::EXPORT_FORMAT_INVALID = -1;
-const int ImageFile::EXPORT_FORMAT_PGM = 0;
-const int ImageFile::EXPORT_FORMAT_PGMASCII = 1;
+const int ImageFile::EXPORT_FORMAT_TEXT = 0;
+const int ImageFile::EXPORT_FORMAT_PGM = 1;
+const int ImageFile::EXPORT_FORMAT_PGMASCII = 2;
 #ifdef HAVE_PNG
-const int ImageFile::EXPORT_FORMAT_PNG = 2;
-const int ImageFile::EXPORT_FORMAT_PNG16 = 3;
+const int ImageFile::EXPORT_FORMAT_PNG = 3;
+const int ImageFile::EXPORT_FORMAT_PNG16 = 4;
 #endif
 #ifdef HAVE_CTN_DICOM
-const int ImageFile::EXPORT_FORMAT_DICOM = 4;
+const int ImageFile::EXPORT_FORMAT_DICOM = 5;
 #endif
 
 const char* ImageFile::s_aszExportFormatName[] = 
 {
+  {"text"},
   {"pgm"},
   {"pgmascii"},
 #ifdef HAVE_PNG
@@ -61,6 +63,7 @@ const char* ImageFile::s_aszExportFormatName[] =
 
 const char* ImageFile::s_aszExportFormatTitle[] = 
 {
+  {"Text"},
   {"PGM"},
   {"PGM ASCII"},
   {"PNG"},
@@ -926,8 +929,8 @@ ImageFile::fftRows (ImageFile& result) const
     
     Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
     for (ix = 0; ix < m_nx; ix++) {
-      vReal[ix][iy] = pcRow[ix].real();
-      vImag[ix][iy] = pcRow[ix].imag();
+      vReal[ix][iy] = pcRow[ix].real() / m_nx;
+      vImag[ix][iy] = pcRow[ix].imag() / m_nx;
     }
   }
   delete [] pcRow;
@@ -994,13 +997,105 @@ ImageFile::ifftRows (ImageFile& result) const
 bool
 ImageFile::fftCols (ImageFile& result) const
 {
-  return false;
+  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_ny];
+  
+  ImageFileArrayConst vReal = getArray();
+  ImageFileArrayConst vImag = getImaginaryArray();
+  
+  fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
+  std::complex<double>* pcCol = new std::complex<double> [m_ny];
+  
+  unsigned int ix, iy;
+  unsigned int iArray = 0;
+  for (ix = 0; ix < m_nx; ix++) {
+    for (iy = 0; iy < m_ny; iy++) {
+      in[iy].re = vReal[ix][iy];
+      if (isComplex())
+        in[iy].im = vImag[ix][iy];
+      else
+        in[iy].im = 0;
+    }
+    
+    fftw_one (plan, in, NULL);
+    
+    for (iy = 0; iy < m_ny; iy++)
+      pcCol[iy] = std::complex<double>(in[iy].re, in[iy].im);
+    
+    Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
+    for (iy = 0; iy < m_ny; iy++) {
+      vReal[ix][iy] = pcCol[iy].real() / m_ny;
+      vImag[ix][iy] = pcCol[iy].imag() / m_ny;
+    }
+  }
+  delete [] pcCol;
+  
+  fftw_destroy_plan (plan);
+  delete in;
+  
+  return true;
 }
 
 bool
 ImageFile::ifftCols (ImageFile& result) const
 {
-  return false;
+  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_ny];
+  
+  ImageFileArrayConst vReal = getArray();
+  ImageFileArrayConst vImag = getImaginaryArray();
+  
+  fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
+  std::complex<double>* pcCol = new std::complex<double> [m_ny];
+  
+  unsigned int ix, iy;
+  unsigned int iArray = 0;
+  for (ix = 0; ix < m_nx; ix++) {
+    for (iy = 0; iy < m_ny; iy++) {
+      double dImag = 0;
+      if (isComplex())
+        dImag = vImag[ix][iy];
+      pcCol[iy] = std::complex<double> (vReal[ix][iy], dImag);
+    }
+    
+    Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
+    
+    for (iy = 0; iy < m_ny; iy++) {
+      in[iy].re = pcCol[iy].real();
+      in[iy].im = pcCol[iy].imag();
+    }
+    
+    fftw_one (plan, in, NULL);
+    
+    for (iy = 0; iy < m_ny; iy++) {
+      vReal[ix][iy] = in[iy].re;
+      vImag[ix][iy] = in[iy].im;
+    }
+  }
+  delete [] pcCol;
+  
+  fftw_destroy_plan (plan);
+  delete in;
+  
+  return true;
 }
 
 #endif // HAVE_FFTW
@@ -1564,10 +1659,14 @@ ImageFile::exportImage (const char* const pszFormat, const char* const pszFilena
     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);
   else if (iFormatID == EXPORT_FORMAT_PGMASCII)
     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);
+  else if (iFormatID == EXPORT_FORMAT_TEXT)
+    return writeImageText (pszFilename);
+#ifdef HAVE_PNG
   else if (iFormatID == EXPORT_FORMAT_PNG)
     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);
   else if (iFormatID == EXPORT_FORMAT_PNG16)
     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);
+#endif
 #ifdef HAVE_CTN_DICOM
   else if (iFormatID == EXPORT_FORMAT_DICOM) {
     DicomExporter dicomExport (this);
@@ -1660,6 +1759,36 @@ ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell
   return true;
 }
 
+bool
+ImageFile::writeImageText (const char* const outfile)
+{
+  FILE *fp;
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  ImageFileArray vImag = getImaginaryArray();
+
+  if ((fp = fopen (outfile, "w")) == NULL)
+    return false;
+  
+  for (int irow = ny - 1; irow >= 0; irow--) {
+    for (int icol = 0; icol < nx; icol++) {
+      if (isComplex()) {
+        if (vImag[icol][irow] >= 0)
+          fprintf (fp, "%.9g+%.9gi ", v[icol][irow], vImag[icol][irow]);
+        else
+          fprintf (fp, "%.9g-%.9gi ", v[icol][irow], -vImag[icol][irow]);
+      } else
+        fprintf (fp, "%12.8g ", v[icol][irow]);
+    }
+    fprintf(fp, "\n");
+  }
+  
+  fclose(fp);
+  
+  return true;
+}
+
 
 #ifdef HAVE_PNG
 bool