-/*****************************************************************************\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
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: fourier.cpp
+** Purpose: Fourier transform functions
+** Programmer: Kevin Rosenberg
+** Date Started: Dec 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2001 Kevin Rosenberg
+**
+** $Id: fourier.cpp,v 1.2 2001/01/02 16:02:13 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
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+
+void
+Fourier::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++) {
+ Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
+ if (im.isComplex())
+ Fourier::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];
+ Fourier::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];
+ Fourier::shuffleFourierToNaturalOrder (pRow, nx);;
+ for (ix = 0; ix < nx; ix++)
+ vImag[ix][iy] = pRow[ix];
+ }
+ }
+ delete pRow;
+}
+
+void
+Fourier::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++) {
+ Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
+ if (im.isComplex())
+ Fourier::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];
+ Fourier::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];
+ Fourier::shuffleNaturalToFourierOrder (pRow, nx);
+ for (ix = 0; ix < nx; ix++)
+ vImag[ix][iy] = pRow[ix];
+ }
+ }
+ delete [] pRow;
+}
+
+
+// Odd Number of Points
+// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
+// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
+// Even Number of Points
+// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
+// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
+
+void
+Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
+{
+ double* pdTemp = new double [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];
+ 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];
+ }
+
+ for (i = 0; i < n; i++)
+ pdVector[i] = pdTemp[i];
+ delete pdTemp;
+}
+
+void
+Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
+{
+ std::complex<double>* pdTemp = new std::complex<double> [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];
+ 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];
+ }
+
+ for (i = 0; i < n; i++)
+ pdVector[i] = pdTemp[i];
+ delete [] pdTemp;
+}
+
+
+void
+Fourier::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];
+ 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];
+ }
+
+ for (i = 0; i < n; i++)
+ pdVector[i] = pdTemp[i];
+ delete pdTemp;
+}
+
+
+
+void
+Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
+{
+ double* pdTemp = new double [n];
+ int i;
+ if (n % 2) { // Odd
+ int iHalfN = (n - 1) / 2;
+
+ pdTemp[iHalfN] = pdVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i] = pdVector[i + iHalfN + 1];
+ } else { // Even
+ int iHalfN = n / 2;
+ pdTemp[iHalfN] = pdVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i] = pdVector[i + iHalfN];
+ for (i = 0; i < iHalfN - 1; i++)
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];
+ }
+
+ for (i = 0; i < n; i++)
+ pdVector[i] = pdTemp[i];
+ delete pdTemp;
+}
+
+
+void
+Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
+{
+ std::complex<double>* pdTemp = new std::complex<double> [n];
+ int i;
+ if (n % 2) { // Odd
+ int iHalfN = (n - 1) / 2;
+
+ pdTemp[iHalfN] = pdVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i] = pdVector[i + iHalfN + 1];
+ } else { // Even
+ int iHalfN = n / 2;
+ pdTemp[iHalfN] = pdVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pdTemp[i] = pdVector[i + iHalfN];
+ for (i = 0; i < iHalfN - 1; i++)
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];
+ }
+
+ for (i = 0; i < n; i++)
+ pdVector[i] = pdTemp[i];
+ delete [] pdTemp;
+}
+
+
+
+
+void
+Fourier::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;
+}
+
+
projects / ctsim.git / blobdiff