libctsim/fourier.cpp

   1 /*****************************************************************************\r
   2 ** FILE IDENTIFICATION\r
   3 **\r
   4 **   Name:          fourier.cpp\r
   5 **   Purpose:       Fourier transform functions\r
   6 **   Programmer:    Kevin Rosenberg\r
   7 **   Date Started:  Dec 2000\r
   8 **\r
   9 **  This is part of the CTSim program\r
  10 **  Copyright (C) 1983-2001 Kevin Rosenberg\r
  11 **\r
  12 **  $Id: fourier.cpp,v 1.1 2001/01/02 06:33:04 kevin Exp $\r
  13 **\r
  14 **  This program is free software; you can redistribute it and/or modify\r
  15 **  it under the terms of the GNU General Public License (version 2) as\r
  16 **  published by the Free Software Foundation.\r
  17 **\r
  18 **  This program is distributed in the hope that it will be useful,\r
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  21 **  GNU General Public License for more details.\r
  22 **\r
  23 **  You should have received a copy of the GNU General Public License\r
  24 **  along with this program; if not, write to the Free Software\r
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA\r
  26 ******************************************************************************/\r
  27 \r
  28 #include "ct.h"\r
  29 \r
  30 \r
  31 \r
  32 void\r
  33 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)\r
  34 {\r
  35   ImageFileArray vReal = im.getArray();\r
  36   ImageFileArray vImag = im.getImaginaryArray();\r
  37   unsigned int ix, iy;\r
  38   unsigned int nx = im.nx();\r
  39   unsigned int ny = im.ny();\r
  40 \r
  41   // shuffle each column\r
  42   for (ix = 0; ix < nx; ix++) {\r
  43     Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);\r
  44     if (im.isComplex())\r
  45       Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);\r
  46   }\r
  47 \r
  48   // shuffle each row\r
  49   float* pRow = new float [nx];\r
  50   for (iy = 0; iy < ny; iy++) {\r
  51     for (ix = 0; ix < nx; ix++)\r
  52       pRow[ix] = vReal[ix][iy];\r
  53     Fourier::shuffleFourierToNaturalOrder (pRow, nx);\r
  54     for (ix = 0; ix < nx; ix++)\r
  55       vReal[ix][iy] = pRow[ix];\r
  56     if (im.isComplex()) {\r
  57       for (ix = 0; ix < nx; ix++)\r
  58         pRow[ix] = vImag[ix][iy];\r
  59       Fourier::shuffleFourierToNaturalOrder (pRow, nx);;\r
  60       for (ix = 0; ix < nx; ix++)\r
  61         vImag[ix][iy] = pRow[ix];\r
  62     }\r
  63   }\r
  64   delete pRow;\r
  65 }\r
  66  \r
  67 void\r
  68 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)\r
  69 {\r
  70   ImageFileArray vReal = im.getArray();\r
  71   ImageFileArray vImag = im.getImaginaryArray();\r
  72   unsigned int ix, iy;\r
  73   unsigned int nx = im.nx();\r
  74   unsigned int ny = im.ny();\r
  75 \r
  76   // shuffle each x column\r
  77   for (ix = 0; ix < nx; ix++) {\r
  78     Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);\r
  79     if (im.isComplex())\r
  80       Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);\r
  81   }\r
  82 \r
  83   // shuffle each y row\r
  84   float* pRow = new float [nx];\r
  85   for (iy = 0; iy < ny; iy++) {\r
  86     for (ix = 0; ix < nx; ix++)\r
  87       pRow[ix] = vReal[ix][iy];\r
  88     Fourier::shuffleNaturalToFourierOrder (pRow, nx);\r
  89     for (ix = 0; ix < nx; ix++)\r
  90       vReal[ix][iy] = pRow[ix];\r
  91     if (im.isComplex()) {\r
  92       for (ix = 0; ix < nx; ix++)\r
  93         pRow[ix] = vImag[ix][iy];\r
  94       Fourier::shuffleNaturalToFourierOrder (pRow, nx);\r
  95       for (ix = 0; ix < nx; ix++)\r
  96         vImag[ix][iy] = pRow[ix];\r
  97     }\r
  98   }\r
  99   delete [] pRow;\r
 100 }\r
 101 \r
 102  \r
 103 // Odd Number of Points\r
 104 //   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2\r
 105 //   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1\r
 106 // Even Number of Points\r
 107 //   Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)\r
 108 //   Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1\r
 109 \r
 110 void\r
 111 Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)\r
 112 {\r
 113   double* pdTemp = new double [n];\r
 114   int i;\r
 115   if (n % 2) { // Odd\r
 116     int iHalfN = (n - 1) / 2;\r
 117     \r
 118     pdTemp[0] = pdVector[iHalfN];\r
 119     for (i = 0; i < iHalfN; i++)\r
 120       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
 121     for (i = 0; i < iHalfN; i++)\r
 122       pdTemp[i + iHalfN + 1] = pdVector[i];\r
 123   } else {     // Even\r
 124     int iHalfN = n / 2;\r
 125     pdTemp[0] = pdVector[iHalfN];\r
 126     for (i = 0; i < iHalfN - 1; i++)\r
 127       pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
 128     for (i = 0; i < iHalfN; i++)\r
 129       pdTemp[i + iHalfN] = pdVector[i];\r
 130   }\r
 131   \r
 132   for (i = 0; i < n; i++)\r
 133     pdVector[i] = pdTemp[i];\r
 134   delete pdTemp;\r
 135 }\r
 136 \r
 137 void\r
 138 Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)\r
 139 {\r
 140   std::complex<double>* pdTemp = new std::complex<double> [n];\r
 141   int i;\r
 142   if (n % 2) { // Odd\r
 143     int iHalfN = (n - 1) / 2;\r
 144     \r
 145     pdTemp[0] = pdVector[iHalfN];\r
 146     for (i = 0; i < iHalfN; i++)\r
 147       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
 148     for (i = 0; i < iHalfN; i++)\r
 149       pdTemp[i + iHalfN + 1] = pdVector[i];\r
 150   } else {     // Even\r
 151     int iHalfN = n / 2;\r
 152     pdTemp[0] = pdVector[iHalfN];\r
 153     for (i = 0; i < iHalfN - 1; i++)\r
 154       pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
 155     for (i = 0; i < iHalfN; i++)\r
 156       pdTemp[i + iHalfN] = pdVector[i];\r
 157   }\r
 158   \r
 159   for (i = 0; i < n; i++)\r
 160     pdVector[i] = pdTemp[i];\r
 161   delete [] pdTemp;\r
 162 }\r
 163 \r
 164 \r
 165 void\r
 166 Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)\r
 167 {\r
 168   float* pdTemp = new float [n];\r
 169   int i;\r
 170   if (n % 2) { // Odd\r
 171     int iHalfN = (n - 1) / 2;\r
 172     \r
 173     pdTemp[0] = pdVector[iHalfN];\r
 174     for (i = 0; i < iHalfN; i++)\r
 175       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
 176     for (i = 0; i < iHalfN; i++)\r
 177       pdTemp[i + iHalfN + 1] = pdVector[i];\r
 178   } else {     // Even\r
 179     int iHalfN = n / 2;\r
 180     pdTemp[0] = pdVector[iHalfN];\r
 181     for (i = 0; i < iHalfN - 1; i++)\r
 182       pdTemp[i + 1] = pdVector[i + iHalfN + 1];\r
 183     for (i = 0; i < iHalfN; i++)\r
 184       pdTemp[i + iHalfN] = pdVector[i];\r
 185   }\r
 186   \r
 187   for (i = 0; i < n; i++)\r
 188     pdVector[i] = pdTemp[i];\r
 189   delete pdTemp;\r
 190 }\r
 191 \r
 192 \r
 193 \r
 194 void\r
 195 Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)\r
 196 {\r
 197   double* pdTemp = new double [n];\r
 198   int i;\r
 199   if (n % 2) { // Odd\r
 200     int iHalfN = (n - 1) / 2;\r
 201     \r
 202     pdTemp[iHalfN] = pdVector[0];\r
 203     for (i = 0; i < iHalfN; i++)\r
 204       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
 205     for (i = 0; i < iHalfN; i++)\r
 206       pdTemp[i] = pdVector[i + iHalfN + 1];\r
 207   } else {     // Even\r
 208     int iHalfN = n / 2;\r
 209     pdTemp[iHalfN] = pdVector[0];\r
 210     for (i = 0; i < iHalfN; i++)\r
 211       pdTemp[i] = pdVector[i + iHalfN];\r
 212     for (i = 0; i < iHalfN - 1; i++)\r
 213       pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
 214   }\r
 215   \r
 216   for (i = 0; i < n; i++)\r
 217     pdVector[i] = pdTemp[i];\r
 218   delete pdTemp;\r
 219 }\r
 220 \r
 221 \r
 222 void\r
 223 Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)\r
 224 {\r
 225   std::complex<double>* pdTemp = new std::complex<double> [n];\r
 226   int i;\r
 227   if (n % 2) { // Odd\r
 228     int iHalfN = (n - 1) / 2;\r
 229     \r
 230     pdTemp[iHalfN] = pdVector[0];\r
 231     for (i = 0; i < iHalfN; i++)\r
 232       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
 233     for (i = 0; i < iHalfN; i++)\r
 234       pdTemp[i] = pdVector[i + iHalfN + 1];\r
 235   } else {     // Even\r
 236     int iHalfN = n / 2;\r
 237     pdTemp[iHalfN] = pdVector[0];\r
 238     for (i = 0; i < iHalfN; i++)\r
 239       pdTemp[i] = pdVector[i + iHalfN];\r
 240     for (i = 0; i < iHalfN - 1; i++)\r
 241       pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
 242   }\r
 243   \r
 244   for (i = 0; i < n; i++)\r
 245     pdVector[i] = pdTemp[i];\r
 246   delete [] pdTemp;\r
 247 }\r
 248 \r
 249 \r
 250 \r
 251 \r
 252 void\r
 253 Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)\r
 254 {\r
 255   float* pTemp = new float [n];\r
 256   int i;\r
 257   if (n % 2) { // Odd\r
 258     int iHalfN = (n - 1) / 2;\r
 259     \r
 260     pTemp[iHalfN] = pVector[0];\r
 261     for (i = 0; i < iHalfN; i++)\r
 262       pTemp[i + 1 + iHalfN] = pVector[i + 1];\r
 263     for (i = 0; i < iHalfN; i++)\r
 264       pTemp[i] = pVector[i + iHalfN + 1];\r
 265   } else {     // Even\r
 266     int iHalfN = n / 2;\r
 267     pTemp[iHalfN] = pVector[0];\r
 268     for (i = 0; i < iHalfN; i++)\r
 269       pTemp[i] = pVector[i + iHalfN];\r
 270     for (i = 0; i < iHalfN - 1; i++)\r
 271       pTemp[i + iHalfN + 1] = pVector[i+1];\r
 272   }\r
 273   \r
 274   for (i = 0; i < n; i++)\r
 275     pVector[i] = pTemp[i];\r
 276   delete [] pTemp;\r
 277 }\r
 278 \r
 279 \r