1 /*****************************************************************************
5 ** Purpose: Fourier transform functions
6 ** Programmer: Kevin Rosenberg
7 ** Date Started: Dec 2000
9 ** This is part of the CTSim program
10 ** Copyright (c) 1983-2001 Kevin Rosenberg
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License (version 2) as
16 ** published by the Free Software Foundation.
18 ** This program is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 ** GNU General Public License for more details.
23 ** You should have received a copy of the GNU General Public License
24 ** along with this program; if not, write to the Free Software
25 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 ******************************************************************************/
33 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)
35 ImageFileArray vReal = im.getArray();
36 ImageFileArray vImag = im.getImaginaryArray();
38 unsigned int nx = im.nx();
39 unsigned int ny = im.ny();
41 // shuffle each column
42 for (ix = 0; ix < nx; ix++) {
43 Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
45 Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);
49 float* pRow = new float [nx];
50 for (iy = 0; iy < ny; iy++) {
51 for (ix = 0; ix < nx; ix++)
52 pRow[ix] = vReal[ix][iy];
53 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
54 for (ix = 0; ix < nx; ix++)
55 vReal[ix][iy] = pRow[ix];
57 for (ix = 0; ix < nx; ix++)
58 pRow[ix] = vImag[ix][iy];
59 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
60 for (ix = 0; ix < nx; ix++)
61 vImag[ix][iy] = pRow[ix];
68 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
70 ImageFileArray vReal = im.getArray();
71 ImageFileArray vImag = im.getImaginaryArray();
73 unsigned int nx = im.nx();
74 unsigned int ny = im.ny();
76 // shuffle each x column
77 for (ix = 0; ix < nx; ix++) {
78 Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
80 Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);
84 float* pRow = new float [nx];
85 for (iy = 0; iy < ny; iy++) {
86 for (ix = 0; ix < nx; ix++)
87 pRow[ix] = vReal[ix][iy];
88 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
89 for (ix = 0; ix < nx; ix++)
90 vReal[ix][iy] = pRow[ix];
92 for (ix = 0; ix < nx; ix++)
93 pRow[ix] = vImag[ix][iy];
94 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
95 for (ix = 0; ix < nx; ix++)
96 vImag[ix][iy] = pRow[ix];
103 void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
105 fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
108 if (isOdd(n)) { // Odd
109 int iHalfN = (n - 1) / 2;
111 pTemp[0][0] = pVector[iHalfN][0];
112 pTemp[0][1] = pVector[iHalfN][1];
113 for (i = 0; i < iHalfN; i++) {
114 pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
115 pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
117 for (i = 0; i < iHalfN; i++) {
118 pTemp[i + iHalfN + 1][0] = pVector[i][0];
119 pTemp[i + iHalfN + 1][1] = pVector[i][1];
123 pTemp[0][0] = pVector[iHalfN][0];
124 pTemp[0][1] = pVector[iHalfN][1];
125 for (i = 0; i < iHalfN - 1; i++) {
126 pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
127 pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
129 for (i = 0; i < iHalfN; i++) {
130 pTemp[i + iHalfN][0] = pVector[i][0];
131 pTemp[i + iHalfN][1] = pVector[i][1];
135 for (i = 0; i < n; i++) {
136 pVector[i][0] = pTemp[i][0];
137 pVector[i][1] = pTemp[i][1];
142 void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
144 fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
146 if (isOdd(n)) { // Odd
147 int iHalfN = (n - 1) / 2;
149 pTemp[iHalfN][0] = pVector[0][0];
150 pTemp[iHalfN][1] = pVector[0][1];
151 for (i = 0; i < iHalfN; i++) {
152 pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
153 pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
155 for (i = 0; i < iHalfN; i++) {
156 pTemp[i][0] = pVector[i + iHalfN + 1][0];
157 pTemp[i][1] = pVector[i + iHalfN + 1][1];
161 pTemp[iHalfN][0] = pVector[0][0];
162 pTemp[iHalfN][1] = pVector[0][1];
163 for (i = 0; i < iHalfN; i++) {
164 pTemp[i][0] = pVector[i + iHalfN][0];
165 pTemp[i][1] = pVector[i + iHalfN][1];
167 for (i = 0; i < iHalfN - 1; i++) {
168 pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
169 pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
173 for (i = 0; i < n; i++) {
174 pVector[i][0] = pTemp[i][0];
175 pVector[i][1] = pTemp[i][1];