libctsim/fourier.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          fourier.cpp
   5 **   Purpose:       Fourier transform functions
   6 **   Programmer:    Kevin Rosenberg
   7 **   Date Started:  Dec 2000
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (c) 1983-2001 Kevin Rosenberg
  11 **
  12 **  $Id: fourier.cpp,v 1.5 2001/03/18 18:08:25 kevin Exp $
  13 **
  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.
  17 **
  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.
  22 **
  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 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30
  31
  32 void
  33 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)
  34 {
  35   ImageFileArray vReal = im.getArray();
  36   ImageFileArray vImag = im.getImaginaryArray();
  37   unsigned int ix, iy;
  38   unsigned int nx = im.nx();
  39   unsigned int ny = im.ny();
  40
  41   // shuffle each column
  42   for (ix = 0; ix < nx; ix++) {
  43     Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
  44     if (im.isComplex())
  45       Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);
  46   }
  47
  48   // shuffle each row
  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];
  56     if (im.isComplex()) {
  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];
  62     }
  63   }
  64   delete pRow;
  65 }
  66
  67 void
  68 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
  69 {
  70   ImageFileArray vReal = im.getArray();
  71   ImageFileArray vImag = im.getImaginaryArray();
  72   unsigned int ix, iy;
  73   unsigned int nx = im.nx();
  74   unsigned int ny = im.ny();
  75
  76   // shuffle each x column
  77   for (ix = 0; ix < nx; ix++) {
  78     Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
  79     if (im.isComplex())
  80       Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);
  81   }
  82
  83   // shuffle each y row
  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];
  91     if (im.isComplex()) {
  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];
  97     }
  98   }
  99   delete [] pRow;
 100 }
 101
 102
 103 // Odd Number of Points
 104 //   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
 105 //   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
 106 // Even Number of Points
 107 //   Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
 108 //   Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
 109
 110 void
 111 Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
 112 {
 113   double* pdTemp = new double [n];
 114   int i;
 115   if (isOdd(n)) { // Odd
 116     int iHalfN = (n - 1) / 2;
 117
 118     pdTemp[0] = pdVector[iHalfN];
 119     for (i = 0; i < iHalfN; i++)
 120       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
 121     for (i = 0; i < iHalfN; i++)
 122       pdTemp[i + iHalfN + 1] = pdVector[i];
 123   } else {     // Even
 124     int iHalfN = n / 2;
 125     pdTemp[0] = pdVector[iHalfN];
 126     for (i = 0; i < iHalfN - 1; i++)
 127       pdTemp[i + 1] = pdVector[i + iHalfN + 1];
 128     for (i = 0; i < iHalfN; i++)
 129       pdTemp[i + iHalfN] = pdVector[i];
 130   }
 131
 132   for (i = 0; i < n; i++)
 133     pdVector[i] = pdTemp[i];
 134   delete pdTemp;
 135 }
 136
 137 void
 138 Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
 139 {
 140   std::complex<double>* pdTemp = new std::complex<double> [n];
 141   int i;
 142   if (isOdd(n)) { // Odd
 143     int iHalfN = (n - 1) / 2;
 144
 145     pdTemp[0] = pdVector[iHalfN];
 146     for (i = 0; i < iHalfN; i++)
 147       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
 148     for (i = 0; i < iHalfN; i++)
 149       pdTemp[i + iHalfN + 1] = pdVector[i];
 150   } else {     // Even
 151     int iHalfN = n / 2;
 152     pdTemp[0] = pdVector[iHalfN];
 153     for (i = 0; i < iHalfN - 1; i++)
 154       pdTemp[i + 1] = pdVector[i + iHalfN + 1];
 155     for (i = 0; i < iHalfN; i++)
 156       pdTemp[i + iHalfN] = pdVector[i];
 157   }
 158
 159   for (i = 0; i < n; i++)
 160     pdVector[i] = pdTemp[i];
 161   delete [] pdTemp;
 162 }
 163
 164
 165 void
 166 Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)
 167 {
 168   float* pdTemp = new float [n];
 169   int i;
 170   if (isOdd (n)) { // Odd
 171     int iHalfN = (n - 1) / 2;
 172
 173     pdTemp[0] = pdVector[iHalfN];
 174     for (i = 0; i < iHalfN; i++)
 175       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
 176     for (i = 0; i < iHalfN; i++)
 177       pdTemp[i + iHalfN + 1] = pdVector[i];
 178   } else {     // Even
 179     int iHalfN = n / 2;
 180     pdTemp[0] = pdVector[iHalfN];
 181     for (i = 0; i < iHalfN - 1; i++)
 182       pdTemp[i + 1] = pdVector[i + iHalfN + 1];
 183     for (i = 0; i < iHalfN; i++)
 184       pdTemp[i + iHalfN] = pdVector[i];
 185   }
 186
 187   for (i = 0; i < n; i++)
 188     pdVector[i] = pdTemp[i];
 189   delete pdTemp;
 190 }
 191
 192
 193
 194 void
 195 Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 196 {
 197   double* pdTemp = new double [n];
 198   int i;
 199   if (isOdd(n)) { // Odd
 200     int iHalfN = (n - 1) / 2;
 201
 202     pdTemp[iHalfN] = pdVector[0];
 203     for (i = 0; i < iHalfN; i++)
 204       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
 205     for (i = 0; i < iHalfN; i++)
 206       pdTemp[i] = pdVector[i + iHalfN + 1];
 207   } else {     // Even
 208     int iHalfN = n / 2;
 209     pdTemp[iHalfN] = pdVector[0];
 210     for (i = 0; i < iHalfN; i++)
 211       pdTemp[i] = pdVector[i + iHalfN];
 212     for (i = 0; i < iHalfN - 1; i++)
 213       pdTemp[i + iHalfN + 1] = pdVector[i+1];
 214   }
 215
 216   for (i = 0; i < n; i++)
 217     pdVector[i] = pdTemp[i];
 218   delete pdTemp;
 219 }
 220
 221
 222 void
 223 Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
 224 {
 225   std::complex<double>* pdTemp = new std::complex<double> [n];
 226   int i;
 227   if (isOdd(n)) { // Odd
 228     int iHalfN = (n - 1) / 2;
 229
 230     pdTemp[iHalfN] = pdVector[0];
 231     for (i = 0; i < iHalfN; i++)
 232       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
 233     for (i = 0; i < iHalfN; i++)
 234       pdTemp[i] = pdVector[i + iHalfN + 1];
 235   } else {     // Even
 236     int iHalfN = n / 2;
 237     pdTemp[iHalfN] = pdVector[0];
 238     for (i = 0; i < iHalfN; i++)
 239       pdTemp[i] = pdVector[i + iHalfN];
 240     for (i = 0; i < iHalfN - 1; i++)
 241       pdTemp[i + iHalfN + 1] = pdVector[i+1];
 242   }
 243
 244   for (i = 0; i < n; i++)
 245     pdVector[i] = pdTemp[i];
 246   delete [] pdTemp;
 247 }
 248
 249
 250
 251
 252 void
 253 Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)
 254 {
 255   float* pTemp = new float [n];
 256   int i;
 257   if (isOdd (n)) { // Odd
 258     int iHalfN = (n - 1) / 2;
 259
 260     pTemp[iHalfN] = pVector[0];
 261     for (i = 0; i < iHalfN; i++)
 262       pTemp[i + 1 + iHalfN] = pVector[i + 1];
 263     for (i = 0; i < iHalfN; i++)
 264       pTemp[i] = pVector[i + iHalfN + 1];
 265   } else {     // Even
 266     int iHalfN = n / 2;
 267     pTemp[iHalfN] = pVector[0];
 268     for (i = 0; i < iHalfN; i++)
 269       pTemp[i] = pVector[i + iHalfN];
 270     for (i = 0; i < iHalfN - 1; i++)
 271       pTemp[i + iHalfN + 1] = pVector[i+1];
 272   }
 273
 274   for (i = 0; i < n; i++)
 275     pVector[i] = pTemp[i];
 276   delete [] pTemp;
 277 }
 278
 279