--- /dev/null
+/*****************************************************************************\r
+** FILE IDENTIFICATION\r
+**\r
+** Name: fourier.h\r
+** Purpose: Header for Fourier transform functions\r
+** Programmer: Kevin Rosenberg\r
+** Date Started: Dec 2000\r
+**\r
+** This is part of the CTSim program\r
+** Copyright (C) 1983-2001 Kevin Rosenberg\r
+**\r
+** $Id: fourier.h,v 1.1 2001/01/02 06:33:04 kevin Exp $\r
+**\r
+** This program is free software; you can redistribute it and/or modify\r
+** it under the terms of the GNU General Public License (version 2) as\r
+** published by the Free Software Foundation.\r
+**\r
+** This program is distributed in the hope that it will be useful,\r
+** but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+** GNU General Public License for more details.\r
+**\r
+** You should have received a copy of the GNU General Public License\r
+** along with this program; if not, write to the Free Software\r
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+******************************************************************************/\r
+\r
+#include <complex>\r
+\r
+class ImageFile;\r
+\r
+class Fourier {\r
+public:\r
+ static void shuffleFourierToNaturalOrder (ImageFile& im);\r
+ static void shuffleNaturalToFourierOrder (ImageFile& im);\r
+\r
+ static void shuffleNaturalToFourierOrder (float* pdVector, const int n);\r
+ static void shuffleNaturalToFourierOrder (double* pdVector, const int n);\r
+ static void shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n);\r
+ static void shuffleFourierToNaturalOrder (float* pdVector, const int n);\r
+ static void shuffleFourierToNaturalOrder (double* pdVector, const int n);\r
+ static void shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n);\r
+\r
+}; // namespace Fourier\r
+\r
--- /dev/null
+/*****************************************************************************\r
+** FILE IDENTIFICATION\r
+**\r
+** Name: fourier.cpp\r
+** Purpose: Fourier transform functions\r
+** Programmer: Kevin Rosenberg\r
+** Date Started: Dec 2000\r
+**\r
+** This is part of the CTSim program\r
+** Copyright (C) 1983-2001 Kevin Rosenberg\r
+**\r
+** $Id: fourier.cpp,v 1.1 2001/01/02 06:33:04 kevin Exp $\r
+**\r
+** This program is free software; you can redistribute it and/or modify\r
+** it under the terms of the GNU General Public License (version 2) as\r
+** published by the Free Software Foundation.\r
+**\r
+** This program is distributed in the hope that it will be useful,\r
+** but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+** GNU General Public License for more details.\r
+**\r
+** You should have received a copy of the GNU General Public License\r
+** along with this program; if not, write to the Free Software\r
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+******************************************************************************/\r
+\r
+#include "ct.h"\r
+\r
+\r
+\r
+void\r
+Fourier::shuffleFourierToNaturalOrder (ImageFile& im)\r
+{\r
+ ImageFileArray vReal = im.getArray();\r
+ ImageFileArray vImag = im.getImaginaryArray();\r
+ unsigned int ix, iy;\r
+ unsigned int nx = im.nx();\r
+ unsigned int ny = im.ny();\r
+\r
+ // shuffle each column\r
+ for (ix = 0; ix < nx; ix++) {\r
+ Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);\r
+ if (im.isComplex())\r
+ Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);\r
+ }\r
+\r
+ // shuffle each row\r
+ float* pRow = new float [nx];\r
+ for (iy = 0; iy < ny; iy++) {\r
+ for (ix = 0; ix < nx; ix++)\r
+ pRow[ix] = vReal[ix][iy];\r
+ Fourier::shuffleFourierToNaturalOrder (pRow, nx);\r
+ for (ix = 0; ix < nx; ix++)\r
+ vReal[ix][iy] = pRow[ix];\r
+ if (im.isComplex()) {\r
+ for (ix = 0; ix < nx; ix++)\r
+ pRow[ix] = vImag[ix][iy];\r
+ Fourier::shuffleFourierToNaturalOrder (pRow, nx);;\r
+ for (ix = 0; ix < nx; ix++)\r
+ vImag[ix][iy] = pRow[ix];\r
+ }\r
+ }\r
+ delete pRow;\r
+}\r
+ \r
+void\r
+Fourier::shuffleNaturalToFourierOrder (ImageFile& im)\r
+{\r
+ ImageFileArray vReal = im.getArray();\r
+ ImageFileArray vImag = im.getImaginaryArray();\r
+ unsigned int ix, iy;\r
+ unsigned int nx = im.nx();\r
+ unsigned int ny = im.ny();\r
+\r
+ // shuffle each x column\r
+ for (ix = 0; ix < nx; ix++) {\r
+ Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);\r
+ if (im.isComplex())\r
+ Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);\r
+ }\r
+\r
+ // shuffle each y row\r
+ float* pRow = new float [nx];\r
+ for (iy = 0; iy < ny; iy++) {\r
+ for (ix = 0; ix < nx; ix++)\r
+ pRow[ix] = vReal[ix][iy];\r
+ Fourier::shuffleNaturalToFourierOrder (pRow, nx);\r
+ for (ix = 0; ix < nx; ix++)\r
+ vReal[ix][iy] = pRow[ix];\r
+ if (im.isComplex()) {\r
+ for (ix = 0; ix < nx; ix++)\r
+ pRow[ix] = vImag[ix][iy];\r
+ Fourier::shuffleNaturalToFourierOrder (pRow, nx);\r
+ for (ix = 0; ix < nx; ix++)\r
+ vImag[ix][iy] = pRow[ix];\r
+ }\r
+ }\r
+ delete [] pRow;\r
+}\r
+\r
+ \r
+// Odd Number of Points\r
+// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2\r
+// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1\r
+// Even Number of Points\r
+// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)\r
+// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1\r
+\r
+void\r
+Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)\r
+{\r
+ double* pdTemp = new double [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN] = pdVector[i];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete pdTemp;\r
+}\r
+\r
+void\r
+Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)\r
+{\r
+ std::complex<double>* pdTemp = new std::complex<double> [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN] = pdVector[i];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete [] pdTemp;\r
+}\r
+\r
+\r
+void\r
+Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)\r
+{\r
+ float* pdTemp = new float [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN] = pdVector[i];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete pdTemp;\r
+}\r
+\r
+\r
+\r
+void\r
+Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)\r
+{\r
+ double* pdTemp = new double [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN + 1];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete pdTemp;\r
+}\r
+\r
+\r
+void\r
+Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)\r
+{\r
+ std::complex<double>* pdTemp = new std::complex<double> [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN + 1];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete [] pdTemp;\r
+}\r
+\r
+\r
+\r
+\r
+void\r
+Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)\r
+{\r
+ float* pTemp = new float [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pTemp[iHalfN] = pVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pTemp[i + 1 + iHalfN] = pVector[i + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pTemp[i] = pVector[i + iHalfN + 1];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pTemp[iHalfN] = pVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pTemp[i] = pVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pTemp[i + iHalfN + 1] = pVector[i+1];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pVector[i] = pTemp[i];\r
+ delete [] pTemp;\r
+}\r
+\r
+\r