From: Kevin M. Rosenberg <kevin@rosenberg.net>
Date: Mon, 1 Jan 2001 10:14:34 +0000 (+0000)
Subject: r326: FFTW additions, filter image generation
X-Git-Tag: debian-4.5.3-3~691
X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=7ec2cd66921180a624813dff9f8bac76c6b268cc

r326: FFTW additions, filter image generation
---

diff --git a/ChangeLog b/ChangeLog
index cf938a0..9decb06 100644
--- 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
 
diff --git a/include/array2dfile.h b/include/array2dfile.h
index 754ad93..22e8b61 100644
--- a/include/array2dfile.h
+++ b/include/array2dfile.h
@@ -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; }
 
+  bool isComplex() const
+  { return m_dataType == DATA_TYPE_COMPLEX; }
+
+  bool isReal() const
+  { return m_dataType == DATA_TYPE_REAL; }
+
   int dataType () const
   { return static_cast<int>(m_dataType); }
 
diff --git a/include/ctsupport.h b/include/ctsupport.h
index 81479c9..84eef25 100644
--- a/include/ctsupport.h
+++ b/include/ctsupport.h
@@ -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)); }
-
+
+inline bool isEven (int n)
+{ return (n % 2) == 0; }
+
+inline bool isOdd (int n)
+{ return (n % 2) != 0; }
+
+#if 0
+inline bool isEven (long n)
+{ return (n % 2) == 0; }
+
+inline bool isOdd (long n)
+{ return (n % 2) != 0; }
+#endif
+
 inline int imax (int a, int b)
 { return (a >= b ? a : b); }
 
diff --git a/include/imagefile.h b/include/imagefile.h
index 57d093f..945bdbf 100644
--- a/include/imagefile.h
+++ b/include/imagefile.h
@@ -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;
-  void statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const;
+  void statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const;
   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;
   bool fourier (ImageFile& result) const;
   bool inverseFourier (ImageFile& result) const;
+#ifdef HAVE_FFTW
+  bool fft (ImageFile& result) const;
+  bool ifft (ImageFile& result) const;
+#endif
   bool magnitude (ImageFile& result) const;
   bool phase (ImageFile& result) const;
 
diff --git a/include/procsignal.h b/include/procsignal.h
index a160154..8eb98c6 100644
--- a/include/procsignal.h
+++ b/include/procsignal.h
@@ -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;
+
+typedef std::complex<double> CTSimComplex;
+
 
 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);
     static void shuffleNaturalToFourierOrder (double* pdVector, const int n);
     static void shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n);
+    static void shuffleFourierToNaturalOrder (float* pdVector, const int n);
     static void shuffleFourierToNaturalOrder (double* pdVector, const int n);
     static void shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n);
+
 
  private:
 	 std::string m_nameFilterMethod;
diff --git a/libctsim/array2dfile.cpp b/libctsim/array2dfile.cpp
index ca45d8b..d0f2d68 100644
--- a/libctsim/array2dfile.cpp
+++ b/libctsim/array2dfile.cpp
@@ -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();
 }
 
 
diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 7926b0d..726e5b8 100644
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -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"},
diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index b6525e1..df647ce 100644
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.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)
 : Array2dFile()
 {
-	setPixelFormat (Array2dFile::PIXEL_FLOAT32);
-	setPixelSize (sizeof(kfloat32));
+  setPixelFormat (Array2dFile::PIXEL_FLOAT32);
+  setPixelSize (sizeof(kfloat32));
   setDataType (Array2dFile::DATA_TYPE_REAL);
 }
 
@@ -49,78 +49,98 @@ F64Image::F64Image (int nx, int ny, int dataType)
 F64Image::F64Image (void)
 : Array2dFile ()
 {
-	setPixelFormat (PIXEL_FLOAT64);
-	setPixelSize (sizeof(kfloat64));
+  setPixelFormat (PIXEL_FLOAT64);
+  setPixelSize (sizeof(kfloat64));
   setDataType (Array2dFile::DATA_TYPE_REAL);
 }
 
 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++) {
+{
+  ImageFileArray v = getArray();
+  SignalFilter filter (filterName, domainName, bw, filt_param);
+
+#if 1
+  int iXCenter, iYCenter;
+  if (isEven (m_nx))
+    iXCenter = m_nx / 2;
+  else
+    iXCenter = (m_nx - 1) / 2;
+  if (isEven (m_ny))
+    iYCenter = m_ny / 2;
+  else
+    iYCenter = (m_ny - 1) / 2;
+
+  for (int ix = 0; ix < m_nx; ix++)
+    for (int iy = 0; iy < m_ny; iy++) {
+      long lD2 = (ix - iXCenter) * (ix - iXCenter) + (iy - iYCenter) * (iy - iYCenter);
+      double r = ::sqrt (static_cast<double>(lD2));
+      v[ix][iy] = filter.response (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);
+    }
+  }
+#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);
-    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
-      os << "Real Component Statistics" << std::endl;
-
-    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;
-
-    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
-      statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
-      os << std::endl << "Imaginary Component Statistics" << std::endl;
+  double min, max, mean, mode, median, stddev;
+  
+  statistics (min, max, mean, mode, median, stddev);
+  if (isComplex())
+    os << "Real Component Statistics" << std::endl;
+  
+  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;
+  
+  if (isComplex()) {
+    statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
+    os << std::endl << "Imaginary Component Statistics" << std::endl;
     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;
-    }
+  }
 }
 
 
@@ -247,45 +267,45 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou
 
 
 void
-ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
+ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 {
-    int nx = m_nx;
-    int ny = m_ny;
-    
-    if (v == NULL || nx == 0 || ny == 0)
-		return;
-
-	std::vector<double> vecImage;
-	int iVec = 0;
-	vecImage.resize (nx * ny);
-    for (int ix = 0; ix < nx; ix++) {
-		for (int iy = 0; iy < ny; iy++)
-			vecImage[iVec++] = v[ix][iy];
-	}
-
-	vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
+  int nx = m_nx;
+  int ny = m_ny;
+  
+  if (v == NULL || nx == 0 || ny == 0)
+    return;
+  
+  std::vector<double> vecImage;
+  int iVec = 0;
+  vecImage.resize (nx * ny);
+  for (int ix = 0; ix < nx; ix++) {
+    for (int iy = 0; iy < ny; iy++)
+      vecImage[iVec++] = v[ix][iy];
+  }
+  
+  vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
 }
 
 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];
     }
+  }
 }
 
 bool
@@ -293,17 +313,17 @@ ImageFile::convertRealToComplex ()
 {
   if (dataType() != Array2dFile::DATA_TYPE_REAL)
     return false;
-
+  
   if (! reallocRealToComplex())
     return false;
-
+  
   ImageFileArray vImag = getImaginaryArray();
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumn vCol = vImag[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
       *vCol++ = 0;
   }
-
+  
   return true;
 }
 
@@ -312,19 +332,19 @@ ImageFile::convertComplexToReal ()
 {
   if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)
     return false;
-
+  
   ImageFileArray vReal = getArray();
   ImageFileArray vImag = getImaginaryArray();
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumn vRealCol = vReal[ix];
     ImageFileColumn vImagCol = vImag[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++) {
-      std::complex<double> c (*vRealCol, *vImagCol);
+      CTSimComplex c (*vRealCol, *vImagCol);
       *vRealCol++ = std::abs (c);
       vImagCol++;
     }
   }
-
+  
   return reallocComplexToReal();
 }
 
@@ -335,20 +355,20 @@ ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
-        *out++ = *in1++ - *in2++;
+      *out++ = *in1++ - *in2++;
   }
-
-    return true;
+  
+  return true;
 }
 
 bool
@@ -358,20 +378,20 @@ ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
-        *out++ = *in1++ + *in2++;
+      *out++ = *in1++ + *in2++;
   }
-
-    return true;
+  
+  return true;
 }
 
 bool
@@ -381,20 +401,20 @@ ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
-        *out++ = *in1++ * *in2++;
+      *out++ = *in1++ * *in2++;
   }
-
-    return true;
+  
+  return true;
 }
 
 bool
@@ -404,11 +424,11 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
@@ -420,8 +440,8 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
         *out++ = 0;
     }
   }
-
-    return true;
+  
+  return true;
 }
 
 
@@ -432,17 +452,17 @@ ImageFile::invertPixelValues (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
-        *out++ = - *in++;
+      *out++ = - *in++;
   }
-
+  
   return true;
 }
 
@@ -453,10 +473,10 @@ ImageFile::sqrt (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
@@ -466,7 +486,7 @@ ImageFile::sqrt (ImageFile& result) const
       else
         *out++ = ::sqrt(*in++);
   }
-
+  
   return true;
 }
 
@@ -477,10 +497,10 @@ ImageFile::log (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
@@ -490,7 +510,7 @@ ImageFile::log (ImageFile& result) const
       else
         *out++ = ::log(*in++);
   }
-
+  
   return true;
 }
 
@@ -501,33 +521,153 @@ ImageFile::exp (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++)
       *out++ = ::exp (*in++);
   }
-
+  
   return true;
 }
 
+namespace ProcessImage {
+
+void
+shuffleFourierToNaturalOrder (ImageFile& im)
+{
+  ImageFileArray vReal = im.getArray();
+  ImageFileArray vImag = im.getImaginaryArray();
+  unsigned int ix, iy;
+  unsigned int nx = im.nx();
+  unsigned int ny = im.ny();
+
+  // shuffle each column
+  for (ix = 0; ix < nx; ix++) {
+    ProcessSignal::shuffleFourierToNaturalOrder (vReal[ix], ny);
+    if (im.isComplex())
+      ProcessSignal::shuffleFourierToNaturalOrder (vImag[ix], ny);
+  }
+
+  // shuffle each row
+  float* pRow = new float [nx];
+  for (iy = 0; iy < ny; iy++) {
+    for (ix = 0; ix < nx; ix++)
+      pRow[ix] = vReal[ix][iy];
+    ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);
+    for (ix = 0; ix < nx; ix++)
+      vReal[ix][iy] = pRow[ix];
+    if (im.isComplex()) {
+      for (ix = 0; ix < nx; ix++)
+        pRow[ix] = vImag[ix][iy];
+      ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);;
+      for (ix = 0; ix < nx; ix++)
+        vImag[ix][iy] = pRow[ix];
+    }
+  }
+  delete pRow;
+}
+ 
+void
+shuffleNaturalToFourierOrder (ImageFile& im)
+{
+  ImageFileArray vReal = im.getArray();
+  ImageFileArray vImag = im.getImaginaryArray();
+  unsigned int ix, iy;
+  unsigned int nx = im.nx();
+  unsigned int ny = im.ny();
+
+  // shuffle each x column
+  for (ix = 0; ix < nx; ix++) {
+    ProcessSignal::shuffleNaturalToFourierOrder (vReal[ix], ny);
+    if (im.isComplex())
+      ProcessSignal::shuffleNaturalToFourierOrder (vImag[ix], ny);
+  }
+
+  // shuffle each y row
+  float* pRow = new float [nx];
+  for (iy = 0; iy < ny; iy++) {
+    for (ix = 0; ix < nx; ix++)
+      pRow[ix] = vReal[ix][iy];
+    ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);
+    for (ix = 0; ix < nx; ix++)
+      vReal[ix][iy] = pRow[ix];
+    if (im.isComplex()) {
+      for (ix = 0; ix < nx; ix++)
+        pRow[ix] = vImag[ix][iy];
+      ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);
+      for (ix = 0; ix < nx; ix++)
+        vImag[ix][iy] = pRow[ix];
+    }
+  }
+  delete [] pRow;
+}
+
+ 
+}; // namespace ProcessIamge
+
+#ifdef HAVE_FFTW
 bool
-ImageFile::inverseFourier (ImageFile& result) const
+ImageFile::fft (ImageFile& result) const
 {
   if (m_nx != result.nx() || m_ny != result.ny()) {
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
 
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.reallocRealToComplex ())
+      return false;
+  }
+
+  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+
+  ImageFileArrayConst vReal = getArray();
+  ImageFileArrayConst vImag = getImaginaryArray();
+
+  unsigned int ix, iy;
+  unsigned int iArray = 0;
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      in[iArray].re = vReal[ix][iy];
+      if (isComplex())
+        in[iArray].im = vImag[ix][iy];
+      else
+        in[iArray].im = 0;
+      iArray++;
+    }
+
+  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
+
+  fftwnd_one (plan, in, NULL);
+
+  ImageFileArray vRealResult = result.getArray();
+  ImageFileArray vImagResult = result.getImaginaryArray();
+  iArray = 0;
+  unsigned int iScale = m_nx * m_ny;
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = in[iArray].re / iScale;
+      vImagResult[ix][iy] = in[iArray].im / iScale;
+      iArray++;
+    }
+
+  fftwnd_destroy_plan (plan);
+  delete in;
+
+
+  ProcessImage::shuffleFourierToNaturalOrder (result);
+
   return true;
 }
 
+
 bool
-ImageFile::fourier (ImageFile& result) const
+ImageFile::ifft (ImageFile& result) const
 {
   if (m_nx != result.nx() || m_ny != result.ny()) {
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
@@ -539,62 +679,230 @@ ImageFile::fourier (ImageFile& result) const
       return false;
   }
 
-  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vReal = getArray();
+  ImageFileArrayConst vImag = getImaginaryArray();
   ImageFileArray vRealResult = result.getArray();
   ImageFileArray vImagResult = result.getImaginaryArray();
-
   unsigned int ix, iy;
-  double* pY = new double [m_ny];
-  std::complex<double>** complexOut = new std::complex<double>* [m_nx];
   for (ix = 0; ix < m_nx; ix++)
-    complexOut[ix] = new std::complex<double> [m_ny];
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = vReal[ix][iy];
+      if (isComplex()) 
+        vImagResult[ix][iy] = vImag[ix][iy];
+      else
+        vImagResult[ix][iy] = 0;
+    }
+
+  ProcessImage::shuffleNaturalToFourierOrder (result);
 
+  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+
+  unsigned int iArray = 0;
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      in[iArray].re = vRealResult[ix][iy];
+      in[iArray].im = vImagResult[ix][iy];
+      iArray++;
+    }
+
+  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
+
+  fftwnd_one (plan, in, NULL);
+
+  iArray = 0;
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = in[iArray].re;
+      vImagResult[ix][iy] = in[iArray].im;
+      iArray++;
+    }
+
+  fftwnd_destroy_plan (plan);
+
+  delete in;
+
+  return true;
+}
+#endif // HAVE_FFTW
+
+
+  
+bool
+ImageFile::fourier (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+    return false;
+  }
+  
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.reallocRealToComplex ())
+      return false;
+  }
+  
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vLHSImag = getImaginaryArray();
+  ImageFileArray vRealResult = result.getArray();
+  ImageFileArray vImagResult = result.getImaginaryArray();
+  
+  unsigned int ix, iy;
+
+  // alloc output matrix
+  CTSimComplex** complexOut = new CTSimComplex* [m_nx];
+  for (ix = 0; ix < m_nx; ix++)
+    complexOut[ix] = new CTSimComplex [m_ny];
+  
+  // fourier each x column
+#if 1
+  CTSimComplex* pY = new CTSimComplex [m_ny];
   for (ix = 0; ix < m_nx; ix++) {
-    for (iy = 0; iy < m_ny; iy++)
+    for (iy = 0; iy < m_ny; iy++) {
+      pY[iy].real (vLHS[ix][iy]);
+      if (isComplex())
+        pY[iy].imag (vLHSImag[ix][iy]);
+      else
+        pY[iy].imag (0);
+    } 
+    ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
+  }
+#else
+  double* pY = new double [m_ny];
+  for (ix = 0; ix < m_nx; ix++) {
+    for (iy = 0; iy < m_ny; iy++) {
       pY[iy] = vLHS[ix][iy];
+    } 
     ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
   }
-  delete pY;
-
-  std::complex<double>* pX = new std::complex<double> [m_nx];
-  std::complex<double>* complexOutCol = new std::complex<double> [m_nx];
+#endif
+  delete [] pY;
+  
+  // fourier each y row
+  CTSimComplex* pX = new CTSimComplex [m_nx];
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
   for (iy = 0; iy < m_ny; iy++) {
     for (ix = 0; ix < m_nx; ix++)
       pX[ix] = complexOut[ix][iy];
-    ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD);
+    ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
     for (ix = 0; ix < m_nx; ix++)
-      complexOut[ix][iy] = complexOutCol[ix];
+      complexOut[ix][iy] = complexOutRow[ix];
   }
   delete [] pX;
-
+  
+  // shuffle each column
   for (ix = 0; ix < m_nx; ix++)
     ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);
+  
+  // shuffle each row
+  for (iy = 0; iy < m_ny; iy++) {
+    for (ix = 0; ix < m_nx; ix++)
+      complexOutRow[ix] = complexOut[ix][iy];
+    ProcessSignal::shuffleFourierToNaturalOrder (complexOutRow, m_nx);;
+    for (ix = 0; ix < m_nx; ix++)
+      complexOut[ix][iy] = complexOutRow[ix];
+    
+  }
+  delete [] complexOutRow;
+  
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = complexOut[ix][iy].real();
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();
+    }
+    
+    for (ix = 0; ix < m_nx; ix++)
+      delete [] complexOut[ix];
+    delete [] complexOut;
+    
+    return true;
+}
 
+bool
+ImageFile::inverseFourier (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+    return false;
+  }
+  
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.reallocRealToComplex ())
+      return false;
+  }
+  
+  ImageFileArrayConst vLHSReal = getArray();
+  ImageFileArrayConst vLHSImag = getImaginaryArray();
   
+  unsigned int ix, iy;
+
+  // alloc 2d complex output matrix
+  CTSimComplex** complexOut = new CTSimComplex* [m_nx];
+  for (ix = 0; ix < m_nx; ix++)
+    complexOut[ix] = new CTSimComplex [m_ny];
+
+  // put input image into complexOut
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      complexOut[ix][iy].real (vLHSReal[ix][iy]);
+      if (isComplex())
+        complexOut[ix][iy].imag (vLHSImag[ix][iy]);
+      else
+        complexOut[ix][iy].imag (0);
+    }
+
+  // shuffle each x column
+  for (ix = 0; ix < m_nx; ix++)
+    ProcessSignal::shuffleNaturalToFourierOrder (complexOut[ix], m_ny);
+
+  // shuffle each y row
+  CTSimComplex* pComplexRow = new CTSimComplex [m_nx];
   for (iy = 0; iy < m_ny; iy++) {
     for (ix = 0; ix < m_nx; ix++)
-      complexOutCol[ix] = complexOut[ix][iy];
-    ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);;
+      pComplexRow[ix] = complexOut[ix][iy];
+    ProcessSignal::shuffleNaturalToFourierOrder (pComplexRow, m_nx);
     for (ix = 0; ix < m_nx; ix++)
-      complexOut[ix][iy] = complexOutCol[ix];
+      complexOut[ix][iy] = pComplexRow[ix];
+  }
+  delete [] pComplexRow;
 
+  // ifourier each x column
+  CTSimComplex* pCol = new CTSimComplex [m_ny];
+  for (ix = 0; ix < m_nx; ix++) {
+    for (iy = 0; iy < m_ny; iy++)
+      pCol[iy] = complexOut[ix][iy];
+    ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);
   }
-  delete [] complexOutCol;
+  delete [] pCol;
 
+  // ifourier each y row
+  CTSimComplex* complexInRow = new CTSimComplex [m_nx];
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+  for (iy = 0; iy < m_ny; iy++) {
+    for (ix = 0; ix < m_nx; ix++)
+      complexInRow[ix] = complexOut[ix][iy];
+    ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
+    for (ix = 0; ix < m_nx; ix++)
+      complexOut[ix][iy] = complexOutRow[ix];
+  }
+  delete [] complexInRow;
+  delete [] complexOutRow;
+
+  ImageFileArray vRealResult = result.getArray();
+  ImageFileArray vImagResult = result.getImaginaryArray();
   for (ix = 0; ix < m_nx; ix++)
     for (iy = 0; iy < m_ny; iy++) {
       vRealResult[ix][iy] = complexOut[ix][iy].real();
       vImagResult[ix][iy] = complexOut[ix][iy].imag();
     }
-
+    
+  // delete complexOut matrix
   for (ix = 0; ix < m_nx; ix++)
     delete [] complexOut[ix];
-
   delete [] complexOut;
-
+    
   return true;
 }
 
+
 bool
 ImageFile::magnitude (ImageFile& result) const
 {
@@ -602,29 +910,23 @@ ImageFile::magnitude (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArray vReal = getArray();
-  ImageFileArray vImag = NULL;
-  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
-    vImag = getImaginaryArray();
-
+  ImageFileArray vImag = getImaginaryArray();
+  
   ImageFileArray vRealResult = result.getArray();
-  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
-    ImageFileArray vImagResult = result.getImaginaryArray();
-    for (unsigned int ix = 0; ix < m_nx; ix++)
-      for (unsigned int iy = 0; iy < m_ny; iy++)
-        vImagResult[ix][iy] = 0;
-  }
-
   for (unsigned int ix = 0; ix < m_nx; ix++)
     for (unsigned int iy = 0; iy < m_ny; iy++) {
-      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) 
+      if (isComplex()) 
         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
       else
         vRealResult[ix][iy] = vReal[ix][iy];
     }
 
-  return true;
+  if (result.isComplex())
+    result.reallocComplexToReal();
+  
+    return true;
 }
 
 bool
@@ -634,28 +936,22 @@ ImageFile::phase (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArray vReal = getArray();
-  ImageFileArray vImag = NULL;
-  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
-    vImag = getImaginaryArray();
-
+  ImageFileArray vImag = getImaginaryArray();
+  
   ImageFileArray vRealResult = result.getArray();
-  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
-    ImageFileArray vImagResult = result.getImaginaryArray();
-    for (unsigned int ix = 0; ix < m_nx; ix++)
-      for (unsigned int iy = 0; iy < m_ny; iy++)
-        vImagResult[ix][iy] = 0;
-  }
-
   for (unsigned int ix = 0; ix < m_nx; ix++)
     for (unsigned int iy = 0; iy < m_ny; iy++) {
-      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) 
-        vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]);
+      if (isComplex())
+        vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
       else
-        vRealResult[ix][iy] = vReal[ix][iy];
+        vRealResult[ix][iy] = 0;
     }
-
+    
+  if (result.isComplex())
+    result.reallocComplexToReal();
+  
   return true;
 }
 
@@ -666,19 +962,19 @@ ImageFile::square (ImageFile& result) const
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
-
+  
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
-
+  
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
     for (unsigned int iy = 0; iy < m_ny; iy++) {
-        *out++ = *in * *in;
-        in++;
+      *out++ = *in * *in;
+      in++;
     }
   }
-
+  
   return true;
 }
 
@@ -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 );
-		}
-	}
-	
-	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 );
+    }
+  }
+  
+  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");
-		}
-	}
-	
-	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");
+    }
+  }
+  
+  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;
-	
-	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;
+  
+  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
 
diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index a7933b2..8c45da4 100644
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -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;
 #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];
-      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
+      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
       delete adSpatialFilter;
       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;
 #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);
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
     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);
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
@@ -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);
+    finiteFourierTransform (inputSignal, fftSignal, FORWARD);
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
@@ -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));
+      sum += input[j] * exponentTerm;
     }
     if (direction < 0) {
       sum /= n;
@@ -878,10 +878,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
   } else {     // Even
     int iHalfN = n / 2;
     pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
     for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1] = pdVector[i + iHalfN];
     for (i = 0; i < iHalfN - 1; i++)
       pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
   }
   
   for (i = 0; i < n; i++)
@@ -905,10 +912,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   } else {     // Even
     int iHalfN = n / 2;
     pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
     for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1] = pdVector[i + iHalfN];
     for (i = 0; i < iHalfN - 1; i++)
       pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
   }
   
   for (i = 0; i < n; i++)
@@ -916,7 +930,43 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;
 }
 
-
+
+void
+ProcessSignal::shuffleNaturalToFourierOrder (float* pdVector, const int n)
+{
+  float* pdTemp = new float [n];
+  int i;
+  if (n % 2) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pdTemp[0] = pdVector[iHalfN];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN + 1] = pdVector[i];
+  } else {     // Even
+    int iHalfN = n / 2;
+    pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN];
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
+  }
+  
+  for (i = 0; i < n; i++)
+    pdVector[i] = pdTemp[i];
+  delete pdTemp;
+}
+
+
+
 void
 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 {
@@ -944,6 +994,7 @@ ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
   delete pdTemp;
 }
 
+
 void
 ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
 {
@@ -971,3 +1022,33 @@ ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;
 }
 
+
+
+
+void
+ProcessSignal::shuffleFourierToNaturalOrder (float* pVector, const int n)
+{
+  float* pTemp = new float [n];
+  int i;
+  if (n % 2) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pTemp[iHalfN] = pVector[0];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i + 1 + iHalfN] = pVector[i + 1];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i] = pVector[i + iHalfN + 1];
+  } else {     // Even
+    int iHalfN = n / 2;
+    pTemp[iHalfN] = pVector[0];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i] = pVector[i + iHalfN];
+    for (i = 0; i < iHalfN - 1; i++)
+      pTemp[i + iHalfN + 1] = pVector[i+1];
+  }
+  
+  for (i = 0; i < n; i++)
+    pVector[i] = pTemp[i];
+  delete [] pTemp;
+}
+
diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp
index 2eaa551..3604b2f 100644
--- a/libctsupport/mathfuncs.cpp
+++ b/libctsupport/mathfuncs.cpp
@@ -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
@@ -22,18 +22,18 @@
 
 
 /* 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 
 vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)
 {
-	if (nPoints <= 0)
-		return;
-
-    mean = 0;
-    min = vec[0];
-    max = vec[0];
-	int i;
-    for (i = 0; i < nPoints; i++) {
-		double v = vec[i];
-		if (v > max)
-			max = v;
-		if (v < min)
-			min = v;
-		mean += v;
-    }
-    mean /= nPoints;
-	
-    static const int nbin = 1024;
-    int hist[ nbin ] = {0};
-    double spread = max - min;
-    mode = 0;
-    stddev = 0;
-    for (i = 0; i < nPoints; i++) {
-		double v = vec[i];
-		int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
-		hist[b]++;
-		double diff = (v - mean);
-		stddev += diff * diff;
+  if (nPoints <= 0)
+    return;
+  
+  mean = 0;
+  min = vec[0];
+  max = vec[0];
+  int i;
+  int iMinPos = 0;
+  for (i = 0; i < nPoints; i++) {
+    double v = vec[i];
+    if (v > max)
+      max = v;
+    if (v < min) {
+      min = v;
+      iMinPos = i;
     }
-    stddev = sqrt (stddev / nPoints);
-	
-    int max_binindex = 0;
-    int max_bin = -1;
-    for (int ibin = 0; ibin < nbin; ibin++) {
-		if (hist[ibin] > max_bin) {
-			max_bin = hist[ibin];
-			max_binindex = ibin;
-		}
+    mean += v;
+  }
+  mean /= nPoints;
+  
+  static const int nbin = 1024;
+  int hist[ nbin ] = {0};
+  double spread = max - min;
+  mode = 0;
+  stddev = 0;
+  for (i = 0; i < nPoints; i++) {
+    double v = vec[i];
+    int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
+    hist[b]++;
+    double diff = (v - mean);
+    stddev += diff * diff;
+  }
+  stddev = sqrt (stddev / nPoints);
+  
+  int max_binindex = 0;
+  int max_bin = -1;
+  for (int ibin = 0; ibin < nbin; ibin++) {
+    if (hist[ibin] > max_bin) {
+      max_bin = hist[ibin];
+      max_binindex = ibin;
     }
-	
-    mode = (max_binindex * spread / (nbin - 1)) + min;
-	
-	std::sort(vec.begin(), vec.end());
-	
-	if (nPoints % 2)  // Odd
-		median = vec[((nPoints - 1) / 2)];
-	else        // Even
-		median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
+  }
+  
+  mode = (max_binindex * spread / (nbin - 1)) + min;
+  
+  std::sort(vec.begin(), vec.end());
+  
+  if (nPoints % 2)  // Odd
+    median = vec[((nPoints - 1) / 2)];
+  else        // Even
+    median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
 }
diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg
index daff73f..fb8f744 100644
--- a/msvc/ctsim/ctsim.plg
+++ b/msvc/ctsim/ctsim.plg
@@ -3,303 +3,78 @@
 <pre>
 <h1>Build Log</h1>
 <h3>
---------------------Configuration: FFTW2st - Win32 Release--------------------
+--------------------Configuration: libctsim - Win32 Debug--------------------
 </h3>
 <h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-FFTW2st.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: FFTW2st - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-FFTW2st.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: RFFTW2st - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-RFFTW2st.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: RFFTW2st - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-RFFTW2st.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: ctsim - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-ctsim.exe - 0 error(s), 0 warning(s)
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29A.tmp" with contents
+[
+/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 
+"C:\ctsim-2.0.6\libctsim\imagefile.cpp"
+]
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29A.tmp" 
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29B.tmp" with contents
+[
+/nologo /out:"Debug\libctsim.lib" 
+".\Debug\array2dfile.obj"
+".\Debug\backprojectors.obj"
+".\Debug\clip.obj"
+".\Debug\consoleio.obj"
+".\Debug\ezplot.obj"
+".\Debug\ezset.obj"
+".\Debug\ezsupport.obj"
+".\Debug\filter.obj"
+".\Debug\fnetorderstream.obj"
+".\Debug\getopt.obj"
+".\Debug\getopt1.obj"
+".\Debug\hashtable.obj"
+".\Debug\imagefile.obj"
+".\Debug\mathfuncs.obj"
+".\Debug\phantom.obj"
+".\Debug\plotfile.obj"
+".\Debug\pol.obj"
+".\Debug\procsignal.obj"
+".\Debug\projections.obj"
+".\Debug\reconstruct.obj"
+".\Debug\scanner.obj"
+".\Debug\sgp.obj"
+".\Debug\strfuncs.obj"
+".\Debug\syserror.obj"
+".\Debug\trace.obj"
+".\Debug\transformmatrix.obj"
+".\Debug\xform.obj"
+]
+Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29B.tmp"
+<h3>Output Window</h3>
+Compiling...
+imagefile.cpp
+Creating library...
 <h3>
 --------------------Configuration: ctsim - Win32 Debug--------------------
 </h3>
 <h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-ctsim.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: if1 - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-if1.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: if1 - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-if1.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: if2 - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-if2.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: if2 - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-if2.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: ifexport - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-ifexport.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: ifexport - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPA0.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29C.tmp" with contents
 [
-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 
-".\Debug\ifexport.obj"
+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" 
+".\Debug\ctsim.obj"
+".\Debug\dialogs.obj"
+".\Debug\dlgprojections.obj"
+".\Debug\dlgreconstruct.obj"
+".\Debug\docs.obj"
+".\Debug\views.obj"
+".\Debug\wx.res"
 "\ctsim-2.0.6\msvc\libctsim\Debug\libctsim.lib"
 "\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib"
 "\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib"
+"\wx2\lib\wxd.lib"
 ]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPA0.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP29C.tmp"
 <h3>Output Window</h3>
 Linking...
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxPen::~wxPen(void)" (??1wxPen@@UAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxBrush *  wxWHITE_BRUSH" (?wxWHITE_BRUSH@@3PAVwxBrush@@A)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "char const * const  wxEmptyString" (?wxEmptyString@@3PBDB)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetWidth(int)" (?SetWidth@wxPen@@QAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxPen::wxPen(void)" (??0wxPen@@QAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxFont  wxNullFont" (?wxNullFont@@3VwxFont@@A)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxBrush::~wxBrush(void)" (??1wxBrush@@UAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetColour(unsigned char,unsigned char,unsigned char)" (?SetColour@wxPen@@QAEXEEE@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "class wxBrush  wxNullBrush" (?wxNullBrush@@3VwxBrush@@A)
-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)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxBrush::wxBrush(void)" (??0wxBrush@@QAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxColour::~wxColour(void)" (??1wxColour@@UAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxColour::wxColour(unsigned char,unsigned char,unsigned char)" (??0wxColour@@QAE@EEE@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetColour(class wxColour const &)" (?SetColour@wxPen@@QAEXABVwxColour@@@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: void __thiscall wxPen::SetStyle(int)" (?SetStyle@wxPen@@QAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "private: void __thiscall wxString::InitWith(char const *,unsigned int,unsigned int)" (?InitWith@wxString@@AAEXPBDII@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetEncoding(enum wxFontEncoding)" (?SetEncoding@wxFont@@UAEXW4wxFontEncoding@@@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetUnderlined(bool)" (?SetUnderlined@wxFont@@UAEX_N@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetFaceName(class wxString const &)" (?SetFaceName@wxFont@@UAEXABVwxString@@@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetWeight(int)" (?SetWeight@wxFont@@UAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetStyle(int)" (?SetStyle@wxFont@@UAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetFamily(int)" (?SetFamily@wxFont@@UAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxFont::SetPointSize(int)" (?SetPointSize@wxFont@@UAEXH@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual enum wxFontEncoding  __thiscall wxFont::GetEncoding(void)const " (?GetEncoding@wxFont@@UBE?AW4wxFontEncoding@@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual class wxString  __thiscall wxFont::GetFaceName(void)const " (?GetFaceName@wxFont@@UBE?AVwxString@@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::GetUnderlined(void)const " (?GetUnderlined@wxFont@@UBE_NXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetWeight(void)const " (?GetWeight@wxFont@@UBEHXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetStyle(void)const " (?GetStyle@wxFont@@UBEHXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetFamily(void)const " (?GetFamily@wxFont@@UBEHXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual int __thiscall wxFont::GetPointSize(void)const " (?GetPointSize@wxFont@@UBEHXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual unsigned long __thiscall wxFont::GetResourceHandle(void)" (?GetResourceHandle@wxFont@@UAEKXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::IsFree(void)const " (?IsFree@wxFont@@UBE_NXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::FreeResource(bool)" (?FreeResource@wxFont@@UAE_N_N@Z)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual bool __thiscall wxFont::RealizeResource(void)" (?RealizeResource@wxFont@@UAE_NXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual void __thiscall wxObject::CopyObject(class wxObject &)const " (?CopyObject@wxObject@@UBEXAAV1@@Z)
-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)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "protected: void __thiscall wxFont::Init(void)" (?Init@wxFont@@IAEXXZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: static class wxClassInfo  wxFont::sm_classwxFont" (?sm_classwxFont@wxFont@@2VwxClassInfo@@A)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxFont::~wxFont(void)" (??1wxFont@@UAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: __thiscall wxObject::wxObject(void)" (??0wxObject@@QAE@XZ)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: static class wxClassInfo  wxGDIObject::sm_classwxGDIObject" (?sm_classwxGDIObject@wxGDIObject@@2VwxClassInfo@@A)
-libctsim.lib(sgp.obj) : error LNK2001: unresolved external symbol "public: virtual __thiscall wxObject::~wxObject(void)" (??1wxObject@@UAE@XZ)
-Debug/ifexport.exe : fatal error LNK1120: 42 unresolved externals
-Error executing link.exe.
 
 
 
 <h3>Results</h3>
-ifexport.exe - 43 error(s), 0 warning(s)
-<h3>
---------------------Configuration: ifinfo - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-ifinfo.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: ifinfo - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-ifinfo.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: libpng - Win32 LIB--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-libpng.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: libpng - Win32 LIB Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-libpng.lib - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: phm2if - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-phm2if.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: phm2if - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-phm2if.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: phm2pj - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-phm2pj.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: phm2pj - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-phm2pj.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pj2if - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pj2if.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pj2if - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pj2if.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pjinfo - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pjinfo.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pjinfo - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pjinfo.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pjrec - Win32 Release--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pjrec.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: pjrec - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-pjrec.exe - 0 error(s), 0 warning(s)
-<h3>
---------------------Configuration: wxvc - Win32 LIB Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-
-
-
-<h3>Results</h3>
-wxd.lib - 0 error(s), 0 warning(s)
+ctsim.exe - 0 error(s), 0 warning(s)
 </pre>
 </body>
 </html>
diff --git a/src/ctsim.cpp b/src/ctsim.cpp
index 8dc7e22..172cda9 100644
--- a/src/ctsim.cpp
+++ b/src/ctsim.cpp
@@ -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
 
-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)
+  EVT_MENU(MAINMENU_FILE_CREATE_FILTER, MainFrame::OnCreateFilter)
   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...");
+    file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");
     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 
+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 
+MainFrame::OnCreateFilter (wxCommandEvent& WXUNUSED(event))
+{
+  DialogGetFilterParameters dialogFilter (this, 256, 256, SignalFilter::FILTER_BANDLIMIT, 1., 10., SignalFilter::DOMAIN_SPATIAL);
+    int dialogReturn = dialogFilter.ShowModal();
+    if (dialogReturn == wxID_OK) {
+      wxString strFilter (dialogFilter.getFilterName());
+      wxString strDomain (dialogFilter.getDomainName());
+      unsigned int nx = dialogFilter.getXSize();
+      unsigned int ny = dialogFilter.getYSize();
+      double dBandwidth = dialogFilter.getBandwidth();
+      double dFilterParam= dialogFilter.getFilterParam();
+      *theApp->getLog() << "Selected filter " << strFilter.c_str() << ", domain " << strDomain.c_str() << ", filterParam " << dFilterParam << ", bandwidth " << dBandwidth << "\n";
+      wxString filename = "untitled.if";
+      ImageFileDocument* pFilterDoc = dynamic_cast<ImageFileDocument*>(theApp->getDocManager()->CreateDocument ("untitled.if", wxDOC_SILENT));
+      if (! pFilterDoc) {
+        sys_error (ERR_SEVERE, "Unable to create filter image");
+        return;
+      }
+      ImageFile& rIF = pFilterDoc->getImageFile();
+      rIF.setArraySize (nx, ny);
+      rIF.filterResponse (strDomain.c_str(), dBandwidth, strFilter.c_str(), dFilterParam);
+      if (theApp->getSetModifyNewDocs())
+        pFilterDoc->Modify (true);
+      pFilterDoc->UpdateAllViews();
+      pFilterDoc->GetFirstView()->OnUpdate (NULL, NULL);
+    }
+}
+
 void
 CTSimApp::getCompatibleImages (const ImageFileDocument* pIFDoc, std::vector<ImageFileDocument*>& vecIF)
 {
diff --git a/src/ctsim.h b/src/ctsim.h
index b944560..2055227 100644
--- a/src/ctsim.h
+++ b/src/ctsim.h
@@ -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);
+    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,
+    MAINMENU_FILE_CREATE_FILTER,
     MAINMENU_FILE_EXIT,
     IFMENU_FILE_PROPERTIES,
     PJMENU_FILE_PROPERTIES,
@@ -154,6 +156,8 @@ enum {
   IFMENU_PROCESS_EXP,
   IFMENU_PROCESS_FOURIER,
   IFMENU_PROCESS_INVERSE_FOURIER,
+  IFMENU_PROCESS_FFT,
+  IFMENU_PROCESS_IFFT,
   IFMENU_PROCESS_MAGNITUDE,
   IFMENU_PROCESS_PHASE,
     PHMMENU_PROCESS_RASTERIZE,
diff --git a/src/dialogs.cpp b/src/dialogs.cpp
index 4a44735..3e2d8fb 100644
--- a/src/dialogs.cpp
+++ b/src/dialogs.cpp
@@ -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)
 }
 
 
+/////////////////////////////////////////////////////////////////////
+// CLASS IDENTIFICATION
+//
+// DialogGetFilterParameters
+/////////////////////////////////////////////////////////////////////
+
+
+DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDefaultXSize, int iDefaultYSize, int iDefaultFilterID, double dDefaultFilterParam,  double dDefaultBandwidth, int iDefaultDomainID)
+: wxDialog (pParent, -1, "Set Filter Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
+{
+  wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);
+  
+  pTopSizer->Add (new wxStaticText (this, -1, "Set Filter Parameters"), 0, wxALIGN_CENTER | wxTOP | wxLEFT | wxRIGHT, 5);
+  
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+  
+  std::ostringstream os;
+  os << iDefaultXSize;
+  m_pTextCtrlXSize = new wxTextCtrl (this, -1, os.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  std::ostringstream osYSize;
+  osYSize << iDefaultYSize;
+  m_pTextCtrlYSize = new wxTextCtrl (this, -1, osYSize.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  std::ostringstream osFilterParam;
+  osFilterParam << dDefaultFilterParam;
+  m_pTextCtrlFilterParam = new wxTextCtrl (this, -1, osFilterParam.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  std::ostringstream osBandwidth;
+  osBandwidth << dDefaultBandwidth;
+  m_pTextCtrlBandwidth = new wxTextCtrl (this, -1, osBandwidth.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  
+  wxFlexGridSizer* pGridSizer = new wxFlexGridSizer (2);
+  pGridSizer->Add (new wxStaticText (this, -1, "Filter"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5);
+  m_pListBoxFilter = new StringValueAndTitleListBox (this, SignalFilter::getFilterCount(), SignalFilter::getFilterTitleArray(), SignalFilter::getFilterNameArray());
+  m_pListBoxFilter->SetSelection (iDefaultFilterID);
+  pGridSizer->Add (m_pListBoxFilter, 0, wxALL | wxALIGN_LEFT | wxEXPAND);
+  
+  m_pListBoxDomain = new StringValueAndTitleListBox (this, SignalFilter::getDomainCount(), SignalFilter::getDomainTitleArray(), SignalFilter::getDomainNameArray());
+  m_pListBoxDomain->SetSelection (iDefaultDomainID);
+  pGridSizer->Add (new wxStaticText (this, -1, "Domain"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5);
+  pGridSizer->Add (m_pListBoxDomain, 0, wxALL | wxALIGN_LEFT | wxEXPAND);
+  
+  pGridSizer->Add (new wxStaticText (this, -1, "X Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlXSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (new wxStaticText (this, -1, "Y Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlYSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (new wxStaticText (this, -1, "Filter Parameter"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlFilterParam, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (new wxStaticText (this, -1, "Bandwidth"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlBandwidth, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  
+  pTopSizer->Add (pGridSizer, 1, wxALL, 3);
+  
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+  
+  wxBoxSizer* pButtonSizer = new wxBoxSizer (wxHORIZONTAL);
+  wxButton* pButtonOk = new wxButton (this, wxID_OK, "Okay");
+  wxButton* pButtonCancel = new wxButton (this, wxID_CANCEL, "Cancel");
+  pButtonSizer->Add (pButtonOk, 0, wxEXPAND | wxALL, 10);
+  pButtonSizer->Add (pButtonCancel, 0, wxEXPAND | wxALL, 10);
+  
+  pTopSizer->Add (pButtonSizer, 0, wxALIGN_CENTER);
+  
+  SetAutoLayout (true);
+  SetSizer (pTopSizer);
+  pTopSizer->Layout();
+  pTopSizer->Fit (this);
+  pTopSizer->SetSizeHints (this);
+}
+
+DialogGetFilterParameters::~DialogGetFilterParameters (void)
+{
+}
+
+
+unsigned int
+DialogGetFilterParameters::getXSize (void)
+{
+  wxString strCtrl = m_pTextCtrlXSize->GetValue();
+  unsigned long lValue;
+  if (strCtrl.ToULong (&lValue))
+    return lValue;
+  else
+    return (m_iDefaultXSize);
+}
+
+unsigned int
+DialogGetFilterParameters::getYSize (void)
+{
+  wxString strCtrl = m_pTextCtrlYSize->GetValue();
+  unsigned long lValue;
+  if (strCtrl.ToULong (&lValue))
+    return lValue;
+  else
+    return (m_iDefaultYSize);
+}
+
+double
+DialogGetFilterParameters::getBandwidth (void)
+{
+  wxString strCtrl = m_pTextCtrlBandwidth->GetValue();
+  double dValue;
+  if (strCtrl.ToDouble (&dValue))
+    return dValue;
+  else
+    return (m_dDefaultBandwidth);
+}
+
+double
+DialogGetFilterParameters::getFilterParam (void)
+{
+  wxString strCtrl = m_pTextCtrlFilterParam->GetValue();
+  double dValue;
+  if (strCtrl.ToDouble (&dValue))
+    return (dValue);
+  else
+    return (m_dDefaultFilterParam);
+}
+
+const char*
+DialogGetFilterParameters::getFilterName (void)
+{
+  return m_pListBoxFilter->getSelectionStringValue();
+}
+
+const char*
+DialogGetFilterParameters::getDomainName (void)
+{
+  return m_pListBoxDomain->getSelectionStringValue();
+}
+
diff --git a/src/dialogs.h b/src/dialogs.h
index 4441ca0..3546021 100644
--- a/src/dialogs.h
+++ b/src/dialogs.h
@@ -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
+{
+ 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 DialogGetFilterParameters : public wxDialog
+{
+ public:
+    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);
+    virtual ~DialogGetFilterParameters ();
+
+    unsigned int getXSize();
+    unsigned int getYSize();
+    const char* getFilterName();
+    const char* getDomainName();
+    double getFilterParam();
+    double getBandwidth();
+
+ private:
+    wxTextCtrl* m_pTextCtrlXSize;
+    wxTextCtrl* m_pTextCtrlYSize;
+    wxTextCtrl* m_pTextCtrlFilterParam;
+    wxTextCtrl* m_pTextCtrlBandwidth;
+
+    StringValueAndTitleListBox* m_pListBoxFilter;
+    StringValueAndTitleListBox* m_pListBoxDomain;
+
+    int m_iDefaultXSize;
+    int m_iDefaultYSize;
+    double m_dDefaultFilterParam;
+    double m_dDefaultBandwidth;
+    int m_iDefaultDomain;
+};
+
 class DialogAutoScaleParameters : public wxDialog
 {
  public:
diff --git a/src/views.cpp b/src/views.cpp
index 86e0cf5..93d49b3 100644
--- a/src/views.cpp
+++ b/src/views.cpp
@@ -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()) {
     if (pt.x >= 0 && pt.x < nx && pt.y >= 0 && pt.y < ny) {
       std::ostringstream os;
-      os << "Image value (" << pt.x << "," << yPt << ") = " << v[pt.x][yPt] << "\n";
+      os << "Image value (" << pt.x << "," << yPt << ") = " << v[pt.x][yPt];
+      if (rIF.isComplex()) {
+        double dImag = rIF.getImaginaryArray()[pt.x][yPt];
+        if (dImag < 0)
+          os << " - " << -dImag;
+        else
+          os << " + " << dImag;
+        os << "i\n";
+      } else
+        os << "\n";
       *theApp->getLog() << os.str().c_str();
     } else
       *theApp->getLog() << "Mouse out of image range (" << pt.x << "," << yPt << ")\n";
@@ -170,6 +179,10 @@ EVT_MENU(IFMENU_PROCESS_LOG, ImageFileView::OnLog)
 EVT_MENU(IFMENU_PROCESS_EXP, ImageFileView::OnExp)
 EVT_MENU(IFMENU_PROCESS_FOURIER, ImageFileView::OnFourier)
 EVT_MENU(IFMENU_PROCESS_INVERSE_FOURIER, ImageFileView::OnInverseFourier)
+#ifdef HAVE_FFTW
+EVT_MENU(IFMENU_PROCESS_FFT, ImageFileView::OnFFT)
+EVT_MENU(IFMENU_PROCESS_IFFT, ImageFileView::OnIFFT)
+#endif
 EVT_MENU(IFMENU_PROCESS_MAGNITUDE, ImageFileView::OnMagnitude)
 EVT_MENU(IFMENU_PROCESS_PHASE, ImageFileView::OnPhase)
 EVT_MENU(IFMENU_PLOT_ROW, ImageFileView::OnPlotRow)
@@ -195,16 +208,16 @@ ImageFileView::OnProperties (wxCommandEvent& event)
     const std::string& rFilename = rIF.getFilename();
     std::ostringstream os;
     double min, max, mean, mode, median, stddev;
-    rIF.statistics (min, max, mean, mode, median, stddev);
+    rIF.statistics (rIF.getArray(), min, max, mean, mode, median, stddev);
     os << "Filename: " << rFilename << "\n";
     os << "Size: (" << rIF.nx() << "," << rIF.ny() << ")\n";
     os << "Data type: ";
-    if (rIF.dataType() == Array2dFile::DATA_TYPE_COMPLEX)
+    if (rIF.isComplex())
       os << "Complex\n";
     else
       os << "Real\n";
     os << "\nMinimum: "<<min<<"\nMaximum: "<<max<<"\nMean: "<<mean<<"\nMedian: "<<median<<"\nMode: "<<mode<<"\nStandard Deviation: "<<stddev << "\n";
-    if (rIF.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
+    if (rIF.isComplex()) {
       rIF.statistics (rIF.getImaginaryArray(), min, max, mean, mode, median, stddev);
       os << "\nImaginary: min: "<<min<<"\nmax: "<<max<<"\nmean: "<<mean<<"\nmedian: "<<median<<"\nmode: "<<mode<<"\nstddev: "<<stddev << "\n";
     }
@@ -279,7 +292,7 @@ ImageFileView::OnCompare (wxCommandEvent& event)
       rIF.statistics (min, max, mean, mode, median, stddev);
       os << rIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
       rCompareIF.statistics (min, max, mean, mode, median, stddev);
-      os << rCompareIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
+      os << pCompareDoc->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
       os << "\n";
       double d, r, e;
       rIF.comparativeStatistics (rCompareIF, d, r, e);
@@ -363,6 +376,7 @@ void
 ImageFileView::OnFourier (wxCommandEvent& event)
 {
   ImageFile& rIF = GetDocument()->getImageFile();
+  wxProgressDialog dlgProgress (wxString("Fourier"), wxString("Fourier Progress"), 1, m_frame, wxPD_APP_MODAL);
   rIF.fourier (rIF);
   m_bMinSpecified = false;
   m_bMaxSpecified = false;
@@ -371,10 +385,39 @@ ImageFileView::OnFourier (wxCommandEvent& event)
   GetDocument()->UpdateAllViews(this);
 }
 
+#ifdef HAVE_FFTW
+void
+ImageFileView::OnFFT (wxCommandEvent& event)
+{
+  ImageFile& rIF = GetDocument()->getImageFile();
+  wxProgressDialog dlgProgress (wxString("FFT"), wxString("FFT Progress"), 1, m_frame, wxPD_APP_MODAL);
+  rIF.fft (rIF);
+  m_bMinSpecified = false;
+  m_bMaxSpecified = false;
+  if (theApp->getSetModifyNewDocs())
+    GetDocument()->Modify(TRUE);
+  GetDocument()->UpdateAllViews(this);
+}
+
+void
+ImageFileView::OnIFFT (wxCommandEvent& event)
+{
+  ImageFile& rIF = GetDocument()->getImageFile();
+  wxProgressDialog dlgProgress (wxString("IFFT"), wxString("IFFT Progress"), 1, m_frame, wxPD_APP_MODAL);
+  rIF.ifft (rIF);
+  m_bMinSpecified = false;
+  m_bMaxSpecified = false;
+  if (theApp->getSetModifyNewDocs())
+    GetDocument()->Modify(TRUE);
+  GetDocument()->UpdateAllViews(this);
+}
+#endif
+
 void
 ImageFileView::OnInverseFourier (wxCommandEvent& event)
 {
   ImageFile& rIF = GetDocument()->getImageFile();
+  wxProgressDialog dlgProgress (wxString("Inverse Fourier"), wxString("Inverse Fourier Progress"), 1, m_frame, wxPD_APP_MODAL);
   rIF.inverseFourier (rIF);
   m_bMinSpecified = false;
   m_bMaxSpecified = false;
@@ -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...");
+  file_menu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...");
   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");
   process_menu->Append (IFMENU_PROCESS_EXP, "&Exp");
   process_menu->AppendSeparator();
+#ifdef HAVE_FFTW
+  process_menu->Append (IFMENU_PROCESS_FFT, "&FFT");
+  process_menu->Append (IFMENU_PROCESS_IFFT, "&IFFT");
+  process_menu->Append (IFMENU_PROCESS_FOURIER, "F&ourier");
+  process_menu->Append (IFMENU_PROCESS_INVERSE_FOURIER, "Inverse Fo&urier");
+#else
   process_menu->Append (IFMENU_PROCESS_FOURIER, "&Fourier");
   process_menu->Append (IFMENU_PROCESS_INVERSE_FOURIER, "&Inverse Fourier");
+#endif
   process_menu->Append (IFMENU_PROCESS_MAGNITUDE, "&Magnitude");
   process_menu->Append (IFMENU_PROCESS_PHASE, "&Phase");
 
@@ -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...");
   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...");
   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...");
   file_menu->Append(wxID_OPEN, "&Open...");
   file_menu->Append(wxID_SAVE, "&Save");
   file_menu->Append(wxID_SAVEAS, "Save &As...");
diff --git a/src/views.h b/src/views.h
index 3512fac..8d46ae3 100644
--- a/src/views.h
+++ b/src/views.h
@@ -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);
     void OnFourier (wxCommandEvent& event);
     void OnInverseFourier (wxCommandEvent& event);
+#ifdef HAVE_FFTW
+    void OnFFT (wxCommandEvent& event);
+    void OnIFFT (wxCommandEvent& event);
+#endif
     void OnMagnitude (wxCommandEvent& event);
     void OnPhase (wxCommandEvent& event);
     void OnScaleAuto (wxCommandEvent& event);