r326: FFTW additions, filter image generation
authorKevin M. Rosenberg <kevin@rosenberg.net>
Mon, 1 Jan 2001 10:14:34 +0000 (10:14 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Mon, 1 Jan 2001 10:14:34 +0000 (10:14 +0000)
17 files changed:
ChangeLog
include/array2dfile.h
include/ctsupport.h
include/imagefile.h
include/procsignal.h
libctsim/array2dfile.cpp
libctsim/backprojectors.cpp
libctsim/imagefile.cpp
libctsim/procsignal.cpp
libctsupport/mathfuncs.cpp
msvc/ctsim/ctsim.plg
src/ctsim.cpp
src/ctsim.h
src/dialogs.cpp
src/dialogs.h
src/views.cpp
src/views.h

index cf938a0..9decb06 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,7 +1,16 @@
-TODO
+       TODO
        Read PlotFile's.
        Add to ctsim export of imagefiles to standard graphic formats.
-       Use FFTW library in imagefile.cpp
+
+
+3.0alpha2 
+       
+       * processsignal.cpp: Fixed "off by one" bug in
+       shuffleNaturalToFourierOrder when n is even.
+
+       * imagefile: Added FFTW library to imagefile processing
+
+       * ctsim: added generation of filter images
        
 3.0alpha1 - Released 12/29/00
 
index 754ad93..22e8b61 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: array2dfile.h,v 1.13 2000/12/29 19:30:08 kevin Exp $
+**  $Id: array2dfile.h,v 1.14 2001/01/01 10:14:34 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
@@ -165,6 +165,12 @@ public:
   kuint32 ny (void) const
       { return m_ny; }
 \r
+  bool isComplex() const\r
+  { return m_dataType == DATA_TYPE_COMPLEX; }\r
+\r
+  bool isReal() const\r
+  { return m_dataType == DATA_TYPE_REAL; }\r
+\r
   int dataType () const\r
   { return static_cast<int>(m_dataType); }\r
 \r
index 81479c9..84eef25 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsupport.h,v 1.19 2000/12/29 20:04:02 kevin Exp $
+**  $Id: ctsupport.h,v 1.20 2001/01/01 10:14:34 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -184,7 +184,21 @@ convertRadiansToDegrees (double x)
 template<class T>
 inline T nearest (double x)
 { return (x > 0 ? static_cast<T>(x+0.5) : static_cast<T>(x-0.5)); }
-
+\r
+inline bool isEven (int n)\r
+{ return (n % 2) == 0; }\r
+\r
+inline bool isOdd (int n)\r
+{ return (n % 2) != 0; }\r
+\r
+#if 0\r
+inline bool isEven (long n)\r
+{ return (n % 2) == 0; }\r
+\r
+inline bool isOdd (long n)\r
+{ return (n % 2) != 0; }\r
+#endif\r
+\r
 inline int imax (int a, int b)\r
 { return (a >= b ? a : b); }\r
 \r
index 57d093f..945bdbf 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.h,v 1.25 2000/12/29 19:30:08 kevin Exp $
+**  $Id: imagefile.h,v 1.26 2001/01/01 10:14:34 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
@@ -135,7 +135,7 @@ class ImageFile : public ImageFileBase
   void filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param);
 
   void statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const;\r
-  void statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const;\r
+  void statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const;\r
   void getMinMax (double& min, double& max) const;
   void printStatistics (std::ostream& os) const;
   bool comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const;
@@ -153,6 +153,10 @@ class ImageFile : public ImageFileBase
   bool exp (ImageFile& result) const;\r
   bool fourier (ImageFile& result) const;\r
   bool inverseFourier (ImageFile& result) const;\r
+#ifdef HAVE_FFTW\r
+  bool fft (ImageFile& result) const;\r
+  bool ifft (ImageFile& result) const;\r
+#endif\r
   bool magnitude (ImageFile& result) const;\r
   bool phase (ImageFile& result) const;\r
 \r
index a160154..8eb98c6 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.h,v 1.9 2000/12/29 15:45:06 kevin Exp $
+**  $Id: procsignal.h,v 1.10 2001/01/01 10:14:34 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
@@ -41,6 +41,9 @@
 
 class SignalFilter;
 class SGP;
+\r
+typedef std::complex<double> CTSimComplex;\r
+\r
 
 class ProcessSignal {
  public:
@@ -104,10 +107,13 @@ class ProcessSignal {
     static void finiteFourierTransform (const std::complex<double> input[], double output[], const int n, const int direction);
 
 
+    static void shuffleNaturalToFourierOrder (float* pdVector, const int n);\r
     static void shuffleNaturalToFourierOrder (double* pdVector, const int n);\r
     static void shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n);\r
+    static void shuffleFourierToNaturalOrder (float* pdVector, const int n);\r
     static void shuffleFourierToNaturalOrder (double* pdVector, const int n);\r
     static void shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n);\r
+\r
 
  private:
         std::string m_nameFilterMethod;
index ca45d8b..d0f2d68 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: array2dfile.cpp,v 1.22 2000/12/29 19:30:08 kevin Exp $
+**  $Id: array2dfile.cpp,v 1.23 2001/01/01 10:14:34 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
@@ -410,7 +410,7 @@ Array2dFile::headerRead (frnetorderstream& fs)
     return false;
   }
   
-  return true;
+  return ! fs.fail();
 }
 
 
@@ -448,7 +448,7 @@ Array2dFile::headerWrite (frnetorderstream& fs)
   fs.seekp (0);
   fs.writeInt16 (m_headersize);
   
-  return true;
+  return ! fs.fail();
 }
 
 
index 7926b0d..726e5b8 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.20 2000/12/18 09:31:26 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.21 2001/01/01 10:14:34 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,7 +84,7 @@ const char* Backprojector::s_aszInterpTitle[] =
 {
   {"Nearest"},
   {"Linear"},
-  {"Frequency Preinterpolationj"},
+  {"Frequency Preinterpolation"},
 #if HAVE_BSPLINE_INTERP
   {"B-Spline"},
   {"B-Spline 1st Order"},
index b6525e1..df647ce 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.27 2000/12/29 20:15:37 kevin Exp $
+**  $Id: imagefile.cpp,v 1.28 2001/01/01 10:14:34 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,8 +36,8 @@ F32Image::F32Image (int nx, int ny, int dataType)
 F32Image::F32Image (void)\r
 : Array2dFile()\r
 {\r
-       setPixelFormat (Array2dFile::PIXEL_FLOAT32);\r
-       setPixelSize (sizeof(kfloat32));\r
+  setPixelFormat (Array2dFile::PIXEL_FLOAT32);\r
+  setPixelSize (sizeof(kfloat32));\r
   setDataType (Array2dFile::DATA_TYPE_REAL);\r
 }\r
 \r
@@ -49,78 +49,98 @@ F64Image::F64Image (int nx, int ny, int dataType)
 F64Image::F64Image (void)\r
 : Array2dFile ()\r
 {\r
-       setPixelFormat (PIXEL_FLOAT64);\r
-       setPixelSize (sizeof(kfloat64));\r
+  setPixelFormat (PIXEL_FLOAT64);\r
+  setPixelSize (sizeof(kfloat64));\r
   setDataType (Array2dFile::DATA_TYPE_REAL);\r
 }\r
 
 void 
 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param)
-{
-       int hx = (m_nx - 1) / 2;
-       int hy = (m_ny - 1) / 2;
-       ImageFileArray v = getArray();
-       SignalFilter filter (filterName, domainName, bw, filt_param);
-       
-       for (int i = -hx; i <= hx; i++) {
-               for (int j = -hy; j <= hy; j++) {
+{\r
+  ImageFileArray v = getArray();\r
+  SignalFilter filter (filterName, domainName, bw, filt_param);\r
+\r
+#if 1\r
+  int iXCenter, iYCenter;\r
+  if (isEven (m_nx))\r
+    iXCenter = m_nx / 2;\r
+  else\r
+    iXCenter = (m_nx - 1) / 2;\r
+  if (isEven (m_ny))\r
+    iYCenter = m_ny / 2;\r
+  else\r
+    iYCenter = (m_ny - 1) / 2;\r
+\r
+  for (int ix = 0; ix < m_nx; ix++)\r
+    for (int iy = 0; iy < m_ny; iy++) {\r
+      long lD2 = (ix - iXCenter) * (ix - iXCenter) + (iy - iYCenter) * (iy - iYCenter);\r
+      double r = ::sqrt (static_cast<double>(lD2));\r
+      v[ix][iy] = filter.response (r);\r
+    }\r
+#else
+  int hx = (m_nx - 1) / 2;
+  int hy = (m_ny - 1) / 2;
+  
+  for (int i = -hx; i <= hx; i++) {
+    for (int j = -hy; j <= hy; j++) {
       double r = ::sqrt (i * i + j * j);
-                       
-                       v[i+hx][j+hy] = filter.response (r);
-               }
-       }
+      
+      v[i+hx][j+hy] = filter.response (r);
+    }
+  }\r
+#endif
 }
 
 int
 ImageFile::display (void) const
 {
-    double pmin, pmax;
-       
-    getMinMax (pmin, pmax);
-       
-    return (displayScaling (1, pmin, pmax));
+  double pmin, pmax;
+  
+  getMinMax (pmin, pmax);
+  
+  return (displayScaling (1, pmin, pmax));
 }
 
 int 
 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
 {
-    int nx = m_nx;
-    int ny = m_ny;
-    ImageFileArrayConst v = getArray();
-    if (v == NULL || nx == 0 || ny == 0)
-               return 0;
-       
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArrayConst v = getArray();
+  if (v == NULL || nx == 0 || ny == 0)
+    return 0;
+  
 #if HAVE_G2_H
-    int* pPens = new int [nx * ny * scale * scale ];
-       
-    double view_scale = 255 / (pmax - pmin);
-    int id_X11 = g2_open_X11 (nx * scale, ny * scale);
-    int grayscale[256];
-    for (int i = 0; i < 256; i++) {
-               double cval = i / 255.;
-               grayscale[i] = g2_ink (id_X11, cval, cval, cval);
-    }
-       
-    for (int iy = ny - 1; iy >= 0; iy--) {
-               int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
-               for (int ix = 0; ix < nx; ix++) {
-                       int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
-                       if (cval < 0)  
-                               cval = 0;
-                       else if (cval > 255) 
-                               cval = 255;
-                       for (int sy = 0; sy < scale; sy++)
-                               for (int sx = 0; sx < scale; sx++)
-                                       pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
-               }
+  int* pPens = new int [nx * ny * scale * scale ];
+  
+  double view_scale = 255 / (pmax - pmin);
+  int id_X11 = g2_open_X11 (nx * scale, ny * scale);
+  int grayscale[256];
+  for (int i = 0; i < 256; i++) {
+    double cval = i / 255.;
+    grayscale[i] = g2_ink (id_X11, cval, cval, cval);
+  }
+  
+  for (int iy = ny - 1; iy >= 0; iy--) {
+    int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
+    for (int ix = 0; ix < nx; ix++) {
+      int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
+      if (cval < 0)  
+        cval = 0;
+      else if (cval > 255) 
+        cval = 255;
+      for (int sy = 0; sy < scale; sy++)
+        for (int sx = 0; sx < scale; sx++)
+          pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
     }
-       
-    g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
-       
-    delete pPens;
-    return (id_X11);
+  }
+  
+  g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
+  
+  delete pPens;
+  return (id_X11);
 #else
-    return 0;
+  return 0;
 #endif
 }
 
@@ -139,102 +159,102 @@ ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const Ima
 bool
 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
 {
-    if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
-               sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
-               return false;
-    }
-    ImageFileArrayConst v = getArray();
-    if (v == NULL || m_nx == 0 || m_ny == 0)
-               return false;
-       
-    ImageFileArrayConst vComp = imComp.getArray();
-       
-    double myMean = 0.;
-    for (unsigned int ix = 0; ix < m_nx; ix++) {
-               for (unsigned int iy = 0; iy < m_ny; iy++) {
-                       myMean += v[ix][iy];
-               }
+  if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
+    sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
+    return false;
+  }
+  ImageFileArrayConst v = getArray();
+  if (v == NULL || m_nx == 0 || m_ny == 0)
+    return false;
+  
+  ImageFileArrayConst vComp = imComp.getArray();
+  
+  double myMean = 0.;
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      myMean += v[ix][iy];
     }
-    myMean /= (m_nx * m_ny);
-       
-    double sqErrorSum = 0.;
-    double absErrorSum = 0.;
-    double sqDiffFromMeanSum = 0.;
-    double absValueSum = 0.;
-    for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
-               for (unsigned int iy = 0; iy < m_ny; iy++) {
-                       double diff = v[ix2][iy] - vComp[ix2][iy];
-                       sqErrorSum += diff * diff;
-                       absErrorSum += fabs(diff);
-                       double diffFromMean = v[ix2][iy] - myMean;
-                       sqDiffFromMeanSum += diffFromMean * diffFromMean;
-                       absValueSum += fabs(v[ix2][iy]);
-               }
+  }
+  myMean /= (m_nx * m_ny);
+  
+  double sqErrorSum = 0.;
+  double absErrorSum = 0.;
+  double sqDiffFromMeanSum = 0.;
+  double absValueSum = 0.;
+  for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      double diff = v[ix2][iy] - vComp[ix2][iy];
+      sqErrorSum += diff * diff;
+      absErrorSum += fabs(diff);
+      double diffFromMean = v[ix2][iy] - myMean;
+      sqDiffFromMeanSum += diffFromMean * diffFromMean;
+      absValueSum += fabs(v[ix2][iy]);
     }
-       
-    d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
-    r = absErrorSum / absValueSum;
-       
-    int hx = m_nx / 2;
-    int hy = m_ny / 2;
-    double eMax = -1;
-    for (int ix3 = 0; ix3 < hx; ix3++) {
-               for (int iy = 0; iy < hy; iy++) {
-                       double avgPixel = 0.25 * (v[2*ix3][2*iy] + v[2*ix3+1][2*iy] + v[2*ix3][2*iy+1] + v[2*ix3+1][2*iy+1]);
-                       double avgPixelComp = 0.25 * (vComp[2*ix3][2*iy] + vComp[2*ix3+1][2*iy] + vComp[2*ix3][2*iy+1] + vComp[2*ix3+1][2*iy+1]);
-                       double error = fabs (avgPixel - avgPixelComp);
-                       if (error > eMax)
-                               eMax = error;
-               }
+  }
+  
+  d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
+  r = absErrorSum / absValueSum;
+  
+  int hx = m_nx / 2;
+  int hy = m_ny / 2;
+  double eMax = -1;
+  for (int ix3 = 0; ix3 < hx; ix3++) {
+    for (int iy = 0; iy < hy; iy++) {
+      double avgPixel = 0.25 * (v[2*ix3][2*iy] + v[2*ix3+1][2*iy] + v[2*ix3][2*iy+1] + v[2*ix3+1][2*iy+1]);
+      double avgPixelComp = 0.25 * (vComp[2*ix3][2*iy] + vComp[2*ix3+1][2*iy] + vComp[2*ix3][2*iy+1] + vComp[2*ix3+1][2*iy+1]);
+      double error = fabs (avgPixel - avgPixelComp);
+      if (error > eMax)
+        eMax = error;
     }
-       
-    e = eMax;
-       
-    return true;
+  }
+  
+  e = eMax;
+  
+  return true;
 }
 
 
 bool
 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
 {
-       double d, r, e;
-       
-       if (comparativeStatistics (imComp, d, r, e)) {
-               os << "  Normalized root mean squared distance (d): " << d << std::endl;
-               os << "      Normalized mean absolute distance (r): " << r << std::endl;
-               os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
-               return true;
-       }
-       return false;
+  double d, r, e;
+  
+  if (comparativeStatistics (imComp, d, r, e)) {
+    os << "  Normalized root mean squared distance (d): " << d << std::endl;
+    os << "      Normalized mean absolute distance (r): " << r << std::endl;
+    os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
+    return true;
+  }
+  return false;
 }
 
 
 void
 ImageFile::printStatistics (std::ostream& os) const
 {
-    double min, max, mean, mode, median, stddev;
-       
-    statistics (min, max, mean, mode, median, stddev);\r
-    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
-      os << "Real Component Statistics" << std::endl;\r
-
-    os << "   min: " << min << std::endl;
-    os << "   max: " << max << std::endl;
-    os << "  mean: " << mean << std::endl;
-    os << "  mode: " << mode << std::endl;
-    os << "median: " << median << std::endl;
-    os << "stddev: " << stddev << std::endl;\r
-\r
-    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
-      statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);\r
-      os << std::endl << "Imaginary Component Statistics" << std::endl;
+  double min, max, mean, mode, median, stddev;
+  
+  statistics (min, max, mean, mode, median, stddev);\r
+  if (isComplex())\r
+    os << "Real Component Statistics" << std::endl;\r
+  
+  os << "   min: " << min << std::endl;
+  os << "   max: " << max << std::endl;
+  os << "  mean: " << mean << std::endl;
+  os << "  mode: " << mode << std::endl;
+  os << "median: " << median << std::endl;
+  os << "stddev: " << stddev << std::endl;\r
+  \r
+  if (isComplex()) {\r
+    statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);\r
+    os << std::endl << "Imaginary Component Statistics" << std::endl;
     os << "   min: " << min << std::endl;\r
     os << "   max: " << max << std::endl;\r
     os << "  mean: " << mean << std::endl;\r
     os << "  mode: " << mode << std::endl;\r
     os << "median: " << median << std::endl;\r
     os << "stddev: " << stddev << std::endl;\r
-    }\r
+  }\r
 }
 
 
@@ -247,45 +267,45 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou
 
 
 void\r
-ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const\r
+ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const\r
 {\r
-    int nx = m_nx;\r
-    int ny = m_ny;\r
-    \r
-    if (v == NULL || nx == 0 || ny == 0)\r
-               return;\r
-\r
-       std::vector<double> vecImage;\r
-       int iVec = 0;\r
-       vecImage.resize (nx * ny);\r
-    for (int ix = 0; ix < nx; ix++) {\r
-               for (int iy = 0; iy < ny; iy++)\r
-                       vecImage[iVec++] = v[ix][iy];\r
-       }\r
-\r
-       vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);\r
+  int nx = m_nx;\r
+  int ny = m_ny;\r
+  \r
+  if (v == NULL || nx == 0 || ny == 0)\r
+    return;\r
+  \r
+  std::vector<double> vecImage;\r
+  int iVec = 0;\r
+  vecImage.resize (nx * ny);\r
+  for (int ix = 0; ix < nx; ix++) {\r
+    for (int iy = 0; iy < ny; iy++)\r
+      vecImage[iVec++] = v[ix][iy];\r
+  }\r
+  \r
+  vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);\r
 }\r
 \r
 void
 ImageFile::getMinMax (double& min, double& max) const
 {
-    int nx = m_nx;
-    int ny = m_ny;
-    ImageFileArrayConst v = getArray();
-    
-    if (v == NULL || nx == 0 || ny == 0)
-               return;
-       
-    min = v[0][0];
-    max = v[0][0];
-    for (int ix = 0; ix < nx; ix++) {
-               for (int iy = 0; iy < ny; iy++) {
-                       if (v[ix][iy] > max)
-                               max = v[ix][iy];
-                       if (v[ix][iy] < min)
-                               min = v[ix][iy];
-               }
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArrayConst v = getArray();
+  
+  if (v == NULL || nx == 0 || ny == 0)
+    return;
+  
+  min = v[0][0];
+  max = v[0][0];
+  for (int ix = 0; ix < nx; ix++) {
+    for (int iy = 0; iy < ny; iy++) {
+      if (v[ix][iy] > max)
+        max = v[ix][iy];
+      if (v[ix][iy] < min)
+        min = v[ix][iy];
     }
+  }
 }
 \r
 bool\r
@@ -293,17 +313,17 @@ ImageFile::convertRealToComplex ()
 {\r
   if (dataType() != Array2dFile::DATA_TYPE_REAL)\r
     return false;\r
-\r
+  \r
   if (! reallocRealToComplex())\r
     return false;\r
-\r
+  \r
   ImageFileArray vImag = getImaginaryArray();\r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumn vCol = vImag[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
       *vCol++ = 0;\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
@@ -312,19 +332,19 @@ ImageFile::convertComplexToReal ()
 {\r
   if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)\r
     return false;\r
-\r
+  \r
   ImageFileArray vReal = getArray();\r
   ImageFileArray vImag = getImaginaryArray();\r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumn vRealCol = vReal[ix];\r
     ImageFileColumn vImagCol = vImag[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
-      std::complex<double> c (*vRealCol, *vImagCol);\r
+      CTSimComplex c (*vRealCol, *vImagCol);\r
       *vRealCol++ = std::abs (c);\r
       vImagCol++;\r
     }\r
   }\r
-\r
+  \r
   return reallocComplexToReal();\r
 }\r
 \r
@@ -335,20 +355,20 @@ ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArrayConst vRHS = rRHS.getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in1 = vLHS[ix];\r
     ImageFileColumnConst in2 = vRHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        *out++ = *in1++ - *in2++;\r
+      *out++ = *in1++ - *in2++;\r
   }\r
-\r
-    return true;\r
+  \r
+  return true;\r
 }\r
 
 bool\r
@@ -358,20 +378,20 @@ ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArrayConst vRHS = rRHS.getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in1 = vLHS[ix];\r
     ImageFileColumnConst in2 = vRHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        *out++ = *in1++ + *in2++;\r
+      *out++ = *in1++ + *in2++;\r
   }\r
-\r
-    return true;\r
+  \r
+  return true;\r
 }\r
 \r
 bool\r
@@ -381,20 +401,20 @@ ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArrayConst vRHS = rRHS.getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in1 = vLHS[ix];\r
     ImageFileColumnConst in2 = vRHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        *out++ = *in1++ * *in2++;\r
+      *out++ = *in1++ * *in2++;\r
   }\r
-\r
-    return true;\r
+  \r
+  return true;\r
 }\r
 \r
 bool\r
@@ -404,11 +424,11 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArrayConst vRHS = rRHS.getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in1 = vLHS[ix];\r
     ImageFileColumnConst in2 = vRHS[ix];\r
@@ -420,8 +440,8 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
         *out++ = 0;\r
     }\r
   }\r
-\r
-    return true;\r
+  \r
+  return true;\r
 }\r
 \r
 \r
@@ -432,17 +452,17 @@ ImageFile::invertPixelValues (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in = vLHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        *out++ = - *in++;\r
+      *out++ = - *in++;\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
@@ -453,10 +473,10 @@ ImageFile::sqrt (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in = vLHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
@@ -466,7 +486,7 @@ ImageFile::sqrt (ImageFile& result) const
       else\r
         *out++ = ::sqrt(*in++);\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
@@ -477,10 +497,10 @@ ImageFile::log (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in = vLHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
@@ -490,7 +510,7 @@ ImageFile::log (ImageFile& result) const
       else\r
         *out++ = ::log(*in++);\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
@@ -501,33 +521,153 @@ ImageFile::exp (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in = vLHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++)\r
       *out++ = ::exp (*in++);\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
+namespace ProcessImage {\r
+\r
+void\r
+shuffleFourierToNaturalOrder (ImageFile& im)\r
+{\r
+  ImageFileArray vReal = im.getArray();\r
+  ImageFileArray vImag = im.getImaginaryArray();\r
+  unsigned int ix, iy;\r
+  unsigned int nx = im.nx();\r
+  unsigned int ny = im.ny();\r
+\r
+  // shuffle each column\r
+  for (ix = 0; ix < nx; ix++) {\r
+    ProcessSignal::shuffleFourierToNaturalOrder (vReal[ix], ny);\r
+    if (im.isComplex())\r
+      ProcessSignal::shuffleFourierToNaturalOrder (vImag[ix], ny);\r
+  }\r
+\r
+  // shuffle each row\r
+  float* pRow = new float [nx];\r
+  for (iy = 0; iy < ny; iy++) {\r
+    for (ix = 0; ix < nx; ix++)\r
+      pRow[ix] = vReal[ix][iy];\r
+    ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);\r
+    for (ix = 0; ix < nx; ix++)\r
+      vReal[ix][iy] = pRow[ix];\r
+    if (im.isComplex()) {\r
+      for (ix = 0; ix < nx; ix++)\r
+        pRow[ix] = vImag[ix][iy];\r
+      ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);;\r
+      for (ix = 0; ix < nx; ix++)\r
+        vImag[ix][iy] = pRow[ix];\r
+    }\r
+  }\r
+  delete pRow;\r
+}\r
\r
+void\r
+shuffleNaturalToFourierOrder (ImageFile& im)\r
+{\r
+  ImageFileArray vReal = im.getArray();\r
+  ImageFileArray vImag = im.getImaginaryArray();\r
+  unsigned int ix, iy;\r
+  unsigned int nx = im.nx();\r
+  unsigned int ny = im.ny();\r
+\r
+  // shuffle each x column\r
+  for (ix = 0; ix < nx; ix++) {\r
+    ProcessSignal::shuffleNaturalToFourierOrder (vReal[ix], ny);\r
+    if (im.isComplex())\r
+      ProcessSignal::shuffleNaturalToFourierOrder (vImag[ix], ny);\r
+  }\r
+\r
+  // shuffle each y row\r
+  float* pRow = new float [nx];\r
+  for (iy = 0; iy < ny; iy++) {\r
+    for (ix = 0; ix < nx; ix++)\r
+      pRow[ix] = vReal[ix][iy];\r
+    ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);\r
+    for (ix = 0; ix < nx; ix++)\r
+      vReal[ix][iy] = pRow[ix];\r
+    if (im.isComplex()) {\r
+      for (ix = 0; ix < nx; ix++)\r
+        pRow[ix] = vImag[ix][iy];\r
+      ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);\r
+      for (ix = 0; ix < nx; ix++)\r
+        vImag[ix][iy] = pRow[ix];\r
+    }\r
+  }\r
+  delete [] pRow;\r
+}\r
+\r
\r
+}; // namespace ProcessIamge\r
+\r
+#ifdef HAVE_FFTW\r
 bool\r
-ImageFile::inverseFourier (ImageFile& result) const\r
+ImageFile::fft (ImageFile& result) const\r
 {\r
   if (m_nx != result.nx() || m_ny != result.ny()) {\r
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
 \r
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
+    if (! result.reallocRealToComplex ())\r
+      return false;\r
+  }\r
+\r
+  fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
+\r
+  ImageFileArrayConst vReal = getArray();\r
+  ImageFileArrayConst vImag = getImaginaryArray();\r
+\r
+  unsigned int ix, iy;\r
+  unsigned int iArray = 0;\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      in[iArray].re = vReal[ix][iy];\r
+      if (isComplex())\r
+        in[iArray].im = vImag[ix][iy];\r
+      else\r
+        in[iArray].im = 0;\r
+      iArray++;\r
+    }\r
+\r
+  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);\r
+\r
+  fftwnd_one (plan, in, NULL);\r
+\r
+  ImageFileArray vRealResult = result.getArray();\r
+  ImageFileArray vImagResult = result.getImaginaryArray();\r
+  iArray = 0;\r
+  unsigned int iScale = m_nx * m_ny;\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      vRealResult[ix][iy] = in[iArray].re / iScale;\r
+      vImagResult[ix][iy] = in[iArray].im / iScale;\r
+      iArray++;\r
+    }\r
+\r
+  fftwnd_destroy_plan (plan);\r
+  delete in;\r
+\r
+\r
+  ProcessImage::shuffleFourierToNaturalOrder (result);\r
+\r
   return true;\r
 }\r
 \r
+\r
 bool\r
-ImageFile::fourier (ImageFile& result) const\r
+ImageFile::ifft (ImageFile& result) const\r
 {\r
   if (m_nx != result.nx() || m_ny != result.ny()) {\r
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
@@ -539,62 +679,230 @@ ImageFile::fourier (ImageFile& result) const
       return false;\r
   }\r
 \r
-  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vReal = getArray();\r
+  ImageFileArrayConst vImag = getImaginaryArray();\r
   ImageFileArray vRealResult = result.getArray();\r
   ImageFileArray vImagResult = result.getImaginaryArray();\r
-\r
   unsigned int ix, iy;\r
-  double* pY = new double [m_ny];\r
-  std::complex<double>** complexOut = new std::complex<double>* [m_nx];\r
   for (ix = 0; ix < m_nx; ix++)\r
-    complexOut[ix] = new std::complex<double> [m_ny];\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      vRealResult[ix][iy] = vReal[ix][iy];\r
+      if (isComplex()) \r
+        vImagResult[ix][iy] = vImag[ix][iy];\r
+      else\r
+        vImagResult[ix][iy] = 0;\r
+    }\r
+\r
+  ProcessImage::shuffleNaturalToFourierOrder (result);\r
 \r
+  fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
+\r
+  unsigned int iArray = 0;\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      in[iArray].re = vRealResult[ix][iy];\r
+      in[iArray].im = vImagResult[ix][iy];\r
+      iArray++;\r
+    }\r
+\r
+  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);\r
+\r
+  fftwnd_one (plan, in, NULL);\r
+\r
+  iArray = 0;\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      vRealResult[ix][iy] = in[iArray].re;\r
+      vImagResult[ix][iy] = in[iArray].im;\r
+      iArray++;\r
+    }\r
+\r
+  fftwnd_destroy_plan (plan);\r
+\r
+  delete in;\r
+\r
+  return true;\r
+}\r
+#endif // HAVE_FFTW\r
+\r
+\r
+  \r
+bool\r
+ImageFile::fourier (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+  \r
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
+    if (! result.reallocRealToComplex ())\r
+      return false;\r
+  }\r
+  \r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vLHSImag = getImaginaryArray();\r
+  ImageFileArray vRealResult = result.getArray();\r
+  ImageFileArray vImagResult = result.getImaginaryArray();\r
+  \r
+  unsigned int ix, iy;\r
+\r
+  // alloc output matrix\r
+  CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    complexOut[ix] = new CTSimComplex [m_ny];\r
+  \r
+  // fourier each x column\r
+#if 1\r
+  CTSimComplex* pY = new CTSimComplex [m_ny];\r
   for (ix = 0; ix < m_nx; ix++) {\r
-    for (iy = 0; iy < m_ny; iy++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      pY[iy].real (vLHS[ix][iy]);\r
+      if (isComplex())\r
+        pY[iy].imag (vLHSImag[ix][iy]);\r
+      else\r
+        pY[iy].imag (0);\r
+    } \r
+    ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
+  }\r
+#else\r
+  double* pY = new double [m_ny];\r
+  for (ix = 0; ix < m_nx; ix++) {\r
+    for (iy = 0; iy < m_ny; iy++) {\r
       pY[iy] = vLHS[ix][iy];\r
+    } \r
     ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
   }\r
-  delete pY;\r
-\r
-  std::complex<double>* pX = new std::complex<double> [m_nx];\r
-  std::complex<double>* complexOutCol = new std::complex<double> [m_nx];\r
+#endif\r
+  delete [] pY;\r
+  \r
+  // fourier each y row\r
+  CTSimComplex* pX = new CTSimComplex [m_nx];\r
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
   for (iy = 0; iy < m_ny; iy++) {\r
     for (ix = 0; ix < m_nx; ix++)\r
       pX[ix] = complexOut[ix][iy];\r
-    ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD);\r
+    ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);\r
     for (ix = 0; ix < m_nx; ix++)\r
-      complexOut[ix][iy] = complexOutCol[ix];\r
+      complexOut[ix][iy] = complexOutRow[ix];\r
   }\r
   delete [] pX;\r
-\r
+  \r
+  // shuffle each column\r
   for (ix = 0; ix < m_nx; ix++)\r
     ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);\r
+  \r
+  // shuffle each row\r
+  for (iy = 0; iy < m_ny; iy++) {\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOutRow[ix] = complexOut[ix][iy];\r
+    ProcessSignal::shuffleFourierToNaturalOrder (complexOutRow, m_nx);;\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOut[ix][iy] = complexOutRow[ix];\r
+    \r
+  }\r
+  delete [] complexOutRow;\r
+  \r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      vRealResult[ix][iy] = complexOut[ix][iy].real();\r
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
+    }\r
+    \r
+    for (ix = 0; ix < m_nx; ix++)\r
+      delete [] complexOut[ix];\r
+    delete [] complexOut;\r
+    \r
+    return true;\r
+}\r
 \r
+bool\r
+ImageFile::inverseFourier (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+  \r
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
+    if (! result.reallocRealToComplex ())\r
+      return false;\r
+  }\r
+  \r
+  ImageFileArrayConst vLHSReal = getArray();\r
+  ImageFileArrayConst vLHSImag = getImaginaryArray();\r
   \r
+  unsigned int ix, iy;\r
+\r
+  // alloc 2d complex output matrix\r
+  CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    complexOut[ix] = new CTSimComplex [m_ny];\r
+\r
+  // put input image into complexOut\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      complexOut[ix][iy].real (vLHSReal[ix][iy]);\r
+      if (isComplex())\r
+        complexOut[ix][iy].imag (vLHSImag[ix][iy]);\r
+      else\r
+        complexOut[ix][iy].imag (0);\r
+    }\r
+\r
+  // shuffle each x column\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    ProcessSignal::shuffleNaturalToFourierOrder (complexOut[ix], m_ny);\r
+\r
+  // shuffle each y row\r
+  CTSimComplex* pComplexRow = new CTSimComplex [m_nx];\r
   for (iy = 0; iy < m_ny; iy++) {\r
     for (ix = 0; ix < m_nx; ix++)\r
-      complexOutCol[ix] = complexOut[ix][iy];\r
-    ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);;\r
+      pComplexRow[ix] = complexOut[ix][iy];\r
+    ProcessSignal::shuffleNaturalToFourierOrder (pComplexRow, m_nx);\r
     for (ix = 0; ix < m_nx; ix++)\r
-      complexOut[ix][iy] = complexOutCol[ix];\r
+      complexOut[ix][iy] = pComplexRow[ix];\r
+  }\r
+  delete [] pComplexRow;\r
 \r
+  // ifourier each x column\r
+  CTSimComplex* pCol = new CTSimComplex [m_ny];\r
+  for (ix = 0; ix < m_nx; ix++) {\r
+    for (iy = 0; iy < m_ny; iy++)\r
+      pCol[iy] = complexOut[ix][iy];\r
+    ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);\r
   }\r
-  delete [] complexOutCol;\r
+  delete [] pCol;\r
 \r
+  // ifourier each y row\r
+  CTSimComplex* complexInRow = new CTSimComplex [m_nx];\r
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
+  for (iy = 0; iy < m_ny; iy++) {\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexInRow[ix] = complexOut[ix][iy];\r
+    ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOut[ix][iy] = complexOutRow[ix];\r
+  }\r
+  delete [] complexInRow;\r
+  delete [] complexOutRow;\r
+\r
+  ImageFileArray vRealResult = result.getArray();\r
+  ImageFileArray vImagResult = result.getImaginaryArray();\r
   for (ix = 0; ix < m_nx; ix++)\r
     for (iy = 0; iy < m_ny; iy++) {\r
       vRealResult[ix][iy] = complexOut[ix][iy].real();\r
       vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
     }\r
-\r
+    \r
+  // delete complexOut matrix\r
   for (ix = 0; ix < m_nx; ix++)\r
     delete [] complexOut[ix];\r
-\r
   delete [] complexOut;\r
-\r
+    \r
   return true;\r
 }\r
 \r
+\r
 bool\r
 ImageFile::magnitude (ImageFile& result) const\r
 {\r
@@ -602,29 +910,23 @@ ImageFile::magnitude (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArray vReal = getArray();\r
-  ImageFileArray vImag = NULL;\r
-  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
-    vImag = getImaginaryArray();\r
-\r
+  ImageFileArray vImag = getImaginaryArray();\r
+  \r
   ImageFileArray vRealResult = result.getArray();\r
-  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
-    ImageFileArray vImagResult = result.getImaginaryArray();\r
-    for (unsigned int ix = 0; ix < m_nx; ix++)\r
-      for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        vImagResult[ix][iy] = 0;\r
-  }\r
-\r
   for (unsigned int ix = 0; ix < m_nx; ix++)\r
     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
-      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX\r
+      if (isComplex()\r
         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
       else\r
         vRealResult[ix][iy] = vReal[ix][iy];\r
     }\r
 \r
-  return true;\r
+  if (result.isComplex())\r
+    result.reallocComplexToReal();\r
+  \r
+    return true;\r
 }\r
 \r
 bool\r
@@ -634,28 +936,22 @@ ImageFile::phase (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArray vReal = getArray();\r
-  ImageFileArray vImag = NULL;\r
-  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
-    vImag = getImaginaryArray();\r
-\r
+  ImageFileArray vImag = getImaginaryArray();\r
+  \r
   ImageFileArray vRealResult = result.getArray();\r
-  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
-    ImageFileArray vImagResult = result.getImaginaryArray();\r
-    for (unsigned int ix = 0; ix < m_nx; ix++)\r
-      for (unsigned int iy = 0; iy < m_ny; iy++)\r
-        vImagResult[ix][iy] = 0;\r
-  }\r
-\r
   for (unsigned int ix = 0; ix < m_nx; ix++)\r
     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
-      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) \r
-        vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]);\r
+      if (isComplex())\r
+        vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);\r
       else\r
-        vRealResult[ix][iy] = vReal[ix][iy];\r
+        vRealResult[ix][iy] = 0;\r
     }\r
-\r
+    \r
+  if (result.isComplex())\r
+    result.reallocComplexToReal();\r
+  \r
   return true;\r
 }\r
 \r
@@ -666,19 +962,19 @@ ImageFile::square (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
     return false;\r
   }\r
-\r
+  \r
   ImageFileArrayConst vLHS = getArray();\r
   ImageFileArray vResult = result.getArray();\r
-\r
+  \r
   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
     ImageFileColumnConst in = vLHS[ix];\r
     ImageFileColumn out = vResult[ix];\r
     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
-        *out++ = *in * *in;\r
-        in++;\r
+      *out++ = *in * *in;\r
+      in++;\r
     }\r
   }\r
-\r
+  \r
   return true;\r
 }\r
 \r
@@ -686,74 +982,74 @@ ImageFile::square (ImageFile& result) const
 void 
 ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 {
-       FILE *fp;
-       int nx = m_nx;
-       int ny = m_ny;
-       ImageFileArray v = getArray();
-       
-       unsigned char* rowp = new unsigned char [nx * nxcell];
-       
-       if ((fp = fopen (outfile, "wb")) == NULL)
-               return;
-       
-       fprintf(fp, "P5\n");
-       fprintf(fp, "%d %d\n", nx, ny);
-       fprintf(fp, "255\n");
-       
-       for (int irow = ny - 1; irow >= 0; irow--) {
-               for (int icol = 0; icol < nx; icol++) {
-                       int pos = icol * nxcell;
-                       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
-                       dens = clamp (dens, 0., 1.);
-                       for (int p = pos; p < pos + nxcell; p++) {
-                               rowp[p] = static_cast<unsigned int> (dens * 255.);
-                       }
-               }
-               for (int ir = 0; ir < nycell; ir++) {
-                       for (int ic = 0; ic < nx * nxcell; ic++) 
-                               fputc( rowp[ic], fp );
-               }
-       }
-       \r
-       delete rowp;
-       fclose(fp);
+  FILE *fp;
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  
+  unsigned char* rowp = new unsigned char [nx * nxcell];
+  
+  if ((fp = fopen (outfile, "wb")) == NULL)
+    return;
+  
+  fprintf(fp, "P5\n");
+  fprintf(fp, "%d %d\n", nx, ny);
+  fprintf(fp, "255\n");
+  
+  for (int irow = ny - 1; irow >= 0; irow--) {
+    for (int icol = 0; icol < nx; icol++) {
+      int pos = icol * nxcell;
+      double dens = (v[icol][irow] - densmin) / (densmax - densmin);
+      dens = clamp (dens, 0., 1.);
+      for (int p = pos; p < pos + nxcell; p++) {
+        rowp[p] = static_cast<unsigned int> (dens * 255.);
+      }
+    }
+    for (int ir = 0; ir < nycell; ir++) {
+      for (int ic = 0; ic < nx * nxcell; ic++) 
+        fputc( rowp[ic], fp );
+    }
+  }
+  \r
+  delete rowp;
+  fclose(fp);
 }
 
 void 
 ImageFile::writeImagePGMASCII (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 {
-       FILE *fp;
-       int nx = m_nx;
-       int ny = m_ny;
-       ImageFileArray v = getArray();
-       
-       unsigned char* rowp = new unsigned char [nx * nxcell];
-       
-       if ((fp = fopen (outfile, "wb")) == NULL)
-               return;
-       
-       fprintf(fp, "P2\n");
-       fprintf(fp, "%d %d\n", nx, ny);
-       fprintf(fp, "255\n");
-       
-       for (int irow = ny - 1; irow >= 0; irow--) {
-               for (int icol = 0; icol < nx; icol++) {
-                       int pos = icol * nxcell;
-                       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
-                       dens = clamp (dens, 0., 1.);
-                       for (int p = pos; p < pos + nxcell; p++) {
-                               rowp[p] = static_cast<unsigned int> (dens * 255.);
-                       }
-               }
-               for (int ir = 0; ir < nycell; ir++) {
-                       for (int ic = 0; ic < nx * nxcell; ic++) 
-                               fprintf(fp, "%d ", rowp[ic]);
-                       fprintf(fp, "\n");
-               }
-       }
-       \r
-       delete rowp;
-       fclose(fp);
+  FILE *fp;
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  
+  unsigned char* rowp = new unsigned char [nx * nxcell];
+  
+  if ((fp = fopen (outfile, "wb")) == NULL)
+    return;
+  
+  fprintf(fp, "P2\n");
+  fprintf(fp, "%d %d\n", nx, ny);
+  fprintf(fp, "255\n");
+  
+  for (int irow = ny - 1; irow >= 0; irow--) {
+    for (int icol = 0; icol < nx; icol++) {
+      int pos = icol * nxcell;
+      double dens = (v[icol][irow] - densmin) / (densmax - densmin);
+      dens = clamp (dens, 0., 1.);
+      for (int p = pos; p < pos + nxcell; p++) {
+        rowp[p] = static_cast<unsigned int> (dens * 255.);
+      }
+    }
+    for (int ir = 0; ir < nycell; ir++) {
+      for (int ic = 0; ic < nx * nxcell; ic++) 
+        fprintf(fp, "%d ", rowp[ic]);
+      fprintf(fp, "\n");
+    }
+  }
+  \r
+  delete rowp;
+  fclose(fp);
 }
 
 
@@ -761,69 +1057,69 @@ ImageFile::writeImagePGMASCII (const char *outfile, int nxcell, int nycell, doub
 void 
 ImageFile::writeImagePNG (const char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
 {
-       FILE *fp;
-       png_structp png_ptr;
-       png_infop info_ptr;
-       double max_out_level = (1 << bitdepth) - 1;
-       int nx = m_nx;
-       int ny = m_ny;
-       ImageFileArray v = getArray();
-       
-       unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
-       
-       if ((fp = fopen (outfile, "wb")) == NULL)
-               return;
-       
-       png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
-       if (! png_ptr)
-               return;
-       
-       info_ptr = png_create_info_struct(png_ptr);
-       if (! info_ptr) {
-               png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
-               fclose(fp);
-               return;
-       }
-       
-       if (setjmp(png_ptr->jmpbuf)) {
-               png_destroy_write_struct(&png_ptr, &info_ptr);
-               fclose(fp);
-               return;
-       }
-       
-       png_init_io(png_ptr, fp);
-       
-       png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT);
-       
-       png_write_info(png_ptr, info_ptr);
-       for (int irow = ny - 1; irow >= 0; irow--) {
-               png_bytep row_pointer = rowp;
-               
-               for (int icol = 0; icol < nx; icol++) {
-                       int pos = icol * nxcell;
-                       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
-                       dens = clamp (dens, 0., 1.);
-                       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
-                       
-                       for (int p = pos; p < pos + nxcell; p++) {
-                               if (bitdepth == 8)
-                                       rowp[p] = outval;
-                               else {
-                                       int rowpos = p * 2;
-                                       rowp[rowpos] = (outval >> 8) & 0xFF;
-                                       rowp[rowpos+1] = (outval & 0xFF);
-                               }
-                       }
-               }
-               for (int ir = 0; ir < nycell; ir++)
-                       png_write_rows (png_ptr, &row_pointer, 1);
-       }
-       
-       png_write_end(png_ptr, info_ptr);
-       png_destroy_write_struct(&png_ptr, &info_ptr);
-       delete rowp;\r
-       
-       fclose(fp);
+  FILE *fp;
+  png_structp png_ptr;
+  png_infop info_ptr;
+  double max_out_level = (1 << bitdepth) - 1;
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  
+  unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
+  
+  if ((fp = fopen (outfile, "wb")) == NULL)
+    return;
+  
+  png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+  if (! png_ptr)
+    return;
+  
+  info_ptr = png_create_info_struct(png_ptr);
+  if (! info_ptr) {
+    png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
+    fclose(fp);
+    return;
+  }
+  
+  if (setjmp(png_ptr->jmpbuf)) {
+    png_destroy_write_struct(&png_ptr, &info_ptr);
+    fclose(fp);
+    return;
+  }
+  
+  png_init_io(png_ptr, fp);
+  
+  png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT);
+  
+  png_write_info(png_ptr, info_ptr);
+  for (int irow = ny - 1; irow >= 0; irow--) {
+    png_bytep row_pointer = rowp;
+    
+    for (int icol = 0; icol < nx; icol++) {
+      int pos = icol * nxcell;
+      double dens = (v[icol][irow] - densmin) / (densmax - densmin);
+      dens = clamp (dens, 0., 1.);
+      unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
+      
+      for (int p = pos; p < pos + nxcell; p++) {
+        if (bitdepth == 8)
+          rowp[p] = outval;
+        else {
+          int rowpos = p * 2;
+          rowp[rowpos] = (outval >> 8) & 0xFF;
+          rowp[rowpos+1] = (outval & 0xFF);
+        }
+      }
+    }
+    for (int ir = 0; ir < nycell; ir++)
+      png_write_rows (png_ptr, &row_pointer, 1);
+  }
+  
+  png_write_end(png_ptr, info_ptr);
+  png_destroy_write_struct(&png_ptr, &info_ptr);
+  delete rowp;\r
+  
+  fclose(fp);
 }
 #endif
 
@@ -834,44 +1130,44 @@ static const int N_GRAYSCALE=256;
 void
 ImageFile::writeImageGIF (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 {
-       gdImagePtr gif;
-       FILE *out;
-       int gs_indices[N_GRAYSCALE];
-       int nx = m_nx;
-       int ny = m_ny;
-       ImageFileArray v = getArray();
-       
-       unsigned char rowp [nx * nxcell];
-       if (rowp == NULL)
-               return;
-       
-       gif = gdImageCreate(nx * nxcell, ny * nycell);
-       for (int i = 0; i < N_GRAYSCALE; i++)
-               gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
-       
-       int lastrow = ny * nycell - 1;
-       for (int irow = 0; irow < ny; irow++) {
-               int rpos = irow * nycell;
-               for (int ir = rpos; ir < rpos + nycell; ir++) {
-                       for (int icol = 0; icol < nx; icol++) {
-                               int cpos = icol * nxcell;
-                               double dens = (v[icol][irow] - densmin) / (densmax - densmin);
-                               dens = clamp(dens, 0., 1.);
-                               for (int ic = cpos; ic < cpos + nxcell; ic++) {
-                                       rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
-                                       gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
-                               }
-                       }
-               }
-       }
-       
-       if ((out = fopen(outfile,"w")) == NULL) {
-               sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
-               return (1);
-       }
-       gdImageGif(gif,out);
-       fclose(out);
-       gdImageDestroy(gif);
+  gdImagePtr gif;
+  FILE *out;
+  int gs_indices[N_GRAYSCALE];
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  
+  unsigned char rowp [nx * nxcell];
+  if (rowp == NULL)
+    return;
+  
+  gif = gdImageCreate(nx * nxcell, ny * nycell);
+  for (int i = 0; i < N_GRAYSCALE; i++)
+    gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
+  
+  int lastrow = ny * nycell - 1;
+  for (int irow = 0; irow < ny; irow++) {
+    int rpos = irow * nycell;
+    for (int ir = rpos; ir < rpos + nycell; ir++) {
+      for (int icol = 0; icol < nx; icol++) {
+        int cpos = icol * nxcell;
+        double dens = (v[icol][irow] - densmin) / (densmax - densmin);
+        dens = clamp(dens, 0., 1.);
+        for (int ic = cpos; ic < cpos + nxcell; ic++) {
+          rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
+          gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
+        }
+      }
+    }
+  }
+  
+  if ((out = fopen(outfile,"w")) == NULL) {
+    sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
+    return (1);
+  }
+  gdImageGif(gif,out);
+  fclose(out);
+  gdImageDestroy(gif);
 }
 #endif
 
index a7933b2..8c45da4 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.cpp,v 1.11 2000/12/29 15:45:06 kevin Exp $
+**  $Id: procsignal.cpp,v 1.12 2001/01/01 10:14:34 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
@@ -199,7 +199,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         pEZPlot->plot (pSGP);
       }
 #endif
-      ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+      ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
       delete adFrequencyFilter;\r
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
@@ -370,10 +370,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       
       m_adFilter = new double [m_nFilterPoints];
       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];\r
-      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
+      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
       delete adSpatialFilter;\r
       for (i = 0; i < m_nFilterPoints; i++)
-        m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc;
+        m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
       delete acInverseFilter;\r
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
@@ -553,12 +553,12 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);\r
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);\r
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);\r
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);\r
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];\r
@@ -570,12 +570,12 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, -1);\r
+    finiteFourierTransform (inputSignal, fftSignal, FORWARD);\r
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, 1);\r
+    finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);\r
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];\r
@@ -733,8 +733,8 @@ ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::
     std::complex<double> sum (0,0);
     for (int j = 0; j < n; j++) {
       double angle = i * j * angleIncrement;
-      std::complex<double> exponentTerm (cos(angle), sin(angle));
-      sum += input[j] * exponentTerm;
+      std::complex<double> exponentTerm (cos(angle), sin(angle));\r
+      sum += input[j] * exponentTerm;\r
     }
     if (direction < 0) {
       sum /= n;
@@ -878,10 +878,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
   } else {     // Even\r
     int iHalfN = n / 2;\r
     pdTemp[0] = pdVector[iHalfN];\r
+#if USE_BROKEN_SHUFFLE\r
     for (i = 0; i < iHalfN; i++)\r
       pdTemp[i + 1] = pdVector[i + iHalfN];\r
     for (i = 0; i < iHalfN - 1; i++)\r
       pdTemp[i + iHalfN + 1] = pdVector[i];\r
+#else\r
+    for (i = 0; i < iHalfN - 1; i++)\r
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + iHalfN] = pdVector[i];\r
+#endif\r
   }\r
   \r
   for (i = 0; i < n; i++)\r
@@ -905,10 +912,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   } else {     // Even\r
     int iHalfN = n / 2;\r
     pdTemp[0] = pdVector[iHalfN];\r
+#if USE_BROKEN_SHUFFLE\r
     for (i = 0; i < iHalfN; i++)\r
       pdTemp[i + 1] = pdVector[i + iHalfN];\r
     for (i = 0; i < iHalfN - 1; i++)\r
       pdTemp[i + iHalfN + 1] = pdVector[i];\r
+#else\r
+    for (i = 0; i < iHalfN - 1; i++)\r
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + iHalfN] = pdVector[i];\r
+#endif\r
   }\r
   \r
   for (i = 0; i < n; i++)\r
@@ -916,7 +930,43 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;\r
 }\r
 \r
-
+\r
+void\r
+ProcessSignal::shuffleNaturalToFourierOrder (float* pdVector, const int n)\r
+{\r
+  float* pdTemp = new float [n];\r
+  int i;\r
+  if (n % 2) { // Odd\r
+    int iHalfN = (n - 1) / 2;\r
+    \r
+    pdTemp[0] = pdVector[iHalfN];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + iHalfN + 1] = pdVector[i];\r
+  } else {     // Even\r
+    int iHalfN = n / 2;\r
+    pdTemp[0] = pdVector[iHalfN];\r
+#if USE_BROKEN_SHUFFLE\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + 1] = pdVector[i + iHalfN];\r
+    for (i = 0; i < iHalfN - 1; i++)\r
+      pdTemp[i + iHalfN + 1] = pdVector[i];\r
+#else\r
+    for (i = 0; i < iHalfN - 1; i++)\r
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pdTemp[i + iHalfN] = pdVector[i];\r
+#endif\r
+  }\r
+  \r
+  for (i = 0; i < n; i++)\r
+    pdVector[i] = pdTemp[i];\r
+  delete pdTemp;\r
+}\r
+\r
+\r
+\r
 void\r
 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)\r
 {\r
@@ -944,6 +994,7 @@ ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
   delete pdTemp;\r
 }\r
 \r
+\r
 void\r
 ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)\r
 {\r
@@ -971,3 +1022,33 @@ ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;\r
 }\r
 \r
+\r
+\r
+\r
+void\r
+ProcessSignal::shuffleFourierToNaturalOrder (float* pVector, const int n)\r
+{\r
+  float* pTemp = new float [n];\r
+  int i;\r
+  if (n % 2) { // Odd\r
+    int iHalfN = (n - 1) / 2;\r
+    \r
+    pTemp[iHalfN] = pVector[0];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pTemp[i + 1 + iHalfN] = pVector[i + 1];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pTemp[i] = pVector[i + iHalfN + 1];\r
+  } else {     // Even\r
+    int iHalfN = n / 2;\r
+    pTemp[iHalfN] = pVector[0];\r
+    for (i = 0; i < iHalfN; i++)\r
+      pTemp[i] = pVector[i + iHalfN];\r
+    for (i = 0; i < iHalfN - 1; i++)\r
+      pTemp[i + iHalfN + 1] = pVector[i+1];\r
+  }\r
+  \r
+  for (i = 0; i < n; i++)\r
+    pVector[i] = pTemp[i];\r
+  delete [] pTemp;\r
+}\r
+\r
index 2eaa551..3604b2f 100644 (file)
@@ -2,7 +2,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: mathfuncs.cpp,v 1.4 2000/12/21 03:40:58 kevin Exp $
+**  $Id: mathfuncs.cpp,v 1.5 2001/01/01 10:14:34 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
 
 
 /* NAME
- *    integrateSimpson         Integrate array of data by Simpson's rule
- *
- * SYNOPSIS
- *    double integrateSimpson (xmin, xmax, y, np)
- *    double xmin, xmax                Extent of integration
- *    double y[]               Function values to be integrated
- *    int np                   number of data points
- *                             (must be an odd number and at least 3)
- *
- * RETURNS
- *    integrand of function
- */
+*    integrateSimpson         Integrate array of data by Simpson's rule
+*
+* SYNOPSIS
+*    double integrateSimpson (xmin, xmax, y, np)
+*    double xmin, xmax         Extent of integration
+*    double y[]                Function values to be integrated
+*    int np                    number of data points
+                             (must be an odd number and at least 3)
+*
+* RETURNS
+*    integrand of function
+*/
 
 double 
 integrateSimpson (const double xmin, const double xmax, const double *y, const int np) 
@@ -42,16 +42,16 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i
     return (0.);
   else if (np == 2)
     return ((xmax - xmin) * (y[0] + y[1]) / 2);
-
+  
   double area = 0;
   int nDiv = (np - 1) / 2;  // number of divisions
   double width = (xmax - xmin) / (double) (np - 1);    // width of cells
-
+  
   for (int i = 1; i <= nDiv; i++) {
     int xr = 2 * i;
     int xl = xr - 2;       // 2 * (i - 1) == 2 * i - 2 == xr - 2
     int xm = xr - 1;       // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
-
+    
     area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
   }
   
@@ -63,13 +63,13 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i
 
 
 /* NAME
- *    normalizeAngle       Normalize angle to 0 to 2 * PI range
- *
- * SYNOPSIS
- *    t = normalizeAngle (theta)
- *    double t        Normalized angle
- *    double theta     Input angle
- */
+*    normalizeAngle       Normalize angle to 0 to 2 * PI range
+*
+* SYNOPSIS
+*    t = normalizeAngle (theta)
+*    double t         Normalized angle
+*    double theta     Input angle
+*/
 
 double 
 normalizeAngle (double theta)
@@ -86,52 +86,55 @@ normalizeAngle (double theta)
 void \r
 vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)\r
 {\r
-       if (nPoints <= 0)\r
-               return;\r
-\r
-    mean = 0;\r
-    min = vec[0];\r
-    max = vec[0];\r
-       int i;\r
-    for (i = 0; i < nPoints; i++) {\r
-               double v = vec[i];\r
-               if (v > max)\r
-                       max = v;\r
-               if (v < min)\r
-                       min = v;\r
-               mean += v;\r
-    }\r
-    mean /= nPoints;\r
-       \r
-    static const int nbin = 1024;\r
-    int hist[ nbin ] = {0};\r
-    double spread = max - min;\r
-    mode = 0;\r
-    stddev = 0;\r
-    for (i = 0; i < nPoints; i++) {\r
-               double v = vec[i];\r
-               int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);\r
-               hist[b]++;\r
-               double diff = (v - mean);\r
-               stddev += diff * diff;\r
+  if (nPoints <= 0)\r
+    return;\r
+  \r
+  mean = 0;\r
+  min = vec[0];\r
+  max = vec[0];\r
+  int i;\r
+  int iMinPos = 0;\r
+  for (i = 0; i < nPoints; i++) {\r
+    double v = vec[i];\r
+    if (v > max)\r
+      max = v;\r
+    if (v < min) {\r
+      min = v;\r
+      iMinPos = i;\r
     }\r
-    stddev = sqrt (stddev / nPoints);\r
-       \r
-    int max_binindex = 0;\r
-    int max_bin = -1;\r
-    for (int ibin = 0; ibin < nbin; ibin++) {\r
-               if (hist[ibin] > max_bin) {\r
-                       max_bin = hist[ibin];\r
-                       max_binindex = ibin;\r
-               }\r
+    mean += v;\r
+  }\r
+  mean /= nPoints;\r
+  \r
+  static const int nbin = 1024;\r
+  int hist[ nbin ] = {0};\r
+  double spread = max - min;\r
+  mode = 0;\r
+  stddev = 0;\r
+  for (i = 0; i < nPoints; i++) {\r
+    double v = vec[i];\r
+    int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);\r
+    hist[b]++;\r
+    double diff = (v - mean);\r
+    stddev += diff * diff;\r
+  }\r
+  stddev = sqrt (stddev / nPoints);\r
+  \r
+  int max_binindex = 0;\r
+  int max_bin = -1;\r
+  for (int ibin = 0; ibin < nbin; ibin++) {\r
+    if (hist[ibin] > max_bin) {\r
+      max_bin = hist[ibin];\r
+      max_binindex = ibin;\r
     }\r
-       \r
-    mode = (max_binindex * spread / (nbin - 1)) + min;\r
-       \r
-       std::sort(vec.begin(), vec.end());\r
-       \r
-       if (nPoints % 2)  // Odd\r
-               median = vec[((nPoints - 1) / 2)];\r
-       else        // Even\r
-               median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;\r
+  }\r
+  \r
+  mode = (max_binindex * spread / (nbin - 1)) + min;\r
+  \r
+  std::sort(vec.begin(), vec.end());\r
+  \r
+  if (nPoints % 2)  // Odd\r
+    median = vec[((nPoints - 1) / 2)];\r
+  else        // Even\r
+    median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;\r
 }\r
index daff73f..fb8f744 100644 (file)
 <pre>\r
 <h1>Build Log</h1>\r
 <h3>\r
---------------------Configuration: FFTW2st - Win32 Release--------------------\r
+--------------------Configuration: libctsim - Win32 Debug--------------------\r
 </h3>\r
 <h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-FFTW2st.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: FFTW2st - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-FFTW2st.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: RFFTW2st - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-RFFTW2st.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: RFFTW2st - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-RFFTW2st.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: ctsim - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-ctsim.exe - 0 error(s), 0 warning(s)\r
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29A.tmp" with contents\r
+[\r
+/nologo /G6 /MTd /W3 /Gm /GR /GX /ZI /Od /I "..\..\..\zlib" /I "..\..\INCLUDE" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2\include" /D "_DEBUG" /D "HAVE_WXWIN" /D "HAVE_STRING_H" /D "HAVE_GETOPT_H" /D "WIN32" /D "_MBCS" /D "_LIB" /D "MSVC" /D "HAVE_FFTW" /D "HAVE_PNG" /D "HAVE_SGP" /D "HAVE_WXWINDOWS" /D "__WXMSW__" /FR"Debug/" /Fp"Debug/libctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c \r
+"C:\ctsim-2.0.6\libctsim\imagefile.cpp"\r
+]\r
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29A.tmp" \r
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29B.tmp" with contents\r
+[\r
+/nologo /out:"Debug\libctsim.lib" \r
+".\Debug\array2dfile.obj"\r
+".\Debug\backprojectors.obj"\r
+".\Debug\clip.obj"\r
+".\Debug\consoleio.obj"\r
+".\Debug\ezplot.obj"\r
+".\Debug\ezset.obj"\r
+".\Debug\ezsupport.obj"\r
+".\Debug\filter.obj"\r
+".\Debug\fnetorderstream.obj"\r
+".\Debug\getopt.obj"\r
+".\Debug\getopt1.obj"\r
+".\Debug\hashtable.obj"\r
+".\Debug\imagefile.obj"\r
+".\Debug\mathfuncs.obj"\r
+".\Debug\phantom.obj"\r
+".\Debug\plotfile.obj"\r
+".\Debug\pol.obj"\r
+".\Debug\procsignal.obj"\r
+".\Debug\projections.obj"\r
+".\Debug\reconstruct.obj"\r
+".\Debug\scanner.obj"\r
+".\Debug\sgp.obj"\r
+".\Debug\strfuncs.obj"\r
+".\Debug\syserror.obj"\r
+".\Debug\trace.obj"\r
+".\Debug\transformmatrix.obj"\r
+".\Debug\xform.obj"\r
+]\r
+Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29B.tmp"\r
+<h3>Output Window</h3>\r
+Compiling...\r
+imagefile.cpp\r
+Creating library...\r
 <h3>\r
 --------------------Configuration: ctsim - Win32 Debug--------------------\r
 </h3>\r
 <h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-ctsim.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: if1 - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-if1.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: if1 - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-if1.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: if2 - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-if2.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: if2 - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-if2.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: ifexport - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-ifexport.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: ifexport - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPA0.tmp" with contents\r
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29C.tmp" with contents\r
 [\r
-rpcrt4.lib comctl32.lib wsock32.lib ../libctsim/Debug/libctsim.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ..\..\..\lpng108\msvc\win32\libpng\lib_dbg\libpng.lib ..\..\..\lpng108\msvc\win32\zlib\lib_dbg\zlib.lib libcpd.lib libcd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib /nologo /subsystem:console /incremental:yes /pdb:"Debug/ifexport.pdb" /debug /machine:I386 /nodefaultlib:"libcd.lib" /nodefaultlib:"libcid.lib" /nodefaultlib:"msvcrt.lib" /nodefaultlib /out:"Debug/ifexport.exe" /pdbtype:sept \r
-".\Debug\ifexport.obj"\r
+comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ../libctsim/Debug/libctsim.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ..\..\..\lpng108\msvc\win32\libpng\lib_dbg\libpng.lib ..\..\..\lpng108\msvc\win32\zlib\lib_dbg\zlib.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib ../../../wx2/lib/wxd.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /nodefaultlib:"libcd.lib" /nodefaultlib:"libcid.lib" /nodefaultlib:"msvcrtd.lib" /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"..\..\..\lpng108\msvc\win32\libpng\lib" /libpath:"..\..\..\lpng108\msvc\win32\zlib\lib" \r
+".\Debug\ctsim.obj"\r
+".\Debug\dialogs.obj"\r
+".\Debug\dlgprojections.obj"\r
+".\Debug\dlgreconstruct.obj"\r
+".\Debug\docs.obj"\r
+".\Debug\views.obj"\r
+".\Debug\wx.res"\r
 "\ctsim-2.0.6\msvc\libctsim\Debug\libctsim.lib"\r
 "\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib"\r
 "\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib"\r
+"\wx2\lib\wxd.lib"\r
 ]\r
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPA0.tmp"\r
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29C.tmp"\r
 <h3>Output Window</h3>\r
 Linking...\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxPen::~wxPen(void)" (??1wxPen@@UAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxBrush *  wxWHITE_BRUSH" (?wxWHITE_BRUSH@@3PAVwxBrush@@A)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "char const * const  wxEmptyString" (?wxEmptyString@@3PBDB)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetWidth(int)" (?SetWidth@wxPen@@QAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxPen::wxPen(void)" (??0wxPen@@QAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxFont  wxNullFont" (?wxNullFont@@3VwxFont@@A)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxBrush::~wxBrush(void)" (??1wxBrush@@UAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetColour(unsigned char,unsigned char,unsigned char)" (?SetColour@wxPen@@QAEXEEE@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxBrush  wxNullBrush" (?wxNullBrush@@3VwxBrush@@A)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxBrush::SetColour(unsigned char,unsigned char,unsigned char)" (?SetColour@wxBrush@@UAEXEEE@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxBrush::wxBrush(void)" (??0wxBrush@@QAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxColour::~wxColour(void)" (??1wxColour@@UAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxColour::wxColour(unsigned char,unsigned char,unsigned char)" (??0wxColour@@QAE@EEE@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetColour(class wxColour const &)" (?SetColour@wxPen@@QAEXABVwxColour@@@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetStyle(int)" (?SetStyle@wxPen@@QAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "private: void __thiscall wxString::InitWith(char const *,unsigned int,unsigned int)" (?InitWith@wxString@@AAEXPBDII@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetEncoding(enum wxFontEncoding)" (?SetEncoding@wxFont@@UAEXW4wxFontEncoding@@@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetUnderlined(bool)" (?SetUnderlined@wxFont@@UAEX_N@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetFaceName(class wxString const &)" (?SetFaceName@wxFont@@UAEXABVwxString@@@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetWeight(int)" (?SetWeight@wxFont@@UAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetStyle(int)" (?SetStyle@wxFont@@UAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetFamily(int)" (?SetFamily@wxFont@@UAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetPointSize(int)" (?SetPointSize@wxFont@@UAEXH@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual enum wxFontEncoding  __thiscall wxFont::GetEncoding(void)const " (?GetEncoding@wxFont@@UBE?AW4wxFontEncoding@@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual class wxString  __thiscall wxFont::GetFaceName(void)const " (?GetFaceName@wxFont@@UBE?AVwxString@@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::GetUnderlined(void)const " (?GetUnderlined@wxFont@@UBE_NXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetWeight(void)const " (?GetWeight@wxFont@@UBEHXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetStyle(void)const " (?GetStyle@wxFont@@UBEHXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetFamily(void)const " (?GetFamily@wxFont@@UBEHXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetPointSize(void)const " (?GetPointSize@wxFont@@UBEHXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual unsigned long __thiscall wxFont::GetResourceHandle(void)" (?GetResourceHandle@wxFont@@UAEKXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::IsFree(void)const " (?IsFree@wxFont@@UBE_NXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::FreeResource(bool)" (?FreeResource@wxFont@@UAE_N_N@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::RealizeResource(void)" (?RealizeResource@wxFont@@UAE_NXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxObject::CopyObject(class wxObject &)const " (?CopyObject@wxObject@@UBEXAAV1@@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: bool __thiscall wxFont::Create(int,int,int,int,bool,class wxString const &,enum wxFontEncoding)" (?Create@wxFont@@QAE_NHHHH_NABVwxString@@W4wxFontEncoding@@@Z)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "protected: void __thiscall wxFont::Init(void)" (?Init@wxFont@@IAEXXZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: static class wxClassInfo  wxFont::sm_classwxFont" (?sm_classwxFont@wxFont@@2VwxClassInfo@@A)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxFont::~wxFont(void)" (??1wxFont@@UAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxObject::wxObject(void)" (??0wxObject@@QAE@XZ)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: static class wxClassInfo  wxGDIObject::sm_classwxGDIObject" (?sm_classwxGDIObject@wxGDIObject@@2VwxClassInfo@@A)\r
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxObject::~wxObject(void)" (??1wxObject@@UAE@XZ)\r
-Debug/ifexport.exe : fatal error LNK1120: 42 unresolved externals\r
-Error executing link.exe.\r
 \r
 \r
 \r
 <h3>Results</h3>\r
-ifexport.exe - 43 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: ifinfo - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-ifinfo.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: ifinfo - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-ifinfo.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: libpng - Win32 LIB--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-libpng.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: libpng - Win32 LIB Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-libpng.lib - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: phm2if - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-phm2if.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: phm2if - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-phm2if.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: phm2pj - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-phm2pj.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: phm2pj - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-phm2pj.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pj2if - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pj2if.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pj2if - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pj2if.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pjinfo - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pjinfo.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pjinfo - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pjinfo.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pjrec - Win32 Release--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pjrec.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: pjrec - Win32 Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-pjrec.exe - 0 error(s), 0 warning(s)\r
-<h3>\r
---------------------Configuration: wxvc - Win32 LIB Debug--------------------\r
-</h3>\r
-<h3>Command Lines</h3>\r
-\r
-\r
-\r
-<h3>Results</h3>\r
-wxd.lib - 0 error(s), 0 warning(s)\r
+ctsim.exe - 0 error(s), 0 warning(s)\r
 </pre>\r
 </body>\r
 </html>\r
index 8dc7e22..172cda9 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsim.cpp,v 1.22 2000/12/29 20:18:59 kevin Exp $
+**  $Id: ctsim.cpp,v 1.23 2001/01/01 10:14:34 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
@@ -61,7 +61,7 @@
 #endif
 #endif
 \r
-static const char* rcsindent = "$Id: ctsim.cpp,v 1.22 2000/12/29 20:18:59 kevin Exp $";
+static const char* rcsindent = "$Id: ctsim.cpp,v 1.23 2001/01/01 10:14:34 kevin Exp $";
 
 class CTSimApp* theApp = NULL;
 
@@ -185,7 +185,8 @@ IMPLEMENT_CLASS(MainFrame, wxDocParentFrame)
 BEGIN_EVENT_TABLE(MainFrame, wxDocParentFrame)
   EVT_MENU(MAINMENU_HELP_ABOUT, MainFrame::OnAbout)
   EVT_MENU(MAINMENU_HELP_CONTENTS, MainFrame::OnHelpContents)
-  EVT_MENU(MAINMENU_FILE_CREATE_PHANTOM, MainFrame::OnCreatePhantom)
+  EVT_MENU(MAINMENU_FILE_CREATE_PHANTOM, MainFrame::OnCreatePhantom)\r
+  EVT_MENU(MAINMENU_FILE_CREATE_FILTER, MainFrame::OnCreateFilter)\r
   EVT_MENU(MAINMENU_FILE_EXIT, MainFrame::OnExit)
   EVT_MENU(MAINMENU_WINDOW_BASE, MainFrame::OnWindowMenu0)
   EVT_MENU(MAINMENU_WINDOW_BASE+1, MainFrame::OnWindowMenu1)
@@ -223,7 +224,8 @@ MainFrame::MainFrame(wxDocManager *manager, wxFrame *frame, wxWindowID id, const
     //// Make a menubar
     wxMenu *file_menu = new wxMenu;
     
-    file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");
+    file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");\r
+    file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");\r
     file_menu->Append(wxID_OPEN, "&Open...");
     
     file_menu->AppendSeparator();
@@ -271,20 +273,49 @@ MainFrame::OnAbout(wxCommandEvent& WXUNUSED(event) )
     wxMessageBox(msg, "About CTSim", wxOK | wxICON_INFORMATION, this);
 }
 
-void 
-MainFrame::OnCreatePhantom(wxCommandEvent& WXUNUSED(event))
-{
-    DialogGetPhantom dialogPhantom (this, Phantom::PHM_HERMAN);
-    int dialogReturn = dialogPhantom.ShowModal();
-    if (dialogReturn == wxID_OK) {
-      wxString selection (dialogPhantom.getPhantom());
-      *theApp->getLog() << "Selected phantom " << selection.c_str() << "\n";
-      wxString filename = selection + ".phm";
-      theApp->getDocManager()->CreateDocument(filename, wxDOC_SILENT);
-    }
-    
-}
-
+void \r
+MainFrame::OnCreatePhantom(wxCommandEvent& WXUNUSED(event))\r
+{\r
+    DialogGetPhantom dialogPhantom (this, Phantom::PHM_HERMAN);\r
+    int dialogReturn = dialogPhantom.ShowModal();\r
+    if (dialogReturn == wxID_OK) {\r
+      wxString selection (dialogPhantom.getPhantom());\r
+      *theApp->getLog() << "Selected phantom " << selection.c_str() << "\n";\r
+      wxString filename = selection + ".phm";\r
+      theApp->getDocManager()->CreateDocument(filename, wxDOC_SILENT);\r
+    }\r
+    \r
+}\r
+\r
+void \r
+MainFrame::OnCreateFilter (wxCommandEvent& WXUNUSED(event))\r
+{\r
+  DialogGetFilterParameters dialogFilter (this, 256, 256, SignalFilter::FILTER_BANDLIMIT, 1., 10., SignalFilter::DOMAIN_SPATIAL);\r
+    int dialogReturn = dialogFilter.ShowModal();\r
+    if (dialogReturn == wxID_OK) {\r
+      wxString strFilter (dialogFilter.getFilterName());\r
+      wxString strDomain (dialogFilter.getDomainName());\r
+      unsigned int nx = dialogFilter.getXSize();\r
+      unsigned int ny = dialogFilter.getYSize();\r
+      double dBandwidth = dialogFilter.getBandwidth();\r
+      double dFilterParam= dialogFilter.getFilterParam();\r
+      *theApp->getLog() << "Selected filter " << strFilter.c_str() << ", domain " << strDomain.c_str() << ", filterParam " << dFilterParam << ", bandwidth " << dBandwidth << "\n";\r
+      wxString filename = "untitled.if";\r
+      ImageFileDocument* pFilterDoc = dynamic_cast<ImageFileDocument*>(theApp->getDocManager()->CreateDocument ("untitled.if", wxDOC_SILENT));\r
+      if (! pFilterDoc) {\r
+        sys_error (ERR_SEVERE, "Unable to create filter image");\r
+        return;\r
+      }\r
+      ImageFile& rIF = pFilterDoc->getImageFile();\r
+      rIF.setArraySize (nx, ny);\r
+      rIF.filterResponse (strDomain.c_str(), dBandwidth, strFilter.c_str(), dFilterParam);\r
+      if (theApp->getSetModifyNewDocs())\r
+        pFilterDoc->Modify (true);\r
+      pFilterDoc->UpdateAllViews();\r
+      pFilterDoc->GetFirstView()->OnUpdate (NULL, NULL);\r
+    }\r
+}\r
+\r
 void\r
 CTSimApp::getCompatibleImages (const ImageFileDocument* pIFDoc, std::vector<ImageFileDocument*>& vecIF)\r
 {\r
index b944560..2055227 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsim.h,v 1.15 2000/12/29 19:30:08 kevin Exp $
+**  $Id: ctsim.h,v 1.16 2001/01/01 10:14:34 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
@@ -59,7 +59,8 @@ public:
     
     void OnAbout (wxCommandEvent& event);
     void OnHelpContents (wxCommandEvent& event);
-    void OnCreatePhantom (wxCommandEvent& event);
+    void OnCreatePhantom (wxCommandEvent& event);\r
+    void OnCreateFilter (wxCommandEvent& event);
     void OnExit (wxCommandEvent& event);
 
     void OnUpdateUI (wxUpdateUIEvent& event);
@@ -134,7 +135,8 @@ extern class CTSimApp* theApp;
 enum {
     MAINMENU_HELP_ABOUT = 500,
     MAINMENU_HELP_CONTENTS,
-    MAINMENU_FILE_CREATE_PHANTOM,
+    MAINMENU_FILE_CREATE_PHANTOM,\r
+    MAINMENU_FILE_CREATE_FILTER,
     MAINMENU_FILE_EXIT,
     IFMENU_FILE_PROPERTIES,
     PJMENU_FILE_PROPERTIES,
@@ -154,6 +156,8 @@ enum {
   IFMENU_PROCESS_EXP,\r
   IFMENU_PROCESS_FOURIER,\r
   IFMENU_PROCESS_INVERSE_FOURIER,\r
+  IFMENU_PROCESS_FFT,\r
+  IFMENU_PROCESS_IFFT,\r
   IFMENU_PROCESS_MAGNITUDE,\r
   IFMENU_PROCESS_PHASE,\r
     PHMMENU_PROCESS_RASTERIZE,
index 4a44735..3e2d8fb 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.cpp,v 1.21 2000/12/27 03:16:02 kevin Exp $
+**  $Id: dialogs.cpp,v 1.22 2001/01/01 10:14:34 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
@@ -814,3 +814,132 @@ DialogGetReconstructionParameters::getFilterGenerationName (void)
 }
 
 \r
+/////////////////////////////////////////////////////////////////////\r
+// CLASS IDENTIFICATION\r
+//\r
+// DialogGetFilterParameters\r
+/////////////////////////////////////////////////////////////////////\r
+\r
+\r
+DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDefaultXSize, int iDefaultYSize, int iDefaultFilterID, double dDefaultFilterParam,  double dDefaultBandwidth, int iDefaultDomainID)\r
+: wxDialog (pParent, -1, "Set Filter Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)\r
+{\r
+  wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);\r
+  \r
+  pTopSizer->Add (new wxStaticText (this, -1, "Set Filter Parameters"), 0, wxALIGN_CENTER | wxTOP | wxLEFT | wxRIGHT, 5);\r
+  \r
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);\r
+  \r
+  std::ostringstream os;\r
+  os << iDefaultXSize;\r
+  m_pTextCtrlXSize = new wxTextCtrl (this, -1, os.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);\r
+  std::ostringstream osYSize;\r
+  osYSize << iDefaultYSize;\r
+  m_pTextCtrlYSize = new wxTextCtrl (this, -1, osYSize.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);\r
+  std::ostringstream osFilterParam;\r
+  osFilterParam << dDefaultFilterParam;\r
+  m_pTextCtrlFilterParam = new wxTextCtrl (this, -1, osFilterParam.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);\r
+  std::ostringstream osBandwidth;\r
+  osBandwidth << dDefaultBandwidth;\r
+  m_pTextCtrlBandwidth = new wxTextCtrl (this, -1, osBandwidth.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);\r
+  \r
+  wxFlexGridSizer* pGridSizer = new wxFlexGridSizer (2);\r
+  pGridSizer->Add (new wxStaticText (this, -1, "Filter"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5);\r
+  m_pListBoxFilter = new StringValueAndTitleListBox (this, SignalFilter::getFilterCount(), SignalFilter::getFilterTitleArray(), SignalFilter::getFilterNameArray());\r
+  m_pListBoxFilter->SetSelection (iDefaultFilterID);\r
+  pGridSizer->Add (m_pListBoxFilter, 0, wxALL | wxALIGN_LEFT | wxEXPAND);\r
+  \r
+  m_pListBoxDomain = new StringValueAndTitleListBox (this, SignalFilter::getDomainCount(), SignalFilter::getDomainTitleArray(), SignalFilter::getDomainNameArray());\r
+  m_pListBoxDomain->SetSelection (iDefaultDomainID);\r
+  pGridSizer->Add (new wxStaticText (this, -1, "Domain"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5);\r
+  pGridSizer->Add (m_pListBoxDomain, 0, wxALL | wxALIGN_LEFT | wxEXPAND);\r
+  \r
+  pGridSizer->Add (new wxStaticText (this, -1, "X Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (m_pTextCtrlXSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (new wxStaticText (this, -1, "Y Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (m_pTextCtrlYSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (new wxStaticText (this, -1, "Filter Parameter"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (m_pTextCtrlFilterParam, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (new wxStaticText (this, -1, "Bandwidth"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);\r
+  pGridSizer->Add (m_pTextCtrlBandwidth, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);\r
+  \r
+  pTopSizer->Add (pGridSizer, 1, wxALL, 3);\r
+  \r
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);\r
+  \r
+  wxBoxSizer* pButtonSizer = new wxBoxSizer (wxHORIZONTAL);\r
+  wxButton* pButtonOk = new wxButton (this, wxID_OK, "Okay");\r
+  wxButton* pButtonCancel = new wxButton (this, wxID_CANCEL, "Cancel");\r
+  pButtonSizer->Add (pButtonOk, 0, wxEXPAND | wxALL, 10);\r
+  pButtonSizer->Add (pButtonCancel, 0, wxEXPAND | wxALL, 10);\r
+  \r
+  pTopSizer->Add (pButtonSizer, 0, wxALIGN_CENTER);\r
+  \r
+  SetAutoLayout (true);\r
+  SetSizer (pTopSizer);\r
+  pTopSizer->Layout();\r
+  pTopSizer->Fit (this);\r
+  pTopSizer->SetSizeHints (this);\r
+}\r
+\r
+DialogGetFilterParameters::~DialogGetFilterParameters (void)\r
+{\r
+}\r
+\r
+\r
+unsigned int\r
+DialogGetFilterParameters::getXSize (void)\r
+{\r
+  wxString strCtrl = m_pTextCtrlXSize->GetValue();\r
+  unsigned long lValue;\r
+  if (strCtrl.ToULong (&lValue))\r
+    return lValue;\r
+  else\r
+    return (m_iDefaultXSize);\r
+}\r
+\r
+unsigned int\r
+DialogGetFilterParameters::getYSize (void)\r
+{\r
+  wxString strCtrl = m_pTextCtrlYSize->GetValue();\r
+  unsigned long lValue;\r
+  if (strCtrl.ToULong (&lValue))\r
+    return lValue;\r
+  else\r
+    return (m_iDefaultYSize);\r
+}\r
+\r
+double\r
+DialogGetFilterParameters::getBandwidth (void)\r
+{\r
+  wxString strCtrl = m_pTextCtrlBandwidth->GetValue();\r
+  double dValue;\r
+  if (strCtrl.ToDouble (&dValue))\r
+    return dValue;\r
+  else\r
+    return (m_dDefaultBandwidth);\r
+}\r
+\r
+double\r
+DialogGetFilterParameters::getFilterParam (void)\r
+{\r
+  wxString strCtrl = m_pTextCtrlFilterParam->GetValue();\r
+  double dValue;\r
+  if (strCtrl.ToDouble (&dValue))\r
+    return (dValue);\r
+  else\r
+    return (m_dDefaultFilterParam);\r
+}\r
+\r
+const char*\r
+DialogGetFilterParameters::getFilterName (void)\r
+{\r
+  return m_pListBoxFilter->getSelectionStringValue();\r
+}\r
+\r
+const char*\r
+DialogGetFilterParameters::getDomainName (void)\r
+{\r
+  return m_pListBoxDomain->getSelectionStringValue();\r
+}\r
+\r
index 4441ca0..3546021 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.h,v 1.16 2000/12/22 04:18:00 kevin Exp $
+**  $Id: dialogs.h,v 1.17 2001/01/01 10:14:34 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
@@ -167,46 +167,76 @@ class DialogGetProjectionParameters : public wxDialog
 
 
 #include "backprojectors.h"
-class DialogGetReconstructionParameters : public wxDialog
-{
- public:
-    DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGeneration = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3, int iDefaultTrace = Trace::TRACE_NONE);
-    virtual ~DialogGetReconstructionParameters ();
-
-    unsigned int getXSize();
-    unsigned int getYSize();
-    const char* getFilterName();
-    double getFilterParam();
-    const char* getFilterMethodName();
-    unsigned int getZeropad();
-    const char* getFilterGenerationName();
-    const char* getInterpName();
-    unsigned int getInterpParam();
-    const char* getBackprojectName();
-    int getTrace ();
-
- private:
-    wxTextCtrl* m_pTextCtrlXSize;
-    wxTextCtrl* m_pTextCtrlYSize;
-    wxTextCtrl* m_pTextCtrlZeropad;
-    wxTextCtrl* m_pTextCtrlFilterParam;
-    wxTextCtrl* m_pTextCtrlInterpParam;
-
-    StringValueAndTitleListBox* m_pListBoxFilter;
-    StringValueAndTitleListBox* m_pListBoxFilterMethod;
-    StringValueAndTitleListBox* m_pListBoxFilterGeneration;
-    StringValueAndTitleListBox* m_pListBoxInterp;
-    StringValueAndTitleListBox* m_pListBoxBackproject;
-    StringValueAndTitleListBox* m_pListBoxTrace;
-
-    int m_iDefaultXSize;
-    int m_iDefaultYSize;
-    double m_dDefaultFilterParam;
-    int m_iDefaultZeropad;
-    int m_iDefaultInterpParam;
-    int m_iDefaultTrace;
-};
+class DialogGetReconstructionParameters : public wxDialog\r
+{\r
+ public:\r
+    DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGeneration = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3, int iDefaultTrace = Trace::TRACE_NONE);\r
+    virtual ~DialogGetReconstructionParameters ();\r
+\r
+    unsigned int getXSize();\r
+    unsigned int getYSize();\r
+    const char* getFilterName();\r
+    double getFilterParam();\r
+    const char* getFilterMethodName();\r
+    unsigned int getZeropad();\r
+    const char* getFilterGenerationName();\r
+    const char* getInterpName();\r
+    unsigned int getInterpParam();\r
+    const char* getBackprojectName();\r
+    int getTrace ();\r
+\r
+ private:\r
+    wxTextCtrl* m_pTextCtrlXSize;\r
+    wxTextCtrl* m_pTextCtrlYSize;\r
+    wxTextCtrl* m_pTextCtrlZeropad;\r
+    wxTextCtrl* m_pTextCtrlFilterParam;\r
+    wxTextCtrl* m_pTextCtrlInterpParam;\r
+\r
+    StringValueAndTitleListBox* m_pListBoxFilter;\r
+    StringValueAndTitleListBox* m_pListBoxFilterMethod;\r
+    StringValueAndTitleListBox* m_pListBoxFilterGeneration;\r
+    StringValueAndTitleListBox* m_pListBoxInterp;\r
+    StringValueAndTitleListBox* m_pListBoxBackproject;\r
+    StringValueAndTitleListBox* m_pListBoxTrace;\r
+\r
+    int m_iDefaultXSize;\r
+    int m_iDefaultYSize;\r
+    double m_dDefaultFilterParam;\r
+    int m_iDefaultZeropad;\r
+    int m_iDefaultInterpParam;\r
+    int m_iDefaultTrace;\r
+};\r
+\r
 
+class DialogGetFilterParameters : public wxDialog\r
+{\r
+ public:\r
+    DialogGetFilterParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_BANDLIMIT, double dDefaultFilterParam = 1., double dDefaultBandwidth = 1., int iDefaultDomain = SignalFilter::DOMAIN_SPATIAL);\r
+    virtual ~DialogGetFilterParameters ();\r
+\r
+    unsigned int getXSize();\r
+    unsigned int getYSize();\r
+    const char* getFilterName();\r
+    const char* getDomainName();\r
+    double getFilterParam();\r
+    double getBandwidth();\r
+\r
+ private:\r
+    wxTextCtrl* m_pTextCtrlXSize;\r
+    wxTextCtrl* m_pTextCtrlYSize;\r
+    wxTextCtrl* m_pTextCtrlFilterParam;\r
+    wxTextCtrl* m_pTextCtrlBandwidth;\r
+\r
+    StringValueAndTitleListBox* m_pListBoxFilter;\r
+    StringValueAndTitleListBox* m_pListBoxDomain;\r
+\r
+    int m_iDefaultXSize;\r
+    int m_iDefaultYSize;\r
+    double m_dDefaultFilterParam;\r
+    double m_dDefaultBandwidth;\r
+    int m_iDefaultDomain;\r
+};\r
+\r
 class DialogAutoScaleParameters : public wxDialog
 {
  public:
index 86e0cf5..93d49b3 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.42 2000/12/29 20:18:59 kevin Exp $
+**  $Id: views.cpp,v 1.43 2001/01/01 10:14:34 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
@@ -129,7 +129,16 @@ ImageFileCanvas::OnMouseEvent(wxMouseEvent& event)
   if (event.RightIsDown()) {\r
     if (pt.x >= 0 && pt.x < nx && pt.y >= 0 && pt.y < ny) {\r
       std::ostringstream os;\r
-      os << "Image value (" << pt.x << "," << yPt << ") = " << v[pt.x][yPt] << "\n";\r
+      os << "Image value (" << pt.x << "," << yPt << ") = " << v[pt.x][yPt];\r
+      if (rIF.isComplex()) {\r
+        double dImag = rIF.getImaginaryArray()[pt.x][yPt];\r
+        if (dImag < 0)\r
+          os << " - " << -dImag;\r
+        else\r
+          os << " + " << dImag;\r
+        os << "i\n";\r
+      } else\r
+        os << "\n";\r
       *theApp->getLog() << os.str().c_str();\r
     } else\r
       *theApp->getLog() << "Mouse out of image range (" << pt.x << "," << yPt << ")\n";\r
@@ -170,6 +179,10 @@ EVT_MENU(IFMENU_PROCESS_LOG, ImageFileView::OnLog)
 EVT_MENU(IFMENU_PROCESS_EXP, ImageFileView::OnExp)\r
 EVT_MENU(IFMENU_PROCESS_FOURIER, ImageFileView::OnFourier)\r
 EVT_MENU(IFMENU_PROCESS_INVERSE_FOURIER, ImageFileView::OnInverseFourier)\r
+#ifdef HAVE_FFTW\r
+EVT_MENU(IFMENU_PROCESS_FFT, ImageFileView::OnFFT)\r
+EVT_MENU(IFMENU_PROCESS_IFFT, ImageFileView::OnIFFT)\r
+#endif\r
 EVT_MENU(IFMENU_PROCESS_MAGNITUDE, ImageFileView::OnMagnitude)\r
 EVT_MENU(IFMENU_PROCESS_PHASE, ImageFileView::OnPhase)\r
 EVT_MENU(IFMENU_PLOT_ROW, ImageFileView::OnPlotRow)\r
@@ -195,16 +208,16 @@ ImageFileView::OnProperties (wxCommandEvent& event)
     const std::string& rFilename = rIF.getFilename();
     std::ostringstream os;
     double min, max, mean, mode, median, stddev;\r
-    rIF.statistics (min, max, mean, mode, median, stddev);\r
+    rIF.statistics (rIF.getArray(), min, max, mean, mode, median, stddev);\r
     os << "Filename: " << rFilename << "\n";\r
     os << "Size: (" << rIF.nx() << "," << rIF.ny() << ")\n";\r
     os << "Data type: ";\r
-    if (rIF.dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
+    if (rIF.isComplex())\r
       os << "Complex\n";\r
     else\r
       os << "Real\n";\r
     os << "\nMinimum: "<<min<<"\nMaximum: "<<max<<"\nMean: "<<mean<<"\nMedian: "<<median<<"\nMode: "<<mode<<"\nStandard Deviation: "<<stddev << "\n";
-    if (rIF.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
+    if (rIF.isComplex()) {\r
       rIF.statistics (rIF.getImaginaryArray(), min, max, mean, mode, median, stddev);\r
       os << "\nImaginary: min: "<<min<<"\nmax: "<<max<<"\nmean: "<<mean<<"\nmedian: "<<median<<"\nmode: "<<mode<<"\nstddev: "<<stddev << "\n";\r
     }\r
@@ -279,7 +292,7 @@ ImageFileView::OnCompare (wxCommandEvent& event)
       rIF.statistics (min, max, mean, mode, median, stddev);\r
       os << rIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";\r
       rCompareIF.statistics (min, max, mean, mode, median, stddev);\r
-      os << rCompareIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";\r
+      os << pCompareDoc->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";\r
       os << "\n";\r
       double d, r, e;\r
       rIF.comparativeStatistics (rCompareIF, d, r, e);\r
@@ -363,6 +376,7 @@ void
 ImageFileView::OnFourier (wxCommandEvent& event)\r
 {\r
   ImageFile& rIF = GetDocument()->getImageFile();\r
+  wxProgressDialog dlgProgress (wxString("Fourier"), wxString("Fourier Progress"), 1, m_frame, wxPD_APP_MODAL);\r
   rIF.fourier (rIF);\r
   m_bMinSpecified = false;\r
   m_bMaxSpecified = false;\r
@@ -371,10 +385,39 @@ ImageFileView::OnFourier (wxCommandEvent& event)
   GetDocument()->UpdateAllViews(this);\r
 }\r
 \r
+#ifdef HAVE_FFTW\r
+void\r
+ImageFileView::OnFFT (wxCommandEvent& event)\r
+{\r
+  ImageFile& rIF = GetDocument()->getImageFile();\r
+  wxProgressDialog dlgProgress (wxString("FFT"), wxString("FFT Progress"), 1, m_frame, wxPD_APP_MODAL);\r
+  rIF.fft (rIF);\r
+  m_bMinSpecified = false;\r
+  m_bMaxSpecified = false;\r
+  if (theApp->getSetModifyNewDocs())\r
+    GetDocument()->Modify(TRUE);\r
+  GetDocument()->UpdateAllViews(this);\r
+}\r
+\r
+void\r
+ImageFileView::OnIFFT (wxCommandEvent& event)\r
+{\r
+  ImageFile& rIF = GetDocument()->getImageFile();\r
+  wxProgressDialog dlgProgress (wxString("IFFT"), wxString("IFFT Progress"), 1, m_frame, wxPD_APP_MODAL);\r
+  rIF.ifft (rIF);\r
+  m_bMinSpecified = false;\r
+  m_bMaxSpecified = false;\r
+  if (theApp->getSetModifyNewDocs())\r
+    GetDocument()->Modify(TRUE);\r
+  GetDocument()->UpdateAllViews(this);\r
+}\r
+#endif\r
+\r
 void\r
 ImageFileView::OnInverseFourier (wxCommandEvent& event)\r
 {\r
   ImageFile& rIF = GetDocument()->getImageFile();\r
+  wxProgressDialog dlgProgress (wxString("Inverse Fourier"), wxString("Inverse Fourier Progress"), 1, m_frame, wxPD_APP_MODAL);\r
   rIF.inverseFourier (rIF);\r
   m_bMinSpecified = false;\r
   m_bMaxSpecified = false;\r
@@ -431,7 +474,8 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   
   wxMenu *file_menu = new wxMenu;
   
-  file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");
+  file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");\r
+  file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");\r
   file_menu->Append(wxID_OPEN, "&Open...");
   file_menu->Append(wxID_SAVE, "&Save");
   file_menu->Append(wxID_SAVEAS, "Save &As...");
@@ -456,8 +500,15 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   process_menu->Append (IFMENU_PROCESS_LOG, "&Log");\r
   process_menu->Append (IFMENU_PROCESS_EXP, "&Exp");\r
   process_menu->AppendSeparator();\r
+#ifdef HAVE_FFTW\r
+  process_menu->Append (IFMENU_PROCESS_FFT, "&FFT");\r
+  process_menu->Append (IFMENU_PROCESS_IFFT, "&IFFT");\r
+  process_menu->Append (IFMENU_PROCESS_FOURIER, "F&ourier");\r
+  process_menu->Append (IFMENU_PROCESS_INVERSE_FOURIER, "Inverse Fo&urier");\r
+#else\r
   process_menu->Append (IFMENU_PROCESS_FOURIER, "&Fourier");\r
   process_menu->Append (IFMENU_PROCESS_INVERSE_FOURIER, "&Inverse Fourier");\r
+#endif\r
   process_menu->Append (IFMENU_PROCESS_MAGNITUDE, "&Magnitude");\r
   process_menu->Append (IFMENU_PROCESS_PHASE, "&Phase");\r
 \r
@@ -1050,6 +1101,7 @@ PhantomView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenu *file_menu = new wxMenu;
   
   file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");
+  file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");\r
   file_menu->Append(wxID_OPEN, "&Open...");
   file_menu->Append(wxID_CLOSE, "&Close");
   
@@ -1339,6 +1391,7 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenu *file_menu = new wxMenu;
   
   file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");
+  file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");\r
   file_menu->Append(wxID_OPEN, "&Open...");
   file_menu->Append(wxID_SAVE, "&Save");
   file_menu->Append(wxID_SAVEAS, "Save &As...");
@@ -1601,6 +1654,7 @@ PlotFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenu *file_menu = new wxMenu;
   
   file_menu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...");
+  file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");\r
   file_menu->Append(wxID_OPEN, "&Open...");
   file_menu->Append(wxID_SAVE, "&Save");
   file_menu->Append(wxID_SAVEAS, "Save &As...");
index 3512fac..8d46ae3 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.h,v 1.19 2000/12/29 20:18:59 kevin Exp $
+**  $Id: views.h,v 1.20 2001/01/01 10:14:34 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
@@ -73,6 +73,10 @@ public:
     void OnExp (wxCommandEvent& event);\r
     void OnFourier (wxCommandEvent& event);\r
     void OnInverseFourier (wxCommandEvent& event);\r
+#ifdef HAVE_FFTW\r
+    void OnFFT (wxCommandEvent& event);\r
+    void OnIFFT (wxCommandEvent& event);\r
+#endif\r
     void OnMagnitude (wxCommandEvent& event);\r
     void OnPhase (wxCommandEvent& event);\r
     void OnScaleAuto (wxCommandEvent& event);