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$
  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 #ifdef HAVE_FFTW
 103 void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
 104 {
 105   fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
 106   int i;
 107
 108   if (isOdd(n)) { // Odd
 109     int iHalfN = (n - 1) / 2;
 110
 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];
 116     }
 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];
 120     }
 121   } else {     // Even
 122     int iHalfN = n / 2;
 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];
 128     }
 129     for (i = 0; i < iHalfN; i++) {
 130       pTemp[i + iHalfN][0] = pVector[i][0];
 131       pTemp[i + iHalfN][1] = pVector[i][1];
 132     }
 133   }
 134
 135   for (i = 0; i < n; i++) {
 136     pVector[i][0] = pTemp[i][0];
 137     pVector[i][1] = pTemp[i][1];
 138   }
 139   fftw_free(pTemp);
 140 }
 141
 142 void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
 143 {
 144   fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
 145   int i;
 146   if (isOdd(n)) { // Odd
 147     int iHalfN = (n - 1) / 2;
 148
 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];
 154     }
 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];
 158     }
 159   } else {     // Even
 160     int iHalfN = n / 2;
 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];
 166     }
 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];
 170     }
 171   }
 172
 173   for (i = 0; i < n; i++) {
 174     pVector[i][0] = pTemp[i][0];
 175     pVector[i][1] = pTemp[i][1];
 176   }
 177
 178   fftw_free(pTemp);
 179 }
 180 #endif
 181